AGB  ·  Datenschutz  ·  Impressum  







Anmelden
Nützliche Links
Registrieren
Zurück Delphi-PRAXiS Sprachen und Entwicklungsumgebungen Sonstige Fragen zu Delphi Delphi Subtraktion mit DECMath bringt Probleme
Thema durchsuchen
Ansicht
Themen-Optionen

Subtraktion mit DECMath bringt Probleme

Ein Thema von Tortus · begonnen am 15. Mär 2006 · letzter Beitrag vom 16. Mär 2006
Antwort Antwort
Seite 2 von 2     12   
Benutzerbild von Luckie
Luckie

Registriert seit: 29. Mai 2002
37.621 Beiträge
 
Delphi 2006 Professional
 
#11

Re: Subtraktion mit DECMath bringt Probleme

  Alt 15. Mär 2006, 21:22
Die aktualisierte Fassung wird gerade bei mir hochgeladen und dürfte in fünf Minuten dort zur Verfügung stehen.
Michael
Ein Teil meines Codes würde euch verunsichern.
  Mit Zitat antworten Zitat
Tortus

Registriert seit: 15. Nov 2003
Ort: Gescher
47 Beiträge
 
Delphi 7 Enterprise
 
#12

Re: Subtraktion mit DECMath bringt Probleme

  Alt 15. Mär 2006, 21:45
Hm,.. kann es sein, das für NSqrt(2) auch nur 1 zurückgegeben wird?
Thorsten Lanfer
  Mit Zitat antworten Zitat
Benutzerbild von negaH
negaH

Registriert seit: 25. Jun 2003
Ort: Thüringen
2.950 Beiträge
 
#13

Re: Subtraktion mit DECMath bringt Probleme

  Alt 15. Mär 2006, 21:49
Mach mich nicht schwach, ich werds die Nacht testen.

Gruß Hagen
  Mit Zitat antworten Zitat
Tortus

Registriert seit: 15. Nov 2003
Ort: Gescher
47 Beiträge
 
Delphi 7 Enterprise
 
#14

Re: Subtraktion mit DECMath bringt Probleme

  Alt 15. Mär 2006, 21:53
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
Thorsten Lanfer
  Mit Zitat antworten Zitat
Benutzerbild von negaH
negaH

Registriert seit: 25. Jun 2003
Ort: Thüringen
2.950 Beiträge
 
#15

Re: Subtraktion mit DECMath bringt Probleme

  Alt 15. Mär 2006, 22:23
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
  Mit Zitat antworten Zitat
Tortus

Registriert seit: 15. Nov 2003
Ort: Gescher
47 Beiträge
 
Delphi 7 Enterprise
 
#16

Re: Subtraktion mit DECMath bringt Probleme

  Alt 15. Mär 2006, 22:57
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
Thorsten Lanfer
  Mit Zitat antworten Zitat
Benutzerbild von negaH
negaH

Registriert seit: 25. Jun 2003
Ort: Thüringen
2.950 Beiträge
 
#17

Re: Subtraktion mit DECMath bringt Probleme

  Alt 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;
  Mit Zitat antworten Zitat
Benutzerbild von negaH
negaH

Registriert seit: 25. Jun 2003
Ort: Thüringen
2.950 Beiträge
 
#18

Re: Subtraktion mit DECMath bringt Probleme

  Alt 16. Mär 2006, 00:25
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
  Mit Zitat antworten Zitat
Antwort Antwort
Seite 2 von 2     12   


Forumregeln

Es ist dir nicht erlaubt, neue Themen zu verfassen.
Es ist dir nicht erlaubt, auf Beiträge zu antworten.
Es ist dir nicht erlaubt, Anhänge hochzuladen.
Es ist dir nicht erlaubt, deine Beiträge zu bearbeiten.

BB-Code ist an.
Smileys sind an.
[IMG] Code ist an.
HTML-Code ist aus.
Trackbacks are an
Pingbacks are an
Refbacks are aus

Gehe zu:

Impressum · AGB · Datenschutz · Nach oben
Alle Zeitangaben in WEZ +1. Es ist jetzt 11:12 Uhr.
Powered by vBulletin® Copyright ©2000 - 2024, Jelsoft Enterprises Ltd.
LinkBacks Enabled by vBSEO © 2011, Crawlability, Inc.
Delphi-PRAXiS (c) 2002 - 2023 by Daniel R. Wolf, 2024 by Thomas Breitkreuz