Delphi-PRAXiS - Delphi Subtraktion mit DECMath bringt Probleme

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;