|
Registriert seit: 25. Jun 2003
Ort: Thüringen
2.950 Beiträge
|
Re: Subtraktion mit DECMath bringt Probleme
15. Mär 2006, 23:41
Schau mal in Unit NInt_1.pas rein. Dort ist mit Sourcen die NExp() Funktion
// A = A * e^(U / V)
procedure NExp(var A: IInteger; const U,V: IInteger);
Allerdings bevor ich die die nächste ZIP fertig mache musst du begreifen das bei Berechnungen mit IRational() und auch zb. NExp() immer mit Ungenauigkeiten zu rechnen sind. Das ist einfach in der Mathematik begründet wenn man mit supergroßen Brüchen rechnen will und nicht exakt deren erforderliche Genauigkeiten vorherberechnen kann. Dh. in deinen Berechnungen musst du auch immer die erforderliche Genauigkeit der Zahlen berechnen und dementsprechend im DEC einstellen.
Bei den IRational gibt es dazu zwei Möglichkeiten
1.) const NRats.DefaultPrecision: Cardinal = 1024; Vorgabegenauigkeit in Bits.
2.) NPrc(var A: IRational; Precision: Cardinal = 0; Base: TBase = 10): Cardinal; liest oder liest und setzt die Genauigkeit einer Bruchvariable A:IRational.
Es gibt zwei Aufrufmöglichkeiten
Delphi-Quellcode:
// 1.) Abfrage der eingestellten Genauigkeit in Bits
if NPrc(A) > 100 then ;
// 2.) setzen der Genauigkeit auf 100 Dezimalstellen
NPrc(A, 100, 10);
// 2a.) setzen auf Genauigkeit von 100 Bits
NPrc(A, 100, 2);
So ich packe nun das ZIP nochmal zusammen. Ich weis nicht wie der Fehler in NSqrt() da reingekommen ist, wahrscheinlich beim "Schönmachen" der Sources.
Bei weiteren Fehler (schon peinlich genug) müssen wir diese aber erstmal auf eine Todo Liste setzen. Ich kann nicht jedesmal Luckie zumailen.
Dann noch zu deine Pi Berechnung:
Ich denke du hast bestimmt schon NPi() in Unit NInt_1.pas als Source angechaut ?? Falls nicht \ DEC\PART_II\NInt_1.pas öffnen.
Und hier eine experimentelle Version die nach dem gleichen Algo. arbeitet aber statt rekursive iterativ arbeitet. Die hatte ich letzte Zeit mal in der Zerre 
Gruß Hagen
Delphi-Quellcode:
procedure NPiX( var A: IInteger; Decimals: Cardinal);
function NSum( var N,D: IInteger; Precision: Cardinal): Cardinal;
var
S: array[0..3*32-1] of IInteger;
I,J,K: Cardinal;
C,E,P,Q,Alpha,Beta,Gamma,Delta,Eps,T,U: IInteger;
begin
Result := 0;
NSet(E, 10939058860032000);
NSet(Alpha, 1);
NSet(Beta, 1);
NSet(Gamma, 5);
NSet(Delta, 53360);
NSet(Eps, 13591409);
NSet(Q, 1);
K := 0;
for I := 1 to (Precision + 197) div 94 do
begin
NMul(S[K], Alpha, Beta);
NMul(S[K], Gamma);
NSwp(S[K +1], Delta);
NInc(Alpha, 2);
NInc(Beta, 6);
NInc(Gamma, 6);
NAdd(P, Q);
NInc(Q, 2);
NAdd(C, E);
NMul(Delta, P, C);
NMul(T, Eps, Delta);
NInc(Eps, 545140134);
NMul(U, S[K], Eps);
NSub(S[K +2], T, U);
NMul(T, Alpha, Beta);
NMul(U, S[K], Gamma);
NMul(S[K], T, U);
NMul(S[K +1], Delta);
NInc(Alpha, 2);
NInc(Beta, 6);
NInc(Gamma, 6);
NAdd(P, Q);
NInc(Q, 2);
NAdd(C, E);
NMul(Delta, P, C);
NInc(Eps, 545140134);
J := I;
while not Odd(J) do
begin
J := J shr 1;
NMul(T, S[K +1], S[K -1]);
NMul(S[K -1], S[K -3], S[K +2]);
NAdd(S[K -1], T); // S[k-3] * S[k +2] + S[k +1] * S[k -1]
NMul(S[K -3], S[K]); // S[k-3] * S[k]
NMul(S[K -2], S[K +1]); // S[k-2] * S[k +1]
Dec(K, 3);
end;
Inc(K, 3);
end;
Dec(K, 3);
while K <> 0 do
begin
Dec(K, 3);
NMul(T, S[K +4], S[K +2]);
NMul(S[K +2], S[K], S[K +5]);
NAdd(S[K +2], T); // S[k] * S[k +5] + S[k +2] * S[k +4]
NMul(S[K +1], S[K +4]);
end;
NSwp(N, S[1]);
NSwp(D, S[2]);
end;
var
N,D: IInteger;
begin
NSet(A, 5);
NPow(A, Decimals);
NSum(N, D, Cardinal(NSize(A)) + Decimals);
NSqr(A);
NMul(A, 640320);
NShl(A, Decimals + Decimals);
NSqrt(A);
NMul(A, N);
NDiv(A, D);
end;
function NPi( var A: IInteger; Decimals: Cardinal; Method: TIIntegerPIMethod): Cardinal;
{ timings on PIV 1500Mhz 512Mb, Win2k in seconds
Decimals : Chud(fast) Chud(iter) AGM Machin(fast) Machin(Iter) CRC
10000 : 0.02 0.05 0.06 0.07 0.23 0022C012
25000 : 0.07 0.33 0.27 0.28 1.44 C95E3B65
50000 : 0.21 1.35 0.79 0.83 5.70 D4509617
100000 : 0.57 5.94 2.32 2.20 22.76 364124EB
250000 : 2.23 44.02 8.08 7.87 145.99 A6211796
500000 : 5.53 200.37 19.49 19.66 783.36 5BCC3354
1000000 : 13.77 841.04 54.99 50.40 3805.81 38845D51
2000000 : 33.26 3384.13 129.67 126.60 - BE0EE27F
2500000 : 46.40 5293.82 182.96 176.08 - 64E69944
4000000 : 82.07 - 354.42 321.70 - A5A8FE00
5000000 : 121.13 - 439.85 459.02 - 672A62A5
8000000 : 198.16 - 825.77 808.90 - FD180141
16000000 : 483.48 - 2196.44 2041.08 - 7F242222
32000000 : 1151.34 - 5061.61 4775.54 - 742593FD
64000000 : 3288.20 - - - - D60B96F3
}
procedure NPi_Iter_Machin( var R: IInteger; Decimals: Cardinal);
// Machin's Arctan series, iterative method
procedure AddArcTan( var R,P: IInteger; X,M: Integer);
var
N,D: IInteger;
K: Cardinal;
begin
NMul(N, P, X * M);
X := X * X +1;
NSet(D, X);
X := X + X;
K := 0;
repeat
NDiv(N, D);
if NSgn(N) = 0 then Break;
Inc(K, 2);
NAdd(R, N);
NMul(N, K);
NAdd(D, X);
until False;
end;
var
I: Cardinal;
P: IInteger;
begin
I := System.Trunc(Ln(Decimals) / Ln(10)) + 2;
NPow(P, 10, Decimals + I);
NSet(R, 0);
// AddArcTan(Self, P, 10, 8 * 4);
// AddArcTan(Self, P, 239, -1 * 4);
// AddArcTan(Self, P, 515, -4 * 4);
AddArcTan(R, P, 18, 12 * 4);
AddArcTan(R, P, 57, 8 * 4);
AddArcTan(R, P, 239, -5 * 4);
NPow(P, 10, I);
NDiv(R, P);
end;
procedure NPi_Fast_Machin( var R: IInteger; Decimals: Integer);
procedure NArcTanSeries( var R: IInteger; const S: array of Integer);
var
P,T: IInteger;
I: Integer;
begin
NShl(R, 32);
I := 0;
repeat
NSet(P, R);
NArcTan(P, S[I +1]);
NMul(P, S[I + 0]);
NAdd(T, P);
Inc(I, 2);
until I >= Length(S);
NShr(R, T, 30);
end;
begin
NPow(R, 5, Decimals);
NShl(R, Decimals);
// pi/4 = 44*arctan(1/57) +7*arctan(1/239) -12*arctan(1/682) +24*arctan(1/12943)
NArcTanSeries(R, [+44, 57, +7, 239, -12, 682, +24, 12943]);
// NArcTanSeries(R, [+1, 2, +1, 3]);
// NArcTanSeries(R, [+12, 18, +8, 57, -5, 239]);
// NArcTanSeries(R, [+4, 5, -1, 239]);
// NArcTanSeries(R, [+6, 8, +2, 57, +1, 239]);
// NArcTanSeries(R, [+88, 111, +7, 239, -44, 515, +32, 682, +24, 12943]);
// NArcTanSeries(R, [+322, 577, +76, 682, +139, 1393, +156, 12943, +132, 32807, +44, 1049433]);
end;
procedure NPi_AGM( var R: IInteger; Decimals: Cardinal);
{AGM start with:
a = 1, b = 1/Sqrt(2), t = 1/2, k = 1
iteration:
s = (a + b) / 2
d = a - s
d = d^2
a = s
s = s^2
t = t - d * 2^k
b = Sqrt(s - d)
k = k +1
final:
pi ~ (a + b)^2 / 2t }
var
A,B,D,T: IInteger;
W: Integer;
begin
(*
---- a formula based on the AGM (Arithmetic-Geometric Mean) ----
c = sqrt(0.125);
a = 1 + 3 * c;
b = sqrt(a);
e = b - 0.625;
b = 2 * b;
c = e - c;
a = a + e;
npow = 4;
do {
npow = 2 * npow;
e = (a + b) / 2;
b = sqrt(a * b);
e = e - b;
b = 2 * b;
c = c - e;
a = e + b;
} while (e > SQRT_SQRT_EPSILON);
e = e * e / 4;
a = a + b;
pi = (a * a - e - e / 2) / (a * c - e) / npow;
*)
Inc(Decimals, 3); // +3 digits reserve
NPow(R, 5, Decimals); // R = 5^Decimals
NShl(A, R, Decimals); // A = 10^Decimals
NShl(D, R, Decimals -1); // D = 10^(Decimals -1)^2
NSqr(D);
NSqr(B, A); // B = (10^Decimals^2 * 2)^0.5 div 2
NShl(B, 1);
NSqrt(B);
NShr(B, 1);
W := 0;
repeat
NAdd(T, A, B); // T = (A + B) div 2, new A
NShr(T, 1);
NMul(B, A); // B = (B * A)^0.5
NSqrt(B);
NSub(A, T); // A = (A - T)^2 * 2^W
NSqr(A);
NShl(A, W);
Inc(W);
NSub(D, A); // D = D - A
NSwp(A, T);
until NCmp(B, A) = 0;
NShr(D, Decimals);
NDiv(D, R);
NMul(D, 1000);
NSqr(R, A);
NDiv(R, D);
end;
procedure NPi_Iter_Chudnovsky( var R: IInteger; Decimals: Cardinal);
// Chudnovsky's brothers Ramanujan iterative computation
// 12 (6n)! * (A + n * B)
// 1/PI = ------- * sum(n = 0..invinite) * (-1)^n ------------------------
// C^(3/2) n!^3 * (3 * n)! * C^3^n
//
const
k_A = 13591409;
k_B = 545140134;
k_C = 53360;
k_1 = 100100025; // k_1 * k_2 = 32817176580096000 = (12 * C)^3 / 8
k_2 = 327843840;
var
N: Integer;
P,S,T: IInteger;
begin
NPow(T, 10, Decimals + 6); // T = 10^Deciamls * 10^6 --> +6 Digit correction
NSqr(R, T); // R = 53360 * Sqrt(640320 * T^2)
NMul(R, k_C * 12);
NSqrt(R);
NMul(R, k_C);
NMul(R, T);
NMul(T, 1000000); // -6 digits correct
NMul(P, T, k_A); // P = T * k_A, P = sum[]
N := 1;
while NSgn(T) > 0 do
begin
NSet(S, 6 * N -5); // T = T * (6n -5)(6n -3)(6n -1)
NMul(S, 6 * N -3);
NMul(S, 6 * N -1);
NMul(T, S);
NPow(S, N, 3); // T = T / n^3 * 32817176580096000
NMul(S, k_1);
NMul(S, k_2);
NDiv(T, S);
NSet(S, k_B); // P = P +- T * (k_A + n * k_B)
NMul(S, N);
NAdd(S, k_A);
NMul(S, T);
NNeg(S, Odd(N)); // (-1)^N --> -1 = odd(N) or 1 = not odd(N)
NAdd(P, S);
Inc(N);
end;
NDiv(R, P); // PI == (R * 10^3) div (P * 10^6)
end;
procedure NPi_Fast_Chudnovsky( var R: IInteger; Decimals: Cardinal);
// Chudnovsky's brothers Ramanujan with binary splitting
// this code is most only 2-3 Times slower as the fastests special PI-programs
// Some optimizations can be done, but here i wanted a readable code.
// One way to made it faster (~2 times) could be to avoid NBinarySplitting() and its
// callback.
// Use same formula as above.
{
1 12 inf (-1)^n (6n)! (A+nB)
-- = ------ SUM: -------------------
pi C^1.5 n=0 (3n)! (n!)^3 C^(3n)
A=13591409 B=545140134 C=640320
a(n) = (A+nB)
p(n) = -1(6n-5)(6n-4)(6n-3)(6n-2)(n6-1)(6n)
b(n) = 1
q(n) = (3n-2)(3n-1)(3n)(n^3)(C^3)
}
const
k_A = 163096908; // 12 * 13591409
k_B = 6541681608; // 12 * 545140134
k_C = 53360; // 640320 / 12
k_D = 1728; // 12^3 = 24 * 72
k_E = 262537412640768000; // 1727 * k_C^3
procedure DoPI(N: Integer; var D: TIIntegerSplitData); register;
// the callback for binary splitting
// max. n < ~306783379 then 6n << MaxInt
// B[1] = 1, P[0] = 1, Q[0] = 1
// don't touch B,P[0],Q[0] because NIL interfaces are here assumed to value +1
begin
if N > 0 then
begin
NSet(D.A, k_B); // a = k_A + n * k_B
NMul(D.A, N);
NAdd(D.A, k_A);
NSet(D.P, 2 * N -1 );
NMul(D.P, 6 * N -1 );
NMul(D.P, -(6 * N -5)); // p = - (6*n -5) * (6*n -1)* (2*n -1)
NSet(D.Q, k_C); // q = 72(n * k_C)^3 // 72 * n^3 * k_C^3
NMul(D.Q, N);
NPow(D.Q, 3);
NMul(D.Q, 72);
end else
NSet(D.A, k_A); // A[0] = k_A
end;
var
P,Q: IInteger;
S: Cardinal;
begin
S := Trunc(Decimals / 14.181) +2;
S := NBinarySplitting(P, Q, S, @DoPI, False); // decimal approximation is very roughly !!
NPow(R, 100, Decimals);
NMul(R, NInt(k_E));
NSqrt(R); // R = (k_C^3 * 12^3 * 100^Decimals)^(1/2)
NMul(R, Q);
NShl(R, S);
NDiv(R, P); // R = R * Q / P = R / S
end;
procedure NPi_Fast_Chudnovsky2( var R: IInteger; Decimals: Cardinal);
// another bin. split. chudnovsky
const
k_A = 13591409;
k_B = 545140134;
k_C = 53360;
procedure DoPI(N: Integer; var D: TIIntegerSplitData); register;
begin
Inc(N);
NSet(D.A, k_B);
NMul(D.A, N);
NAdd(D.A, k_A);
NNeg(D.A, Odd(N));
NSet(D.Q, N);
NMul(D.Q, N);
NMul(D.Q, N);
NMul(D.Q, (k_C div 2) * (k_C div 2));
NMul(D.Q, k_C * 12 * 12 * 2);
NSet(D.P, N + N -1);
N := N * 6;
NMul(D.P, N -1);
NMul(D.P, N -5);
end;
var
P,Q: IInteger;
S: Cardinal;
begin
Inc(Decimals, 3);
S := NBinarySplitting(P, Q, Trunc(Decimals / 14.181) +2, @DoPI, False); // decimal approximation is very roughly !!
// R := (10^(Decimals*2) * 12 * k_C)^0.5 * k_C * Q / (P + Q * k_A)
NMul(R, Q, k_A);
NShl(R, S);
NAdd(P, R);
NPow(R, 5, Decimals + Decimals);
NMul(R, k_C * 12);
NShl(R, Decimals + Decimals);
NSqrt(R);
NMul(R, k_C);
NMul(R, Q);
NShl(R, S);
NMul(P, 1000);
NDiv(R, P);
end;
begin
Result := 1;
case Method of
piFastChudnovsky : NPi_Fast_Chudnovsky(A, Decimals);
piIterChudnovsky : NPi_Iter_Chudnovsky(A, Decimals);
piAGM : NPi_AGM(A, Decimals);
piIterMachin : NPi_Iter_Machin(A, Decimals);
piFastMachin : NPi_Fast_Machin(A, Decimals);
end;
end;
|