Program Integrals;
{$E-,N+,G+,D-}
Uses CRT;
var
a, k : Extended;
a0, b0, hx, hz, z0, ze, x0, xe : Extended;
i : Longint;
Pi : Extended;
Function Pi_ : Extended;Assembler;
asm
fldpi
end;
Procedure Nops;
inline($90/$90/$90/$90/$90);
Function Ro(Theta : Extended) : Extended;
begin
Ro := a*sqr(cos(Theta));
end;
Function Roasm(Theta : Extended) : Extended;Assembler;
asm
fld a {a}
fld Theta {Theta\a}
db 0D9h, 0FFh{fcos : 80387} {cos(Theta)\a}
fld st(0) {cos\cos\a}
fmul {eq fmul st(0), st(1) | sqr(cos)\a}
fmul {a*sqr(cos(Theta))}
end;
Function koeff : Extended;Assembler;
asm
fld a {a}
fld1 {1\a}
fld1 {1\1\a}
fld1 {1\1\1\a}
fadd {2\1\a}
fdiv st(2), st(0) {2\1\(a/2)}
fadd st(1), st(0) {2\3\(a/2)}
fxch {3\2\(a/2)}
fdivp st(1), st(0) {(2/3)\(a/2)}
f2xm1 {2^(2/3)-1\(a/2)}
fsqrt {sqrt(2^(2/3)-1)\(a/2)}
fmul st(0), st(1) {sqrt(2^(2/3)-1)*(a/2)}
end;
Function sqr_r(z : Extended) : Extended;
begin
sqr_r := exp( (1/3)*ln(sqr(a*z*z)) ) - z*z;
end;
Function r(z : Extended) : Extended;
begin
r := sqrt( exp( (1/3)*ln(sqr(a*z*z)) ) - z*z );
end;
Function y(x, z : Extended) : Extended;
begin
y := sqrt(exp((1/3)*ln(sqr(a*z*z)))-x*x-z*z);
end;
Function Trapezoids(N : Longint): Extended;
var
S, z : Extended;
i, j : Longint;
begin
hz := (ze-z0)/N;{hz := a/N}
{ xe := r(z0);{}
{ x0 := -xe;}
S := y(0, ze){0};
for i := 1 to N-1 do
begin
z := z0+hz*i;
xe := r(z);{}
x0 := -xe;
{Ничего не прибавляем - это нули}
hx := (xe-x0)/N;{2*r(z0+hz*i)/N}
for j := 1 to N-1 do
S := S + 2*y(x0+j*hx, z)*hz*hx;{2 - две поверхности}
end;
Trapezoids := S;
end;
Function Simpson(N : Longint): Extended;
var
S, z : Extended;
i, j : Longint;
begin
S := 0;
hz := (ze-z0)/N;
for i := 0 to (N div 2 - 1) do
begin
z := z0+hz*(2*i+1);
xe := r(z);{}
x0 := -xe;
{Ничего не прибавляем - это нули}
hx := (xe-x0)/N;{2*r(z0+hz*i)/N}
for j := 0 to (N div 2 - 1) do
S := S + 8*y(x0+(2*j+1)*hx, z)*hz*hx/3;{2 - две поверхности}
for j := 1 to (N div 2 - 1) do
S := S + 4*y(x0+2*j*hx, z)*hz*hx/3;{2 - две поверхности}
end;
for i := 1 to (N div 2 - 1) do
begin
z := z0+hz*2*i;
xe := r(z);{}
x0 := -xe;
{Ничего не прибавляем - это нули}
hx := (xe-x0)/N;{2*r(z0+hz*i)/N}
for j := 0 to (N div 2 - 1) do
S := S + 8*y(x0+(2*j+1)*hx, z)*hz*hx/3;{2 - две поверхности}
for j := 1 to (N div 2 - 1) do
S := S + 4*y(x0+2*j*hx, z)*hz*hx/3;{2 - две поверхности}
end;
Simpson := S;
end;
Function MonteCarlo(N : Longint): Extended;
var
i : Longint;
x, y, z : Extended;
finds : Longint;
begin
finds := 0;
for i := 1 to N do
begin
x := k*(2*Random-1.0);
y := k*(2*Random-1.0);
z := a*Random;
if sqr(y)<=(exp((1/3)*ln(sqr(a*z*z)))-x*x-z*z) then inc(finds);
end;
MonteCarlo := 4*k*k*a*finds/N;
end;
Function MonteCarlosSantana(N : Longint): Extended;
var
i : Longint;
x, y, z : Extended;
finds : Longint;
begin
finds := 0;
for i := 1 to N do
begin
x := k*(2*Random-1.0);
y := k*(2*Random-1.0);
z := a*Random;
if sqrt(x*x+y*y)<=r(z) then inc(finds);
end;
MonteCarlosSantana := 4*k*k*a*finds/N;
end;
Function Trapezza(N : Longint): Extended;
var
S : Extended;
i : Longint;
begin
hz := (ze-z0)/N;
S := sqr_r(ze);
for i := 1 to N-1 do
S := S + sqr_r(z0+i*hz)*hz;
Trapezza := Pi*S;
end;
Function TheSimpsons(N : Longint): Extended;
var
S : Extended;
i : Longint;
begin
hz := (ze-z0)/N;
S := sqr_r(ze)*hz/3;
for i := 0 to (N div 2-1) do
S := S + 4*sqr_r(z0+(2*i+1)*hz)*hz/3;
for i := 1 to (N div 2-1) do
S := S + 2*sqr_r(z0+2*i*hz)*hz/3;
TheSimpsons := Pi*S;
end;
{Function CMonteCarlosSantana(N : Longint): Extended;
var
i : Longint;
rr, z : Extended;
finds : Longint;
begin
finds := 0;
for i := 1 to N do
begin
z := a*Random;
rr := k*Random;
if rr*rr<=sqr(r(z)) then inc(finds);
end;
CMonteCarlosSantana := finds/N;
end;}
Function Counter : Longint;Assembler;
asm
push ds
xor ax, ax
mov ds, ax
mov bx, 046Ch
mov ax, word ptr [bx]
@q:
cmp ax, word ptr [bx]
je @q
mov ax, word ptr [bx]
mov dx, word ptr [bx+2]
pop ds
end;
Function PowerIncrement : Longint;Assembler;
asm
push ds
xor ax, ax
mov ds, ax
mov bx, 046Ch
xor dx, dx
xor cx, cx
mov ax, word ptr [bx]
@q:
cmp ax, word ptr [bx]
je @q
mov ax, word ptr [bx]
@q1:
inc cx
jnz @q1
inc dx
cmp ax, word ptr [bx]
je @q1
mov ax, cx
pop ds
end;
Function PowerAddition : Longint;Assembler;
asm
push ds
xor ax, ax
mov ds, ax
xor dx, dx
xor cx, cx
mov bx, 046Ch
mov ax, word ptr [bx]
@q:
cmp ax, word ptr [bx]
je @q
mov ax, word ptr [bx]
@q1:
add cx, 1
adc dx, 0
cmp ax, word ptr [bx]
je @q1
mov ax, cx
pop ds
end;
Function PowerFPUAdd : Extended;Assembler;
asm
push ds
xor ax, ax
mov ds, ax
mov bx, 046Ch
mov ax, word ptr [bx]
@q:
cmp ax, word ptr [bx]
je @q
fld1
fld1
mov ax, word ptr [bx]
@q1:
fadd st(0), st(1)
cmp ax, word ptr [bx]
je @q1
faddp st(1), st(0)
fld1
fsubp st(1), st(0)
pop ds
end;
var
Ctr : Longint;
begin
ClrScr;
Randomize;
Pi := 0.0;
Pi := Pi_;
a := 1.0;
Nops;
k := koeff;{a*sqrt(2^(2/3)-1)/2}
WriteLn('k= ', k:1:20);
WriteLn;
z0 := 0;
ze := a;
for i := 1 to 5 do
WriteLn('Команды "add\adc": ', (1000/55)*PowerAddition:1:0, ' операций в секунду');
for i := 1 to 5 do
WriteLn('Команда FPU "fadd": ', (1000/55)*PowerFPUAdd:1:0, ' операций в секунду');
for i := 1 to 5 do
WriteLn('Команды "inc\inc": ', (1000/55)*PowerIncrement:1:0, ' операций в секунду');
WriteLn;
WriteLn('Теоретическое значение: ', '':0, (2/21)*Pi*a*a*a:1:20);
WriteLn;
{ Ctr := Counter;
Write('Метод трапеций: ', '':8, Trapezoids(5000):1:20);
WriteLn('':5, (Counter-Ctr)*55/1000:1:3, ' секунд');
Ctr := Counter;
Write('Метод Симпсона: ', '':8, Simpson(5000):1:20);
WriteLn('':5, (Counter-Ctr)*55/1000:1:3, ' секунд');
}
Ctr := Counter;
Write('Метод трапеций: ', '':8, Trapezoids(2000):1:20);
WriteLn('':5, (Counter-Ctr)*55/1000:1:3, ' секунд');
Ctr := Counter;
Write('Метод Симпсона: ', '':8, Simpson(2000):1:20);
WriteLn('':5, (Counter-Ctr)*55/1000:1:3, ' секунд');
Ctr := Counter;
Write('Метод Монте-Карло: ', '':5, MonteCarlo(10000000):1:20);
WriteLn('':5, (Counter-Ctr)*55/1000:1:3, ' секунд');
WriteLn;
Ctr := Counter;
Write('Метод трапецца: ', '':8, Trapezza(10000000):1:20);
WriteLn('':5, (Counter-Ctr)*55/1000:1:3, ' секунд');
Ctr := Counter;
Write('Метод Симпсонов: ', '':7, TheSimpsons(10000000):1:20);
WriteLn('':5, (Counter-Ctr)*55/1000:1:3, ' секунд');
Ctr := Counter;
Write('Метод Карлоса: ', '':9, MonteCarlosSantana(10000000):5:20);
WriteLn('':5, (Counter-Ctr)*55/1000:1:3, ' секунд');
{ WriteLn('Метод цилиндрч Карлоса: ', '':0, CMonteCarlosSantana(100000):5:20);}
ReadKey;
end.
< -->