Subtraktion mit DECMath bringt Probleme

[Luckie]

Die aktualisierte Fassung wird gerade bei mir hochgeladen und dürfte in fünf Minuten dort zur Verfügung stehen.

[Tortus]

Hm,.. kann es sein, das für NSqrt(2) auch nur 1 zurückgegeben wird?

[negaH]

Mach mich nicht schwach, ich werds die Nacht testen.

Gruß Hagen

[Tortus]

Zitat von negaH:

Mach mich nicht schwach, ich werds die Nacht testen.

Gruß Hagen

Ich wunder mich ja auch schon, is ja schließlich nicht so, dass ich jetzt der Profiprogrammierer bin

[negaH]

Du könntest mir einen Gefallen machen, und ich werde die Wurzel noch zum laufen bringen.

Da du ja schon am testen bist, eventuell kannst du in der Zwischenzeit schnell noch die anderen Funktionen testen. So brauche ich nicht jedesmal 6 Mb an Luckie mailen, falls du noch andere Fehler findest.

Gruß Hagen

[Tortus]

Hm,.. wenn du schon dabei bist, gibts evtl eine Möglichkeit für Exponentialfunktionen mit Nicht-Integer-Exponenten?

/E: Ansonsten sieht eigentlich bislang alles recht vielversprechend aus

[negaH]

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;

[negaH]

Zitat:

Hm,.. wenn du schon dabei bist, gibts evtl eine Möglichkeit für Exponentialfunktionen mit Nicht-Integer-Exponenten?

was ist mit

Delphi-Quellcode:

procedure NExp(var A: IRational; const B: IRational);
// A = exp(B)

procedure NExp(var A: IRational);
// A = exp(A)

in Unit NRats ? Öffne mal \DEC\Part_II\LibIntf\NRats.pas

Gruß Hagen