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:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:
221:
222:
223:
224:
225:
226:
227:

  function TVarInt.Sqroot: boolean;
  var qRest: TVarInt;
  begin
    qRest:= TVarInt.Create;
    try
      self.Sqroot(qRest);
      result:=(qRest.RLen = 0);
    finally
      qRest.Free;
    end;
  end;

  function TVarInt.Sqroot(var qREM: TVarInt; mitRest: boolean=false): boolean;
  var
    xNorm, xBit, x: integer;
    q4, {a3,} aL, xrem, w: int64;
    a3Ext, q4Ext: Extended;
    j, hQ: integer;
    i, n, m, q: integer;
    qminREM: TDig_V;
    aC, qDivisor, spSelf: TVarInt;
  //-----------
      procedure product(var x, y: TDig_V; f, a: integer); //x = y*f
      var i,j: integer;
        temp, carry: int64;
      begin
        carry:=0; j:=a + a;
        for i:=a to a + m -1 do begin
          temp:=int64(y[i])*f + carry;
          carry:=temp shr LOGBASE_V;
          x[j]:=temp and MODMASK_V; //temp - (carry shl LOGBASE_V); //
          j:=j+1;
        end;
         x[j]:=carry;
         x[j+1]:=0;
      end;
  //-----------
      procedure differenz(var x, y: TDig_V; k, a: integer); //x = x-y
      var i: integer;                                       //x > y
      begin
        for i:=k to k + a do begin
          if(x[i] >= y[i])then
            x[i]:=x[i] - y[i]
          else begin
            x[i]:=x[i] or BASE_V - y[i];
            y[i+1]:=y[i+1] + 1;
          end;
        end;
      end;
  //-----------
      procedure add(var x: TDig_V; k, a: integer); //x = x+k
      begin
        x[a]:=x[a] + k;
        if(x[a] and BASE_V <> 0)then begin
          while(x[a] and BASE_V <> 0)do begin
            x[a]:=x[a] - BASE_V; //and MODMASK_V; //
            a:=a + 1;
            x[a]:=x[a] + 1;
          end;
        end;
      end;
  //-----------
      function Compare(var i1, i2: TDig_V; a, b:integer): boolean;
      var i: integer;                                       //true, wenn i1 > i2
      begin
        result:=false;
        for i:=a + b downto a do begin
          if i1[i] < i2[i] then
            exit
          else
            if i1[i] > i2[i] then begin
              result:=true; exit;
            end;
        end;
      end;
  //------------
    procedure anpassenW;
    begin
      qDivisor.RNum[hQ shr 1]:=w;
      product(qminREM, qDivisor.RNum, w, hQ shr 1);//qminREM = qDivisor.RNum * w
    end;
  //------------
  begin
    if not(RSig)then begin    //wie Wurzelziehen von Hand
     // raise ERangeError.Create(' Radikant darf nicht negativ sein! (TVarInt.Sqroot)');
     // Erg.Fehler(111); 12060478383186831971436118892604006569990026212581001918540970352126506291069610672437
      result:=false;
    end
    else begin
      IF(self.isInt64)THEN BEGIN //if(RLen < 3)
        self.ConvertToInt64(aL);
        w:=trunc(Sqrt(aL/1));
        xrem:=aL - int64(w*w);   //qREM
        self.ToVar(w);
        qREM.ToVar(xrem);
        result:=(xrem = 0);
      END
      ELSE BEGIN
        a3Ext:=0; // a3:=0;
       // uTest:=0;
        xNorm:=0;
        qDivisor:=TVarInt.Create;
        qREM.Assign(self);  //direkt Weitergeben an qREM

        if(qREM.RLen mod 2 = 1)then begin //wenn self.RLen ungerade
          qREM.LeftShiftbyBASE(1);
          xNorm:=xNorm + 15;             //qREM.RLen immer gerade machen
        end;

        xBit:=GetHighestBitIndex30(qREM.RNum[qREM.RLen-1]);  //highest digit
        case xBit of
           0..25: begin
                    x:=(27 - xBit) shr 1;
                    xNorm:=xNorm + x;
                    qREM.LeftShift(x shl 1);
                  end;
          26, 27: ;
          else
            xNorm:=xNorm + 29;     //ab (2^28) wird qDivisor > BASE_V
            qREM.LeftShift(58);
            //qREM.RightShift(2);
            //xNorm:=xNorm -1;   //schwieriger zu realisieren
        end;
  //jetzt ist highest digit 2^26 .. (2^28)-1  und qREM.RLen gerade
                    //67 108 864 .. 268 435 455  optimal; wenige Korrekturen

        hQ:=qREM.RLen -1; //high(FNumber);                  //[hQ --> ungerade]
        aL:=int64(qREM.RNum[hQ]);    //die 2 ersten self-Ziffern convertieren !!
        aL:=aL shl LOGBASE_V + int64(qREM.RNum[hQ-1]);

        w:=trunc(Sqrt(aL/1));     //1.Ergebnis-FNum der Länge 1
        qDivisor.ToVar(w + w);//doppeltes, 1.Ergebnis muss ebenso 1 lang sein

        xrem:=aL - int64(w*w);      //1.qREM
//in sehr seltenen Fällen wird xrem < 0, was eigentlich nicht sein dürfte.
//Ursache ist eine fehlerhafte Berechnung von trunc(Sqrt(aL/1)): trunced- Fehler
//bei z.B. sqrt(288230376151711743) = 536870911,999999999068677
//                                --> 536870912, anstatt 536870911
//dieses Phänomen tritt nur bei bestimmten Aktionen auf wie z.B. Memo1.MouseDown
if(xrem < 0)then begin
 // showMessage('Rundungs-FEHLER bei Sqrt-Berechnung wird behoben: bitte OK drücken.');
  w:=w -1;
  qDivisor.RNum[0]:=w + w;
  qDivisor.RNum[1]:=0;
  xrem:=aL - int64(w*w);
end;

        aC:=TVarInt.Create(xrem);  //aC --> qREM
        SetLength(qminREM, hQ+1);


//Aufbereitung von qREM:
        for j:=hQ-1 to hQ do            //die 2 ersten Ziffern in qREM -->
          qREM.RNum[j]:=aC.RNum[j+1-hQ];       //--> ersetzen von 2 Ziffern aC

        q:=(hQ + 1)shr 1;
        qDivisor.LenTest(q);
        qDivisor.LeftShiftbyBASE(q -1);

spSelf:= TVarInt.Create(self);
        RLen:=q;
        RNum[q-1]:=w; // shr LOGBASE; //div BASE;  //1.Ergebnis-Ziffer
        TRY
          dec(hQ, 2);
          n:=hQ + 1;  //Length(self.FNum); //n immer gerade
          m:=1; //Length(qDivisor.FNum);  //m anfangs immer Länge 1, wegen *****
          WHILE(hQ > 0)DO BEGIN   //Berechnung
            n:=n - 2;
            m:=m + 1;
            if(m < 4)then begin   //nur die ersten 2 Durchgänge
                      //a3:=qDivisor.RNum[q-1];
              a3Ext:=qDivisor.RNum[q-1] + qDivisor.RNum[q-2]/BASE_V
            end;
            q4:=0;
            for i:=n + m downto n + m -1 do
              q4:=q4 shl LOGBASE_V + qREM.RNum[i];

            q4Ext:=q4 + qREM.RNum[n+m-2]/BASE_V;
                      //w:=q4 div a3; //Probedivision
            w:=trunc(q4Ext/a3Ext);
            if(w = 0)then begin
              RNum[hQ shr 1]:=0;
              hQ:=hQ - 2;
              continue;
            end;
            if(w >= BASE_V)then     //+1198898461484465545599392243
              w:=high(TBase_V);
            anpassenW;
            while(Compare(qminREM, qREM.RNum, n, m))do begin //Korrektur nötig
              w:= w -1; //inc(uTest);
              anpassenW;
            end;
            differenz(qREM.RNum, qminREM, n, m);
            RNum[hQ shr 1]:=w;              //weitere Ergebnis-RNums
            hQ:=hQ - 2;
            add(qDivisor.RNum, w, q - m);  //oder doppeltes, bisheriges Ergebnis
          END;

          if(xNorm > 0)then
            self.RightShift(xNorm);

          qREM.Trim;
          result:=(qREM.RLen = 0);  //max 2*Ergebnis

          if(not result)and(mitRest)then //qRest braucht nur dann berechnet zu werden,
            if(xNorm > 0)then begin //...wenn er tatsächlich verlangt wird(mitRest=true)
              xNorm:=xNorm shl 1;
              if(qREM.isSelf_mod_2k_ZERO(xNorm))then begin  //qREM durch xNorm teilbar ?
                qREM.RightShift(xNorm);
              end
              else begin
                aC.Assign(self);
                aC.Square;
                qREM.Sub(spSelf, aC);
              end;
            end;

        FINALLY
          aC.Free;
          qDivisor.Free;
         // if(assigned(spSelf))then
            spSelf.Free;
//hier tritt der Fehler auf !!!!!!!!!!!!!!!!!
        END;
      END;
    end;
  end;