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:

//gespeichert unter sqrt2gmp.dpr
program sqrt2gmp;
{$IFdef fpc}
  {$MODE Delphi}
{$ELSE}
  {$APPTYPE CONSOLE}
{$ENDIF}

// Wurzel 2 im Heronverfahren für cDigits Dezimalstellen
// Ausgabe mp_printf funktionierte unter linux64 mit 32bitlibs nur bis 8192 Zeichen
// Unter Windows kein Problem
uses
  sysutils,
{$IFdef fpc}
    gmp;
{$ELSE}
  gmp_lib;
{$ENDIF}

const
 cDigits =100000;
var
  T2,T1,T0: TDateTime;
  n,cnt : integer;
  a,Naeherung,erg,die2: mpq_t;
  f : mpf_t;
  buf : string;
begin
  mpq_init(die2);
  mpq_init(a);
  mpq_init(Naeherung);
  mpq_init(erg);
  T0:= Time;
try
  mpq_set_str(die2,'10',2); // 2 im Dualsystem
  mpq_set(a,die2);
  mpq_set(Naeherung,die2);
  n := round(ln(cDigits*(2-sqrt(2))/sqrt(2))/ln(2)+1)+1;
  cnt := n;
  repeat
//   Naeherung := ( a / Naeherung+Naeherung) / die2 ;
    mpq_div(erg,a,Naeherung);           //  a/x
    mpq_add(Naeherung,erg,Naeherung);   //  a/x+x
    //Statt Mul shiften um n=1 Bit, ob das wirklich heutzutage noch schneller ist?
    mpq_div_2exp(Naeherung,Naeherung,1);// x = (a/x+x)/2
    dec(cnt);
  until cnt<0;
  //Jetzt in Fliesskommazahl umwandeln
  mpf_init2(f,round(cDigits/0.30103));// cDigits dezimalen  = 100000/ log10(2) = 33219?.???

  mpf_set_q(f,erg);

  T1:= Time;
  setlength(buf,cDigits+7);
  {$IFDEF FPC}
    mp_snprintf(PAnsiChar(buf),cDigits,'%.*Ff',cdigits,@f);
  {$else}
    gmp_printf(PAnsiChar(buf),'%.*Ff',cdigits,@f);
//    gmp_printf(PAnsiChar(buf), '%.*Ff ',@f,cDigits);
  {$endif}
  T2:= Time;
  writeln;
  writeln(buf);
  writeln( n,' Durchlaeufe');
  Writeln('Rechenzeit  ',FormatDateTime('NN:SS.ZZZ' ,T1-T0));
  Writeln('Ausgabezeit ',FormatDateTime('NN:SS.ZZZ' ,T2-T1));
finally
  mpf_clear(f);
  mpq_clear(erg);
  mpq_clear(Naeherung);
  mpq_clear(a);
  mpq_clear(die2);
  readln
end;
end.                    

1:
2:
3:
4:

...4012197700603576101524970318470624937068864601669535559159840183770081805610147523
17 Durchlaeufe
Rechenzeit  00:00.493
Ausgabezeit 00:00.158                    

1:
2:
3:
4:
5:

..
  mpf_init2(f,round(cDigits/0.30103));
  mpfset_ui(f,2);
  mpf_sqrt(f,f);
..                    

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:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:

Procedure TWurzForm.StartClick(Sender: TObject);
Var I, N, P, W: LongInt;
Begin

  Start.Enabled := False;
  ZeitX := GetTickCount;
  ZeitR.Text := '0,000 sec.';
  WurzEdit.Text := '';
  StelleX.Text := '0';
  Go := True;

  SetCursor(LoadCursor(0, PChar(idc_Wait)));
  Y := 'Wurzel.Doc';
  AssignFile(Fil, Y);

  Z := 0;
  Zahl    := 0.0;
  Stellen := 1000.0;
  Summe1  := 0;  Summe2   := 0;
  Q_Test  := 0;  Q_Wurzel := 0.0;

  For I := -1 To 2000000 Do Begin
    A[I] := 0;
  End; {* For *}

  For I:=-2 To 4012000 Do Begin
    S[I] := 0;
    K[I] := 0;
  End; {* For *}

  Val(Wurzel.Text, Zahl, R);
  Val(Stelle.Text, Stellen, R);

  N := Trunc(Stellen);

  If (Zahl <= 0.0) Or (Zahl >= 1000) Then Begin
   MessageDlg('Ungültige Zahl', mtError, [mbOk], 0);
   Stop.Enabled := False;
   Clear.Enabled := True;
   Exit;
  End;

  If (Stellen < 1) Or (Stellen > 1000000) Then Begin
   MessageDlg('Ungültige Anzahl Stellen', mtError, [mbOk], 0);
   Stop.Enabled := False;
   Clear.Enabled := True;
   Exit;
  End;

  P := Trunc(Sqrt(Zahl));

  K[-2] := Trunc((P * P)/100);
  K[-1] := Trunc((P * P)/10) - 10 * Trunc((P * P)/100);
  K[0]  := P * P - 10 * Trunc((P * P)/10);
  A[-1] := Trunc(P/10);
  A[0]  := P - 10 * Trunc(P/10);
  S[-2] := K[-2];
  S[-1] := K[-1];
  S[0]  := K[0];

  Q_Test   := 0;
  Q_Wurzel := 0.0;
  Q_Test   := Round(10000.0 * Zahl);

  If (Sqrt(Q_Test) = Trunc(Sqrt(Q_Test))) Then Begin
     Q_Wurzel := Sqrt(Q_Test)/100;
     Str(Q_Wurzel:2:4, StrZ);
     WurzEdit.Lines[0] := 'Unsinnige Vorgabe, QWurzel = ' +
     StrZ;
     Stop.Enabled := False;
     Clear.Enabled := True;
     Exit;
  End; {* IF *}

  I := 1;  {* Wichtig: I Ist nach den obigen Schleifen nicht definiert *}


  {* Zentrale Schleife, die die Nachkommastellen mit Vektoroperationen *}
  {* berechnet. Jede Komponente ist eine Ziffer                        *}

  While (I <= N) And Go Do Begin

     Application.ProcessMessages;
     If (I Mod 100 = 0) Then Begin
       ZeitY := GetTickCount - ZeitX;
       ZeitReal := ZeitY / 1000.0;
       ZeitR.Text := FloatToStrF(ZeitReal, ffFixed, 6, 3) + ' sec.';
       Str(I, AnzStel);
       StelleX.Text := AnzStel;
     End;

     Stop.Enabled := True;
     Clear.Enabled := False;
     {* Ausmultiplizieren in Vektorschreibweise *}

     For W := -1 To I Do Begin
         Inc(K[W + I], 2 * A[W]);
     End;  {* For *}

     K[2 * I] := K[2 * I] + 1;


     {* Stellenübertrag *}

     For W := 2 * I DownTo -1 Do Begin
         Z := 0;
         While (K[W] > 9) Do Begin
           Dec(K[W], 10);
           Inc(Z);
         End;
         Inc(K[W - 1], Z);
     End;  {* For *}


     {* 10'000 fache Eingabe und 10'000 fache Multiplikationssumme *}

     Summe1 := Round(1.0E06 * K[-2] + 1.0E05 * K[-1] + 1.0E04 * K[0] +
               1000.0 * K[1]);
     Summe1 := Summe1 + Round(100.0 * K[2] + 10.0 * K[3] + 1.0 * K[4]);
     Summe2 := Round(10000.0 * ZAHL);


     {* Grösse der Eingabe erreicht / Oder letzte Ziffer = 9 *}

     If (Summe1 >= Summe2) or (A[I] = 9) Then Begin
       Inc(I);
       CopyMemory(@K[-2], @S[-2], 2 * I + 3);
     End {* If (Summe1 ... Then *}

     {* Sonst: Letzte Ziffer um 1 erhöhen *}

     Else Begin
       Inc(A[I]);
       CopyMemory(@S[-2], @K[-2], 2 * I + 3);
     End;  {* IF (Summe1 ... Else *}

  End;  {* While I <= N *}

  If (Not Go) Then N := I;
  ZeitY := GetTickCount - ZeitX;
  ZeitReal := ZeitY / 1000.0;

  If (Q_Wurzel = 0.0) Then Begin
   ReWrite(Fil);
   Writeln(Fil,'Protokoll:');
   Writeln(Fil,'__________');
   Writeln(Fil,'          ');
   Writeln(Fil,'Quadratwurzel von ',Zahl:3:4,' auf ',N,' Stellen ist: ');
   Write(Fil, A[-1], A[0], '.');

   For I := 1 To N Do Begin
       If (I/72) = Trunc(I/72) Then Writeln(Fil, ' ');
       Write(Fil, A[I]);
   End; {* For *}

   Writeln(Fil, ' '); Writeln(Fil, ' ');
   Writeln(Fil, 'Ende des Protokolles.');
   CloseFile(Fil);
   ZeitR.Text := FloatToStrF(ZeitReal, ffFixed, 6, 3) + ' sec.';
   WurzEdit.Lines.LoadFromFile(Y);

   Stop.Enabled := False;
   Clear.Enabled := True;

   SetCursor(LoadCursor(0, PChar(idc_Arrow)));
  End;
End;