src/integer/algebraic/cl_I_rootp_aux.cc

   1 // cl_rootp_aux().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/integer.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_I.h"
  13 #include "cl_DS.h"
  14 #include "cl_2D.h"
  15 #include "cl_2DS.h"
  16
  17 namespace cln {
  18
  19 // Stellt fest, ob ein Integer >=0 eine n-te Potenz ist.
  20 // rootp(x,n,&w)
  21 // > x: ein Integer >=0
  22 // > n: ein Integer >0
  23 // > Annahme: x > 1 und n < (integer-length x).
  24 // < w: Integer (expt x (/ n)) falls x eine n-te Potenz
  25 // < ergebnis: cl_true         ........................, cl_false sonst
  26
  27 cl_boolean cl_rootp_aux (cl_I x, uintL n, cl_I* w)
  28 {
  29
  30 // Methode:
  31 // Falls x=0 oder x=1: x = x^n -> JA, x als Ergebnis.
  32 // Hier also x>1. Suche ein Integer y > 1 mit x=y^n.
  33 // Falls n >= integer_length(x): NEIN. (Da y>=2, müßte x>=2^n gelten.)
  34 // Hier also n>0 klein.
  35 // Solange n gerade ist: x := (sqrt x), n := (/ n 2). x nicht ganz -> NEIN.
  36 // Hier also n>0 ungerade.
  37 // Falls n=1: x = x^n -> JA, x als Ergebnis.
  38 // Falls o := ord2(x) nicht durch n teilbar ist: NEIN.
  39 // Sonst dividiere x durch 2^o, am Schluß y mit 2^(o/n) multiplizieren.
  40 // Hier also n>0 ungerade, x ungerade.
  41 // beta := 2^intDsize, m := ceiling(integer_length(x)/(intDsize*n)).
  42 // Suche ein y mit y>=0, y<beta^m mit  x == y^n mod beta^m :
  43 //   Mit dem Hensel-Lemma. Das Polynom f(X) = X^n-x hat die Diskriminante
  44 //   (-1)^((n-1)*n/2) * Res(X^n-x,n*X^(n-1)) = +- n^n * x^(n-1), und diese ist
  45 //   nicht durch p=2 teilbar. Daher ist das Hensel-Lemma mit p=2 anwendbar.
  46 //   Verwende quadratisches Hensel-Lifting, bei linearem Hensel-Lifting der
  47 //   der Verwaltungsaufwand vergleichsweise größer ist und die schnelle
  48 //   Multiplikation nicht zum Zuge kommt.
  49 //   Sei  y0 mod beta^k  mit  y0^n == x mod beta^k  bekannt. k=m -> fertig.
  50 //   Setze  y == y0 + beta^k*y1 mod beta^2k  an, wobei 2k := min(2*k,m).
  51 //   Dabei wird y1 mod beta^(2k-k) so gewählt, daß mod beta^2k
  52 //   x == y^n == y0^n + n * y0^(n-1) * beta^k*y1. Dazu wird
  53 //   (x - y0^n) mod beta^2k errechnet, durch beta^k dividiert (die Division
  54 //   muß nach Voraussetzung an y0 aufgehen) und
  55 //   y1 := ((x-y0^n)/beta^k) / (n*y0^(n-1)) mod beta^(2k-k) gebildet.
  56 //   Damit hat man  (y0 + beta^k*y1)^n == x mod beta^2k . 2k=m -> fertig.
  57 //   Den Anfang (k=1) bekommt man analog, mit beta:=2 und k=1,k=2,k=4,...
  58 // Dann testet man, ob wirklich x = y^n, und ist fertig.
  59
  60       while ((n % 2) == 0) // n gerade?
  61         { if (!sqrtp(x,&x)) // Quadratwurzel ziehen versuchen
  62             { return cl_false; } // nicht ganzzahlig -> fertig
  63           n = n >> 1; // n := (/ n 2)
  64         }
  65       // Nun ist n ungerade.
  66       if (n==1) { *w = x; return cl_true; } // n=1 -> x als Ergebnis
  67       var uintL oq = 0; // Shift von y am Schluß
  68       {var uintL o = ord2(x);
  69        if (!(o==0))
  70          {var uintL o_r; divu_3232_3232(o,n, oq=,o_r=); // o_r = o mod n
  71           if (!(o_r==0)) { return cl_false; } // o nicht durch n teilbar -> fertig
  72           // oq = o/n.
  73           // dividiere x durch 2^o:
  74           x = ash(x,-(sintL)o);
  75       }  }
  76       // Nun ist n ungerade, x ungerade.
  77       CL_ALLOCA_STACK;
  78       var uintC n_len;
  79       var uintD* n_LSDptr;
  80       var uintC x_len;
  81       var const uintD* x_LSDptr;
  82       // UDS zu n bilden, 0<n_len<=ceiling(32/intDsize):
  83       var uintD n_UDS[ceiling(32,intDsize)];
  84       #if (intDsize==64)
  85       arrayLSref(n_UDS,1,0) = n; n_LSDptr = arrayLSDptr(n_UDS,1); n_len = 1;
  86       #else // (intDsize<=32)
  87       {var uintD* n_MSDptr = arrayMSDptr(n_UDS,32/intDsize);
  88        set_32_Dptr(n_MSDptr,n); n_LSDptr = arrayLSDptr(n_UDS,32/intDsize);
  89        n_len = 32/intDsize; // und (zwecks Effizienz) normieren:
  90        doconsttimes(32/intDsize-1,
  91          { if (!(msprefnext(n_MSDptr) == 0)) goto n_UDS_ok; n_len--; }
  92          );
  93        n_UDS_ok: ; // n_MSDptr/n_len/n_LSDptr ist NUDS zu n.
  94       }
  95       #endif
  96       I_to_NDS_nocopy(x, ,x_len=,x_LSDptr=,cl_false,); // UDS zu x bilden, x_len>0
  97       var uintD x_lsd = lspref(x_LSDptr,0); // letztes Digit von x
  98       var uintD y_lsd; // n-te Wurzel von x_lsd mod 2^intDsize
  99       y_lsd = 1; // Wurzel mod 2^1
 100       // Für k=1,2,4,...:
 101       // y_lsd := y_lsd + 2^k * (x_lsd-y_lsd^n)/2^k / (n*y_lsd^(n-1))
 102       //        = y_lsd + (x_lsd-y_lsd^n) / (n*y_lsd^(n-1))
 103       doconsttimes(log2_intDsize, // log2(intDsize) Iterationen reichen aus
 104         { var uintD y_lsd_n1 = expt_pos(y_lsd,n-1); // y_lsd^(n-1)
 105           var uintD y_lsd_n = mul2adic(y_lsd_n1,y_lsd); // y_lsd^n
 106           var uintD delta = x_lsd-y_lsd_n; // x_lsd - y_lsd^n
 107           if (delta==0) goto y_lsd_ok;
 108           y_lsd = y_lsd + div2adic(delta,mul2adic((uintD)n,y_lsd_n1));
 109         });
 110       y_lsd_ok:
 111       ASSERT(expt_pos(y_lsd,n)==x_lsd);
 112       // Nun ist y_lsd^n == x_lsd mod beta=2^intDsize.
 113       { var uintL m = ceiling((uintL)x_len,n); // für y nötige Länge, >0, <=x_len
 114         var uintD* y_LSDptr;
 115         { var uintD* z1_LSDptr;
 116           var uintD* z2_LSDptr;
 117           var uintD* z3_LSDptr;
 118           num_stack_alloc_1(m, ,y_LSDptr=); // Platz für y
 119           {var uintL need = 2*m+(32/intDsize-1); // >= max(2*m,m+32/intDsize)
 120            num_stack_alloc(need, ,z1_LSDptr=); // Platz für Rechenregister 1
 121            num_stack_alloc(need, ,z2_LSDptr=); // Platz für Rechenregister 2
 122            num_stack_alloc(need, ,z3_LSDptr=); // Platz für Rechenregister 3
 123           }
 124          {var uintL k = 1; // y ist bisher mod beta^k bekannt
 125           lspref(y_LSDptr,0) = y_lsd; // Startwert von y
 126           until (k==m)
 127             { var uintL k2 = 2*k; if (k2>m) { k2=m; } // k2 = min(2*k,m) > k
 128               // bisheriges y mod beta^k2 mit n-1 potenzieren:
 129               // Methode für z := y^(n-1) :
 130               //   zz:=y, e:=n-1.
 131               //   Solange e gerade, setze zz:=zz*zz, e:=e/2.
 132               //   z:=zz.
 133               //   Solange (e:=floor(e/2)) >0,
 134               //     setze zz:=zz*zz, und falls e ungerade, setze z:=z*zz.
 135               var uintL e = n-1; // e:=n-1
 136               var uintD* free_LSDptr = z1_LSDptr;
 137               var uintD* zz_LSDptr = z2_LSDptr;
 138               var uintD* z_LSDptr;
 139               // Ab jetzt {zz_LSDptr,free_LSDptr} = {z1_LSDptr,z2_LSDptr}.
 140               clear_loop_lsp(y_LSDptr lspop k,k2-k); // y auf k2 Digits erweitern
 141               copy_loop_lsp(y_LSDptr,zz_LSDptr,k2); // zz:=y
 142               do { var uintD* new_zz_LSDptr = free_LSDptr;
 143                    cl_UDS_mul(zz_LSDptr,k2,zz_LSDptr,k2,new_zz_LSDptr); // zz:=zz*zz
 144                    free_LSDptr = zz_LSDptr; zz_LSDptr = new_zz_LSDptr;
 145                    e = e>>1; // e:=e/2
 146                  }
 147                  while ((e & bit(0)) ==0); // solange e gerade
 148               z_LSDptr = zz_LSDptr; // z:=zz
 149               // (zz nicht kopieren; ab der nächsten Veränderung von zz wird
 150               // {zz_LSDptr,z_LSDptr,free_LSDptr} = {z1_LSDptr,z2_LSDptr,z3_LSDptr}
 151               // gelten.)
 152               until ((e = e>>1) == 0)
 153                 { {var uintD* new_zz_LSDptr = free_LSDptr;
 154                    cl_UDS_mul(zz_LSDptr,k2,zz_LSDptr,k2,new_zz_LSDptr); // zz:=zz*zz
 155                    free_LSDptr = (z_LSDptr==zz_LSDptr ? z3_LSDptr : zz_LSDptr);
 156                    zz_LSDptr = new_zz_LSDptr;
 157                   }
 158                   if (!((e & bit(0)) == 0))
 159                     {var uintD* new_z_LSDptr = free_LSDptr;
 160                      cl_UDS_mul(z_LSDptr,k2,zz_LSDptr,k2,new_z_LSDptr); // z:=z*zz
 161                      free_LSDptr = z_LSDptr; z_LSDptr = new_z_LSDptr;
 162                 }   }
 163               // z = y^(n-1) mod beta^k2 ist fertig.
 164               if (z_LSDptr==zz_LSDptr) { zz_LSDptr = z3_LSDptr; } // zz ist jetzt auch frei
 165               cl_UDS_mul(z_LSDptr,k2,y_LSDptr,k2,free_LSDptr); // y^n
 166               sub_loop_lsp(x_LSDptr,free_LSDptr,zz_LSDptr,k2); // zz:=x-y^n
 167               ASSERT(!DS_test_loop(zz_LSDptr lspop k,k,zz_LSDptr)); // zz == 0 mod beta^k
 168               cl_UDS_mul(z_LSDptr,k2-k,n_LSDptr,n_len,free_LSDptr); // n*y^(n-1)
 169               // Quotienten mod beta^(k2-k) bilden und an y mod beta^k ankleben:
 170               div2adic(k2-k,zz_LSDptr lspop k,free_LSDptr,y_LSDptr lspop k);
 171               k = k2; // jetzt gilt y^n == x sogar mod beta^k2.
 172         }}  }
 173         // y mit y^n == x mod beta^m ist gefunden.
 174         var cl_I y = UDS_to_I(y_LSDptr lspop m,m); // y als Integer >=0
 175         // y^n (mit n ungerade) bilden:
 176         //   c:=a:=y, b:=n.
 177         //   Solange b:=floor(b/2) >0 ist,
 178         //     setze a:=a*a, und falls b ungerade, setze c:=a*c.
 179         //   Liefere c.
 180         { var cl_I c = y;
 181           var cl_I a = y;
 182           until ((n = n>>1) == 0)
 183             { a = square(a); if (!((n & bit(0)) == 0)) { c = a * c; } }
 184           // c = y^n
 185           // mit x vergleichen:
 186           if (!(x == c))
 187             // Die ganze Rechnung war umsonst.
 188             { return cl_false; }
 189         }
 190         // y ist tatsächlich n-te Wurzel von x.
 191         // Noch mit 2^oq multiplizieren:
 192         if (oq==0) // kein Shift nötig?
 193           { *w = y; }
 194           else
 195           { *w = ash(y,oq); }
 196         return cl_true;
 197       }
 198 }
 199
 200 }  // namespace cln