home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / intsqrt.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  28.1 KB  |  621 lines
  1. # Wurzel aus ganzen Zahlen
  2.  
  3. # Zieht die Ganzzahl-Wurzel aus einer 32-Bit-Zahl und
  4. # liefert eine 16-Bit-Wurzel.
  5. # UL_sqrt_UW(x)
  6. # > uintL x : Radikand, >=0, <2^32
  7. # < uintWL ergebnis : Wurzel, >=0, <2^16
  8.   local uintWL UL_sqrt_UW (uintL x);
  9.   # Methode:
  10.   # x=0 -> y=0, fertig.
  11.   # y := 2^k als Anfangswert, wobei k>0, k<=16 mit 2^(2k-2) <= x < 2^(2k) sei.
  12.   # y := floor((y + floor(x/y))/2) als nächster Wert,
  13.   # solange z := floor(x/y) < y, setze y := floor((y+z)/2).
  14.   # y ist fertig.
  15.   # (Beweis:
  16.   #  1. Die Folge der y ist streng monoton fallend.
  17.   #  2. Stets gilt y >= floor(sqrt(x)) (denn für alle y>0 ist
  18.   #     y + x/y >= 2*sqrt(x) und daher  floor((y + floor(x/y))/2) =
  19.   #     floor(y/2 + x/(2*y)) >= floor(sqrt(x)) ).
  20.   #  3. Am Schluß gilt x >= y^2.
  21.   # )
  22.   local uintWL UL_sqrt_UW(x)
  23.     var reg3 uintL x;
  24.     { if (x==0) return 0; # x=0 -> y=0
  25.      { var reg6 uintC k2; integerlength32(x,k2=); # 2^(k2-1) <= x < 2^k2
  26.       {var reg5 uintC k1 = floor(k2-1,2); # k1 = k-1, k wie oben
  27.        if (k1 < 16-1)
  28.          # k < 16
  29.          { var reg1 uintWL y = (x >> (k1+2)) | bit(k1); # stets 2^(k-1) <= y < 2^k
  30.            loop
  31.              { var reg2 uintWL z;
  32.                divu_3216_1616(x,y, z=,_EMA_); # Dividiere x/y (geht, da x/y < 2^(2k)/2^(k-1) = 2^(k+1) <= 2^16)
  33.                if (z >= y) break;
  34.                y = floor(z+y,2); # geht, da z+y < 2*y < 2^(k+1) <= 2^16
  35.              }
  36.            return y;
  37.          }
  38.          else
  39.          # k = 16, Vorsicht!
  40.          { var reg4 uintWL x1 = high16(x);
  41.            var reg1 uintWL y = (x >> (16+1)) | bit(16-1); # stets 2^(k-1) <= y < 2^k
  42.            loop
  43.              { var reg2 uintWL z;
  44.                if (x1 >= y) break; # Division x/y ergäbe Überlauf -> z > y
  45.                divu_3216_1616(x,y, z=,_EMA_); # Dividiere x/y
  46.                if (z >= y) break;
  47.                y = floor(z+y,2);
  48.                if (intWLsize<=16) { y |= bit(16-1); } # y muß >= 2^(k-1) bleiben
  49.              }
  50.            return y;
  51.          }
  52.     }}}
  53.  
  54. #if 0 # unbenutzt
  55. # Zieht die Ganzzahl-Wurzel aus einer 64-Bit-Zahl und
  56. # liefert eine 32-Bit-Wurzel.
  57. # UL2_sqrt_UL(x1,x0)
  58. # > uintL2 x = x1*2^32+x0 : Radikand, >=0, <2^64
  59. # < uintL ergebnis : Wurzel, >=0, <2^32
  60.   local uintL UL2_sqrt_UL (uintL x1, uintL x0);
  61.   # Methode:
  62.   # x=0 -> y=0, fertig.
  63.   # y := 2^k als Anfangswert, wobei k>0, k<=32 mit 2^(2k-2) <= x < 2^(2k) sei.
  64.   # y := floor((y + floor(x/y))/2) als nächster Wert,
  65.   # solange z := floor(x/y) < y, setze y := floor((y+z)/2).
  66.   # y ist fertig.
  67.   # (Beweis:
  68.   #  1. Die Folge der y ist streng monoton fallend.
  69.   #  2. Stets gilt y >= floor(sqrt(x)) (denn für alle y>0 ist
  70.   #     y + x/y >= 2*sqrt(x) und daher  floor((y + floor(x/y))/2) =
  71.   #     floor(y/2 + x/(2*y)) >= floor(sqrt(x)) ).
  72.   #  3. Am Schluß gilt x >= y^2.
  73.   # )
  74.   local uintL UL2_sqrt_UL(x1,x0)
  75.     var reg3 uintL x1;
  76.     var reg4 uintL x0;
  77.     { if (x1==0) { return UL_sqrt_UW(x0); } # x klein?
  78.      { var reg6 uintC k2; integerlength32(x1,k2=); # 2^(k2+32-1) <= x < 2^(k2+32)
  79.       {var reg5 uintC k = ceiling(k2+32,2); # k wie oben
  80.        if (k < 32)
  81.          # k < 32
  82.          { var reg1 uintL y = ((x1 << (32-k)) | (x0 >> k) | bit(k)) >> 1; # stets 2^(k-1) <= y < 2^k
  83.            loop
  84.              { var reg2 uintL z;
  85.                divu_6432_3232(x1,x0,y, z=,_EMA_); # Dividiere x/y (geht, da x/y < 2^(2k)/2^(k-1) = 2^(k+1) <= 2^32)
  86.                if (z >= y) break;
  87.                y = floor(z+y,2); # geht, da z+y < 2*y < 2^(k+1) <= 2^32
  88.              }
  89.            return y;
  90.          }
  91.          else
  92.          # k = 32, Vorsicht!
  93.          { var reg1 uintL y = (x1 >> 1) | bit(32-1); # stets 2^(k-1) <= y < 2^k
  94.            loop
  95.              { var reg2 uintL z;
  96.                if (x1 >= y) break; # Division x/y ergäbe Überlauf -> z > y
  97.                divu_6432_3232(x1,x0,y, z=,_EMA_); # Dividiere x/y
  98.                if (z >= y) break;
  99.                y = floor(z+y,2) | bit(32-1); # y muß >= 2^(k-1) bleiben
  100.              }
  101.            return y;
  102.          }
  103.     }}}
  104. #endif
  105.  
  106. # Bildet zu einer Unsigned Digit sequence a die Wurzel
  107. # (genauer: Gaußklammer aus Wurzel aus a).
  108. # UDS_sqrt(a_MSDptr,a_len,a_LSDptr, &b, squarep=)
  109. # > a_MSDptr/a_len/a_LSDptr: eine UDS
  110. # < NUDS b: Gaußklammer der Wurzel aus a
  111. # < squarep: TRUE falls a = b^2, FALSE falls b^2 < a < (b+1)^2.
  112. # a wird nicht modifiziert.
  113. # Vorzeichenerweiterung von b ist erlaubt.
  114. # num_stack wird erniedrigt.
  115.   #define UDS_sqrt(a_MSDptr,a_len,a_LSDptr,b_,squarep_zuweisung)  \
  116.     { # ceiling(a_len,2) Digits Platz fürs Ergebnis machen:     \
  117.       var reg1 uintC _a_len = (a_len);                          \
  118.       num_stack_need_1(ceiling(_a_len,2),(b_)->MSDptr=,_EMA_);       \
  119.       squarep_zuweisung UDS_sqrt_(a_MSDptr,_a_len,a_LSDptr,b_); \
  120.     }
  121.   local boolean UDS_sqrt_ (uintD* a_MSDptr, uintC a_len, uintD* a_LSDptr, DS* b_);
  122. # Methode:
  123. # erst A normalisieren. A=0 --> B=0, fertig.
  124. # Wähle n so, daß beta^(2n-2) <= A < beta^(2n).
  125. # Wähle s (0<=s<16) so, daß beta^(2n)/4 <= A*2^(2s) < beta^(2n).
  126. # Setze A:=A*2^(2s) und kopiere dabei A. Suche B=floor(sqrt(A)).
  127. # Mache Platz für B=[0,b[n-1],...,b[0]], (mit einem Nulldigit Platz davor,
  128. # da dort nicht B, sondern 2*B abgespeichert werden wird).
  129. # Auf den Plätzen [a[2n-1],...,a[2n-2j]] wird die Differenz
  130. # [a[2n-1],...,a[2n-2j]] - [b[n-1],...,b[n-j]] ^ 2 abgespeichert.
  131. # Bestimme b[n-1] = floor(sqrt(a[2n-1]*beta+a[2n-2])) mit Heron/Newton:
  132. #   {x:=beta als vorheriger Anfangswert, dann:}
  133. #   x := floor((beta+a[2n-1])/2)
  134. #   wiederhole: d:=floor((a[2n-1]*beta+a[2n-2])/x).
  135. #               Falls d<beta (kein Überlauf) und d<x,
  136. #                 setze x:=floor((x+d)/2), nochmals.
  137. #   b[n-1]:=x. In B um ein Bit nach links verschoben abspeichern.
  138. # {Wegen a[2n-1]>=beta/4 ist b[n-1]>=beta/2.}
  139. # Erniedrige [a[2n-1],a[2n-2]] um b[n-1]^2.
  140. # Für j=1,...,n:
  141. #   {Hier [b[n-1],...,b[n-j]] = floor(sqrt(altes [a[2n-1],...,a[2n-2j]])),
  142. #     in [a[2n-1],...,a[2n-2j]] steht jetzt der Rest
  143. #     [a[2n-1],...,a[2n-2j]] - [b[n-1],...,b[n-j]]^2, er ist >=0 und
  144. #     und <= 2 * [b[n-1],...,b[n-j]], belegt daher höchstens j Digits und 1 Bit.
  145. #     Daher sind nur [a[2n-j],...,a[2n-2j]] von Belang.}
  146. #   Für j<n: Bestimme die nächste Ziffer:
  147. #     b* := min(beta-1,floor([a[2n-j],...,a[2n-2j-1]]/(2*[b[n-1],...,b[n-j]]))).
  148. #     und [a[2n-j],...,a[2n-2j-1]] :=
  149. #         [a[2n-j],...,a[2n-2j-1]] - b* * 2 * [b[n-1],...,b[n-j]] (>= 0).
  150. #     Im einzelnen:
  151. #       b* := min(beta-1,floor([a[2n-j],a[2n-j-1],a[2n-j-2]]/(2*b[n-1]))),
  152. #       [a[2n-j],...,a[2n-2j-1]] wie angegeben erniedigen.
  153. #       Solange die Differenz <0 ist, setze b* := b* - 1 und
  154. #         erhöhe [a[2n-j],...,a[2n-2j-1]] um 2 * [b[n-1],...,b[n-j]].
  155. #     Erniedrige [a[2n-j],...,a[2n-2j-2]] um b* ^ 2.
  156. #     Tritt dabei ein negativer Carry auf,
  157. #       so setze b* := b* - 1,
  158. #          setze b[n-j-1] := b* (im Speicher um 1 Bit nach links verschoben),
  159. #          erhöhe [a[2n-j],...,a[2n-2j-2]] um 2*[b[n-1],...,b[n-j-1]]+1.
  160. #       Sonst setze b[n-j-1] := b* (im Speicher um 1 Bit nach links verschoben).
  161. #     Nächstes j.
  162. #   Für j=n:
  163. #     Falls [a[n],...,a[0]] = [0,...,0], ist die Wurzel exakt, sonst nicht.
  164. #     Ergebnis ist [b[n-1],...,b[0]] * 2^(-s), schiebe also im Speicher
  165. #       [b[n],...,b[0]] um s+1 Bits nach rechts.
  166. #     Das Ergebnis ist eine NUDS der Länge n.
  167.   local boolean UDS_sqrt_(a_MSDptr,a_len,a_LSDptr,b_)
  168.     var reg10 uintD* a_MSDptr;
  169.     var reg10 uintC a_len;
  170.     var reg10 uintD* a_LSDptr;
  171.     var reg10 DS* b_;
  172.     { # A normalisieren:
  173.       while ((a_len>0) && (a_MSDptr[0]==0)) { a_MSDptr++; a_len--; }
  174.       if (a_len==0) # A=0 -> B := NUDS 0
  175.         { b_->LSDptr = b_->MSDptr; b_->len = 0; return TRUE; }
  176.      {SAVE_NUM_STACK # num_stack retten
  177.       # n und s bestimmen:
  178.       var reg10 uintC n = ceiling(a_len,2); # a_len = 2n oder 2n-1, n>0.
  179.       var reg10 uintL s;
  180.       { var reg1 uintD msd = a_MSDptr[0]; # a[2n] bzw. a[2n-1]
  181.         #if 0
  182.         s = 0;
  183.         while /* ((msd & (bit(intDsize-1)|bit(intDsize-2))) ==0) */
  184.               (((sintD)msd >= 0) && ((sintD)(msd<<1) >= 0))
  185.           { msd = msd<<2; s++; }
  186.         #else
  187.         integerlengthD(msd, s = intDsize - ); s = s>>1;
  188.         #endif
  189.       }
  190.       # Noch ist s nur modulo intDsize/2 bestimmt.
  191.       # A um 2s Bits nach links verschoben kopieren:
  192.       { var reg1 uintD* new_a_LSDptr;
  193.         num_stack_need(2*(uintL)n,a_MSDptr=,new_a_LSDptr=); # 2n Digits Platz belegen
  194.        {var reg2 uintL shiftcount = 2*s;
  195.         if (!((a_len & bit(0)) ==0)) # a_len ungerade?
  196.           { s += intDsize/2; *--new_a_LSDptr = 0; } # ja -> ein Nulldigit einschieben
  197.         if (shiftcount==0)
  198.           { copy_loop_down(a_LSDptr,new_a_LSDptr,a_len); }
  199.           else
  200.           { shiftleftcopy_loop_down(a_LSDptr,new_a_LSDptr,a_len,shiftcount); }
  201.       }}
  202.       # Nun ist A = a_MSDptr/2n/..
  203.       # Platz für B belegen:
  204.       { var reg10 uintD* b_MSDptr = &(b_->MSDptr)[-1]; # ab hier n+1 Digits Platz
  205.         var reg9 uintD b_msd;
  206.         # B = [0,b[n-1],...,b[0]] = b_MSDptr/n+1/..
  207.         # Bestimmung von b[n-1]:
  208.         { var reg3 uintD a_msd = a_MSDptr[0]; # a[2n-1]
  209.           var reg4 uintD a_2msd = a_MSDptr[1]; # a[2n-2]
  210.           #if HAVE_DD
  211.           var reg5 uintDD a_msdd = highlowDD(a_msd,a_2msd); # a[2n-1]*beta+a[2n-2]
  212.           #endif
  213.           # Anfangswert: x := floor((beta + a[2n-1])/2)
  214.           var reg1 uintD x = floor(a_msd,2) | bit(intDsize-1);
  215.           loop # Heron-Iterationsschleife
  216.             { var reg2 uintD d;
  217.               # Dividiere d := floor((a[2n-1]*beta+a[2n-2])/x) :
  218.               if (a_msd>=x) break; # Überlauf -> d>=beta -> fertig
  219.               #if HAVE_DD
  220.                 divuD(a_msdd,x, d=,_EMA_);
  221.               #else
  222.                 divuD(a_msd,a_2msd,x, d=,_EMA_);
  223.               #endif
  224.               if (d >= x) break; # d>=x -> fertig
  225.               # Nächste Iteration: x := floor((x+d)/2)
  226.               # (Da die Folge der x bekanntlich monoton fallend ist
  227.               # und bei b[n-1] >= beta/2 endet, muß x >= beta/2 werden,
  228.               # d.h. x+d>=beta.)
  229.               #if HAVE_DD
  230.                 x = (uintD)(floor((uintDD)x + (uintDD)d, 2));
  231.               #else
  232.                 x = floor((uintD)(x+d),2) | bit(intDsize-1);
  233.               #endif
  234.             }
  235.           # x = b[n-1] fertig berechnet.
  236.           b_msd = x;
  237.           # Quadrieren und von [a[2n-1],a[2n-2]] abziehen:
  238.           #if HAVE_DD
  239.             a_msdd -= muluD(x,x);
  240.             a_MSDptr[0] = highD(a_msdd); a_MSDptr[1] = lowD(a_msdd);
  241.           #else
  242.             {var reg5 uintD x2hi;
  243.              var reg2 uintD x2lo;
  244.              muluD(x,x, x2hi=,x2lo=);
  245.              a_MSDptr[0] = a_msd - x2hi;
  246.              if (a_2msd < x2lo) { a_MSDptr[0] -= 1; }
  247.              a_MSDptr[1] = a_2msd - x2lo;
  248.             }
  249.           #endif
  250.           b_MSDptr[0] = 1; b_MSDptr[1] = x<<1; # b[n-1] ablegen
  251.         }
  252.        {var reg8 uintC j = 0;
  253.         var reg3 uintD* a_mptr = &a_MSDptr[0];
  254.         var reg1 uintD* a_lptr = &a_MSDptr[2];
  255.         var reg2 uintD* b_ptr = &b_MSDptr[2];
  256.         # Wurzel-Hauptschleife
  257.         until (++j == n) # j=1,...,n
  258.           { # b_MSDptr = Pointer auf b[n], b_ptr = Pointer hinter b[n-j].
  259.             # a_mptr = Pointer auf a[2n-j], a_lptr = Pointer hinter a[2n-2j].
  260.             # Bestimme b* :
  261.             var reg4 uintD b_stern;
  262.             { var reg7 uintD a_1d = a_mptr[0]; # a[2n-j], =0 oder =1
  263.               var reg5 uintD a_2d = a_mptr[1]; # a[2n-j-1]
  264.               var reg6 uintD a_3d = a_mptr[2]; # a[2n-j-2]
  265.               # a[2n-j]*beta^2+a[2n-j-1]*beta+a[2n-j-2] durch 2 dividieren,
  266.               # dann durch b_msd = b[n-1] dividieren:
  267.               #if HAVE_DD
  268.                 var reg5 uintDD a_123dd = highlowDD(a_2d,a_3d);
  269.                 a_123dd = a_123dd>>1; if (!(a_1d==0)) { a_123dd |= bit(2*intDsize-1); }
  270.                 if (highD(a_123dd) >= b_msd)
  271.                   { b_stern = bitm(intDsize)-1; } # bei Überlauf: beta-1
  272.                   else
  273.                   { divuD(a_123dd,b_msd, b_stern=,_EMA_); }
  274.               #else
  275.                 a_3d = a_3d>>1; if (!((a_2d & bit(0)) ==0)) { a_3d |= bit(intDsize-1); }
  276.                 a_2d = a_2d>>1; if (!(a_1d==0)) { a_2d |= bit(intDsize-1); }
  277.                 if (a_2d >= b_msd)
  278.                   { b_stern = bitm(intDsize)-1; } # bei Überlauf: beta-1
  279.                   else
  280.                   { divuD(a_2d,a_3d,b_msd, b_stern=,_EMA_); }
  281.               #endif
  282.             }
  283.             # b_stern = b* in der ersten Schätzung.
  284.             a_lptr++; # Pointer hinter a[2n-2j-1]
  285.             # Subtraktion [a[2n-j],...,a[2n-2j-1]] -= b* * [b[n],b[n-1],...,b[n-j]] :
  286.             { var reg5 uintD carry = mulusub_loop_down(b_stern,b_ptr,a_lptr,j+1);
  287.               if (a_mptr[0] >= carry)
  288.                 { a_mptr[0] -= carry; }
  289.                 else
  290.                 { a_mptr[0] -= carry; # a[2n-j] wird <0
  291.                   # negativer Übertrag -> b* nach unten korrigieren:
  292.                   loop
  293.                     { b_stern = b_stern-1; # b* := b* - 1
  294.                       # erhöhe [a[2n-j],...,a[2n-2j-1]] um [b[n],...,b[n-j]]:
  295.                       if (!(( addto_loop_down(b_ptr,a_lptr,j+1) ==0)))
  296.                         if ((a_mptr[0] += 1) ==0) # Übertrag zu a[2n-j]
  297.                           break; # macht a[2n-j] wieder >=0 -> Subtraktionsergebnis >=0
  298.             }   }   }
  299.             # b_stern = b* in der zweiten Schätzung.
  300.             a_mptr++; # Pointer auf a[2n-j-1]
  301.             a_lptr++; # Pointer hinter a[2n-2j-2]
  302.             # Ziehe b* ^ 2 von [a[2n-j],...,a[2n-2j-2]] ab:
  303.             #if HAVE_DD
  304.             { var reg7 uintDD b_stern_2 = muluD(b_stern,b_stern);
  305.               var reg6 uintDD a_12dd = highlowDD(a_lptr[-2],a_lptr[-1]); # a[2n-2j-1]*beta+a[2n-2j-2]
  306.               var reg5 uintDD a_12dd_new = a_12dd - b_stern_2;
  307.               a_lptr[-2] = highD(a_12dd_new); a_lptr[-1] = lowD(a_12dd_new);
  308.               if (a_12dd >= b_stern_2) goto b_stern_ok;
  309.             }
  310.             #else
  311.             { var reg5 uintD b_stern_2_hi;
  312.               var reg6 uintD b_stern_2_lo;
  313.               muluD(b_stern,b_stern, b_stern_2_hi=,b_stern_2_lo=);
  314.              {var reg5 uintD a_1d = a_lptr[-2]; # a[2n-2j-1]
  315.               var reg6 uintD a_2d = a_lptr[-1]; # a[2n-2j-2]
  316.               var reg7 uintD a_1d_new = a_1d - b_stern_2_hi;
  317.               var reg7 uintD a_2d_new = a_2d - b_stern_2_lo;
  318.               if (a_2d < b_stern_2_lo) { a_1d_new -= 1; }
  319.               a_lptr[-2] = a_1d_new; a_lptr[-1] = a_2d_new;
  320.               if ((a_1d > b_stern_2_hi)
  321.                   || ((a_1d == b_stern_2_hi) && (a_2d >= b_stern_2_lo))
  322.                  )
  323.                 goto b_stern_ok;
  324.             }}
  325.             #endif
  326.             if (TRUE)
  327.               { # muß noch [a[2n-j],...,a[2n-2j]] um 1 erniedrigen:
  328.                 if ( dec_loop_down(&a_lptr[-2],j+1) ==0) goto b_stern_ok;
  329.                 # Subtraktion von b*^2 lieferte negativen Carry
  330.                 b_stern = b_stern-1; # b* := b* - 1
  331.                 # erhöhe [a[2n-j-1],...,a[2n-2j-2]] um [b[n],...,b[n-j],0] + 2 * b* + 1
  332.                 if ((sintD)b_stern < 0) { b_ptr[-1] |= bit(0); } # höchstes Bit von b* in b[n-j] ablegen
  333.                 b_ptr[0] = (uintD)(b_stern<<1)+1; # niedrige Bits von b* und eine 1 als b[n-j-1] ablegen
  334.                 addto_loop_down(&b_ptr[1],a_lptr,j+2);
  335.                 # (a[2n-j] wird nicht mehr gebraucht.)
  336.                 b_ptr[0] -= 1; # niedrige Bits von b* in b[n-j-1] ablegen
  337.                 b_ptr++;
  338.               }
  339.               else
  340.               b_stern_ok:
  341.               { # b* als b[n-j-1] ablegen:
  342.                 if ((sintD)b_stern < 0) { b_ptr[-1] |= bit(0); } # höchstes Bit von b* in b[n-j] ablegen
  343.                 b_ptr[0] = (uintD)(b_stern<<1); # niedrige Bits von b* als b[n-j-1] ablegen
  344.                 b_ptr++;
  345.               }
  346.           }
  347.         # b_MSDptr = Pointer auf b[n], b_ptr = Pointer hinter b[0].
  348.         # a_mptr = Pointer auf a[n].
  349.         # Schiebe [b[n],...,b[0]] um s+1 Bits nach rechts:
  350.         if (s == intDsize-1)
  351.           { b_ptr--; }
  352.           else
  353.           { shiftright_loop_up(b_MSDptr,n+1,s+1); b_MSDptr++; }
  354.         # b = b_MSDptr/n/b_ptr ist fertig, eine NUDS.
  355.         b_->MSDptr = b_MSDptr; b_->len = n; b_->LSDptr = b_ptr;
  356.         RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  357.         # Teste, ob alle a[n],...,a[0]=0 sind:
  358.         if (test_loop_up(a_mptr,n+1))
  359.           { return FALSE; }
  360.           else
  361.           { return TRUE; } # ja -> Wurzel exakt
  362.     }}}}
  363.  
  364. # Zieht die Wurzel (ISQRT x) aus einem Integer.
  365. # I_isqrt_I(x)
  366. # > x: Integer (sollte >=0 sein)
  367. # < STACK_0: (isqrt x)
  368. # < ergebnis: TRUE falls x Quadratzahl, FALSE sonst
  369. # erniedrigt STACK um 1
  370. # kann GC auslösen
  371.   local boolean I_isqrt_I (object x);
  372.   local boolean I_isqrt_I(x)
  373.     var reg1 object x;
  374.     { if (R_minusp(x))
  375.         { pushSTACK(x); # Wert für Slot DATUM von TYPE-ERROR
  376.           pushSTACK(O(type_posinteger)); # Wert für Slot EXPECTED-TYPE von TYPE-ERROR
  377.           pushSTACK(x);
  378.           pushSTACK(S(isqrt));
  379.           //: DEUTSCH "~ auf negative Zahl ~ angewandt"
  380.           //: ENGLISH "~ applied to negative number ~"
  381.           //: FRANCAIS "~ appliqué au nombre négatif ~"
  382.           fehler(type_error, GETTEXT("~ applied to negative number ~"));
  383.         }
  384.       { SAVE_NUM_STACK # num_stack retten
  385.         var reg2 uintD* x_MSDptr;
  386.         var reg4 uintC x_len;
  387.         var reg3 uintD* x_LSDptr;
  388.         I_to_NDS_nocopy(x, x_MSDptr=,x_len=,x_LSDptr=); # Digit sequence >=0 zu x
  389.        {var DS y;
  390.         var reg6 boolean squarep;
  391.         UDS_sqrt(x_MSDptr,x_len,x_LSDptr, &y, squarep=); # Wurzel ziehen
  392.         RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  393.         pushSTACK(NUDS_to_I(y.MSDptr,y.len)); # als Integer
  394.         return squarep;
  395.     } }}
  396.  
  397. # Stellt fest, ob ein Integer >=0 eine Quadratzahl ist.
  398. # I_sqrtp(x)
  399. # > x: ein Integer >=0
  400. # < ergebnis: Integer (sqrt x) falls x Quadratzahl, nullobj sonst
  401. # kann GC auslösen
  402.   local object I_sqrtp (object x);
  403. # Methode:
  404. # Damit x eine Quadratzahl ist, muß es ==0,1 mod 4 sein, und
  405. # bei ISQRT muß ein Rest 0 herauskommen.
  406.   local object I_sqrtp(x)
  407.     var reg1 object x;
  408.     { if (I_I_logbitp(Fixnum_1,x)) # Bit 1 von x gesetzt?
  409.         { return nullobj; } # ja -> x==2,3 mod 4, also kein Quadrat
  410.       { SAVE_NUM_STACK # num_stack retten
  411.         var reg2 uintD* x_MSDptr;
  412.         var reg4 uintC x_len;
  413.         var reg3 uintD* x_LSDptr;
  414.         I_to_NDS_nocopy(x, x_MSDptr=,x_len=,x_LSDptr=); # Digit sequence >=0 zu x
  415.        {var DS y;
  416.         var reg6 boolean squarep;
  417.         UDS_sqrt(x_MSDptr,x_len,x_LSDptr, &y, squarep=); # Wurzel ziehen
  418.         RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  419.         if (squarep)
  420.           { return NUDS_to_I(y.MSDptr,y.len); } # als Integer
  421.           else
  422.           { return nullobj; }
  423.     } }}
  424.  
  425. # Stellt fest, ob ein Integer >=0 eine n-te Potenz ist.
  426. # I_rootp(x,n)
  427. # > x: ein Integer >=0
  428. # > n: ein Integer >0
  429. # < ergebnis: Integer (expt x (/ n)) falls x eine n-te Potenz, nullobj sonst
  430. # kann GC auslösen
  431.   local object I_rootp (object x, object n);
  432. # Methode:
  433. # Falls x=0 oder x=1: x = x^n -> JA, x als Ergebnis.
  434. # Hier also x>1. Suche ein Integer y > 1 mit x=y^n.
  435. # Falls n >= integer_length(x): NEIN. (Da y>=2, müßte x>=2^n gelten.)
  436. # Hier also n>0 klein.
  437. # Solange n gerade ist: x := (sqrt x), n := (/ n 2). x nicht ganz -> NEIN.
  438. # Hier also n>0 ungerade.
  439. # Falls n=1: x = x^n -> JA, x als Ergebnis.
  440. # Falls o := ord2(x) nicht durch n teilbar ist: NEIN.
  441. # Sonst dividiere x durch 2^o, am Schluß y mit 2^(o/n) multiplizieren.
  442. # Hier also n>0 ungerade, x ungerade.
  443. # beta := 2^intDsize, m := ceiling(integer_length(x)/(intDsize*n)).
  444. # Suche ein y mit y>=0, y<beta^m mit  x == y^n mod beta^m :
  445. #   Mit dem Hensel-Lemma. Das Polynom f(X) = X^n-x hat die Diskriminante
  446. #   (-1)^((n-1)*n/2) * Res(X^n-x,n*X^(n-1)) = +- n^n * x^(n-1), und diese ist
  447. #   nicht durch p=2 teilbar. Daher ist das Hensel-Lemma mit p=2 anwendbar.
  448. #   Verwende quadratisches Hensel-Lifting, bei linearem Hensel-Lifting der
  449. #   der Verwaltungsaufwand vergleichsweise größer ist und die schnelle
  450. #   Multiplikation nicht zum Zuge kommt.
  451. #   Sei  y0 mod beta^k  mit  y0^n == x mod beta^k  bekannt. k=m -> fertig.
  452. #   Setze  y == y0 + beta^k*y1 mod beta^2k  an, wobei 2k := min(2*k,m).
  453. #   Dabei wird y1 mod beta^(2k-k) so gewählt, daß mod beta^2k
  454. #   x == y^n == y0^n + n * y0^(n-1) * beta^k*y1. Dazu wird
  455. #   (x - y0^n) mod beta^2k errechnet, durch beta^k dividiert (die Division
  456. #   muß nach Voraussetzung an y0 aufgehen) und
  457. #   y1 := ((x-y0^n)/beta^k) / (n*y0^(n-1)) mod beta^(2k-k) gebildet.
  458. #   Damit hat man  (y0 + beta^k*y1)^n == x mod beta^2k . 2k=m -> fertig.
  459. #   Den Anfang (k=1) bekommt man analog, mit beta:=2 und k=1,k=2,k=4,...
  460. # Dann testet man, ob wirklich x = y^n, und ist fertig.
  461.   local object I_rootp(x,n)
  462.     var reg1 object x;
  463.     var reg2 object n;
  464.     { if (eq(x,Fixnum_0) || eq(x,Fixnum_1)) # x=0 oder x=1 ?
  465.         { return x; } # ja -> x als Ergebnis
  466.       # Nun ist x>1.
  467.       pushSTACK(x); pushSTACK(n);
  468.       {var reg1 object l = I_integer_length_I(x); # (integer-length x)
  469.        if (I_I_comp(STACK_0,l) >= 0) # n >= (integer-length x) ?
  470.          { skipSTACK(2); return nullobj; }
  471.       }
  472.       # Nun ist n < (integer-length x). Also paßt n in ein uintL.
  473.      {var reg2 uintL n = I_to_UL(popSTACK());
  474.       var reg1 object x = popSTACK();
  475.       while ((n % 2) == 0) # n gerade?
  476.         { x = I_sqrtp(x); # Quadratwurzel ziehen versuchen
  477.           if (eq(x,nullobj)) { return nullobj; } # nicht ganzzahlig -> fertig
  478.           n = n >> 1; # n := (/ n 2)
  479.         }
  480.       # Nun ist n ungerade.
  481.       if (n==1) { return x; } # n=1 -> x als Ergebnis
  482.       {var reg10 uintL oq = 0; # Shift von y am Schluß
  483.        {var reg4 uintL o = I_ord2(x);
  484.         if (!(o==0))
  485.           {var reg3 uintL or; divu_3232_3232(o,n, oq=,or=); # or = o mod n
  486.            if (!(or==0)) { return nullobj; } # o nicht durch n teilbar -> fertig
  487.            # oq = o/n.
  488.            # dividiere x durch 2^o:
  489.            {var reg5 object temp;
  490.             pushSTACK(x); temp = L_to_I(-o); x = I_I_ash_I(popSTACK(),temp);
  491.        }  }}
  492.        # Nun ist n ungerade, x ungerade.
  493.        pushSTACK(x); # x retten
  494.        { SAVE_NUM_STACK # num_stack retten
  495.          var reg10 uintC n_len;
  496.          var reg10 uintD* n_LSDptr;
  497.          var reg10 uintC x_len;
  498.          var reg10 uintD* x_LSDptr;
  499.          # UDS zu n bilden, 0<n_len<=32/intDsize:
  500.          var uintD n_UDS[32/intDsize];
  501.          {var reg1 uintD* n_MSDptr = &n_UDS[0];
  502.           set_32_Dptr(n_MSDptr,n); n_LSDptr = &n_UDS[32/intDsize];
  503.           n_len = 32/intDsize; # und (zwecks Effizienz) normieren:
  504.           doconsttimes(32/intDsize-1,
  505.             { if (!(*n_MSDptr++ == 0)) goto n_UDS_ok; n_len--; }
  506.             );
  507.           n_UDS_ok: ; # n_MSDptr/n_len/n_LSDptr ist NUDS zu n.
  508.          }
  509.          I_to_NDS_nocopy(x,_EMA_,x_len=,x_LSDptr=); # UDS zu x bilden, x_len>0
  510.         {var reg3 uintD x_lsd = x_LSDptr[-1]; # letztes Digit von x
  511.          var reg2 uintD y_lsd; # n-te Wurzel von x_lsd mod 2^intDsize
  512.          y_lsd = 1; # Wurzel mod 2^1
  513.          # Für k=1,2,4,...:
  514.          # y_lsd := y_lsd + 2^k * (x_lsd-y_lsd^n)/2^k / (n*y_lsd^(n-1))
  515.          #        = y_lsd + (x_lsd-y_lsd^n) / (n*y_lsd^(n-1))
  516.          doconsttimes(log2_intDsize, # log2(intDsize) Iterationen reichen aus
  517.            { var reg1 uintD y_lsd_n1 = D_UL_expt_D(y_lsd,n-1); # y_lsd^(n-1)
  518.              var reg1 uintD y_lsd_n = D_D_mal2adic_D(y_lsd_n1,y_lsd); # y_lsd^n
  519.              var reg1 uintD delta = x_lsd-y_lsd_n; # x_lsd - y_lsd^n
  520.              if (delta==0) goto y_lsd_ok;
  521.              y_lsd = y_lsd + D_D_durch2adic_D(delta,D_D_mal2adic_D((uintD)n,y_lsd_n1));
  522.            });
  523.          y_lsd_ok:
  524.          ASSERT(D_UL_expt_D(y_lsd,n)==x_lsd);
  525.          # Nun ist y_lsd^n == x_lsd mod beta=2^intDsize.
  526.          { var reg9 uintL m = ceiling((uintL)x_len,n); # für y nötige Länge, >0, <=x_len
  527.            var reg6 uintD* y_LSDptr;
  528.            { var reg10 uintD* z1_LSDptr;
  529.              var reg10 uintD* z2_LSDptr;
  530.              var reg10 uintD* z3_LSDptr;
  531.              num_stack_need_1(m,_EMA_,y_LSDptr=); # Platz für y
  532.              {var reg1 uintL need = 2*m+(32/intDsize-1); # >= max(2*m,m+32/intDsize)
  533.               num_stack_need(need,_EMA_,z1_LSDptr=); # Platz für Rechenregister 1
  534.               num_stack_need(need,_EMA_,z2_LSDptr=); # Platz für Rechenregister 2
  535.               num_stack_need(need,_EMA_,z3_LSDptr=); # Platz für Rechenregister 3
  536.              }
  537.             {var reg8 uintL k = 1; # y ist bisher mod beta^k bekannt
  538.              y_LSDptr[-1] = y_lsd; # Startwert von y
  539.              until (k==m)
  540.                { var reg5 uintL k2 = 2*k; if (k2>m) { k2=m; } # k2 = min(2*k,m) > k
  541.                  # bisheriges y mod beta^k2 mit n-1 potenzieren:
  542.                  # Methode für z := y^(n-1) :
  543.                  #   zz:=y, e:=n-1.
  544.                  #   Solange e gerade, setze zz:=zz*zz, e:=e/2.
  545.                  #   z:=zz.
  546.                  #   Solange (e:=floor(e/2)) >0,
  547.                  #     setze zz:=zz*zz, und falls e ungerade, setze z:=z*zz.
  548.                 {var reg2 uintL e = n-1; # e:=n-1
  549.                  var reg7 uintD* free_LSDptr = z1_LSDptr;
  550.                  var reg3 uintD* zz_LSDptr = z2_LSDptr;
  551.                  var reg4 uintD* z_LSDptr;
  552.                  # Ab jetzt {zz_LSDptr,free_LSDptr} = {z1_LSDptr,z2_LSDptr}.
  553.                  clear_loop_down(&y_LSDptr[-(uintP)k],k2-k); # y auf k2 Digits erweitern
  554.                  copy_loop_down(y_LSDptr,zz_LSDptr,k2); # zz:=y
  555.                  do { var reg1 uintD* new_zz_LSDptr = free_LSDptr;
  556.                       mulu_2loop_down(zz_LSDptr,k2,zz_LSDptr,k2,new_zz_LSDptr); # zz:=zz*zz
  557.                       free_LSDptr = zz_LSDptr; zz_LSDptr = new_zz_LSDptr;
  558.                       e = e>>1; # e:=e/2
  559.                     }
  560.                     while ((e & bit(0)) ==0); # solange e gerade
  561.                  z_LSDptr = zz_LSDptr; # z:=zz
  562.                  # (zz nicht kopieren; ab der nächsten Veränderung von zz wird
  563.                  # {zz_LSDptr,z_LSDptr,free_LSDptr} = {z1_LSDptr,z2_LSDptr,z3_LSDptr}
  564.                  # gelten.)
  565.                  until ((e = e>>1) == 0)
  566.                    { {var reg1 uintD* new_zz_LSDptr = free_LSDptr;
  567.                       mulu_2loop_down(zz_LSDptr,k2,zz_LSDptr,k2,new_zz_LSDptr); # zz:=zz*zz
  568.                       free_LSDptr = (z_LSDptr==zz_LSDptr ? z3_LSDptr : zz_LSDptr);
  569.                       zz_LSDptr = new_zz_LSDptr;
  570.                      }
  571.                      if (!((e & bit(0)) == 0))
  572.                        {var reg1 uintD* new_z_LSDptr = free_LSDptr;
  573.                         mulu_2loop_down(z_LSDptr,k2,zz_LSDptr,k2,new_z_LSDptr); # z:=z*zz
  574.                         free_LSDptr = z_LSDptr; z_LSDptr = new_z_LSDptr;
  575.                    }   }
  576.                  # z = y^(n-1) mod beta^k2 ist fertig.
  577.                  if (z_LSDptr==zz_LSDptr) { zz_LSDptr = z3_LSDptr; } # zz ist jetzt auch frei
  578.                  mulu_2loop_down(z_LSDptr,k2,y_LSDptr,k2,free_LSDptr); # y^n
  579.                  sub_loop_down(x_LSDptr,free_LSDptr,zz_LSDptr,k2); # zz:=x-y^n
  580.                  ASSERT(!test_loop_up(&zz_LSDptr[-(uintP)k],k)); # zz == 0 mod beta^k
  581.                  mulu_2loop_down(z_LSDptr,k2-k,n_LSDptr,n_len,free_LSDptr); # n*y^(n-1)
  582.                  # Quotienten mod beta^(k2-k) bilden und an y mod beta^k ankleben:
  583.                  UDS_UDS_durch2adic_UDS(k2-k,&zz_LSDptr[-(uintP)k],free_LSDptr,&y_LSDptr[-(uintP)k]);
  584.                  k = k2; # jetzt gilt y^n == x sogar mod beta^k2.
  585.            }}  }}
  586.            # y mit y^n == x mod beta^m ist gefunden.
  587.            RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  588.            {var reg1 object y = UDS_to_I(&y_LSDptr[-(uintP)m],m); # y als Integer >=0
  589.             pushSTACK(y);
  590.        # Stackaufbau: x, y.
  591.        # y^n (mit n ungerade) bilden:
  592.        #   c:=a:=y, b:=n.
  593.        #   Solange b:=floor(b/2) >0 ist,
  594.        #     setze a:=a*a, und falls b ungerade, setze c:=a*c.
  595.        #   Liefere c.
  596.             pushSTACK(y); pushSTACK(y);
  597.        }}} } # Stackaufbau: x, y, c, a.
  598.        until ((n = n>>1) == 0)
  599.          { var reg1 object a = STACK_0; STACK_0 = a = I_I_mal_I(a,a);
  600.            if (!((n & bit(0)) == 0)) { STACK_1 = I_I_mal_I(a,STACK_1); }
  601.          }
  602.        {var reg1 object temp = STACK_1; skipSTACK(2); # temp = y^n
  603.         # mit x vergleichen:
  604.         if (!(I_I_comp(STACK_1,temp)==0))
  605.           # Die ganze Rechnung war umsonst.
  606.           { skipSTACK(2); return nullobj; }
  607.        }
  608.        # y ist tatsächlich n-te Wurzel von x.
  609.        # Noch mit 2^oq multiplizieren:
  610.        if (oq==0) # kein Shift nötig?
  611.          { var reg1 object temp = STACK_0;
  612.            skipSTACK(2); return temp;
  613.          }
  614.          else
  615.          { var reg1 object temp = UL_to_I(oq);
  616.            temp = I_I_ash_I(STACK_0,temp);
  617.            skipSTACK(2); return temp;
  618.          }
  619.     }}}
  620.  
  621.