Autor |
Beitrag |
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
Verfasst: Mi 19.11.14 23:02
Hallo,
unter de.wikipedia.org/wik...C3%B6pfelalgorithmus wurde wieder einmal ein interessanter Beitrag veröffentlicht.
Natürlich muss ich es sofort testen, zuerst den Tröpfelalgorithmus für die Eulersche Zahl.
Es funktioniert auch wie gewünscht, nur ist es mir, wie immer , zu langsam. 10000 Ziffern erhalte ich in 310 ms, 20000 in 1,1 s.
Sieht jemand eine Möglichkeit, wie man vor allem die Schleife
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:
| for j:=1 to n do begin ueber:=0; for i:=m downto 1 do begin z:=10*c[i]+ueber; c[i]:=z mod (i+1); ueber:=z div (i+1); end; k:=k+chr(48+ueber); end; |
beschleunigen kann? i,j,m,n,z sind integer, c ein Array of shortint.
Wahrscheinlich wird's mit Assembler schneller, aber da muss ich passen.
Das Schöne an diesem Tröpfelalgorithmus ist, dass man nur ein Array benötigt. Die Eulersche Zahl e zu berechnen ist ja nicht weiter dramatisch. Interessant wird es für Pi. Der entsprechende Algorithmus wird bei Wikipedia ja auch genannt.
Danke für jeden Hinweis.
Mathematiker
Edit: Der Anhang enthält die geänderte Version mit ASM und besserer Stringverarbeitung.
Einloggen, um Attachments anzusehen!
_________________ Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein
Zuletzt bearbeitet von Mathematiker am Do 20.11.14 09:21, insgesamt 1-mal bearbeitet
|
|
jfheins
Beiträge: 918
Erhaltene Danke: 158
Win 10
VS 2013, VS2015
|
Verfasst: Mi 19.11.14 23:36
Eine Sache fällt mir auf: Du kannst
Delphi-Quelltext 1: 2:
| c[i]:=z mod (i+1); ueber:=z div (i+1); |
zu einer Assembler-Instruktion (divl) zusammenfassen. Ich tippe mal darauf, dass Delphi das nicht macht (bitte prüfen). Da die Division eine vergleichsweise komplexe Operation ist, lässt sich da vielleicht etwas herausholen.
Konnte man inline Assembler-Blöcke einbauen?
Und was ist k? Sag bitte nicht "Ein string, der jede Iteration um 1 Stelle verlängert wird"...
Ach, und ich erinnere mich an den Thread hier: www.entwickler-ecke....ewtopic.php?t=111154
Mitgenommene Lehren: Die innere Schleife (die multipliziert) lässt sich vielleicht parallelisieren, und nicht immer nur mit 10 sondern mit 100 oder so multiplizieren und damit direkt 2-4 Ziffern der Zahl auf einmal herausziehen.
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
Verfasst: Do 20.11.14 00:17
Hallo,
jfheins hat folgendes geschrieben : | Konnte man inline Assembler-Blöcke einbauen? |
inline Assembler-Blöcke gehen bei Delphi 5. Mal sehen, ob ich was "bauen" kann.
Ich habe:
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:
| procedure DivMod32(Dividend, Divisor: Cardinal; var Quotient, Remainder: Cardinal); asm PUSH EBX MOV EBX,EDX XOR EDX,EDX DIV EBX MOV [ECX],EAX MOV EBX,Remainder MOV [EBX],EDX POP EBX end; | gefunden. Damit wird es etwas schneller. Allerdings brauche ich jetzt ein Array of longword.
jfheins hat folgendes geschrieben : | Und was ist k? Sag bitte nicht "Ein string, der jede Iteration um 1 Stelle verlängert wird"... |
Doch, es ist ein String. Aber eine ziffernweise Ausgabe in die Memo oder das Speichern in einem weiteren Feld dauern noch länger. Und es muss je Schleifendurchlauf ausgewertet werden, da in der nächsten Runde die Ziffer wieder "verschwindet".
jfheins hat folgendes geschrieben : | und nicht immer nur mit 10 sondern mit 100 oder so multiplizieren und damit direkt 2-4 Ziffern der Zahl auf einmal herausziehen. |
Das ist eine Idee. Ob das bei diesem Algorithmus funktioniert, muss ich mir mal überlegen.
Beste Grüße
Mathematiker
Nachtrag: Mit
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9:
| function DivModI32(Dividend: Integer; Divisor: Integer; var Remainder: Integer): Integer; register; asm push ebx mov ebx, edx cdq idiv ebx mov [ecx], edx pop ebx end; |
wird es deutlich schneller. Für 10000 Ziffern 255 ms, für 20000 nun 920 ms.
Bleibt noch das Problem mit dem String.
Der Anhang im 1.Beitrag enthält den geänderten Quelltext.
_________________ Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein
Zuletzt bearbeitet von Mathematiker am Do 20.11.14 09:22, insgesamt 2-mal bearbeitet
|
|
C#
Beiträge: 561
Erhaltene Danke: 65
Windows 10, Kubuntu, Android
Visual Studio 2017, C#, C++/CLI, C++/CX, C++, F#, R, Python
|
Verfasst: Do 20.11.14 00:48
Ich sage nur Duff's Device
_________________ Der längste Typ-Name im .NET-Framework ist: ListViewVirtualItemsSelectionRangeChangedEventHandler
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
Verfasst: Do 20.11.14 00:58
Hallo,
Ich habe
Delphi-Quelltext
durch
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12:
| case ueber of 0 : k:=k+'0'; 1 : k:=k+'1'; 2 : k:=k+'2'; 3 : k:=k+'3'; 4 : k:=k+'4'; 5 : k:=k+'5'; 6 : k:=k+'6'; 7 : k:=k+'7'; 8 : k:=k+'8'; 9 : k:=k+'9'; end; |
ersetzt, was ein paar Millisekunden bringt.
Beste Grüße
Mathematiker
_________________ Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein
|
|
C#
Beiträge: 561
Erhaltene Danke: 65
Windows 10, Kubuntu, Android
Visual Studio 2017, C#, C++/CLI, C++/CX, C++, F#, R, Python
|
Verfasst: Do 20.11.14 01:58
Hallo
Also vorab: ich kann kein Delphi und bin auch kein Profi. Ich bin noch ein Erstsemestler, also bitte ich um Nachsicht wenn ich hier was falsch deute
Bei Duff's Device geht's doch letztendlich um das Ausrollen von Schleifen um weniger Iterationen zu verwenden. In deinem Fall müsste das dann so aussehen:
Delphi-Quelltext 1: 2: 3: 4: 5: 6:
| for i:=m downto 1 do begin z:=10*c[i]+ueber; c[i]:=z mod (i+1); ueber:=z div (i+1); end; |
Würde zu
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9:
| for i:=m downto 1 step 2 do begin z:=10*c[i]+ueber; c[i]:=z mod (i+1); ueber:=z div (i+1); z:=10*c[i-1]+ueber; c[i-1]:=z mod (i); ueber:=z div (i); end; |
Je weiter die Schleife ausgerollt wird, desto schneller sollte die Funktion theoretisch sein.
Die Modulo Rechnung ist ja nur da, um den Rest abzuarbeiten, falls der Ausrollfaktor (hier 2) teilerfremd zur Wertespanne (hier m) ist.
_________________ Der längste Typ-Name im .NET-Framework ist: ListViewVirtualItemsSelectionRangeChangedEventHandler
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Narses
Beiträge: 10182
Erhaltene Danke: 1255
W10ent
TP3 .. D7pro .. D10.2CE
|
Verfasst: Do 20.11.14 05:40
Moin!
Ich hab da jetzt nicht soo intensiv reingesehen, aber mal grundsätzlich was dazu:
Mathematiker hat folgendes geschrieben : | Ich habe
Delphi-Quelltext |
NIEMALS (in egal welcher Sprache) Stringoperationen dieser Art durchführen, wenn es auf Performance ankommt. Das ist zwar sehr "anschaulich" (füge ein Zeichen an), aber die Compilermagic kriegt bei sowas immer Probleme, denn Strukturen mit impliziter Größe (Strings verwaltet die Compiler-Magic ja intern) können einfach nicht "schnell" sein. Man stelle sich vor, was da passiert: - k sei zu Beginn der Leerstring
- k := k +Char(x); Speicherblock der Länge 1 allozieren (OK, vereinfacht, so kleine Blöcke verwaltet der Heap gar nicht, aber das ist konzeptionell egal) und das Zeichen reinschreiben
- k := k +Char(x); neuen Speicherblock der Länge 2 allozieren, den alten Block dahin kopieren und das neue Zeichen anfügen
- ...
- k := k +Char(x); neuen Speicherblock der Länge n+1 allozieren, den alten Block dahin kopieren und das neue Zeichen anfügen
Das kann weder schnell noch speicherplatzschonend ablaufen...
Lösung: "Raten", wie groß der Speicherblock vermutlich werden wird (irgendwie abschätzen halt ) und dann mit einem Zeiger (der gleichzeitig Position und Länge hergibt) reinschreiben. Das mit dem "Raten" ist zugegeben schwer, man muss im Zweifel den Speicherbedarf anschauen. Hat man sich verschätzt und muss korrigieren, dann wieder in größeren Blöcken "nachlegen".
cu
Narses
_________________ There are 10 types of people - those who understand binary and those who don´t.
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Horst_H
Beiträge: 1653
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Do 20.11.14 08:32
Hallo,
Ich habe das Programm laufen lassen ( es hat noch keiner heruntergeladen gehabt ??? )
Die letzten Ziffern sind nicht überzeugend...
Ab etwa n = 8645 enstehen falsche Zeichen am Ende.
Quelltext 1:
| 1867079953012004064493547882649367381033176214569394018075890.* |
EDIT:
Zurück marsch marsch....
In der zweiten Version funktioniert es...
Bitte alle Revisionen nach ganz oben.... gähn...
Gruß Horst
Schnelles Edit:
Delphi kennt DivMod schon ewig....
Quelltext 1: 2: 3: 4:
| DivMod procedure DivMod(Dividend: Integer; Divisor: Word; var Result, Remainder: Word);
Also see: MulDiv SysUtils (D3) |
Noch ein Edit:
Die Stringbehadlung kommt doch erst nach n Rechenschritte vor. Da spart man nichts...
HIer mal als
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54: 55: 56: 57: 58: 59: 60: 61: 62: 63: 64: 65:
| program utropf; {$IFDEF FPC} {$MODE DELPHI} {$ELSE} {$APPTYPE CONSOLE} {$ENDIF} uses sysutils; function DivModI32(Dividend: Integer; Divisor: Integer; var Remainder: Integer): Integer;assembler; asm push ebx mov ebx, edx cdq idiv ebx mov [ecx], edx pop ebx end; const digit : array[0..9] of char = ('0','1','2','3','4','5','6','7','8','9');
procedure Button1Click(n: integer); var m,i,j,kIdx : integer; st:double; ueber:integer; c:array of integer; k:string; Time1, Time2, Time3: TDateTime;
begin Time1:= time; st:=ln(2)/ln(10); i:=2; while st<n do begin inc(i); st:=st+ln(i)/ln(10); end; m:=i;
setlength(c,m+1); setlength(k,n+3); for i:=1 to m do c[i]:=1;
k[1]:='2'; k[2]:=','; kIdx := 3;
for j:=1 to n do begin ueber:=0; for i:=m downto 1 do ueber:=DivModI32(10*c[i]+ueber,i+1,c[i]); k[kidx] := digit[ueber]; inc(kidx); end; setlength(k,kidx-1);
Time2:= time; writeln(k); Time3:= time; setlength(c,0); writeln('Rechenzeit :',FormatDatetime('NN:SS.ZZZ',time2-time1)); writeln('Ausgabezeit :',FormatDatetime('NN:SS.ZZZ',time3-time2)); end; |
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Mathematiker
Beiträge: 2622
Erhaltene Danke: 1447
Win 7, 8.1, 10
Delphi 5, 7, 10.1
|
Verfasst: Do 20.11.14 09:19
Hallo
und Danke für die vielen Hinweise.
Narses hat folgendes geschrieben : | NIEMALS (in egal welcher Sprache) Stringoperationen dieser Art durchführen, wenn es auf Performance ankommt. ... "Raten", wie groß der Speicherblock vermutlich werden wird ... |
Verstehe ich. Das Raten ist auch kein Problem, da der String eben die gewünschte Anzahl Stellen hat. Werde ich ändern.
Horst_H hat folgendes geschrieben : | Die Stringbehadlung kommt doch erst nach n Rechenschritte vor. Da spart man nichts... |
Habe ich gemerkt, als ich den String testweise ganz herausgenommen habe.
Deinen Vorschlag zur Verarbeitung des Strings habe ich übernommen. Danke.
Gestern Abend ist mir aber noch klar geworden, dass nicht viel mehr herauszuholen zu ist, denn die eigentliche Tücke bei extrem vielen Stellen ist der dann wahrscheinlich auftretende arithmetische Überlauf der Zwischenwerte. Deshalb musste ich ja von shortint auf integer wechseln.
Im Moment sind es gute 18 Sekunden für 100000 Stellen. Das geht schon.
Der Algorithmus ist theoretisch interessant, für praktisches Verarbeiten aber weniger geeignet. Leider.
Beste Grüße
Mathematiker
_________________ Töten im Krieg ist nach meiner Auffassung um nichts besser als gewöhnlicher Mord. Albert Einstein
|
|
Horst_H
Beiträge: 1653
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Do 20.11.14 09:35
Hallo,
einen Vorteil hat er, er bleibt dezimal..
Gruß Horst
Edit:
ich habe mich gefragt:"was macht der denn da?"
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9:
| st:=ln(2); i:=2; dezis = ln(10)*n; while st<dezis do begin inc(i); st:=st+ln(i); end; m:=i/ln(10); |
das ist ln(n!)/ln(10). Es wird also die Falkultät gesucht die n Dezimalen entspricht.
Da gibt es doch sicher was für ln(n!):
de.wikipedia.org/wiki/Stirlingformel
Das lohnt nicht, per Intervallschachtelung i zu suchen.
Edit2:
Gibt es eine schnellere Annäherung für e.
Bisher habe ich nur
e =Summe(i= 0..unendlich) (2*i+2)/(2*i+1)!
Der Zähler wächst linear ~2i Der Nenner wird mit ~(2*i)³ größer 1!->4!->7!->10!->13! Also quadratische Genauigkeitszunahme.
Bei i= 5 also 2 Dezimalen pro Schritt, bei i= 50 4 Dezimalen bei 500 schon 6 Dezimalen. pro Schritt
|
|
Horst_H
Beiträge: 1653
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Sa 22.11.14 12:24
Hallo,
ein kleiner Nachtrag.
Natürlich kann man die Rechenzeit verringern, indem man nur die noch benötigten Stellen bearbeitet.
Wenn ich 9998 von 10000 Stellen schon ausgegeben habe, muss ich ja nicht mehr bis m = 450 Stellen dahinter rechnen, die ja ohnehin falsch werden.( Wie bei Die Kreiszahl Pi auch )
Damit halbiert sich in etwa die Laufzeit.
Ich benutze das schon genutzte st.
Delphi-Quelltext 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12:
| m0 := m; st := ln(m+5*ln(10)/ln(m)); ... For i := m0 downto 1 do ... st := st-ln(10); IF st < 0 then begin dec(m0); st := st+ln(m0); end; |
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54: 55: 56: 57: 58: 59: 60: 61: 62: 63: 64: 65: 66: 67: 68: 69: 70: 71: 72: 73: 74: 75: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105:
| program utropf; {$IFDEF FPC} {$MODE DELPHI} {$ENDIF} uses sysutils;
function OneLoop(c: pInteger; m0: Integer): Integer;assembler;
asm PUSH ESI; LEA ESI,[EAX+4*EDX] XOR EAX,EAX; TEST EDX,EDX JBE @Nothingtodo PUSH ECX; MOV ECX,EDX; INC ECX XOR EDX,EDX @Loop: MOV EAX,DWORD PTR[ESI] LEA EAX,[EAX+4*EAX] LEA EAX,[EDX+2*EAX] XOR EDX,EDX DIV ECX DEC ECX MOV DWORD PTR[ESI],EDX MOV EDX,EAX
SUB ESI,4 CMP ECX,2 JAE @Loop;
MOV EAX,EDX POP ECX;
@Nothingtodo: POP ESI;
end;
function DivModI32(Dividend: Integer; Divisor: Integer; var Remainder: Integer): Integer;assembler; asm push ebx mov ebx, edx xor edx,edx; idiv ebx mov [ecx], edx pop ebx end;
const digit : array[0..9] of char = ('0','1','2','3','4','5','6','7','8','9');
procedure Button1Click(n: integer); var m,i,j,m0 : integer; dezis,st:double; ueber:integer; c:array of integer; k:string; Time1, Time2, Time3: TDateTime;
begin Time1:= time; st:=ln(2); i:=2; dezis := n*ln(10); while st<dezis do begin inc(i); st:=st+ln(i); end; m:=i; writeln(' Stellen ',i); setlength(c,m+1); for i:=1 to m do c[i]:=1;
k:='2,'; m0 := m; st := ln(m+5*ln(10)/ln(m)); for j:= 1 to n do begin ueber:= OneLoop(@c[0],m0); k:= k+digit[ueber]; st := st-ln(10); IF st < 0 then begin dec(m0); st := st+ln(m0); end; end;
Time2:= time; writeln(k); Time3:= time; setlength(c,0); writeln('Rechenzeit :',FormatDatetime('NN:SS.ZZZ',time2-time1)); writeln('Ausgabezeit :',FormatDatetime('NN:SS.ZZZ',time3-time2)); end;
BEGIN Button1Click(100000); end. |
Gruß Horst
EDIT:
Die Hauptschleife ganz in Assembler.Das bringt keine 10 %.
ob die Lea schneller als MUL sind, weiß ich nicht, aber ich muss den Uebertrag der Division nicht in ein extra Register packen.
Zuletzt bearbeitet von Horst_H am Mo 24.11.14 21:19, insgesamt 1-mal bearbeitet
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Fiete
Beiträge: 611
Erhaltene Danke: 347
W7
Delphi 6 pro
|
Verfasst: Mo 24.11.14 18:58
Hier meine Version zur e - Berechnung.
Bestimmung von e mittels Reihenentwicklung
e=1/0!+1/1!+1/2!+1/3!+1/4!+1/5!+... oder
e=1/0!+1/1!+1/1!/2+1/2!/3+1/3!/4+1/4!/5+...
Mein Delphi 6 arbeitet bei Extended mit 80 Bit!
Viel Spaß beim Testen!
Gruß Fiete
Edit1: korrigierte Version - Berechnung von Digits
Einloggen, um Attachments anzusehen!
_________________ Fietes Gesetz: use your brain (THINK)
Zuletzt bearbeitet von Fiete am Mi 26.11.14 12:23, insgesamt 1-mal bearbeitet
Für diesen Beitrag haben gedankt: Mathematiker
|
|
Horst_H
Beiträge: 1653
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Mo 24.11.14 21:15
Hallo,
Hier kann man mit 1 Mio stellen von e vergleichen.
apod.nasa.gov/htmltest/gifcity/e.1mil
Was soll ich sagen, bei 1. Mio Stellen funktioniert es nicht richtig. Bei der 892407 Nachkommastelle stimmt es noch, dann nicht mehr.
Dummerweise ist der Tröpfelalgorithmus schneller und stimmt auch. ( 12 zu 11 min )
Gruß Horst
EDIT: wieder weg...
|
|
Fiete
Beiträge: 611
Erhaltene Danke: 347
W7
Delphi 6 pro
|
Verfasst: Di 25.11.14 16:45
Moin,
die Berechnung der Digits stimmte nicht ganz.
Delphi-Quelltext 1: 2: 3: 4: 5: 6:
| ZiffernAnzahl:=SpinEditN.Value; Ausgabe.Clear;Ausgabe.Repaint; Digits:=trunc(19-ln(ZiffernAnzahl)/ln(10)); if Digits>14 then Digits:=14; Basis:=1; for K:=1 to Digits do Basis:=Basis*10; |
Vorher stand hier:
Delphi-Quelltext 1: 2:
| Digits:=14; Basis:=100000000000000; |
Jetzt stimmen die Ziffern mit de.wikipedia.org/wik...C3%B6pfelalgorithmus überein.
Das Tröpfelprogramm ist bei mir langsamer
Gruß Fiete
_________________ Fietes Gesetz: use your brain (THINK)
|
|
Horst_H
Beiträge: 1653
Erhaltene Danke: 243
WIN10,PuppyLinux
FreePascal,Lazarus
|
Verfasst: Di 25.11.14 18:10
Hallo,
wer sucht denn auch bei der 892407.ten Nachkommastelle....
Das der Tröpfelalgo langsamer ist kann auch an der CPU liegen.Haswell kann nunmal sehr schnell dividieren. ( 4 Bit pro Takt, AMD 1 Bit/Takt ) .
Wahrscheinlicher: Wieder an der Memo.Ausgabe mit zig tausend Zeilen.
Ich versuche as jetzt mit GMP (oder mparith ab 65000 Zeichen geht es nicht oder nur über Umwege, müßte ich genauer suchen )
e = Summe i= 0..n (2*i+2)/(2*i+1)!
Das sah immer so gut aus, aber irgendwann stimmten die Stellen nicht mehr.
Von der Genauigkeit der Fliesskommazahl, die man setzen muss, wusste ich nicht.
Es klappt, ich bin erlöst....
1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38: 39: 40: 41: 42: 43: 44: 45: 46: 47: 48: 49: 50: 51: 52: 53: 54: 55: 56: 57: 58: 59: 60: 61: 62: 63: 64: 65: 66: 67: 68: 69: 70: 71: 72: 73: 74: 75: 76: 77: 78: 79: 80: 81: 82: 83: 84: 85: 86: 87: 88: 89: 90: 91: 92: 93: 94: 95: 96: 97: 98: 99: 100: 101: 102: 103: 104: 105: 106: 107: 108: 109: 110: 111: 112: 113: 114: 115: 116: 117: 118: 119: 120:
| program e_testgmp; {$IFdef fpc} {$MODE Delphi} {$ELSE} {$APPTYPE CONSOLE} {$ENDIF}
uses sysutils,classes, {$IFdef fpc} gmp; {$ELSE} gmp_lib; {$ENDIF}
const rowlen = 80; cDigits =1000*1000;
var Zaehler,Nenner,j: mpz_t; f,g : mpf_t;
buf : ansistring; teil : string[rowlen];
digCnt : double; n2,delta : integer;
procedure Vergleich; var se,so : TSTringlist; i: longint; begin se:= TStringlist.create; so:= TStringlist.create; se.loadfromfile('e_test.txt'); so.loadfromfile('out.txt');
For i := 1 to se.count do if se[i] <> so[i] then begin writeln('Abweichung in Zeile ',i+1); writeln(se[i]); writeln(so[i]); BREAK; end; so.free; se.free; end;
procedure Ausgabe(const f:mpf_t); var i,n : integer; txt : text; begin {$IFDEF FPC} mp_snprintf(PAnsiChar(buf),cDigits,'%.*Ff',cdigits,@f); {$else} gmp_printf(PAnsiChar(buf),'%.*Ff',cdigits,@f); {$endif} Assign(txt,'out.txt'); rewrite(txt); n := length(buf); i := 1; setlength(teil,rowlen); fillchar(teil[1],rowlen,' '); move(buf[1],teil[3],rowlen-2); writeln(txt,teil); inc(i,rowlen-2); dec(n,rowlen-2); repeat move(buf[i],teil[1],rowlen); writeln(txt,teil); inc(i,rowlen); dec(n,rowlen); until n < rowlen; Close(txt); end;
begin mpz_init(Zaehler); mpz_init(Nenner); mpz_init(j); setlength(buf,cDigits+10); try digCnt := 0.0; mpz_set_ui(Zaehler,0); mpz_set_ui(Nenner,1); n2 := 2; delta := 14; mpz_set_ui(j,6); repeat mpz_mul(Zaehler,j,Zaehler); mpz_add_ui(Zaehler,Zaehler,n2+2); mpz_mul(Nenner,j,Nenner); mpz_add_ui(j,j,delta); inc(delta,8); digCnt := digCnt+2*ln(n2); inc(n2,2); until digCnt >=ln(10)*cDigits; f_set_default_prec(Nenner.Size*32+32); mpf_init(f); mpf_init(g); mpf_set_z(f,Zaehler); mpf_set_z(g,Nenner); mpf_div(f,f,g); Ausgabe(f); finally mpf_clear(f); mpf_clear(g); mpz_clear(Zaehler); mpz_clear(Nenner); mpz_clear(j); end; Vergleich end. |
Quelltext 1: 2: 3: 4: 5:
| Abweichung in Zeile 12501 8188837471151 56239682713 8188837471151 3676005641232409486509687019237400185250412670296895354569172379931
real 0m14.526s |
Wie man ausrechnen kann 12500*80 = 1 Mio und schön schnell!
Gruß Horst
Für diesen Beitrag haben gedankt: MK2291
|
|
|