|
Registriert seit: 25. Jun 2003
Ort: Thüringen
2.950 Beiträge
|
Re: Unbegrenzt viele Nachkommastellen
16. Jan 2004, 14:47
Hi
Mehrere Optimierungsmöglichkeiten fallen mir da ein.
1.) es gibt bessere Typen des verwendeten Chudnovsky Algorithmus, sprich höherdimensionale Formeln. Statt 15 Stellen auf einmal zu berechnen können diese Chudnovsky Formeln z.B. 31 Stellen auf einmal berechnen. Allerdings sind diese Formeln viel komplizierter als die von mir benutzte
2.) die benutzte Formel basiert in weiten Teilen auf dem sogenannten "Binary Splitting". Man kann sich das so vorstellen das das Ziel eine Zahl X mit Y Stellen ist. Nun zerlegt man diese Y Stellen in zwei Teile einen Linken und Rechten Teil. Man versucht nun eine Formel zu finden die den Linken und Rechten Teil separat berechnen kann und dann nur noch den Linken und Rechten Teil miteinander multipliziert. Somit haben wir zwei kleine Berechnungen mit halb so großen Zahlen und deren Resultate werden multipliziert. Man kann nun Rekursiv wiederum den Linken und Rechten Teil teilen, usw. usw. Somit reduziert man die Zahlengrößen in den Berechnungen immer weiter. Wichtig ist es dabei das die beiden Zahlen der Linken und Rechten Seite der Berechnung möglichst gleichgroß sind. Denn dadurch wird der maximale Durchsatz, eg. Geschwindigkeit der Berechnungen erzielt. Es ist nämlich so das man zwei gleichgroße Zahlen am schnellsten Multiplizieren kann, d.h. zB. zwei Zahlen bei denen die eine Zahl 2 mal so lange wie die andere ist, wird immer langsammer pro Bit sein als wenn die beiden Zahlen etwa gleichgroß sind. Nun, exakt hier könnte man ebenfalls ansetzen, denn die von mir benutze Formel zur Pi Berechnung erzeugt NICHT immer gute Rechte + Linke Zahlen. Ziel der Optimierung ist es also die Berechnungen besser auszugleichen.
3.) Der von mir benutzte Binary Splitting Algorithmus ist ein Algo der auf Universalität codiert wurde. Im Falle der Pi Berechnung schleppt dieser Algo. z.B. eine Variable zusätzlich sinnlos mit. Das kostet Zeit. Ziel der Optimierung wäre es also einen eigenen Binary Splitting Algo. zu coden der besser auf die Anforderungen der Pi Formel angepasst ist. Unten habe ich mal den derzeitigen Bin. Splitting Algo. gepostet.
4.) Die Bin. Splitting Callback der Pi Berechnung berechnet viele succzessive Variablen. Im Grunde sind diese voneinander abhängig in ihren Zahlenwerten. Aber! meine Callback berechnet diese Werte immer absolut. Ziel der Optimierung wäre es nun diese Zahlenwerte ebenfalls inkrementell zu erzeugen, sprich die vielen kleinen Multiplikationen durch schnellere inkrementelle Additionen zu ersetzen. Dies ist durchaus möglich, man muß nur die besten Formeln dafür finden.
5.) meine Math. Library nutzt ganz verschiedene Multiplikations-verfahren abhänig von den Zahlengrößen. Für die längst-andauerenen Berechnungen = größten Zahlen wird der asymptotisch schnellste Algorithmus der heutzutage bekannt ist verwendet. Es ist die Modulare Fermat Fast Fourier Transformation von Schönhage und Strassen. Wichtig bei solchen FFT-Multiplikationen ist das Wörtchen "asympthotisch", denn es bedeutet das die Performance-Kurve im Verhältnis zu den Zahlengrößen eben NICHT linear sondern sprunghaft ist. D.h. die FFT berechnet immer in weiten Schranken unterschiedlich große Zahlen mit gleicher Geschwindigkeit. Erst wenn bestimmte Zahlkengrößen überschrittten werden steigt der Zeitaufwand der FFT sprunghaft an. Sozusagen entsteht eine Treppenfunktion in der Performance beim Ansteigen der Zahlengrößen. Nun, Ziel der Optimierung wäre es die Formeln so umzustellen das sie bei der Multiplizierung eben diese Treppenfunktion=asymptotisches Verhalten der FFT berücksichtigt. Kurz gesgt: es ist für die meisten FFT's am besten wenn die Zahlengröße in Bits ein Vielfaches von 2 ist.
6.) An vielen Stellen mit sehr wichtigen math. Funktionen habe ich in Assembler optimiert. Dabei wurden auch die neueren SSE2 Features der CPU's berücksichtigt. Beim SSE2 kann man sagen das dies ein mindestens 2 mal höhere Performance erzeugt. Allerdings es gibt noch genügend Source-Stückchen die noch nicht mit solchem SSE2 Code programmiert sind. Besodners die Funktionen für die FFT hätten das noch Potential. Ziel der Optimierung wäre es also diese wichtigen Sources in SSE2 Assembler zu coden.
7.) ein weiteres Beispiel für eine solche Optimierung ist der verwendete Division Algortihmus. Auch dieser basiert auf dem "Divide and Conquer" Verfahren. Der benutzte Algorithums ist der zur Zeit schnellste bekannte, die "Fast Recursive Karatsuba Division" von Burnikel & Ziegler. Auch dieser Algo. zerlegt das Gesamtproblem der Division in viel immer kleiner werdende Divisionen + Multiplikationen. Die kleinest Einheit von Divisions-Algo. ist aber eine normale Assemblercodierte Division nach Donald E. Knuth. Dieser Code ist aber in i486 Code programmiert. Ziel der Optimierung wäre es diesen Algo. auf SSE2 zu portieren. Dadurch düfte die Division im gesammten ca. 2 mal schneller werden. In einigen Performance-Vergleichen zu anderen Pi-Berechnugs-Programmen zeigte es sich das gerade meine Division um den Faktor 2 bis 3 langsammer ist als die der Konkurrenten. Dies liegt hauptsächlich daran das dieses Programme eine komplett andere Division von großen Zahlen benutzten. Die meisten nutzen eine Reziprokale Division die auf einer Newton Iteration basiert. Normalerweise ist dieser Algo. ineffizienter als der von mir benutzte, mit einer einzigsten Ausnahme. Eben der Divsion von sehr sehr großen Zahlen die die besonderen Eigenschaften besitzen wie sie bei der Pi-Formel von Natur aus entstehen. Somit wäre ein zweite Ziel der Oprimierung diese Reziprokale Division zu codieren.
8.) mein Speichermanagement für die IInteger ist darauf ausgelegt möglichst uuniversell flexibel zu sein. Dies äußerst sich auch darin das im Zahlenbereich bis zu ca. 500.000 Dezimalstellen von Pi mein jetziger Code eben schon ca. 2 mal schneller ist als die anderen Programme. D.h. mein Speichermanagement ist enorm effizient für kleinere Zahlen für kryptographische Algortihmen. Will man aber sehr große Zahlen berechnen, wie eben Pi, so ist es sinnvoller eine komplett andere Strategie im Speichermanagement zu verfolgen. Statt die viele kleinen Speicherzugriffe zu optimieren wird es bei großen Zahlen wichtig den Speicher möglichst linear wiederzuverwenden. Ziel der Optimierung wäre es also alle FFT Operationen innerhalb eines einmalig sehr groß allozierten Speicherbereiches inplaced durchzuführen. Dies bedingt aber eine komplette Neuprogrammierung der internen Zahlendarstellungen und Speicherzugriffe.
9.) der von mir benutzte Binary Splitting Algorithmus enthält eine partielle Formel-Optimierung. D.h. er optimiert bis zu einer Tiefe von 4, bestimmmte Multiplikationen weg. Theoretisch könnte man diese Optimierungstiefe aber viel stärker ausbauen. Dazu müssten alle Zwischenresulate der Callback bekannt sein. Nun könnte man das Bin. Splitting so optimieren das es mit der geringsten Anzahl an nötigen Multiplikationen arbeitet. Allein dies dürfte das Vrfahren mindestens 2 bis 3 mal beschleunigen. Allerdings steigt der Speicherverbrauch enorm an.
Alles in allem würde ich schätzen das alle obigen Optimierungen meine Pi Berechnungen ca. 3 bis 5 mal schneller machen könnten. Damit wäre diese in der Performance ca. 2 mal schneller als das schnellste Pi Programm das es zur Zeit für PC's gibt.
Allerdings, jede der obigen Optimierungen ist wesentlich komplizierter als es die jetzigen Algorithmen eh schon sind. D.h. es ist nicht so einfach wie es sich anhört. Alle bisherigen Algos. sind von mir schon bis zur max. Leistungsgrenze hin optimiert worden.
Gruß Hagen
Hier mal der Binary Spliting Code von mir:
Delphi-Quellcode:
type
TIIntegerSplitData = packed record
P,Q,A,B: IInteger;
end;
TIIntegerBinarySplittingCallback = procedure(N: Cardinal; var P: TIIntegerSplitData); register;
function NBinarySplitting( var P,Q: IInteger; Count: Integer;
Callback: TIIntegerBinarySplittingCallback; ImplicitShift: Boolean = True): Cardinal;
type
TSplitParam = packed record
P,Q,A,B: IInteger;
Callback: TIIntegerBinarySplittingCallback;
CallerEBP: Cardinal;
CalcP: Boolean;
end;
function BinarySplitting(n1,n2: Cardinal; var R: TSplitParam): Cardinal;
procedure Split_C(n: Cardinal; var R: TSplitParam);
asm
PUSH [EDX].TSplitParam.CallerEBP
CALL [EDX].TSplitParam.CallBack
POP EAX
end;
function Split_1(n: Cardinal; var R: TSplitParam): Cardinal;
begin
Split_C(n, R);
if R.P <> nil then
if R.A <> nil then NMul(R.A, R.P)
else R.A := R.P;
if R.Q <> nil then Result := NShr(R.Q)
else Result := 0;
end;
function Split_2(n: Cardinal; var R: TSplitParam): Cardinal;
var
L: TSplitParam;
rQs,lQs: Cardinal;
begin
rQs := 0;
lQs := 0;
L.Callback := R.Callback;
L.CallerEBP := R.CallerEBP;
L.CalcP := True;
Split_C(n +0, L);
Split_C(n +1, R);
// L.A := (L.A * R.B * R.Q * L.P) shl rQs;
// R.P = R.P * L.P;
// R.Q = R.Q * L.Q;
// R.B = R.B * L.B;
// R.A = R.A * L.B * R.P + L.A;
if R.Q <> nil then
begin
rQs := NShr(R.Q);
if L.A <> nil then NMul(L.A, R.Q)
else NSet(L.A, R.Q);
end else
if L.A = nil then NSet(L.A, 1);
if L.Q <> nil then
begin
lQs := NShr(L.Q);
if R.Q <> nil then NMul(R.Q, L.Q)
else NSwp(R.Q, L.Q);
end;
if L.B <> nil then
if R.A <> nil then NMul(R.A, L.B)
else R.A := L.B;
if R.B <> nil then
begin
NMul(L.A, R.B);
if L.B <> nil then NMul(R.B, L.B);
end else NSwp(R.B, L.B);
if L.P <> nil then
begin
NMul(L.A, L.P);
if R.P <> nil then NMul(R.P, L.P)
else NSwp(R.P, L.P);
end;
if R.P <> nil then
if R.A <> nil then NMul(R.A, R.P)
else R.A := R.P;
if rQs > 0 then NShl(L.A, rQs);
if R.A = nil then
begin
NSwp(R.A, L.A);
NInc(R.A);
end else NAdd(R.A, L.A);
Result := rQs + lQs;
end;
function Split_3(n: Cardinal; var R: TSplitParam): Cardinal;
var
L,M: TSplitParam;
rQs,lQs,mQs: Cardinal;
begin
lQs := 0;
mQs := 0;
rQs := 0;
L.Callback := R.Callback;
L.CallerEBP := R.CallerEBP;
M.Callback := R.Callback;
M.CallerEBP := R.CallerEBP;
Split_C(n +0, L);
Split_C(n +1, M);
Split_C(n +2, R);
if L.Q <> nil then lQs := NShr(L.Q);
if M.Q <> nil then mQs := NShr(M.Q);
if R.Q <> nil then rQs := NShr(R.Q);
// P := [R.P * [M.P * L.P]]
// Q := [R.Q * M.Q] * L.Q
// B := [R.B * M.B] * L.B
// R.A := ([R.B * M.B] * [R.Q * M.Q] * L.A * L.P) shl (rQs + mQs) +
// (([M.P * L.P] * R.B * R.Q * M.A) shl rQs) * L.B +
// [R.P * M.P * L.P] * M.B * R.A
// R.P := P
// R.Q := Q
// R.B := B
if L.P <> nil then
if M.P <> nil then NMul(M.P, L.P)
else M.P := L.P;
if M.P <> nil then
if R.P <> nil then NMul(R.P, M.P)
else R.P := M.P;
if M.A = nil then NSet(M.A, 1);
if M.P <> nil then NMul(M.A, M.P);
if R.B <> nil then NMul(M.A, R.B);
if R.Q <> nil then NMul(M.A, R.Q);
if rQs > 0 then NShl(M.A, rQs);
if R.A = nil then NSet(R.A, 1);
if R.P <> nil then NMul(R.A, R.P);
if M.B <> nil then NMul(R.A, M.B);
NAdd(R.A, M.A);
if L.B <> nil then NMul(R.A, L.B);
if M.Q <> nil then
if R.Q <> nil then NMul(R.Q, M.Q)
else R.Q := M.Q;
if M.B <> nil then
if R.B <> nil then NMul(R.B, M.B)
else R.B := M.B;
if L.A = nil then NSet(L.A, 1);
if L.P <> nil then NMul(L.A, L.P);
if R.B <> nil then NMul(L.A, R.B);
if R.Q <> nil then NMul(L.A, R.Q);
Inc(rQs, mQs);
if rQs > 0 then NShl(L.A, rQs);
NAdd(R.A, L.A);
if L.Q <> nil then
if R.Q <> nil then NMul(R.Q, L.Q)
else R.Q := L.Q;
if L.B <> nil then
if R.B <> nil then NMul(R.B, L.B)
else R.B := L.B;
Result := lQs + rQs;
end;
function Split_4(n: Cardinal; var R: TSplitParam): Cardinal;
var
L,M,K: TSplitParam;
rQs,mQs,kQs,lQs: Cardinal;
begin
lQs := 0;
mQs := 0;
kQs := 0;
rQs := 0;
L.Callback := R.Callback;
L.CallerEBP := R.CallerEBP;
M.Callback := R.Callback;
M.CallerEBP := R.CallerEBP;
K.Callback := R.Callback;
K.CallerEBP := R.CallerEBP;
Split_C(n +0, L);
Split_C(n +1, M);
Split_C(n +2, K);
Split_C(n +3, R);
if L.Q <> nil then lQs := NShr(L.Q);
if M.Q <> nil then mQs := NShr(M.Q);
if K.Q <> nil then kQs := NShr(K.Q);
if R.Q <> nil then rQs := NShr(R.Q);
if L.P <> nil then
if M.P <> nil then NMul(M.P, L.P)
else M.P := L.P;
if M.P <> nil then
if K.P <> nil then NMul(K.P, M.P)
else K.P := M.P;
if K.P <> nil then
if R.P <> nil then NMul(R.P, K.P)
else R.P := K.P;
if M.A = nil then NSet(M.A, 1);
if M.P <> nil then NMul(M.A, M.P);
if L.B <> nil then NMul(M.A, L.B);
if R.Q <> nil then
if K.Q <> nil then NMul(K.Q, R.Q)
else K.Q := R.Q;
if K.Q <> nil then NMul(M.A, K.Q);
if M.Q <> nil then NMul(K.Q, M.Q);
if L.A = nil then NSet(L.A, 1);
if L.P <> nil then NMul(L.A, L.P);
if K.Q <> nil then NMul(L.A, K.Q);
if M.B <> nil then NMul(L.A, M.B);
if mQs > 0 then NShl(L.A, mQs);
NAdd(L.A, M.A);
if R.B <> nil then
if K.B <> nil then NMul(K.B, R.B)
else K.B := R.B;
if K.B <> nil then NMul(L.A, K.B);
Inc(kQs, rQs);
if kQs > 0 then NShl(L.A, kQs);
if K.A = nil then NSet(K.A, 1);
if K.P <> nil then NMul(K.A, K.P);
if R.Q <> nil then NMul(K.A, R.Q);
if R.B <> nil then NMul(K.A, R.B);
if rQs > 0 then NShl(K.A, rQs);
if R.A = nil then NSet(R.A, 1);
if R.P <> nil then NMul(R.A, R.P);
if K.B <> nil then NMul(R.A, K.B);
NAdd(R.A, K.A);
if M.B <> nil then
if L.B <> nil then NMul(L.B, M.B)
else L.B := M.B;
if L.B <> nil then NMul(R.A, L.B);
NAdd(R.A, L.A);
NSwp(R.Q, L.Q);
if K.Q <> nil then
if R.Q <> nil then NMul(R.Q, K.Q)
else NSwp(R.Q, K.Q);
NSwp(R.B, L.B);
if K.B <> nil then
if R.B <> nil then NMul(R.B, K.B)
else NSwp(R.B, K.B);
Result := lQs + kQs + mQs;
end;
function Split_D(n1, n2: Cardinal; var R: TSplitParam): Cardinal;
var
L: TSplitParam;
nM,lQs,rQs: Cardinal;
begin
nM := (n1 + n2) shr 1;
L.Callback := R.Callback;
L.CallerEBP := R.CallerEBP;
L.CalcP := True;
lQs := BinarySplitting(n1, nM, L);
rQs := BinarySplitting(nM, n2, R);
// R.A := (R.A * L.B * L.P) + (L.A * R.B * R.Q) shl rQs
// R.P := R.P * L.P
// R.Q := R.Q * L.Q
// R.B := R.B * L.B
if R.B <> nil then
begin
if R.Q <> nil then
begin
if L.A <> nil then
begin
NMul(L.A, R.B);
NMul(L.A, R.Q);
if rQs > 0 then NShl(L.A, rQs);
end else
begin
NMul(L.A, R.B, R.Q);
if rQs > 0 then NShl(L.A, rQs);
end;
end else L.A := R.B;
end else
if R.Q <> nil then
if L.A <> nil then
begin
NMul(L.A, R.Q);
if rQs > 0 then NShl(L.A, rQs);
end else
if rQs > 0 then NShl(L.A, R.Q, rQs)
else L.A := R.Q;
if L.B <> nil then
begin
if R.A <> nil then
begin
NMul(R.A, L.B);
if L.P <> nil then NMul(R.A, L.P);
end else
if L.P <> nil then NMul(R.A, L.B, L.P);
end else
if L.P <> nil then
if R.A <> nil then NMul(R.A, L.P)
else NSet(R.A, L.P);
if R.A = nil then
begin
NSwp(R.A, L.A);
NInc(R.A);
end else
if L.A <> nil then NAdd(R.A, L.A)
else NInc(R.A);
if R.CalcP then
if L.P <> nil then
if R.P <> nil then NMul(R.P, L.P)
else NSwp(R.P, L.P);
if L.Q <> nil then
if R.Q <> nil then NMul(R.Q, L.Q)
else NSwp(R.Q, L.B);
if L.B <> nil then
if R.B <> nil then NMul(R.B, L.B)
else NSwp(R.B, L.B);
Result := lQs + rQs;
end;
resourcestring
sNBinarySplitting = ' NBinarySplitting(), internal error index = 0';
begin
case n2 - n1 of
0: begin
NRaise(@sNBinarySplitting);
Result := 0;
end;
1: Result := Split_1(n1, R);
2: Result := Split_2(n1, R);
3: Result := Split_3(n1, R);
4: Result := Split_4(n1, R);
else
begin
Result := Split_D(n1, n2, R);
end;
end;
end;
// apply a custom binary splitting computation for linear convergent series
//
// p[0]...p[n] a[n]
// S = P / Q = sum(n=0..Count -1) ----------- ----
// q[0]...q[n] b[n]
//
// TI97-7.ps, Fast multiprecision evaluation of series of rational numbers
// Bruno Haible & Thomas Papanikolaou
// for each coefficient p[n], q[n], a[n], b[n] are callback(N, P, Q, A, B) called.
// The custom callback procedure can access as nested procedure to the owner stack.
// If the IIntegers P,Q,A,B of the callback are unchanged = not assigned = NIL then
// there assumed to a value of +1. That reduce memory, interfaces and speedup.
// Means we can implicate that the callback is called by Callback(n, 1, 1, 1, 1);
resourcestring
sNBinarySplitting = ' NBinarySplitting(), Count must be >= 0 and Callback <> nil';
var
D: TSplitParam;
begin
if (Count < 0) or not Assigned(Callback) then NRaise(@sNBinarySplitting);
if Count = 0 then
begin
NSet(P, 0);
NSet(Q, 1);
Result := 0;
end else
begin
asm
MOV EAX,[EBP]
MOV D.CallerEBP,EAX
end;
D.Callback := Callback;
D.CalcP := False;
Result := BinarySplitting(0, Count, D);
if D.B <> nil then
if D.Q <> nil then NMul(D.Q, D.B)
else NSwp(D.Q, D.B);
if ImplicitShift then
begin
NShl(D.Q, Result);
Result := 0;
end;
NSwp(P, D.A);
NSwp(Q, D.Q);
end;
end;
|