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