Branch data Line data Source code
1 : : /* Copyright (C) 2000 The PARI group.
2 : :
3 : : This file is part of the PARI/GP package.
4 : :
5 : : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : : terms of the GNU General Public License as published by the Free Software
7 : : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
8 : : ANY WARRANTY WHATSOEVER.
9 : :
10 : : Check the License for details. You should have received a copy of it, along
11 : : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 : :
14 : : /*********************************************************************/
15 : : /** **/
16 : : /** ARITHMETIC FUNCTIONS **/
17 : : /** (first part) **/
18 : : /** **/
19 : : /*********************************************************************/
20 : : #include "pari.h"
21 : : #include "paripriv.h"
22 : :
23 : : /******************************************************************/
24 : : /* */
25 : : /* GENERATOR of (Z/mZ)* */
26 : : /* */
27 : : /******************************************************************/
28 : : static GEN
29 [ + - ]: 35 : remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
30 : : static ulong
31 : 20494 : u_remove2(ulong q) { return q >> vals(q); }
32 : : GEN
33 : 35 : odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
34 : : static GEN
35 : 20494 : u_odd_prime_divisors(ulong q) { return gel(factoru(u_remove2(q)), 1); }
36 : : /* p odd prime, q=(p-1)/2; L0 list of (some) divisors of q = (p-1)/2 or NULL
37 : : * (all prime divisors of q); return the q/l, l in L0 */
38 : : static GEN
39 : 33 : is_gener_expo(GEN p, GEN L0)
40 : : {
41 : 33 : GEN L, q = shifti(p,-1);
42 : : long i, l;
43 [ + + ]: 33 : if (L0) {
44 : 14 : l = lg(L0);
45 : 14 : L = cgetg(l, t_VEC);
46 : : } else {
47 : 19 : L0 = L = odd_prime_divisors(q);
48 : 19 : l = lg(L);
49 : : }
50 [ + + ]: 102 : for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
51 : 33 : return L;
52 : : }
53 : : static GEN
54 : 21789 : u_is_gener_expo(ulong p, GEN L0)
55 : : {
56 : 21789 : const ulong q = p >> 1;
57 : : long i, l;
58 : : GEN L;
59 [ + + ]: 21789 : if (L0) {
60 : 2163 : l = lg(L0);
61 : 2163 : L = cgetg(l, t_VECSMALL);
62 : : } else {
63 : 19626 : L0 = L = u_odd_prime_divisors(q);
64 : 19626 : l = lg(L);
65 : : }
66 [ + + ]: 61048 : for (i=1; i<l; i++) L[i] = q / uel(L0,i);
67 : 21789 : return L;
68 : : }
69 : :
70 : : int
71 : 53086 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
72 : : {
73 : : long i;
74 [ + + ]: 53086 : if (krouu(x, p) >= 0) return 0;
75 [ + + ]: 68972 : for (i=lg(L)-1; i; i--)
76 : : {
77 : 45841 : ulong t = Fl_powu(x, uel(L,i), p);
78 [ + + ][ + + ]: 45841 : if (t == p_1 || t == 1) return 0;
79 : : }
80 : 53086 : return 1;
81 : : }
82 : : /* assume p prime */
83 : : ulong
84 : 50810 : pgener_Fl_local(ulong p, GEN L0)
85 : : {
86 : 50810 : const pari_sp av = avma;
87 : 50810 : const ulong p_1 = p-1;
88 : : long x;
89 : : GEN L;
90 [ + + ]: 50810 : if (p <= 19) switch(p)
[ + + + ]
91 : : { /* quick trivial cases */
92 : 42 : case 2: return 1;
93 : : case 7:
94 : 3537 : case 17: return 3;
95 : 25466 : default: return 2;
96 : : }
97 : 21765 : L = u_is_gener_expo(p,L0);
98 [ + + ]: 50993 : for (x=2;;x++) { if (is_gener_Fl(x,p,p_1,L)) break; }
99 : 50810 : avma = av; return x;
100 : : }
101 : : ulong
102 : 40780 : pgener_Fl(ulong p) { return pgener_Fl_local(p, NULL); }
103 : :
104 : : /* L[i] = set of (p-1)/2l, l ODD prime divisor of p-1 (l=2 can be included,
105 : : * but wasteful) */
106 : : int
107 : 127 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
108 : : {
109 [ + + ]: 127 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
110 [ + + ]: 127 : if (t >= 0) return 0;
111 [ + + ]: 283 : for (i = lg(L)-1; i; i--)
112 : : {
113 : 199 : GEN t = Fp_pow(x, gel(L,i), p);
114 [ + + ][ + + ]: 199 : if (equalii(t, p_1) || equali1(t)) return 0;
115 : : }
116 : 127 : return 1;
117 : : }
118 : :
119 : : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
120 : : GEN
121 : 10122 : pgener_Fp_local(GEN p, GEN L0)
122 : : {
123 : 10122 : pari_sp av0 = avma;
124 : : GEN x, p_1, L;
125 [ + + ]: 10122 : if (lgefint(p) == 3)
126 : : {
127 : : ulong z;
128 [ + + ]: 10093 : if (p[2] == 2) return gen_1;
129 [ + + ]: 7678 : if (L0) L0 = ZV_to_nv(L0);
130 : 7678 : z = pgener_Fl_local(uel(p,2), L0);
131 : 7678 : avma = av0; return utoipos(z);
132 : : }
133 : 29 : p_1 = subis(p,1); L = is_gener_expo(p, L0);
134 : 29 : x = utoipos(2);
135 [ + + ]: 32 : for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
136 : 10122 : avma = av0; return utoipos(uel(x,2));
137 : : }
138 : :
139 : : GEN
140 : 9618 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
141 : :
142 : : ulong
143 : 52 : pgener_Zl(ulong p)
144 : : {
145 [ - + ]: 52 : if (p == 2) pari_err_DOMAIN("pgener_Zl","p","=",gen_2,gen_2);
146 : : /* only p < 2^32 such that znprimroot(p) != znprimroot(p^2) */
147 [ - + ]: 52 : if (p == 40487) return 10;
148 : : #ifndef LONG_IS_64BIT
149 : 4 : return pgener_Fl(p);
150 : : #else
151 [ + + ]: 48 : if (p < (1UL<<32)) return pgener_Fl(p);
152 : : else
153 : : {
154 : 24 : const pari_sp av = avma;
155 : 24 : const ulong p_1 = p-1;
156 : : long x ;
157 : 24 : GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
158 : 24 : for (x=2;;x++)
159 [ + + ][ + - ]: 96 : if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2))) break;
160 : 120 : avma = av; return x;
161 : : }
162 : : #endif
163 : : }
164 : :
165 : : /* p prime. Return a primitive root modulo p^e, e > 1 */
166 : : GEN
167 : 56 : pgener_Zp(GEN p)
168 : : {
169 [ + + ]: 56 : if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
170 : : else
171 : : {
172 : 4 : const pari_sp av = avma;
173 : 4 : GEN p_1 = subis(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
174 : 4 : GEN x = utoipos(2);
175 : 12 : for (;; x[2]++)
176 [ + + ][ + - ]: 16 : if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
177 : 68 : avma = av; return utoipos(uel(x,2));
178 : : }
179 : : }
180 : :
181 : : static GEN
182 : 224 : gener_Zp(GEN q, GEN F)
183 : : {
184 : 224 : GEN p = NULL;
185 : 224 : long e = 0;
186 [ + + ]: 224 : if (F)
187 : : {
188 : 14 : GEN P = gel(F,1), E = gel(F,2);
189 : 14 : long i, l = lg(P);
190 [ + + ]: 42 : for (i = 1; i < l; i++)
191 : : {
192 : 28 : p = gel(P,i);
193 [ + + ]: 28 : if (equaliu(p, 2)) continue;
194 [ - + ]: 14 : if (i < l-1) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
195 : 14 : e = itos(gel(E,i));
196 : : }
197 [ - + ]: 14 : if (!p) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
198 : : }
199 : : else
200 : 210 : e = Z_isanypower(q, &p);
201 [ + + ]: 224 : return e > 1? pgener_Zp(p): pgener_Fp(q);
202 : : }
203 : :
204 : : GEN
205 : 294 : znprimroot(GEN N)
206 : : {
207 : 294 : pari_sp av = avma;
208 : : GEN x, n, F;
209 : :
210 [ + + ]: 294 : if ((F = check_arith_non0(N,"znprimroot")))
211 : : {
212 : 14 : F = clean_Z_factor(F);
213 [ + + ]: 14 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
214 : : }
215 [ - + ]: 287 : if (signe(N) < 0) N = absi(N);
216 [ + + ]: 287 : if (cmpiu(N, 4) <= 0) { avma = av; return mkintmodu(N[2]-1,N[2]); }
217 [ + + + ]: 238 : switch(mod4(N))
218 : : {
219 : : case 0: /* N = 0 mod 4 */
220 : 14 : pari_err_DOMAIN("znprimroot", "argument","=",N,N);
221 : 0 : x = NULL; break;
222 : : case 2: /* N = 2 mod 4 */
223 : 21 : n = shifti(N,-1); /* becomes odd */
224 [ - + ]: 21 : x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
225 : 21 : break;
226 : : default: /* N odd */
227 : 203 : x = gener_Zp(N,F);
228 : 203 : break;
229 : : }
230 : 273 : return gerepilecopy(av, mkintmod(x, N));
231 : : }
232 : :
233 : : /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
234 : : GEN
235 : 0 : rootsof1_Fp(GEN n, GEN p)
236 : : {
237 : 0 : pari_sp av = avma;
238 : 0 : GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
239 : 0 : GEN z = pgener_Fp_local(p, L);
240 : 0 : z = Fp_pow(z, diviiexact(subis(p,1), n), p); /* prim. n-th root of 1 */
241 : 0 : return gerepileuptoint(av, z);
242 : : }
243 : :
244 : : GEN
245 : 322 : rootsof1u_Fp(ulong n, GEN p)
246 : : {
247 : 322 : pari_sp av = avma;
248 : 322 : GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
249 : 322 : z = pgener_Fp_local(p, Flv_to_ZV(L));
250 : 322 : z = Fp_pow(z, diviuexact(subis(p,1), n), p); /* prim. n-th root of 1 */
251 : 322 : return gerepileuptoint(av, z);
252 : : }
253 : :
254 : : ulong
255 : 546 : rootsof1_Fl(ulong n, ulong p)
256 : : {
257 : 546 : pari_sp av = avma;
258 : 546 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
259 : 546 : ulong z = pgener_Fl_local(p, L);
260 : 546 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
261 : 546 : avma = av; return z;
262 : : }
263 : :
264 : : GEN
265 : 3598 : znstar(GEN N)
266 : : {
267 : 3598 : GEN F = NULL, P, E, cyc, gen, mod;
268 : : long i, j, nbp, sizeh;
269 : 3598 : pari_sp av = avma;
270 : :
271 [ + + ]: 3598 : if ((F = check_arith_all(N,"znstar")))
272 : : {
273 : 14 : F = clean_Z_factor(F);
274 [ - + ]: 14 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
275 : : }
276 [ + + ]: 3598 : if (!signe(N))
277 : : {
278 : 14 : avma = av;
279 : 14 : retmkvec3(gen_2, mkvec(gen_2), mkvec(gen_m1));
280 : : }
281 [ + + ]: 3584 : if (signe(N) < 0) N = absi(N);
282 [ + + ]: 3584 : if (cmpiu(N,2) <= 0)
283 : : {
284 : 28 : avma = av;
285 : 28 : retmkvec3(gen_1, cgetg(1,t_VEC), cgetg(1,t_VEC));
286 : : }
287 [ + + ]: 3556 : if (!F) F = Z_factor(N);
288 : 3556 : P = gel(F,1);
289 : 3556 : E = gel(F,2); nbp = lg(P)-1;
290 : 3556 : cyc = cgetg(nbp+2,t_VEC);
291 : 3556 : gen = cgetg(nbp+2,t_VEC);
292 : 3556 : mod = cgetg(nbp+2,t_VEC);
293 [ + + + + ]: 3556 : switch(mod8(N))
294 : : {
295 : : case 0: {
296 : 140 : long v2 = itos(gel(E,1));
297 : 140 : gel(cyc,1) = int2n(v2-2);
298 : 140 : gel(cyc,2) = gen_2;
299 : 140 : gel(gen,1) = utoipos(5);
300 : 140 : gel(gen,2) = addis(int2n(v2-1), -1);
301 : 140 : gel(mod,1) = gel(mod,2) = int2n(v2);
302 : 140 : sizeh = nbp+1; i = 3; j = 2; break;
303 : : }
304 : : case 4:
305 : 217 : gel(cyc,1) = gen_2;
306 : 217 : gel(gen,1) = utoipos(3);
307 : 217 : gel(mod,1) = utoipos(4);
308 : 217 : sizeh = nbp; i = j = 2; break;
309 : : case 2: case 6:
310 : 21 : sizeh = nbp-1; i=1; j=2; break;
311 : : default: /* 1, 3, 5, 7 */
312 : 3178 : sizeh = nbp; i = j = 1;
313 : : }
314 [ + + ]: 6832 : for ( ; j<=nbp; i++,j++)
315 : : {
316 : 3276 : long e = itos(gel(E,j));
317 : 3276 : GEN p = gel(P,j), q = powiu(p, e-1), Q = mulii(p, q);
318 : 3276 : gel(cyc,i) = subii(Q, q); /* phi(p^e) */
319 [ + + ]: 3276 : gel(gen,i) = e > 1? pgener_Zp(p): pgener_Fp(p);
320 : 3276 : gel(mod,i) = Q;
321 : : }
322 : : /* gen[i] has order cyc[i] and generates (Z/mod[i]Z)^* */
323 : 3556 : setlg(gen, sizeh+1);
324 : 3556 : setlg(cyc, sizeh+1);
325 [ + + ]: 3556 : if (nbp > 1)
326 [ + + ]: 238 : for (i=1; i<=sizeh; i++)
327 : : {
328 : 168 : GEN Q = gel(mod,i), g = gel(gen,i), qinv = Fp_inv(Q, diviiexact(N,Q));
329 : 168 : g = addii(g, mulii(mulii(subsi(1,g),qinv),Q));
330 : 168 : gel(gen,i) = modii(g, N);
331 : : }
332 : :
333 : : /*The cyc[i] are > 1 and remain so in the loop, gen[i] = 1 mod (N/mod[i]) */
334 [ + + ]: 3773 : for (i=sizeh; i>=2; i--)
335 : : {
336 : 217 : GEN ci = gel(cyc,i), gi = gel(gen,i);
337 [ + + ]: 497 : for (j=i-1; j>=1; j--) /* we want cyc[i] | cyc[j] */
338 : : {
339 : 280 : GEN cj = gel(cyc,j), gj, qj, v, d;
340 : :
341 : 280 : d = bezout(ci,cj,NULL,&v); /* > 1 */
342 [ + + ]: 280 : if (absi_equal(ci, d)) continue; /* ci | cj */
343 [ + + ]: 105 : if (absi_equal(cj, d)) { /* cj | ci */
344 : 84 : swap(gel(gen,j),gel(gen,i)); gi = gel(gen,i);
345 : 84 : swap(gel(cyc,j),gel(cyc,i)); ci = gel(cyc,i); continue;
346 : : }
347 : :
348 : 21 : gj = gel(gen,j);
349 : 21 : qj = diviiexact(cj,d);
350 : 21 : gel(cyc,j) = mulii(ci,qj);
351 : 21 : gel(cyc,i) = d;
352 : : /* [1,v*cj/d; 0,1]*[1,0;-1,1]*diag(cj,ci)*[ci/d,-v; cj/d,u]
353 : : * = diag(lcm,gcd), with u ci + v cj = d */
354 : :
355 : : /* (gj, gi) *= [1,0; -1,1]^-1 */
356 : 21 : gj = Fp_mul(gj, gi, N); /* order ci*qj = lcm(ci,cj) */
357 : : /* (gj,gi) *= [1,v*qj; 0,1]^-1 */
358 : 21 : togglesign_safe(&v);
359 [ - + ]: 21 : if (signe(v) < 0) v = modii(v,ci); /* >= 0 to avoid inversions */
360 : 21 : gi = Fp_mul(gi, Fp_pow(gj, mulii(qj, v), N), N);
361 : 21 : gel(gen,i) = gi;
362 : 21 : gel(gen,j) = gj;
363 [ + - ]: 280 : ci = d; if (equaliu(ci, 2)) break;
364 : : }
365 : : }
366 : 3598 : return gerepilecopy(av, mkvec3(ZV_prod(cyc), cyc, FpV_to_mod(gen,N)));
367 : : }
368 : :
369 : : /*********************************************************************/
370 : : /** **/
371 : : /** INVERSE TOTIENT FUNCTION **/
372 : : /** **/
373 : : /*********************************************************************/
374 : : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
375 : : * primes. Return factor(N) */
376 : : GEN
377 : 350651 : Z_factor_listP(GEN N, GEN L)
378 : : {
379 : 350651 : long i, k, l = lg(L);
380 : 350651 : GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
381 [ + + ]: 1346688 : for (i = k = 1; i < l; i++)
382 : : {
383 : 996037 : GEN p = gel(L,i);
384 : 996037 : long v = Z_pvalrem(N, p, &N);
385 [ + + ]: 996037 : if (v)
386 : : {
387 : 792176 : gel(P,k) = p;
388 : 792176 : gel(E,k) = utoipos(v);
389 : 792176 : k++;
390 : : }
391 : : }
392 : 350651 : setlg(P, k);
393 : 350651 : setlg(E, k); return mkmat2(P,E);
394 : : }
395 : :
396 : : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
397 : : * L is a list of primes containing all prime divisors of n. */
398 : : static long
399 : 621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
400 : : {
401 : 621565 : pari_sp av = avma, av2;
402 : : GEN k, D;
403 : : long i, v;
404 [ + + ][ + + ]: 621565 : if (m && mod2(n))
405 : : {
406 [ + + ]: 270914 : if (!equali1(n)) return 0;
407 [ + - ]: 69986 : if (px) *px = gen_1;
408 : 69986 : return 1;
409 : : }
410 : 350651 : D = divisors(Z_factor_listP(shifti(n, -1), L));
411 : : /* loop through primes p > m, d = p-1 | n */
412 : 350651 : av2 = avma;
413 [ + + ]: 350651 : if (!m)
414 : : { /* special case p = 2, d = 1 */
415 : 69986 : k = n;
416 : 69986 : for (v = 1;; v++) {
417 [ + - ]: 69986 : if (istotient_i(k, gen_2, L, px)) {
418 [ + - ]: 69986 : if (px) *px = shifti(*px, v);
419 : 69986 : return 1;
420 : : }
421 [ # # ]: 0 : if (mod2(k)) break;
422 : 0 : k = shifti(k,-1);
423 : 0 : }
424 : 0 : avma = av2;
425 : : }
426 [ + + ]: 1099462 : for (i = 1; i < lg(D); ++i)
427 : : {
428 : 1001588 : GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
429 [ + - ][ + + ]: 1001588 : if (m && cmpii(d, m) < 0) continue;
430 : 677782 : p = addiu(d, 1);
431 [ + + ]: 677782 : if (!isprime(p)) continue;
432 : 442064 : k = diviiexact(n, d);
433 : 442064 : for (v = 1;; v++) {
434 : : GEN r;
435 [ + + ]: 481593 : if (istotient_i(k, p, L, px)) {
436 [ + - ]: 182791 : if (px) *px = mulii(*px, powiu(p, v));
437 : 182791 : return 1;
438 : : }
439 : 298802 : k = dvmdii(k, p, &r);
440 [ + + ]: 298802 : if (r != gen_0) break;
441 : 39529 : }
442 : 259273 : avma = av2;
443 : : }
444 : 621565 : avma = av; return 0;
445 : : }
446 : :
447 : : /* find x such that phi(x) = n */
448 : : long
449 : 70000 : istotient(GEN n, GEN *px)
450 : : {
451 : 70000 : pari_sp av = avma;
452 [ - + ]: 70000 : if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
453 [ - + ]: 70000 : if (signe(n) < 1) return 0;
454 [ + + ]: 70000 : if (mod2(n))
455 : : {
456 [ - + ]: 14 : if (!equali1(n)) return 0;
457 [ + - ]: 14 : if (px) *px = gen_1;
458 : 14 : return 1;
459 : : }
460 [ + - ]: 69986 : if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
461 : : {
462 [ - + ]: 69986 : if (!px) avma = av;
463 : : else
464 : 69986 : *px = gerepileuptoint(av, *px);
465 : 69986 : return 1;
466 : : }
467 : 70000 : avma = av; return 0;
468 : : }
469 : :
470 : : /*********************************************************************/
471 : : /** **/
472 : : /** INTEGRAL LOGARITHM **/
473 : : /** **/
474 : : /*********************************************************************/
475 : :
476 : : /* y > 1 and B > 0 integers. Return e such that y^(e-1) <= B < y^e, i.e
477 : : * e = 1 + floor(log_y B). Set *ptq = y^e if non-NULL */
478 : : long
479 : 214713 : logint(GEN B, GEN y, GEN *ptq)
480 : : {
481 : 214713 : pari_sp av = avma, av2;
482 : : long eB, ey, e, i, fl;
483 : 214713 : GEN q,pow2, r = y;
484 : :
485 [ - + ]: 214713 : if (typ(B) != t_INT) B = ceil_safe(B);
486 : 214713 : eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
487 : 214713 : ey = expi(y); /* result e satisfies e > eB / (ey+1) */
488 [ + + ]: 214713 : if (eB <= 13 * ey) /* e small, be naive */
489 : : {
490 : 192941 : for (e=1;; e++)
491 : : { /* here, r = y^e */
492 : 583593 : fl = cmpii(r, B);
493 [ + + ]: 583593 : if (fl > 0) goto END;
494 : 390652 : r = mulii(r,y);
495 : 390652 : }
496 : : }
497 : : /* e > 13 ey / (ey + 1) >= 6.5 */
498 : :
499 : : /* binary splitting: compute bits of e one by one */
500 : : /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
501 : 21772 : pow2 = new_chunk((long)log2(eB)+2);
502 : 21772 : gel(pow2,0) = y;
503 : 21772 : for (i=0,q=r;; )
504 : : { /* r = y^2^i */
505 : 118285 : fl = cmpii(r,B);
506 [ - + ]: 118285 : if (!fl) { e = 1L<<i; e++; r = mulii(r,y); goto END; }
507 [ + + ]: 118285 : if (fl > 0) break;
508 : 96513 : q = r; r = sqri(q);
509 : 96513 : i++; gel(pow2,i) = r;
510 : 96513 : }
511 : :
512 : 21772 : av2 = avma;
513 : 21772 : for (i--, e=1L<<i;;)
514 : : { /* y^e = q < B <= r = q * y^(2^i) */
515 [ + + ]: 96499 : if (!fl) break; /* B = r */
516 : : /* q < B < r */
517 [ + + ][ + + ]: 96492 : if (--i < 0) { if (fl > 0) e++; break; }
518 : 74727 : r = mulii(q, gel(pow2,i));
519 : 74727 : fl = cmpii(r, B);
520 [ + + ]: 74727 : if (fl <= 0) { e += (1L<<i); q = r = gerepileuptoint(av2, r); }
521 : 74727 : }
522 [ + + ]: 21772 : if (fl <= 0) { e++; r = mulii(r,y); }
523 : : END:
524 [ + + ]: 214713 : if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
525 : 214713 : return e;
526 : : }
527 : :
528 : : long
529 : 42 : logint0(GEN B, GEN y, GEN *ptq)
530 : : {
531 : : long e;
532 [ - + ]: 42 : if (typ(B) != t_INT) pari_err_TYPE("logint",B);
533 [ - + ]: 42 : if (signe(B)<=0)
534 : 0 : pari_err_DOMAIN("logint", "x" ,"<=", gen_0, B);
535 [ - + ]: 42 : if (typ(y) != t_INT) pari_err_TYPE("logint",y);
536 [ + - ][ - + ]: 42 : if (signe(y)<=0 || equali1(y))
537 : 0 : pari_err_DOMAIN("logint", "b" ,"<=", gen_1, y);
538 [ + + ]: 42 : if (equaliu(y, 2))
539 : : {
540 : 21 : e = expi(B);
541 [ + + ]: 21 : if (ptq) *ptq = int2n(e);
542 : 21 : return e;
543 : : }
544 : 21 : e = logint(B,y,ptq)-1;
545 [ - + ]: 21 : if (ptq)
546 : : {
547 : 0 : pari_sp av = avma;
548 : 0 : *ptq = gerepileuptoint(av, diviiexact(*ptq, y));
549 : : }
550 : 42 : return e;
551 : : }
552 : :
553 : : /*********************************************************************/
554 : : /** **/
555 : : /** INTEGRAL SQUARE ROOT **/
556 : : /** **/
557 : : /*********************************************************************/
558 : : GEN
559 : 17214 : sqrtint(GEN a)
560 : : {
561 [ - + ]: 17214 : if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
562 [ + + + ]: 17214 : switch (signe(a))
563 : : {
564 : 17200 : case 1: return sqrti(a);
565 : 7 : case 0: return gen_0;
566 : 7 : default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
567 : : }
568 : 17207 : return NULL; /* not reached */
569 : : }
570 : :
571 : : /*********************************************************************/
572 : : /** **/
573 : : /** PERFECT SQUARE **/
574 : : /** **/
575 : : /*********************************************************************/
576 : : static int
577 : 8661933 : carremod(ulong A)
578 : : {
579 : 8661933 : const int carresmod64[]={
580 : : 1,1,0,0,1,0,0,0,0,1, 0,0,0,0,0,0,1,1,0,0, 0,0,0,0,0,1,0,0,0,0,
581 : : 0,0,0,1,0,0,1,0,0,0, 0,1,0,0,0,0,0,0,0,1, 0,0,0,0,0,0,0,1,0,0, 0,0,0,0};
582 : 8661933 : const int carresmod63[]={
583 : : 1,1,0,0,1,0,0,1,0,1, 0,0,0,0,0,0,1,0,1,0, 0,0,1,0,0,1,0,0,1,0,
584 : : 0,0,0,0,0,0,1,1,0,0, 0,0,0,1,0,0,1,0,0,1, 0,0,0,0,0,0,0,0,1,0, 0,0,0};
585 : 8661933 : const int carresmod65[]={
586 : : 1,1,0,0,1,0,0,0,0,1, 1,0,0,0,1,0,1,0,0,0, 0,0,0,0,0,1,1,0,0,1,
587 : : 1,0,0,0,0,1,1,0,0,1, 1,0,0,0,0,0,0,0,0,1, 0,1,0,0,0,1,1,0,0,0, 0,1,0,0,1};
588 : 8661933 : const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
589 : 8661933 : return (carresmod64[A & 0x3fUL]
590 [ + + ]: 3176156 : && carresmod63[A % 63UL]
591 [ + + ]: 1896709 : && carresmod65[A % 65UL]
592 [ + + ][ + + ]: 11838089 : && carresmod11[A % 11UL]);
593 : : }
594 : :
595 : : /* emulate Z_issquareall on single-word integers */
596 : : long
597 : 8496249 : uissquareall(ulong A, ulong *sqrtA)
598 : : {
599 [ - + ]: 8496249 : if (!A) { *sqrtA = 0; return 1; }
600 [ + + ]: 8496249 : if (carremod(A))
601 : : {
602 : 1457792 : ulong a = usqrt(A);
603 [ + + ]: 1457792 : if (a * a == A) { *sqrtA = a; return 1; }
604 : : }
605 : 8496249 : return 0;
606 : : }
607 : : long
608 : 115438 : uissquare(ulong A)
609 : : {
610 [ - + ]: 115438 : if (!A) return 1;
611 [ + + ]: 115438 : if (carremod(A))
612 : : {
613 : 2369 : ulong a = usqrt(A);
614 [ + + ]: 2369 : if (a * a == A) return 1;
615 : : }
616 : 115438 : return 0;
617 : : }
618 : :
619 : : long
620 : 2850231 : Z_issquareall(GEN x, GEN *pt)
621 : : {
622 : : pari_sp av;
623 : : GEN y, r;
624 : :
625 [ + + + ]: 2850231 : switch(signe(x))
626 : : {
627 : 112357 : case -1: return 0;
628 [ - + ]: 35 : case 0: if (pt) *pt=gen_0; return 1;
629 : : }
630 [ + + ]: 2737839 : if (lgefint(x) == 3)
631 : : {
632 : 2687593 : ulong u = uel(x,2), a;
633 [ + + ]: 2687593 : if (!pt) return uissquare(u);
634 [ + + ]: 2572155 : if (!uissquareall(u, &a)) return 0;
635 : 2687593 : *pt = utoipos(a); return 1;
636 : : }
637 [ + + ]: 50246 : if (!carremod(umodiu(x, 64*63*65*11))) return 0;
638 : 5936 : av = avma; y = sqrtremi(x, &r);
639 [ + + ]: 5936 : if (r != gen_0) { avma = av; return 0; }
640 [ + + ]: 5350 : if (pt) { *pt = y; avma = (pari_sp)y; } else avma = av;
641 : 2850231 : return 1;
642 : : }
643 : :
644 : : /* a t_INT, p prime */
645 : : long
646 : 0 : Zp_issquare(GEN a, GEN p)
647 : : {
648 : : long v;
649 : : GEN ap;
650 : :
651 [ # # ][ # # ]: 0 : if (!signe(a) || gequal1(a)) return 1;
652 : 0 : v = Z_pvalrem(a, p, &ap);
653 [ # # ]: 0 : if (v&1) return 0;
654 [ # # ]: 0 : return equaliu(p, 2)? umodiu(ap, 8) == 1
655 : 0 : : kronecker(ap,p) == 1;
656 : : }
657 : :
658 : : static long
659 : 1260 : polissquareall(GEN x, GEN *pt)
660 : : {
661 : : pari_sp av;
662 : : long v;
663 : : GEN y, a, b, p;
664 : :
665 [ + + ]: 1260 : if (!signe(x))
666 : : {
667 [ + - ]: 7 : if (pt) *pt = gcopy(x);
668 : 7 : return 1;
669 : : }
670 [ + + ]: 1253 : if (odd(degpol(x))) return 0; /* odd degree */
671 : 1008 : av = avma;
672 : 1008 : v = RgX_valrem(x, &x);
673 [ + + ]: 1008 : if (v & 1) { avma = av; return 0; }
674 : 1001 : a = gel(x,2); /* test constant coeff */
675 [ + + ]: 1001 : if (!pt)
676 [ - + ]: 63 : { if (!issquare(a)) { avma = av; return 0; } }
677 : : else
678 [ - + ]: 938 : { if (!issquareall(a,&b)) { avma = av; return 0; } }
679 [ + + ]: 1001 : if (!degpol(x)) { /* constant polynomial */
680 [ + + ]: 77 : if (!pt) { avma = av; return 1; }
681 : 35 : y = scalarpol(b, varn(x)); goto END;
682 : : }
683 : 924 : p = characteristic(x);
684 [ + + ][ + + ]: 924 : if (signe(p) && !mod2(p))
685 : : {
686 : : long i, lx;
687 [ + + ]: 35 : if (!equaliu(p,2)) pari_err_IMPL("issquare for even characteristic != 2");
688 : 28 : x = gmul(x, mkintmod(gen_1, gen_2));
689 : 28 : lx = lg(x);
690 [ + + ]: 28 : if ((lx-3) & 1) { avma = av; return 0; }
691 [ + + ]: 49 : for (i = 3; i < lx; i+=2)
692 [ - + ]: 28 : if (!gequal0(gel(x,i))) { avma = av; return 0; }
693 [ + + ]: 21 : if (pt) {
694 : 14 : y = cgetg((lx+3) / 2, t_POL);
695 [ + + ]: 49 : for (i = 2; i < lx; i+=2)
696 [ - + ]: 35 : if (!issquareall(gel(x,i), &gel(y, (i+2)>>1))) { avma = av; return 0; }
697 : 14 : y[1] = evalsigne(1) | evalvarn(varn(x));
698 : 14 : goto END;
699 : : } else {
700 [ + + ]: 21 : for (i = 2; i < lx; i+=2)
701 [ - + ]: 14 : if (!issquare(gel(x,i))) { avma = av; return 0; }
702 : 7 : avma = av; return 1;
703 : : }
704 : : }
705 : : else
706 : : {
707 : 889 : long m = 1;
708 : 889 : x = RgX_Rg_div(x,a);
709 : : /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
710 [ + + ]: 889 : if (!signe(p)) x = RgX_deflate_max(x,&m);
711 : 889 : y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
712 [ + + ]: 889 : if (!RgX_equal(RgX_sqr(y), x)) { avma = av; return 0; }
713 [ - + ]: 644 : if (!pt) { avma = av; return 1; }
714 [ + + ]: 644 : if (!gequal1(a)) y = gmul(b, y);
715 [ + + ]: 889 : if (m != 1) y = RgX_inflate(y,m);
716 : : }
717 : : END:
718 [ + + ]: 693 : if (v) y = RgX_shift_shallow(y, v>>1);
719 : 1253 : *pt = gerepilecopy(av, y); return 1;
720 : : }
721 : :
722 : : /* b unit mod p */
723 : : static int
724 : 252 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
725 : : {
726 [ + + ]: 252 : if (d == 1)
727 : : { /* mod p: faster */
728 [ - + ]: 168 : if (!Fp_ispower(b, K, p)) return 0;
729 [ + + ]: 168 : if (pt) *pt = Fp_sqrtn(b, K, p, NULL);
730 : : }
731 : : else
732 : : { /* mod p^{2 +} */
733 [ + + ]: 84 : if (!ispower(cvtop(b, p, d), K, pt)) return 0;
734 [ + + ]: 63 : if (pt) *pt = gtrunc(*pt);
735 : : }
736 : 252 : return 1;
737 : : }
738 : :
739 : : /* We're studying whether a mod (q*p^e) is a K-th power, (q,p) = 1.
740 : : * Decide mod p^e, then reduce a mod q unless q = NULL. */
741 : : static int
742 : 392 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
743 : : {
744 : : GEN t, A;
745 : 392 : long v = Z_pvalrem(*pa, p, &A), d = e - v;
746 [ + + ]: 392 : if (d <= 0) t = gen_0;
747 : : else
748 : : {
749 : : ulong r;
750 : 336 : v = udivui_rem(v, K, &r);
751 [ + + ][ + + ]: 336 : if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
[ + + ]
752 [ + + ][ + + ]: 231 : if (L && v) t = mulii(t, powiu(p, v));
753 : : }
754 [ + - ]: 287 : if (q) *pa = modii(*pa, q);
755 [ + + ]: 287 : if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
756 : 392 : return 1;
757 : : }
758 : : long
759 : 280 : Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
760 : : {
761 : : GEN L, N;
762 : : pari_sp av;
763 : : long e, i, l;
764 : : ulong pp;
765 : : forprime_t S;
766 : :
767 [ + + ]: 280 : if (!signe(a))
768 : : {
769 [ + - ]: 7 : if (pt) {
770 : 7 : GEN t = cgetg(3, t_INTMOD);
771 : 7 : gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
772 : : }
773 : 7 : return 1;
774 : : }
775 : : /* a != 0 */
776 : 273 : av = avma;
777 : :
778 [ - + ]: 273 : if (typ(q) != t_INT) /* integer factorization */
779 : : {
780 : 0 : GEN P = gel(q,1), E = gel(q,2);
781 : 0 : l = lg(P);
782 [ # # ]: 0 : L = pt? vectrunc_init(l): NULL;
783 [ # # ]: 0 : for (i = 1; i < l; i++)
784 : : {
785 : 0 : GEN p = gel(P,i);
786 : 0 : long e = itos(gel(E,i));
787 [ # # ]: 0 : if (!handle_pe(&a, NULL, L, K, p, e)) { avma = av; return 0; }
788 : : }
789 : 0 : goto END;
790 : : }
791 [ + + ]: 273 : if (!mod2(K)
792 [ + + ]: 182 : && kronecker(a, shifti(q,-vali(q))) == -1) { avma = av; return 0; }
793 [ + + ]: 266 : L = pt? vectrunc_init(expi(q)+1): NULL;
794 : 266 : u_forprime_init(&S, 2, tridiv_bound(q));
795 [ + + ]: 883316 : while ((pp = u_forprime_next(&S)))
796 : : {
797 : : int stop;
798 : 883162 : e = Z_lvalrem_stop(&q, pp, &stop);
799 [ + + ]: 883162 : if (!e) continue;
800 [ + + ]: 147 : if (!handle_pe(&a, q, L, K, utoipos(pp), e)) { avma = av; return 0; }
801 [ + + ]: 126 : if (stop)
802 : : {
803 [ + + ][ - + ]: 147 : if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
804 : : goto END;
805 : : }
806 : : }
807 : 154 : l = lg(primetab);
808 [ - + ]: 154 : for (i = 1; i < l; i++)
809 : : {
810 : 0 : GEN p = gel(primetab,i);
811 : 0 : e = Z_pvalrem(q, p, &q);
812 [ # # ]: 0 : if (!e) continue;
813 [ # # ]: 0 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
814 [ # # ]: 0 : if (is_pm1(q)) goto END;
815 : : }
816 : 154 : N = gcdii(a,q);
817 [ + + ]: 154 : if (!is_pm1(N))
818 : : {
819 [ + + ]: 112 : if (ifac_isprime(N))
820 : : {
821 : 70 : e = Z_pvalrem(q, N, &q);
822 [ + - ]: 70 : if (!handle_pe(&a, q, L, K, N, e)) { avma = av; return 0; }
823 : : }
824 : : else
825 : : {
826 : 42 : GEN part = ifac_start(N, 0);
827 : : for(;;)
828 : : {
829 : : long e;
830 : : GEN p;
831 [ + + ]: 84 : if (!ifac_next(&part, &p, &e)) break;
832 : 42 : e = Z_pvalrem(q, p, &q);
833 [ - + ]: 42 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
834 : 84 : }
835 : : }
836 : : }
837 [ + - ]: 84 : if (!is_pm1(q))
838 : : {
839 [ + + ]: 84 : if (ifac_isprime(q))
840 : : {
841 [ - + ]: 28 : if (!handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
842 : : }
843 : : else
844 : : {
845 : 56 : GEN part = ifac_start(q, 0);
846 : : for(;;)
847 : : {
848 : : long e;
849 : : GEN p;
850 [ + + ]: 140 : if (!ifac_next(&part, &p, &e)) break;
851 [ + + ]: 98 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
852 : 126 : }
853 : : }
854 : : }
855 : : END:
856 [ + + ]: 161 : if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
857 : 280 : return 1;
858 : : }
859 : :
860 : : long
861 : 145481 : issquareall(GEN x, GEN *pt)
862 : : {
863 : 145481 : long tx = typ(x);
864 : : GEN F;
865 : : pari_sp av;
866 : :
867 [ + + ]: 145481 : if (!pt) return issquare(x);
868 [ + + + + : 2660 : switch(tx)
+ + + - ]
869 : : {
870 : 1176 : case t_INT: return Z_issquareall(x, pt);
871 : 182 : case t_FRAC: av = avma;
872 : 182 : F = cgetg(3, t_FRAC);
873 [ + + ]: 182 : if ( !Z_issquareall(gel(x,1), &gel(F,1))
874 [ - + ]: 182 : || !Z_issquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
875 : 105 : *pt = F; return 1;
876 : :
877 : 1183 : case t_POL: return polissquareall(x,pt);
878 : 7 : case t_RFRAC: av = avma;
879 : 7 : F = cgetg(3, t_RFRAC);
880 [ + - ]: 7 : if ( !issquareall(gel(x,1), &gel(F,1))
881 [ - + ]: 7 : || !polissquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
882 : 7 : *pt = F; return 1;
883 : :
884 : : case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
885 [ - + ]: 7 : if (!issquare(x)) return 0;
886 : 7 : *pt = gsqrt(x, DEFAULTPREC); return 1;
887 : :
888 : : case t_INTMOD:
889 : 63 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
890 : :
891 : 42 : case t_FFELT: return FF_issquareall(x, pt);
892 : :
893 : : }
894 : 0 : pari_err_TYPE("issquareall",x);
895 : 145474 : return 0; /* not reached */
896 : : }
897 : :
898 : : long
899 : 142996 : issquare(GEN x)
900 : : {
901 : : pari_sp av;
902 : : GEN a, p;
903 : : long i, v;
904 : :
905 [ + + + + : 142996 : switch(typ(x))
+ + + + +
+ - ]
906 : : {
907 : : case t_INT:
908 : 142716 : return Z_issquare(x);
909 : :
910 : : case t_REAL:
911 : 14 : return (signe(x)>=0);
912 : :
913 : : case t_INTMOD:
914 : 70 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, NULL);
915 : :
916 : : case t_FRAC:
917 [ + - ][ + + ]: 21 : return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
918 : :
919 : 7 : case t_FFELT: return FF_issquareall(x, NULL);
920 : :
921 : : case t_COMPLEX:
922 : 7 : return 1;
923 : :
924 : : case t_PADIC:
925 [ - + ]: 63 : a = gel(x,4); if (!signe(a)) return 1;
926 [ + + ]: 63 : if (valp(x)&1) return 0;
927 : 49 : p = gel(x,2);
928 [ + + ]: 49 : if (!equaliu(p, 2)) return (kronecker(a,p) != -1);
929 : :
930 : 14 : v = precp(x); /* here p=2, a is odd */
931 [ + - ][ + + ]: 14 : if ((v>=3 && mod8(a) != 1 ) ||
[ - + ]
932 [ # # ]: 7 : (v==2 && mod4(a) != 1)) return 0;
933 : 7 : return 1;
934 : :
935 : : case t_POL:
936 : 70 : return polissquareall(x,NULL);
937 : :
938 : : case t_SER:
939 [ + + ]: 21 : if (!signe(x)) return 1;
940 [ + + ]: 14 : if (valp(x)&1) return 0;
941 : 7 : return issquare(gel(x,2));
942 : :
943 : : case t_RFRAC:
944 : 7 : av = avma; i = issquare(gmul(gel(x,1),gel(x,2)));
945 : 7 : avma = av; return i;
946 : : }
947 : 0 : pari_err_TYPE("issquare",x);
948 : 142989 : return 0; /* not reached */
949 : : }
950 [ # # ]: 0 : GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
951 [ # # ]: 0 : GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
952 : :
953 : : long
954 : 1386 : ispolygonal(GEN x, GEN S, GEN *N)
955 : : {
956 : 1386 : pari_sp av = avma;
957 : : GEN D, d, n;
958 [ - + ]: 1386 : if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
959 [ - + ]: 1386 : if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
960 [ - + ]: 1386 : if (cmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
961 [ - + ]: 1386 : if (signe(x) < 0) return 0;
962 [ + + ][ + - ]: 1386 : if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
963 [ + + ][ + - ]: 1260 : if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
964 : : /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
965 [ + + ]: 1134 : if (cmpiu(S, 1<<16) < 0) /* common case ! */
966 : : {
967 : 441 : ulong s = S[2], r;
968 [ + + ]: 441 : if (s == 4) return Z_issquareall(x, N);
969 [ - + ]: 378 : if (s == 3)
970 : 0 : D = addiu(shifti(x, 3), 1);
971 : : else
972 : 378 : D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
973 [ - + ]: 378 : if (!Z_issquareall(D, &d)) { avma = av; return 0; }
974 [ - + ]: 378 : if (s == 3)
975 : 0 : d = subiu(d, 1);
976 : : else
977 : 378 : d = addiu(d, s - 4);
978 : 378 : n = diviu_rem(d, 2*s - 4, &r);
979 [ - + ]: 441 : if (r) { avma = av; return 0; }
980 : : }
981 : : else
982 : : {
983 : 693 : GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
984 : 693 : D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
985 [ - + ]: 693 : if (!Z_issquareall(D, &d)) { avma = av; return 0; }
986 : 693 : d = addii(d, S_4);
987 : 693 : n = dvmdii(shifti(d,-1), S_2, &r);
988 [ - + ]: 693 : if (r != gen_0) { avma = av; return 0; }
989 : : }
990 [ + - ]: 1071 : if (N) *N = gerepileuptoint(av, n); else avma = av;
991 : 1386 : return 1;
992 : : }
993 : :
994 : : /*********************************************************************/
995 : : /** **/
996 : : /** PERFECT POWER **/
997 : : /** **/
998 : : /*********************************************************************/
999 : : static long
1000 : 329 : polispower(GEN x, GEN K, GEN *pt)
1001 : : {
1002 : : pari_sp av;
1003 : 329 : long v, k = itos(K);
1004 : : GEN y, a, b;
1005 : :
1006 [ + + ]: 329 : if (!signe(x))
1007 : : {
1008 [ - + ]: 7 : if (pt) *pt = gcopy(x);
1009 : 7 : return 1;
1010 : : }
1011 [ + + ]: 322 : if (degpol(x) % k) return 0; /* degree not multiple of k */
1012 : 315 : av = avma;
1013 : 315 : v = RgX_valrem(x, &x);
1014 [ + + ]: 315 : if (v % k) return 0;
1015 : 308 : a = gel(x,2); b = NULL;
1016 [ + + ]: 308 : if (!ispower(a, K, &b)) { avma = av; return 0; }
1017 [ + + ]: 294 : if (degpol(x))
1018 : : {
1019 : 154 : x = RgX_Rg_div(x,a);
1020 : 154 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
1021 [ - + ]: 154 : if (!RgX_equal(powgi(y, K), x)) { avma = av; return 0; }
1022 : : }
1023 : : else
1024 : 140 : y = pol_1(varn(x));
1025 [ + - ]: 294 : if (pt)
1026 : : {
1027 [ + + ]: 294 : if (!gequal1(a))
1028 : : {
1029 [ - + ]: 91 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
1030 : 91 : y = gmul(b,y);
1031 : : }
1032 [ + + ]: 294 : *pt = v? gerepilecopy(av, RgX_shift_shallow(y, v/k)): gerepileupto(av, y);
1033 : : }
1034 : 0 : else avma = av;
1035 : 329 : return 1;
1036 : : }
1037 : :
1038 : : long
1039 : 1085 : Z_ispowerall(GEN x, ulong k, GEN *pt)
1040 : : {
1041 : 1085 : long s = signe(x);
1042 : : ulong mask;
1043 [ - + ][ # # ]: 1085 : if (!s) { if (pt) *pt = gen_0; return 1; }
1044 [ + + ]: 1085 : if (s > 0) {
1045 [ + + ]: 980 : if (k == 2) return Z_issquareall(x, pt);
1046 [ + + ]: 721 : if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
1047 [ + + ]: 161 : if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
1048 [ + + ]: 147 : if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
1049 : 140 : return is_kth_power(x, k, pt);
1050 : : }
1051 [ + + ]: 105 : if (!odd(k)) return 0;
1052 [ + - ]: 91 : if (Z_ispowerall(absi(x), k, pt))
1053 : : {
1054 [ + + ]: 91 : if (pt) *pt = negi(*pt);
1055 : 91 : return 1;
1056 : : };
1057 : 1085 : return 0;
1058 : : }
1059 : :
1060 : : /* is x a K-th power mod p ? Assume p prime. */
1061 : : int
1062 : 168 : Fp_ispower(GEN x, GEN K, GEN p)
1063 : : {
1064 : 168 : pari_sp av = avma;
1065 : : GEN p_1;
1066 : : long r;
1067 : 168 : x = modii(x, p);
1068 [ + - ][ + + ]: 168 : if (!signe(x) || equali1(x)) { avma = av; return 1; }
1069 : : /* implies p > 2 */
1070 : 91 : p_1 = subiu(p,1);
1071 : 91 : K = gcdii(K, p_1);
1072 [ + + ]: 91 : if (equaliu(K, 2)) { r = kronecker(x,p); avma = av; return (r > 0); }
1073 : 35 : x = Fp_pow(x, diviiexact(p_1,K), p);
1074 : 168 : avma = av; return equali1(x);
1075 : : }
1076 : :
1077 : : /* x unit defined modulo 2^e, e > 0, p prime */
1078 : : static int
1079 : 2373 : U2_issquare(GEN x, long e)
1080 : : {
1081 [ - + ]: 2373 : long r = signe(x)>=0?mod8(x):8-mod8(x);
1082 [ - + ]: 2373 : if (e==1) return 1;
1083 [ + + ]: 2373 : if (e==2) return (r&3L) == 1;
1084 : 2373 : return r == 1;
1085 : : }
1086 : : /* x unit defined modulo p^e, e > 0, p prime */
1087 : : static int
1088 : 4690 : Up_issquare(GEN x, GEN p, long e)
1089 [ + + ]: 4690 : { return (equaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
1090 : :
1091 : : long
1092 : 2548 : Zn_issquare(GEN d, GEN fn)
1093 : : {
1094 : : long j, np;
1095 [ - + ]: 2548 : if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
1096 [ - + ]: 2548 : if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
1097 : : /* integer factorization */
1098 : 2548 : np = nbrows(fn);
1099 [ + + ]: 5320 : for (j = 1; j <= np; ++j)
1100 : : {
1101 : 4970 : GEN r, p = gcoeff(fn, j, 1);
1102 : 4970 : long e = itos(gcoeff(fn, j, 2));
1103 : 4970 : long v = Z_pvalrem(d,p,&r);
1104 [ + + ][ + - ]: 4970 : if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
[ + + ]
1105 : : }
1106 : 2548 : return 1;
1107 : : }
1108 : :
1109 : : long
1110 : 7002191 : ispower(GEN x, GEN K, GEN *pt)
1111 : : {
1112 : : GEN z;
1113 : :
1114 [ + + ]: 7002191 : if (!K) return gisanypower(x, pt);
1115 [ - + ]: 2037 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
1116 [ - + ]: 2037 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
1117 [ + + ][ + - ]: 2037 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
1118 [ + + + + : 2009 : switch(typ(x)) {
+ + + + +
+ - ]
1119 : : case t_INT:
1120 : 350 : return Z_ispowerall(x, itou(K), pt);
1121 : : case t_FRAC:
1122 : : {
1123 : 21 : GEN a = gel(x,1), b = gel(x,2);
1124 : 21 : ulong k = itou(K);
1125 [ + + ]: 21 : if (pt) {
1126 : 14 : z = cgetg(3, t_FRAC);
1127 [ + - ][ + - ]: 14 : if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
1128 : 14 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1129 : : }
1130 : 0 : avma = (pari_sp)(z + 3); return 0;
1131 : : }
1132 [ + - ][ + - ]: 21 : return Z_ispower(a, k) && Z_ispower(b, k);
1133 : : }
1134 : : case t_INTMOD:
1135 : 147 : return Zn_ispower(gel(x,2), gel(x,1), K, pt);
1136 : : case t_FFELT:
1137 : 28 : return FF_ispower(x, K, pt);
1138 : :
1139 : : case t_PADIC:
1140 : 1113 : z = Qp_sqrtn(x, K, NULL);
1141 [ + + ]: 1113 : if (!z) return 0;
1142 [ + + ]: 819 : if (pt) *pt = z;
1143 : 819 : return 1;
1144 : :
1145 : : case t_POL:
1146 : 322 : return polispower(x, K, pt);
1147 : : case t_RFRAC: {
1148 : 7 : GEN a = gel(x,1), b = gel(x,2);
1149 [ + - ]: 7 : if (pt) {
1150 : 7 : z = cgetg(3, t_RFRAC);
1151 [ + - ][ + - ]: 7 : if (ispower(a, K, &a) && polispower(b, K, &b)) {
1152 : 7 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1153 : : }
1154 : 0 : avma = (pari_sp)(z + 3); return 0;
1155 : : }
1156 [ # # ][ # # ]: 7 : return (ispower(a, K, NULL) && polispower(b, K, NULL));
1157 : : }
1158 : : case t_REAL:
1159 [ + - ][ - + ]: 7 : if (signe(x) < 0 && !mpodd(K)) return 0;
1160 : : case t_COMPLEX:
1161 [ - + ]: 14 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1162 : 14 : return 1;
1163 : :
1164 : : case t_SER:
1165 [ + - ][ + - ]: 7 : if (signe(x) && (!dvdsi(valp(x), K) || !ispower(gel(x,2), K, NULL)))
[ - + ]
1166 : 0 : return 0;
1167 [ + - ]: 7 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1168 : 7 : return 1;
1169 : :
1170 : 0 : default: pari_err_TYPE("ispower",x);
1171 : 7002191 : return 0; /* not reached */
1172 : : }
1173 : : }
1174 : :
1175 : : long
1176 : 7000154 : gisanypower(GEN x, GEN *pty)
1177 : : {
1178 : 7000154 : long tx = typ(x);
1179 : : ulong k, h;
1180 [ + + ]: 7000154 : if (tx == t_INT) return Z_isanypower(x, pty);
1181 [ + - ]: 14 : if (tx == t_FRAC)
1182 : : {
1183 : 14 : pari_sp av = avma;
1184 : 14 : GEN fa, P, E, a = gel(x,1), b = gel(x,2);
1185 : : long i, j, p, e;
1186 : 14 : int sw = (absi_cmp(a, b) > 0);
1187 : :
1188 [ - + ]: 14 : if (sw) swap(a, b);
1189 [ + + ]: 14 : k = Z_isanypower(a, pty? &a: NULL);
1190 [ + + ]: 14 : if (!k)
1191 : : { /* a = -1,1 or not a pure power */
1192 [ - + ]: 7 : if (!is_pm1(a)) { avma = av; return 0; }
1193 [ + - ]: 7 : if (signe(a) < 0) b = negi(b);
1194 [ + - ]: 7 : k = Z_isanypower(b, pty? &b: NULL);
1195 [ + - ][ - + ]: 7 : if (!k || !pty) { avma = av; return k; }
1196 : 7 : *pty = gerepileupto(av, ginv(b));
1197 : 7 : return k;
1198 : : }
1199 : 7 : fa = factoru(k);
1200 : 7 : P = gel(fa,1);
1201 : 7 : E = gel(fa,2); h = k;
1202 [ + + ]: 14 : for (i = lg(P) - 1; i > 0; i--)
1203 : : {
1204 : 7 : p = P[i];
1205 : 7 : e = E[i];
1206 [ + + ]: 21 : for (j = 0; j < e; j++)
1207 [ - + ]: 14 : if (!is_kth_power(b, p, &b)) break;
1208 [ - + ]: 7 : if (j < e) k /= upowuu(p, e - j);
1209 : : }
1210 [ - + ]: 7 : if (k == 1) { avma = av; return 0; }
1211 [ + - ]: 7 : if (!pty) { avma = av; return k; }
1212 [ # # ]: 0 : if (k != h) a = powiu(a, h/k);
1213 : 0 : *pty = gerepilecopy(av, mkfrac(a, b));
1214 : 14 : return k;
1215 : : }
1216 : 0 : pari_err_TYPE("gisanypower", x);
1217 : 7000154 : return 0; /* not reached */
1218 : : }
1219 : :
1220 : : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
1221 : : * No need to optimize for 2,3,5,7 powers (done before) */
1222 : : static long
1223 : 505708 : split_exponent(ulong e, GEN *x)
1224 : : {
1225 : : GEN fa, P, E;
1226 : 505708 : long i, j, l, k = 1;
1227 [ + + ]: 505708 : if (e == 1) return 1;
1228 : 7 : fa = factoru(e);
1229 : 7 : P = gel(fa,1);
1230 : 7 : E = gel(fa,2); l = lg(P);
1231 [ + + ]: 14 : for (i = 1; i < l; i++)
1232 : : {
1233 : 7 : ulong p = P[i];
1234 [ + + ]: 14 : for (j = 0; j < E[i]; j++)
1235 : : {
1236 : : GEN y;
1237 [ + - ]: 7 : if (!is_kth_power(*x, p, &y)) break;
1238 : 7 : k *= p; *x = y;
1239 : : }
1240 : : }
1241 : 505708 : return k;
1242 : : }
1243 : :
1244 : : static long
1245 : 864302 : Z_isanypower_nosmalldiv(GEN *px)
1246 : : { /* any prime divisor of x is > 102 */
1247 : 864302 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
1248 : 864302 : const double LOG103 = 4.6347; /* lower bound for log(103) */
1249 : : forprime_t T;
1250 : 864302 : ulong mask = 7, e2;
1251 : : long k, ex;
1252 : 864302 : GEN y, x = *px;
1253 : :
1254 : 864302 : k = 1;
1255 [ + + ]: 865492 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
1256 [ + + ]: 864408 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
1257 : 864302 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
1258 [ + + ]: 864302 : if (u_forprime_init(&T, 11, e2))
1259 : : {
1260 : 16800 : GEN logx = NULL;
1261 : 16800 : const ulong Q = 30011; /* prime */
1262 : : ulong p, xmodQ;
1263 : 16800 : double dlogx = 0;
1264 : : /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
1265 : : * for large p the modular checks are no longer competitively fast */
1266 [ + + ]: 16821 : while ( (ex = is_pth_power(x, &y, &T, 30)) )
1267 : : {
1268 : 21 : k *= ex; x = y;
1269 : 21 : e2 = (ulong)((expi(x) + 1) / LOG2_103);
1270 : 21 : u_forprime_restrict(&T, e2);
1271 : : }
1272 [ - + ]: 16800 : if (DEBUGLEVEL>4) err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x));
1273 : 16800 : xmodQ = umodiu(x, Q);
1274 : : /* test Q | x, just in case */
1275 [ + + ]: 16800 : if (!xmodQ) return k * split_exponent(Z_lval(x,Q), px);
1276 : : /* x^(1/p) < 2^31 */
1277 : 16793 : p = T.p;
1278 [ + + ]: 16793 : if (p <= e2)
1279 : : {
1280 : 16779 : logx = logr_abs( itor(x, DEFAULTPREC) );
1281 : 16779 : dlogx = rtodbl(logx);
1282 : 16779 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1283 : : }
1284 [ + + ][ + + ]: 134659 : while (p && p <= e2)
1285 : : { /* is x a p-th power ? By computing y = round(x^(1/p)).
1286 : : * Check whether y^p = x, first mod Q, then exactly. */
1287 : 117866 : pari_sp av = avma;
1288 : : long e;
1289 : 117866 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
1290 : 117866 : ulong ymodQ = umodiu(y,Q);
1291 [ + + ][ + + ]: 117866 : if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
1292 [ - + ]: 21 : || !equalii(powiu(y, p), x)) avma = av;
1293 : : else
1294 : : {
1295 : 21 : k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
1296 : 21 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1297 : 21 : u_forprime_restrict(&T, e2);
1298 : 21 : continue; /* if success, retry same p */
1299 : : }
1300 : 117845 : p = u_forprime_next(&T);
1301 : : }
1302 : : }
1303 : 864302 : *px = x; return k;
1304 : : }
1305 : :
1306 : : static ulong tinyprimes[] = {
1307 : : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
1308 : : 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
1309 : : 157, 163, 167, 173, 179, 181, 191, 193, 197, 199
1310 : : };
1311 : :
1312 : : /* disregard the sign of x, caller will take care of x < 0 */
1313 : : static long
1314 : 7000560 : Z_isanypower_aux(GEN x, GEN *pty)
1315 : : {
1316 : : long ex, v, i, l, k;
1317 : : GEN y, P, E;
1318 : 7000560 : ulong mask, e = 0;
1319 : :
1320 [ + + ]: 7000560 : if (absi_cmp(x, gen_2) < 0) return 0; /* -1,0,1 */
1321 : :
1322 : 7000546 : k = l = 1;
1323 : 7000546 : P = cgetg(26 + 1, t_VECSMALL);
1324 : 7000546 : E = cgetg(26 + 1, t_VECSMALL);
1325 : : /* trial division */
1326 [ + + ]: 67128845 : for(i = 0; i < 26; i++)
1327 : : {
1328 : 60128299 : ulong p = tinyprimes[i];
1329 : : int stop;
1330 : 60128299 : v = Z_lvalrem_stop(&x, p, &stop);
1331 [ + + ]: 60128299 : if (v)
1332 : : {
1333 : 7922411 : P[l] = p;
1334 : 7922411 : E[l] = v; l++;
1335 [ + + ]: 7922411 : e = ugcd(e, v); if (e == 1) goto END;
1336 : : }
1337 [ + + ]: 54709984 : if (stop) {
1338 [ + + ]: 248052 : if (is_pm1(x)) k = e;
1339 : : goto END;
1340 : : }
1341 : : }
1342 : :
1343 [ + + ]: 1334179 : if (e)
1344 : : { /* Bingo. Result divides e */
1345 : : long v3, v5, v7;
1346 : 505701 : ulong e2 = e;
1347 : 505701 : v = u_lvalrem(e2, 2, &e2);
1348 [ + + ]: 505701 : if (v)
1349 : : {
1350 [ + + ]: 375249 : for (i = 0; i < v; i++)
1351 : : {
1352 [ + + ]: 374171 : if (!Z_issquareall(x, &y)) break;
1353 : 1288 : k <<= 1; x = y;
1354 : : }
1355 : : }
1356 : 505701 : mask = 0;
1357 [ + + ]: 505701 : v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
1358 [ + + ]: 505701 : v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
1359 [ + + ]: 505701 : v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
1360 [ + + ]: 505778 : while ( (ex = is_357_power(x, &y, &mask)) ) {
1361 : 77 : x = y;
1362 [ + + + - ]: 77 : switch(ex)
1363 : : {
1364 [ + - ]: 28 : case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
1365 [ + - ]: 28 : case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
1366 [ + - ]: 21 : case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
1367 : : }
1368 : : }
1369 : 505701 : k *= split_exponent(e2, &x);
1370 : : }
1371 : : else
1372 : 828478 : k = Z_isanypower_nosmalldiv(&x);
1373 : : END:
1374 [ + + ][ + + ]: 7000546 : if (pty && k != 1)
1375 : : {
1376 [ + + ]: 7945 : if (e)
1377 : : { /* add missing small factors */
1378 : 6867 : y = powuu(P[1], E[1] / k);
1379 [ + + ]: 14021 : for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
1380 [ + + ]: 6867 : x = equali1(x)? y: mulii(x,y);
1381 : : }
1382 : 7945 : *pty = x;
1383 : : }
1384 [ + + ]: 7000560 : return k == 1? 0: k;
1385 : : }
1386 : :
1387 : : long
1388 : 7000560 : Z_isanypower(GEN x, GEN *pty)
1389 : : {
1390 : 7000560 : pari_sp av = avma;
1391 : 7000560 : long k = Z_isanypower_aux(x, pty);
1392 [ + + ]: 7000560 : if (!k) { avma = av; return 0; }
1393 [ + + ]: 7980 : if (signe(x) < 0)
1394 : : {
1395 : 35 : long v = vals(k);
1396 [ + + ]: 35 : if (v)
1397 : : {
1398 : : GEN y;
1399 : 21 : k >>= v;
1400 [ + + ]: 21 : if (k == 1) { avma = av; return 0; }
1401 [ - + ]: 14 : if (!pty) { avma = av; return k; }
1402 : 14 : y = *pty;
1403 : 14 : y = powiu(y, 1<<v);
1404 : 14 : togglesign(y);
1405 : 14 : *pty = gerepileuptoint(av, y);
1406 : 14 : return k;
1407 : : }
1408 [ + - ]: 14 : if (pty) togglesign_safe(pty);
1409 : : }
1410 [ + + ]: 7959 : if (pty) *pty = gerepilecopy(av, *pty); else avma = av;
1411 : 7000560 : return k;
1412 : : }
1413 : :
1414 : : /* Faster than */
1415 : : /* !cmpii(n, int2n(vali(n))) */
1416 : : /* !cmpis(shifti(n, -vali(n)), 1) */
1417 : : /* expi(n) == vali(n) */
1418 : : /* hamming(n) == 1 */
1419 : : /* even for single-word values, and much faster for multiword values. */
1420 : : /* If all you have is a word, you can just use n & !(n & (n-1)). */
1421 : : long
1422 : 22444 : Z_ispow2(GEN n)
1423 : : {
1424 : : GEN xp;
1425 : : long i, lx;
1426 : : ulong u;
1427 [ - + ]: 22444 : if (signe(n) != 1) return 0;
1428 : 22444 : xp = int_LSW(n);
1429 : 22444 : lx = lgefint(n);
1430 : 22444 : u = *xp;
1431 [ + + ]: 22476 : for (i = 3; i < lx; ++i)
1432 : : {
1433 [ + + ]: 20034 : if (u) return 0;
1434 : 32 : xp = int_nextW(xp);
1435 : 32 : u = *xp;
1436 : : }
1437 : 22444 : return !(u & (u-1)); /* faster than hamming_word(u) == 1 */
1438 : : }
1439 : :
1440 : : static long
1441 : 841524 : isprimepower_i(GEN n, GEN *pt, long flag)
1442 : : {
1443 : 841524 : pari_sp av = avma;
1444 : : long i, v;
1445 : :
1446 [ - + ]: 841524 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1447 [ - + ]: 841524 : if (signe(n) <= 0) return 0;
1448 : :
1449 [ + + ]: 841524 : if (lgefint(n) == 3)
1450 : : {
1451 : : ulong p;
1452 : 540963 : v = uisprimepower(n[2], &p);
1453 [ + + ]: 540963 : if (v)
1454 : : {
1455 [ + + ]: 54751 : if (pt) *pt = utoipos(p);
1456 : 54751 : return v;
1457 : : }
1458 : 540963 : return 0;
1459 : : }
1460 [ + + ]: 1661657 : for (i = 0; i < 26; i++)
1461 : : {
1462 : 1625833 : ulong p = tinyprimes[i];
1463 : 1625833 : v = Z_lvalrem(n, p, &n);
1464 [ + + ]: 1625834 : if (v)
1465 : : {
1466 : 264738 : avma = av;
1467 [ + + ]: 264738 : if (!is_pm1(n)) return 0;
1468 [ + + ]: 342 : if (pt) *pt = utoipos(p);
1469 : 340 : return v;
1470 : : }
1471 : : }
1472 : : /* p | n => p >= 103 */
1473 : 35824 : v = Z_isanypower_nosmalldiv(&n); /* expensive */
1474 [ + + ][ + + ]: 35824 : if (!(flag? isprime(n): BPSW_psp(n))) { avma = av; return 0; }
1475 [ + + ]: 5512 : if (pt) *pt = gerepilecopy(av, n); else avma = av;
1476 : 841523 : return v;
1477 : : }
1478 : : long
1479 : 840091 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
1480 : : long
1481 : 1433 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1482 : :
1483 : : long
1484 : 541572 : uisprimepower(ulong n, ulong *pp)
1485 : : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1486 : : * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1487 : 541572 : const ulong CUTOFF = 200UL;
1488 : 541572 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1489 : 541572 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1490 : : #ifdef LONG_IS_64BIT
1491 : : /* primes preceeding the appropriate root of ULONG_MAX. */
1492 : 481356 : const ulong ROOT9 = 137;
1493 : 481356 : const ulong ROOT8 = 251;
1494 : 481356 : const ulong ROOT7 = 563;
1495 : 481356 : const ulong ROOT5 = 7129;
1496 : 481356 : const ulong ROOT4 = 65521;
1497 : : #else
1498 : 60216 : const ulong ROOT9 = 11;
1499 : 60216 : const ulong ROOT8 = 13;
1500 : 60216 : const ulong ROOT7 = 23;
1501 : 60216 : const ulong ROOT5 = 83;
1502 : 60216 : const ulong ROOT4 = 251;
1503 : : #endif
1504 : : ulong mask;
1505 : : long v, i;
1506 : : int e;
1507 [ + + ]: 541572 : if (n < 2) return 0;
1508 [ + + ]: 541558 : if (!odd(n)) {
1509 [ + + ]: 270648 : if (n & (n-1)) return 0;
1510 : 851 : *pp = 2; return vals(n);
1511 : : }
1512 [ + + ]: 270910 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1513 [ + + ]: 3653526 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1514 : : {
1515 : 3594457 : ulong p = tinyprimes[i];
1516 [ + + ]: 3594457 : if (n % p == 0)
1517 : : {
1518 : 211470 : v = u_lvalrem(n, p, &n);
1519 [ + + ]: 211470 : if (n == 1) { *pp = p; return v; }
1520 : 209438 : return 0;
1521 : : }
1522 : : }
1523 : : /* p | n => p >= CUTOFF */
1524 : :
1525 [ + + ]: 59069 : if (n < CUTOFF3)
1526 : : {
1527 [ + + ][ + - ]: 46354 : if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1528 [ # # ]: 0 : if (uissquareall(n, &n)) { *pp = n; return 2; }
1529 : 0 : return 0;
1530 : : }
1531 : :
1532 : : /* Check for squares, fourth powers, and eighth powers as appropriate. */
1533 : 12715 : v = 1;
1534 [ - + ]: 12715 : if (uissquareall(n, &n)) {
1535 : 0 : v <<= 1;
1536 [ # # ][ # # ]: 0 : if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1537 : 0 : v <<= 1;
1538 [ # # ][ # # ]: 0 : if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1539 : : }
1540 : : }
1541 : :
1542 [ + + ]: 12715 : if (CUTOFF > ROOT5) mask = 1;
1543 : : else
1544 : : {
1545 : 12714 : const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1546 [ + + ]: 12714 : if (n < CUTOFF5) mask = 1; else mask = 3;
1547 [ + - ]: 12714 : if (CUTOFF <= ROOT7)
1548 : : {
1549 : 12714 : const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1550 [ - + ]: 12714 : if (n >= CUTOFF7) mask = 7;
1551 : : }
1552 : : }
1553 : :
1554 [ - + ][ # # ]: 12715 : if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1555 [ + + ]: 12715 : if ((e = uis_357_power(n, &n, &mask))) v *= e;
1556 : :
1557 [ + + ]: 12715 : if (uisprime_101(n)) { *pp = n; return v; }
1558 : 541572 : return 0;
1559 : : }
1560 : :
1561 : : /*********************************************************************/
1562 : : /** **/
1563 : : /** KRONECKER SYMBOL **/
1564 : : /** **/
1565 : : /*********************************************************************/
1566 : : /* t = 3,5 mod 8 ? (= 2 not a square mod t) */
1567 : : static int
1568 : 37044452 : ome(long t)
1569 : : {
1570 [ + + ]: 37044452 : switch(t & 7)
1571 : : {
1572 : : case 3:
1573 : 21151538 : case 5: return 1;
1574 : 37044452 : default: return 0;
1575 : : }
1576 : : }
1577 : : /* t a t_INT, is t = 3,5 mod 8 ? */
1578 : : static int
1579 : 2620530 : gome(GEN t)
1580 [ + - ]: 2620530 : { return signe(t)? ome( mod2BIL(t) ): 0; }
1581 : :
1582 : : /* t a t_INT, return 1 if t = 3 (mod 4), 0 otherwise */
1583 : : static int
1584 : 18431 : eps(GEN t)
1585 : : {
1586 [ + + - ]: 18431 : switch(signe(t))
1587 : : {
1588 : 4977 : case -1: return mod4(t) == 1;
1589 : 13454 : case 1: return mod4(t) == 3;
1590 : 18431 : default: return 0;
1591 : : }
1592 : : }
1593 : :
1594 : : /* assume y odd, return kronecker(x,y) * s */
1595 : : static long
1596 : 29543059 : krouu_s(ulong x, ulong y, long s)
1597 : : {
1598 : 29543059 : ulong x1 = x, y1 = y, z;
1599 [ + + ]: 120667926 : while (x1)
1600 : : {
1601 : 91124867 : long r = vals(x1);
1602 [ + + ]: 91124867 : if (r)
1603 : : {
1604 [ + + ][ + + ]: 48705783 : if (odd(r) && ome(y1)) s = -s;
1605 : 48705783 : x1 >>= r;
1606 : : }
1607 [ + + ]: 91124867 : if (x1 & y1 & 2) s = -s;
1608 : 91124867 : z = y1 % x1; y1 = x1; x1 = z;
1609 : : }
1610 [ + + ]: 29543059 : return (y1 == 1)? s: 0;
1611 : : }
1612 : :
1613 : : long
1614 : 4966935 : kronecker(GEN x, GEN y)
1615 : : {
1616 : 4966935 : pari_sp av = avma;
1617 : 4966935 : long s = 1, r;
1618 : : ulong xu, yu;
1619 : :
1620 [ - + ]: 4966935 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
1621 [ - + ]: 4966935 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
1622 [ - - + ]: 4966935 : switch (signe(y))
1623 : : {
1624 [ # # ]: 0 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
1625 : 0 : case 0: return is_pm1(x);
1626 : : }
1627 : 4966935 : r = vali(y);
1628 [ + + ]: 4966935 : if (r)
1629 : : {
1630 [ + + ]: 9317 : if (!mpodd(x)) { avma = av; return 0; }
1631 [ + - ][ + + ]: 9121 : if (odd(r) && gome(x)) s = -s;
1632 : 9121 : y = shifti(y,-r);
1633 : : }
1634 : 4966739 : x = modii(x,y);
1635 [ + + ]: 5031128 : while (lgefint(x) > 3) /* x < y */
1636 : : {
1637 : : GEN z;
1638 : 64389 : r = vali(x);
1639 [ + + ]: 64389 : if (r)
1640 : : {
1641 [ + + ][ + + ]: 35224 : if (odd(r) && gome(y)) s = -s;
1642 : 35224 : x = shifti(x,-r);
1643 : : }
1644 : : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1645 [ + + ]: 64389 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
1646 : 64389 : z = remii(y,x); y = x; x = z;
1647 [ - + ]: 64389 : if (gc_needed(av,2))
1648 : : {
1649 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
1650 : 0 : gerepileall(av, 2, &x, &y);
1651 : : }
1652 : : }
1653 : 4966739 : xu = itou(x);
1654 [ + + ][ + + ]: 4966739 : if (!xu) return is_pm1(y)? s: 0;
1655 : 4947664 : r = vals(xu);
1656 [ + + ]: 4947664 : if (r)
1657 : : {
1658 [ + + ][ + + ]: 3462277 : if (odd(r) && gome(y)) s = -s;
1659 : 3462277 : xu >>= r;
1660 : : }
1661 : : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1662 [ + + ]: 4947664 : if (xu & mod2BIL(y) & 2) s = -s;
1663 : 4947664 : yu = umodiu(y, xu);
1664 : 4966935 : avma = av; return krouu_s(yu, xu, s);
1665 : : }
1666 : :
1667 : : long
1668 : 528 : krois(GEN x, long y)
1669 : : {
1670 : : ulong yu;
1671 : 528 : long s = 1;
1672 : :
1673 [ - + ]: 528 : if (y <= 0)
1674 : : {
1675 [ # # ]: 0 : if (y == 0) return is_pm1(x);
1676 [ # # ]: 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
1677 : : }
1678 : : else
1679 : 528 : yu = (ulong)y;
1680 [ - + ]: 528 : if (!odd(yu))
1681 : : {
1682 : : long r;
1683 [ # # ]: 0 : if (!mpodd(x)) return 0;
1684 : 0 : r = vals(yu); yu >>= r;
1685 [ # # ][ # # ]: 0 : if (odd(r) && gome(x)) s = -s;
1686 : : }
1687 : 528 : return krouu_s(umodiu(x, yu), yu, s);
1688 : : }
1689 : : /* assume y != 0 */
1690 : : long
1691 : 1212815 : kroiu(GEN x, ulong y)
1692 : : {
1693 : : long r;
1694 [ + + ]: 1212815 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
1695 [ + + ]: 10010 : if (!mpodd(x)) return 0;
1696 : 6727 : r = vals(y); y >>= r;
1697 [ + + ][ + + ]: 1212815 : return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
1698 : : }
1699 : :
1700 : : /* assume y > 0, odd, return s * kronecker(x,y) */
1701 : : static long
1702 : 28587 : krouodd(ulong x, GEN y, long s)
1703 : : {
1704 : : long r;
1705 [ + + ]: 28587 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
1706 [ - + ]: 9586 : if (!x) return 0; /* y != 1 */
1707 : 9586 : r = vals(x);
1708 [ + + ]: 9586 : if (r)
1709 : : {
1710 [ + + ][ + + ]: 6251 : if (odd(r) && gome(y)) s = -s;
1711 : 6251 : x >>= r;
1712 : : }
1713 : : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1714 [ + + ]: 9586 : if (x & mod2BIL(y) & 2) s = -s;
1715 : 28587 : return krouu_s(umodiu(y,x), x, s);
1716 : : }
1717 : :
1718 : : long
1719 : 28465 : krosi(long x, GEN y)
1720 : : {
1721 : 28465 : const pari_sp av = avma;
1722 : 28465 : long s = 1, r;
1723 [ - - + ]: 28465 : switch (signe(y))
1724 : : {
1725 [ # # ]: 0 : case -1: y = negi(y); if (x < 0) s = -1; break;
1726 [ # # ][ # # ]: 0 : case 0: return (x==1 || x==-1);
1727 : : }
1728 : 28465 : r = vali(y);
1729 [ - + ]: 28465 : if (r)
1730 : : {
1731 [ # # ]: 0 : if (!odd(x)) { avma = av; return 0; }
1732 [ # # ][ # # ]: 0 : if (odd(r) && ome(x)) s = -s;
1733 : 0 : y = shifti(y,-r);
1734 : : }
1735 [ + + ][ + + ]: 28465 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
1736 : 28465 : s = krouodd((ulong)x, y, s);
1737 : 28465 : avma = av; return s;
1738 : : }
1739 : :
1740 : : long
1741 : 122 : kroui(ulong x, GEN y)
1742 : : {
1743 : 122 : const pari_sp av = avma;
1744 : 122 : long s = 1, r;
1745 [ - - + ]: 122 : switch (signe(y))
1746 : : {
1747 : 0 : case -1: y = negi(y); break;
1748 : 0 : case 0: return x==1UL;
1749 : : }
1750 : 122 : r = vali(y);
1751 [ - + ]: 122 : if (r)
1752 : : {
1753 [ # # ]: 0 : if (!odd(x)) { avma = av; return 0; }
1754 [ # # ][ # # ]: 0 : if (odd(r) && ome(x)) s = -s;
1755 : 0 : y = shifti(y,-r);
1756 : : }
1757 : 122 : s = krouodd(x, y, s);
1758 : 122 : avma = av; return s;
1759 : : }
1760 : :
1761 : : long
1762 : 65632 : kross(long x, long y)
1763 : : {
1764 : : ulong yu;
1765 : 65632 : long s = 1;
1766 : :
1767 [ + + ]: 65632 : if (y <= 0)
1768 : : {
1769 [ - + ]: 385 : if (y == 0) return (labs(x)==1);
1770 [ - + ]: 385 : yu = (ulong)-y; if (x < 0) s = -1;
1771 : : }
1772 : : else
1773 : 65247 : yu = (ulong)y;
1774 [ + + ]: 65632 : if (!odd(yu))
1775 : : {
1776 : : long r;
1777 [ + + ]: 455 : if (!odd(x)) return 0;
1778 : 301 : r = vals(yu); yu >>= r;
1779 [ + + ][ + + ]: 301 : if (odd(r) && ome(x)) s = -s;
1780 : : }
1781 [ + + ]: 65478 : x %= (long)yu; if (x < 0) x += yu;
1782 : 65632 : return krouu_s((ulong)x, yu, s);
1783 : : }
1784 : :
1785 : : long
1786 : 23291270 : krouu(ulong x, ulong y)
1787 : : {
1788 : : long r;
1789 [ + + ]: 23291270 : if (odd(y)) return krouu_s(x, y, 1);
1790 [ - + ]: 462519 : if (!odd(x)) return 0;
1791 : 462519 : r = vals(y); y >>= r;
1792 [ + + ][ - + ]: 23291270 : return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
1793 : : }
1794 : :
1795 : : /*********************************************************************/
1796 : : /** **/
1797 : : /** HILBERT SYMBOL **/
1798 : : /** **/
1799 : : /*********************************************************************/
1800 : : /* x,y are t_INT or t_REAL */
1801 : : static long
1802 : 9968 : mphilbertoo(GEN x, GEN y)
1803 : : {
1804 : 9968 : long sx = signe(x), sy = signe(y);
1805 [ + - ][ - + ]: 9968 : if (!sx || !sy) return 0;
1806 [ + + ][ + + ]: 9968 : return (sx < 0 && sy < 0)? -1: 1;
1807 : : }
1808 : :
1809 : : long
1810 : 53095 : hilbertii(GEN x, GEN y, GEN p)
1811 : : {
1812 : : pari_sp av;
1813 : : long oddvx, oddvy, z;
1814 : :
1815 [ + + ]: 53095 : if (!p) return mphilbertoo(x,y);
1816 [ + - ][ - + ]: 43148 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
1817 [ + + ][ + + ]: 43148 : if (!signe(x) || !signe(y)) return 0;
1818 : 43127 : av = avma;
1819 : 43127 : oddvx = odd(Z_pvalrem(x,p,&x));
1820 : 43127 : oddvy = odd(Z_pvalrem(y,p,&y));
1821 : : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
1822 [ + + ]: 43127 : if (equaliu(p, 2))
1823 : : {
1824 [ + + ][ + + ]: 10675 : z = (eps(x) && eps(y))? -1: 1;
1825 [ + + ][ + + ]: 10675 : if (oddvx && gome(y)) z = -z;
1826 [ + + ][ + + ]: 10675 : if (oddvy && gome(x)) z = -z;
1827 : : }
1828 : : else
1829 : : {
1830 [ + + ][ + + ]: 32452 : z = (oddvx && oddvy && eps(p))? -1: 1;
[ + + ]
1831 [ + + ][ + + ]: 32452 : if (oddvx && kronecker(y,p) < 0) z = -z;
1832 [ + + ][ + + ]: 32452 : if (oddvy && kronecker(x,p) < 0) z = -z;
1833 : : }
1834 : 53095 : avma = av; return z;
1835 : : }
1836 : :
1837 : : static void
1838 : 196 : err_prec(void) { pari_err_PREC("hilbert"); }
1839 : : static void
1840 : 161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
1841 : : static void
1842 : 56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
1843 : :
1844 : : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
1845 : : * Return lift(x) provided it's p-adic accuracy is large enough to decide
1846 : : * hilbert()'s value [ problem at p = 2 ] */
1847 : : static GEN
1848 : 420 : lift_intmod(GEN x, GEN *pp)
1849 : : {
1850 : 420 : GEN p = *pp, N = gel(x,1);
1851 : 420 : x = gel(x,2);
1852 [ + + ]: 420 : if (!p)
1853 : : {
1854 : 266 : *pp = p = N;
1855 [ + + ]: 266 : switch(itos_or_0(p))
1856 : : {
1857 : : case 2:
1858 : 126 : case 4: err_prec();
1859 : : }
1860 : 140 : return x;
1861 : : }
1862 [ + + ]: 154 : if (!signe(p)) err_oo(N);
1863 [ + + ]: 112 : if (equaliu(p,2))
1864 [ + + ]: 42 : { if (vali(N) <= 2) err_prec(); }
1865 : : else
1866 [ + + ]: 70 : { if (!dvdii(N,p)) err_p(N,p); }
1867 [ + + ]: 28 : if (!signe(x)) err_prec();
1868 : 161 : return x;
1869 : : }
1870 : : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
1871 : : * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
1872 : : * to decide hilbert()'s value [ problem at p = 2 ]*/
1873 : : static GEN
1874 : 210 : lift_padic(GEN x, GEN *pp)
1875 : : {
1876 : 210 : GEN p = *pp, q = gel(x,2), y = gel(x,4);
1877 [ + + ]: 210 : if (!p) *pp = p = q;
1878 [ + + ]: 147 : else if (!equalii(p,q)) err_p(p, q);
1879 [ + + ][ + + ]: 105 : if (equaliu(p,2) && precp(x) <= 2) err_prec();
1880 [ - + ]: 70 : if (!signe(y)) err_prec();
1881 [ - + ]: 70 : return odd(valp(x))? mulii(p,y): y;
1882 : : }
1883 : :
1884 : : long
1885 : 658 : hilbert(GEN x, GEN y, GEN p)
1886 : : {
1887 : 658 : pari_sp av = avma;
1888 : 658 : long tx = typ(x), ty = typ(y), z;
1889 : :
1890 [ + + ][ - + ]: 658 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
1891 [ + + ]: 658 : if (tx == t_REAL)
1892 : : {
1893 [ + + ][ + + ]: 77 : if (p && signe(p)) err_oo(p);
1894 [ + - + ]: 63 : switch (ty)
1895 : : {
1896 : : case t_INT:
1897 : 7 : case t_REAL: return mphilbertoo(x,y);
1898 : 0 : case t_FRAC: return mphilbertoo(x,gel(y,1));
1899 : 56 : default: pari_err_TYPE2("hilbert",x,y);
1900 : : }
1901 : : }
1902 [ + + ]: 581 : if (ty == t_REAL)
1903 : : {
1904 [ - + ][ # # ]: 14 : if (p && signe(p)) err_oo(p);
1905 [ + - - ]: 14 : switch (tx)
1906 : : {
1907 : : case t_INT:
1908 : 14 : case t_REAL: return mphilbertoo(x,y);
1909 : 0 : case t_FRAC: return mphilbertoo(gel(x,1),y);
1910 : 0 : default: pari_err_TYPE2("hilbert",x,y);
1911 : : }
1912 : : }
1913 [ + + ]: 567 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
1914 [ + + ]: 364 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
1915 : :
1916 [ + + ]: 308 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
1917 [ + + ]: 245 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
1918 : :
1919 [ + + ][ + - ]: 168 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
1920 [ + + ][ + - ]: 168 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
1921 : :
1922 [ + - ][ - + ]: 168 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
1923 [ + + ][ + + ]: 168 : if (p && !signe(p)) p = NULL;
1924 : 189 : z = hilbertii(x,y,p); avma = av; return z;
1925 : : }
1926 : :
1927 : : /*******************************************************************/
1928 : : /* */
1929 : : /* SQUARE ROOT MODULO p */
1930 : : /* */
1931 : : /*******************************************************************/
1932 : :
1933 : : static ulong
1934 : 3776551 : Fl_2gener_pre_all(long e, ulong p, ulong pi, ulong *pt_m)
1935 : : {
1936 : : ulong y, m;
1937 : : long k, i;
1938 : 3776551 : ulong q = (p-1) >> e; /* q = (p-1)/2^oo is odd */
1939 : 3776551 : for (k=2; ; k++)
1940 : : { /* loop terminates for k < p (even if p composite) */
1941 : 8964590 : i = krouu(k, p);
1942 [ + + ]: 8964590 : if (i >= 0)
1943 : : {
1944 [ + + ]: 5188046 : if (i) continue;
1945 : 7 : pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1946 : : }
1947 : 3776544 : y = m = Fl_powu_pre(k, q, p, pi);
1948 [ + + ]: 11265412 : for (i=1; i<e; i++)
1949 [ - + ]: 7488868 : if ((m = Fl_sqr_pre(m, p, pi)) == 1) break;
1950 [ + - ]: 3776544 : if (i == e) break; /* success */
1951 : 5188039 : }
1952 : 3776544 : *pt_m = m;
1953 : 3776544 : return y;
1954 : : }
1955 : :
1956 : : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
1957 : : static ulong
1958 : 7874567 : Fl_sqrt_i(ulong a, ulong p, ulong pi, ulong y, ulong m)
1959 : : {
1960 : : long i, e, k;
1961 : : ulong p1, q, v, w;
1962 : :
1963 [ + + ]: 7874567 : if (!a) return 0;
1964 : 7808114 : p1 = p - 1; e = vals(p1);
1965 [ + + ]: 7808114 : if (e == 0) /* p = 2 */
1966 : : {
1967 [ + + ]: 7133 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1968 : 7126 : return ((a & 1) == 0)? 0: 1;
1969 : : }
1970 : 7800981 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
1971 [ + + ]: 7800981 : if (e == 1) y = p1;
1972 [ + - ]: 3776551 : else if (y==0) y = Fl_2gener_pre_all(e, p, pi, &m);
1973 : 7800974 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
1974 [ - + ]: 7800974 : if (!p1) return 0;
1975 : 7800974 : v = Fl_mul_pre(a, p1, p, pi);
1976 : 7800974 : w = Fl_mul_pre(v, p1, p, pi);
1977 [ + + ]: 10903007 : while (w != 1)
1978 : : { /* a*w = v^2, y primitive 2^e-th root of 1
1979 : : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1980 : 3113375 : p1 = Fl_sqr_pre(w,p,pi);
1981 [ + + ][ + + ]: 6159281 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
1982 [ + + ]: 3113375 : if (k == e) return ~0UL;
1983 : : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1984 : 3102033 : p1 = y;
1985 [ + + ]: 4566873 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
1986 : 3102033 : y = Fl_sqr_pre(p1, p, pi); e = k;
1987 : 3102033 : w = Fl_mul_pre(y, w, p, pi);
1988 : 3102033 : v = Fl_mul_pre(v, p1, p, pi);
1989 : : }
1990 [ + + ]: 7789632 : p1 = p - v; if (v > p1) v = p1;
1991 : 7874553 : return v;
1992 : : }
1993 : :
1994 : : ulong
1995 : 7002148 : Fl_sqrt(ulong a, ulong p)
1996 : : {
1997 : 7002148 : ulong pi = get_Fl_red(p);
1998 : 7002148 : return Fl_sqrt_i(a, p, pi, 0, 0);
1999 : : }
2000 : :
2001 : : ulong
2002 : 872419 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
2003 : : {
2004 : 872419 : return Fl_sqrt_i(a, p, pi, 0, 0);
2005 : : }
2006 : :
2007 : : static ulong
2008 : 8239 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
2009 : : {
2010 : : ulong x, y, m;
2011 : 8239 : ulong le1 = upowuu(l, e-1);
2012 : 8239 : for (x = 2; ; x++)
2013 : : {
2014 : 14582 : y = Fl_powu_pre(x, r, p, pi);
2015 [ + + ]: 14582 : if (y==1) continue;
2016 : 10459 : m = Fl_powu_pre(y, le1, p, pi);
2017 [ + + ]: 10459 : if (m != 1) break;
2018 : 6343 : }
2019 : 8239 : *pt_m = m;
2020 : 8239 : return y;
2021 : : }
2022 : :
2023 : : /* solve x^l = a , l prime in G of order q.
2024 : : *
2025 : : * q = (l^e)*r, e >= 1, (r,l) = 1
2026 : : * y generates the l-Sylow of G
2027 : : * m = y^(l^(e-1)) != 1 */
2028 : : static ulong
2029 : 16080 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
2030 : : {
2031 : : ulong p1, v, w, z, dl, zm;
2032 : : ulong r, e, u2;
2033 [ + + ]: 16080 : if (a==0) return a;
2034 : 16074 : e = u_lvalrem(p-1, l, &r);
2035 : 16074 : u2 = Fl_inv(l%r, r);
2036 : 16074 : v = Fl_powu_pre(a, u2, p,pi);
2037 : 16074 : w = Fl_powu_pre(v, l, p,pi);
2038 : 16074 : w = Fl_mul_pre(w, Fl_inv(a, p),p,pi);
2039 [ + + ]: 16074 : if (w==1) return v;
2040 [ + - ]: 8239 : if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
2041 [ + + ]: 11124 : while (w!=1)
2042 : : {
2043 : 8999 : ulong k = 0;
2044 : 8999 : p1 = w;
2045 : : do
2046 : : {
2047 : 12942 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
2048 : 12942 : k++;
2049 [ + + ]: 12942 : } while (p1!=1);
2050 [ + + ]: 8999 : if (k==e) return ~0UL;
2051 : 2885 : dl = 0; zm = 1;
2052 [ + + ]: 7191 : while (z!=zm)
2053 : : {
2054 : 4306 : zm = Fl_mul_pre(zm, m, p, pi); dl++;
2055 : : }
2056 : 2885 : dl = Fl_neg(dl, l);
2057 : 2885 : p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
2058 : 2885 : m = Fl_powu_pre(m, dl, p, pi);
2059 : 2885 : e = k;
2060 : 2885 : v = Fl_mul_pre(p1,v,p,pi);
2061 : 2885 : y = Fl_powu_pre(p1,l,p,pi);
2062 : 2885 : w = Fl_mul_pre(y,w,p,pi);
2063 : : }
2064 : 16080 : return v;
2065 : : }
2066 : :
2067 : : ulong
2068 : 16080 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
2069 : : {
2070 : 16080 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2071 : : }
2072 : :
2073 : : ulong
2074 : 0 : Fl_sqrtl(ulong a, ulong l, ulong p)
2075 : : {
2076 : 0 : ulong pi = get_Fl_red(p);
2077 : 0 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2078 : : }
2079 : :
2080 : : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
2081 : : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
2082 : : * and in average is about 2 or 2.5 times worse. But try both algorithms for
2083 : : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
2084 : : *
2085 : : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
2086 : : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
2087 : : * If (a|p)=1, then sqrt(a) is in F_p.
2088 : : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
2089 : :
2090 : : /* compute y^2, y = y[1] + y[2] X */
2091 : : static GEN
2092 : 449 : sqrt_Cipolla_sqr(void *data, GEN y)
2093 : : {
2094 : 449 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
2095 : 449 : GEN u2 = sqri(u), v2 = sqri(v);
2096 : 449 : v = subii(sqri(addii(v,u)), addii(u2,v2));
2097 : 449 : u = addii(u2, mulii(v2,n));
2098 : : /* NOT mkvec2: must be gerepileupto-able */
2099 : 449 : retmkvec2(modii(u,p), modii(v,p));
2100 : : }
2101 : : /* compute (t+X) y^2 */
2102 : : static GEN
2103 : 23 : sqrt_Cipolla_msqr(void *data, GEN y)
2104 : : {
2105 : 23 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
2106 : 23 : ulong t = gt[2];
2107 : 23 : GEN d = addii(u, mului(t,v)), d2= sqri(d);
2108 : 23 : GEN b = remii(mulii(a,v), p);
2109 : 23 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
2110 : 23 : v = subii(d2, mulii(b,v));
2111 : : /* NOT mkvec2: must be gerepileupto-able */
2112 : 23 : retmkvec2(modii(u,p), modii(v,p));
2113 : : }
2114 : : /* assume a reduced mod p [ otherwise correct but inefficient ] */
2115 : : static GEN
2116 : 8 : sqrt_Cipolla(GEN a, GEN p)
2117 : : {
2118 : : pari_sp av1;
2119 : : GEN u, v, n, y, pov2;
2120 : : ulong t;
2121 : :
2122 [ - + ]: 8 : if (kronecker(a, p) < 0) return NULL;
2123 : 8 : pov2 = shifti(p,-1);
2124 [ + - ]: 8 : if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
2125 : :
2126 : 8 : av1 = avma;
2127 : 8 : for(t=1; ; t++)
2128 : : {
2129 : 41 : n = subsi((long)(t*t), a);
2130 [ + + ]: 41 : if (kronecker(n, p) < 0) break;
2131 : 33 : avma = av1;
2132 : 33 : }
2133 : :
2134 : : /* compute (t+X)^((p-1)/2) =: u+vX */
2135 : 8 : u = utoipos(t);
2136 : 8 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
2137 : : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
2138 : : /* Now u+vX = (t+X)^((p-1)/2); thus
2139 : : * (u+vX)(t+X) = sqrt(a) + 0 X
2140 : : * Whence,
2141 : : * sqrt(a) = (u+vt)t - v*a
2142 : : * 0 = (u+vt)
2143 : : * Thus a square root is v*a */
2144 : :
2145 : 8 : v = Fp_mul(gel(y, 2), a, p);
2146 [ + + ]: 8 : if (cmpii(v,pov2) > 0) v = subii(p,v);
2147 : 8 : return v;
2148 : : }
2149 : :
2150 : : #define sqrmod(x,p) (remii(sqri(x),p))
2151 : :
2152 : : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
2153 : : GEN
2154 : 1625447 : Fp_sqrt(GEN a, GEN p)
2155 : : {
2156 : 1625447 : pari_sp av = avma, av1;
2157 : : long i, k, e;
2158 : : GEN p1, q, v, y, w, m;
2159 : :
2160 [ - + ]: 1625447 : if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
2161 [ - + ]: 1625447 : if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
2162 [ + - ][ - + ]: 1625447 : if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
2163 [ + + ]: 1625447 : if (lgefint(p) == 3)
2164 : : {
2165 : 1619090 : ulong pp = uel(p,2), u = Fl_sqrt(umodiu(a, pp), pp);
2166 [ + + ]: 1619076 : if (u == ~0UL) return NULL;
2167 : 1619034 : return utoi(u);
2168 : : }
2169 : :
2170 : 6357 : p1 = addsi(-1,p); e = vali(p1);
2171 : 6357 : a = modii(a, p);
2172 : :
2173 : : /* On average, the algorithm of Cipolla is better than the algorithm of
2174 : : * Tonelli and Shanks if and only if e(e-1)>8*log2(n)+20
2175 : : * see LNCS 2286 pp 430 [GTL] */
2176 [ + + ]: 6357 : if (e*(e-1) > 20 + 8 * expi(p))
2177 : : {
2178 : 8 : v = sqrt_Cipolla(a,p);
2179 [ - + ]: 8 : if (!v) { avma = av; return NULL; }
2180 : 8 : return gerepileuptoint(av,v);
2181 : : }
2182 : :
2183 [ - + ]: 6349 : if (e == 0) /* p = 2 */
2184 : : {
2185 : 0 : avma = av;
2186 [ # # ]: 0 : if (!equaliu(p,2)) pari_err_PRIME("Fp_sqrt [modulus]",p);
2187 [ # # ][ # # ]: 0 : if (!signe(a) || !mod2(a)) return gen_0;
2188 : 0 : return gen_1;
2189 : : }
2190 : 6349 : q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
2191 [ + + ]: 6349 : if (e == 1) y = p1;
2192 : : else /* look for an odd power of a primitive root */
2193 : 3264 : for (k=2; ; k++)
2194 : : { /* loop terminates for k < p (even if p composite) */
2195 : :
2196 : 6058 : i = krosi(k,p);
2197 [ + + ]: 6058 : if (i >= 0)
2198 : : {
2199 [ + - ]: 2794 : if (i) continue;
2200 : 0 : pari_err_PRIME("Fp_sqrt [modulus]",p);
2201 : : }
2202 : 3264 : av1 = avma;
2203 : 3264 : y = m = Fp_pow(utoipos((ulong)k),q,p);
2204 [ + + ]: 8617 : for (i=1; i<e; i++)
2205 [ - + ]: 5353 : if (gequal1(m = sqrmod(m,p))) break;
2206 [ + - ]: 3264 : if (i == e) break; /* success */
2207 : 0 : avma = av1;
2208 : 2794 : }
2209 : :
2210 : 6349 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ [(q-1)/2] */
2211 [ + + ]: 6349 : if (!signe(p1)) { avma=av; return gen_0; }
2212 : 6342 : v = Fp_mul(a, p1, p);
2213 : 6342 : w = Fp_mul(v, p1, p);
2214 [ + + ]: 9014 : while (!equali1(w))
2215 : : { /* a*w = v^2, y primitive 2^e-th root of 1
2216 : : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
2217 : 2675 : p1 = sqrmod(w,p);
2218 [ + + ][ + + ]: 4734 : for (k=1; !equali1(p1) && k < e; k++) p1 = sqrmod(p1,p);
2219 [ + + ]: 2675 : if (k == e) { avma=av; return NULL; } /* p composite or (a/p) != 1 */
2220 : : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
2221 : 2672 : p1 = y;
2222 [ + + ]: 3410 : for (i=1; i < e-k; i++) p1 = sqrmod(p1,p);
2223 : 2672 : y = sqrmod(p1, p); e = k;
2224 : 2672 : w = Fp_mul(y, w, p);
2225 : 2672 : v = Fp_mul(v, p1, p);
2226 [ - + ]: 2672 : if (gc_needed(av,1))
2227 : : {
2228 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
2229 : 0 : gerepileall(av,3, &y,&w,&v);
2230 : : }
2231 : : }
2232 : 6339 : av1 = avma;
2233 [ + + ]: 6339 : p1 = subii(p,v); if (cmpii(v,p1) > 0) v = p1; else avma = av1;
2234 : 1625433 : return gerepileuptoint(av, v);
2235 : : }
2236 : :
2237 : : /*********************************************************************/
2238 : : /** **/
2239 : : /** GCD & BEZOUT **/
2240 : : /** **/
2241 : : /*********************************************************************/
2242 : :
2243 : : GEN
2244 : 2694591 : lcmii(GEN x, GEN y)
2245 : : {
2246 : : pari_sp av;
2247 : : GEN a, b;
2248 [ + - ][ - + ]: 2694591 : if (!signe(x) || !signe(y)) return gen_0;
2249 : 2694591 : av = avma;
2250 [ + + ]: 2694591 : a = gcdii(x,y); if (!equali1(a)) y = diviiexact(y,a);
2251 : 2694591 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
2252 : : }
2253 : :
2254 : : /*********************************************************************/
2255 : : /** **/
2256 : : /** CHINESE REMAINDERS **/
2257 : : /** **/
2258 : : /*********************************************************************/
2259 : :
2260 : : /* P.M. & M.H.
2261 : : *
2262 : : * Chinese Remainder Theorem. x and y must have the same type (integermod,
2263 : : * polymod, or polynomial/vector/matrix recursively constructed with these
2264 : : * as coefficients). Creates (with the same type) a z in the same residue
2265 : : * class as x and the same residue class as y, if it is possible.
2266 : : *
2267 : : * We also allow (during recursion) two identical objects even if they are
2268 : : * not integermod or polymod. For example, if
2269 : : *
2270 : : * x = [1. mod(5, 11), mod(X + mod(2, 7), X^2 + 1)]
2271 : : * y = [1, mod(7, 17), mod(X + mod(0, 3), X^2 + 1)],
2272 : : *
2273 : : * then chinese(x, y) returns
2274 : : *
2275 : : * [1, mod(16, 187), mod(X + mod(9, 21), X^2 + 1)]
2276 : : *
2277 : : * Someone else may want to allow power series, complex numbers, and
2278 : : * quadratic numbers.
2279 : : */
2280 : :
2281 : : GEN
2282 : 42 : chinese1(GEN x) { return gassoc_proto(chinese,x,NULL); }
2283 : :
2284 : : GEN
2285 : 16443 : chinese(GEN x, GEN y)
2286 : : {
2287 : : pari_sp av,tetpil;
2288 : 16443 : long tx = typ(x);
2289 : : GEN z,p1,p2,d,u,v;
2290 : :
2291 [ + + ]: 16443 : if (!y) return chinese1(x);
2292 [ + + ]: 16401 : if (gequal(x,y)) return gcopy(x);
2293 [ + - ][ + + : 16394 : if (tx == typ(y)) switch(tx)
+ + - ]
2294 : : {
2295 : : case t_POLMOD:
2296 : : {
2297 : 7 : GEN A = gel(x,1), B = gel(y,1);
2298 : 7 : GEN a = gel(x,2), b = gel(y,2);
2299 : 7 : z = cgetg(3, t_POLMOD);
2300 [ - + ]: 7 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
2301 [ - + ]: 7 : if (RgX_equal(A,B)) /* same modulus */
2302 : : {
2303 : 0 : gel(z,1) = gcopy(A);
2304 : 0 : gel(z,2) = chinese(a,b);
2305 : 0 : return z;
2306 : : }
2307 : 7 : av = avma;
2308 : 7 : d = RgX_extgcd(A,B,&u,&v);
2309 : 7 : p2 = gsub(b, a);
2310 [ - + ]: 7 : if (!gequal0(gmod(p2, d))) break;
2311 : 7 : p1 = gdiv(A,d);
2312 : 7 : p2 = gadd(a, gmul(gmul(u,p1), p2));
2313 : :
2314 : 7 : tetpil = avma;
2315 : 7 : gel(z,1) = gmul(p1,B);
2316 : 7 : gel(z,2) = gmod(p2,gel(z,1));
2317 : 7 : gerepilecoeffssp(av,tetpil,z+1,2); return z;
2318 : : }
2319 : : case t_INTMOD:
2320 : : {
2321 : 16373 : GEN A = gel(x,1), B = gel(y,1);
2322 : 16373 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
2323 : 16373 : z = cgetg(3,t_INTMOD);
2324 : 16373 : Z_chinese_pre(A, B, &C, &U, &d);
2325 : 16373 : c = Z_chinese_post(a, b, C, U, d);
2326 [ - + ]: 16373 : if (!c) pari_err_OP("chinese", x,y);
2327 : 16373 : gel(z,1) = icopy_avma(C, (pari_sp)z);
2328 : 16373 : gel(z,2) = icopy_avma(c, (pari_sp)gel(z,1));
2329 : 16373 : avma = (pari_sp)gel(z,2); return z;
2330 : : }
2331 : : case t_POL:
2332 : : {
2333 : 7 : long i, lx = lg(x), ly = lg(y);
2334 [ - + ]: 7 : if (varn(x) != varn(y)) break;
2335 [ + - ]: 7 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
2336 : 7 : z = cgetg(lx, t_POL); z[1] = x[1];
2337 [ + + ]: 21 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2338 [ + + ]: 14 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
2339 : 7 : return z;
2340 : : }
2341 : :
2342 : : case t_VEC: case t_COL: case t_MAT:
2343 : : {
2344 : : long i, lx;
2345 [ + - ]: 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
2346 [ + + ]: 21 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2347 : 7 : return z;
2348 : : }
2349 : : }
2350 : 0 : pari_err_OP("chinese",x,y);
2351 : 16436 : return NULL; /* not reached */
2352 : : }
2353 : :
2354 : : /* init chinese(Mod(.,A), Mod(.,B)) */
2355 : : void
2356 : 126577 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
2357 : : {
2358 : 126577 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
2359 : 126577 : GEN t = diviiexact(A,d);
2360 : 126577 : *pU = mulii(u, t);
2361 : 126577 : *pC = mulii(t, B);
2362 [ + + ]: 126577 : if (pd) *pd = d;
2363 : 126577 : }
2364 : : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
2365 : : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
2366 : : * If d not NULL, check wether a = b mod d. */
2367 : : GEN
2368 : 616557 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
2369 : : {
2370 : : GEN b_a;
2371 [ + + ]: 616557 : if (!signe(a))
2372 : : {
2373 [ + + ][ - + ]: 133735 : if (d && remii(b, d) != gen_0) return NULL;
2374 : 133735 : return Fp_mul(b, U, C);
2375 : : }
2376 : 482822 : b_a = subii(b,a);
2377 [ + + ][ - + ]: 482822 : if (d && remii(b_a, d) != gen_0) return NULL;
2378 : 616557 : return modii(addii(a, mulii(U, b_a)), C);
2379 : : }
2380 : : GEN
2381 : 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
2382 : : {
2383 : 2142 : pari_sp av = avma;
2384 : 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
2385 : 2142 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
2386 : : }
2387 : : GEN
2388 : 108055 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
2389 : : {
2390 : 108055 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
2391 : 108055 : return Z_chinese_post(a,b, *pC, U, NULL);
2392 : : }
2393 : :
2394 : : /* return lift(chinese(a mod A, b mod B))
2395 : : * assume(A,B)=1, a,b,A,B integers and C = A*B */
2396 : : GEN
2397 : 424305 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
2398 : : {
2399 : 424305 : pari_sp av = avma;
2400 : 424305 : GEN U = mulii(Fp_inv(A,B), A);
2401 : 424305 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2402 : : }
2403 : :
2404 : : /* chinese1 for coprime moduli in Z */
2405 : : static GEN
2406 : 65668 : chinese1_coprime_Z_aux(GEN x, GEN y)
2407 : : {
2408 : 65668 : GEN z = cgetg(3, t_INTMOD);
2409 : 65668 : GEN A = gel(x,1), a = gel(x, 2);
2410 : 65668 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
2411 : 65668 : pari_sp av = avma;
2412 : 65668 : GEN U = mulii(Fp_inv(A,B), A);
2413 : 65668 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2414 : 65668 : gel(z,1) = C; return z;
2415 : : }
2416 : : GEN
2417 : 75685 : chinese1_coprime_Z(GEN x) {return gassoc_proto(chinese1_coprime_Z_aux,x,NULL);}
2418 : :
2419 : : /*********************************************************************/
2420 : : /** **/
2421 : : /** MODULAR EXPONENTIATION **/
2422 : : /** **/
2423 : : /*********************************************************************/
2424 : :
2425 : : /* xa, ya = t_VECSMALL */
2426 : : GEN
2427 : 145100 : ZV_producttree(GEN xa)
2428 : : {
2429 : 145100 : long n = lg(xa)-1;
2430 [ + - ]: 145100 : long m = n==1 ? 1: expu(n-1)+1;
2431 : 145100 : GEN T = cgetg(m+1, t_VEC), t;
2432 : : long i, j, k;
2433 : 145100 : t = cgetg(((n+1)>>1)+1, t_VEC);
2434 [ + + ]: 145100 : if (typ(xa)==t_VECSMALL)
2435 : : {
2436 [ + + ]: 116850 : for (j=1, k=1; k<n; j++, k+=2)
2437 : 70962 : gel(t, j) = muluu(xa[k], xa[k+1]);
2438 [ + + ]: 45888 : if (k==n) gel(t, j) = utoi(xa[k]);
2439 : : } else {
2440 [ + + ]: 267620 : for (j=1, k=1; k<n; j++, k+=2)
2441 : 168408 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
2442 [ + + ]: 99212 : if (k==n) gel(t, j) = icopy(gel(xa,k));
2443 : : }
2444 : 145100 : gel(T,1) = t;
2445 [ + + ]: 249912 : for (i=2; i<=m; i++)
2446 : : {
2447 : 104812 : GEN u = gel(T, i-1);
2448 : 104812 : long n = lg(u)-1;
2449 : 104812 : t = cgetg(((n+1)>>1)+1, t_VEC);
2450 [ + + ]: 245879 : for (j=1, k=1; k<n; j++, k+=2)
2451 : 141067 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
2452 [ + + ]: 104812 : if (k==n) gel(t, j) = gel(u, k);
2453 : 104812 : gel(T, i) = t;
2454 : : }
2455 : 145100 : return T;
2456 : : }
2457 : :
2458 : : static GEN
2459 : 765652 : Z_ZV_mod_tree(GEN P, GEN xa, GEN T)
2460 : : {
2461 : : long i,j,k;
2462 : 765652 : long m = lg(T)-1, n = lg(xa)-1;
2463 : : GEN t;
2464 : 765652 : GEN Tp = cgetg(m+1, t_VEC);
2465 : 765652 : gel(Tp, m) = mkvec(P);
2466 [ + + ]: 1300738 : for (i=m-1; i>=1; i--)
2467 : : {
2468 : 535086 : GEN u = gel(T, i);
2469 : 535086 : GEN v = gel(Tp, i+1);
2470 : 535086 : long n = lg(u)-1;
2471 : 535086 : t = cgetg(n+1, t_VEC);
2472 [ + + ]: 1247997 : for (j=1, k=1; k<n; j++, k+=2)
2473 : : {
2474 : 712911 : gel(t, k) = modii(gel(v, j), gel(u, k));
2475 : 712911 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
2476 : : }
2477 [ + + ]: 535086 : if (k==n) gel(t, k) = gel(v, j);
2478 : 535086 : gel(Tp, i) = t;
2479 : : }
2480 : : {
2481 : 765652 : GEN u = gel(T, i+1);
2482 : 765652 : GEN v = gel(Tp, i+1);
2483 : 765652 : long l = lg(u)-1;
2484 [ + + ]: 765652 : if (typ(xa)==t_VECSMALL)
2485 : : {
2486 : 620552 : GEN R = cgetg(n+1, t_VECSMALL);
2487 [ + + ]: 1812948 : for (j=1, k=1; j<=l; j++, k+=2)
2488 : : {
2489 : 1192396 : uel(R,k) = umodiu(gel(v, j), xa[k]);
2490 [ + + ]: 1192396 : if (k < n)
2491 : 1003898 : uel(R,k+1) = umodiu(gel(v, j), xa[k+1]);
2492 : : }
2493 : 620552 : return R;
2494 : : }
2495 : : else
2496 : : {
2497 : 145100 : GEN R = cgetg(n+1, t_VEC);
2498 [ + + ]: 431267 : for (j=1, k=1; j<=l; j++, k+=2)
2499 : : {
2500 : 286167 : gel(R,k) = modii(gel(v, j), gel(xa,k));
2501 [ + + ]: 286167 : if (k < n)
2502 : 239370 : gel(R,k+1) = modii(gel(v, j), gel(xa,k+1));
2503 : : }
2504 : 765652 : return R;
2505 : : }
2506 : : }
2507 : : }
2508 : :
2509 : : static GEN
2510 : 145100 : ZV_polint_tree(GEN T, GEN R, GEN xa, GEN ya)
2511 : : {
2512 : 145100 : long m = lg(T)-1, n = lg(ya)-1;
2513 : : long i,j,k;
2514 : 145100 : GEN Tp = cgetg(m+1, t_VEC);
2515 : 145100 : GEN M = gel(T, 1);
2516 : 145100 : GEN t = cgetg(lg(M), t_VEC);
2517 [ + + ]: 145100 : if (typ(xa)==t_VECSMALL)
2518 : : {
2519 [ + + ]: 116850 : for (j=1, k=1; k<n; j++, k+=2)
2520 : : {
2521 : 70962 : pari_sp av = avma;
2522 : 70962 : GEN a = mului(ya[k], gel(R,k)), b = mului(ya[k+1], gel(R,k+1));
2523 : 70962 : GEN tj = modii(addii(mului(xa[k],b), mului(xa[k+1],a)), gel(M,j));
2524 : 70962 : gel(t, j) = gerepileuptoint(av, tj);
2525 : : }
2526 [ + + ]: 45888 : if (k==n) gel(t, j) = modii(mului(ya[k], gel(R,k)), gel(M, j));
2527 : : } else
2528 : : {
2529 [ + + ]: 267620 : for (j=1, k=1; k<n; j++, k+=2)
2530 : : {
2531 : 168408 : pari_sp av = avma;
2532 : 168408 : GEN a = mulii(gel(ya,k), gel(R,k)), b = mulii(gel(ya,k+1), gel(R,k+1));
2533 : 168408 : GEN tj = modii(addii(mulii(gel(xa,k),b), mulii(gel(xa,k+1),a)), gel(M,j));
2534 : 168408 : gel(t, j) = gerepileuptoint(av, tj);
2535 : : }
2536 [ + + ]: 99212 : if (k==n) gel(t, j) = modii(mulii(gel(ya,k), gel(R,k)), gel(M, j));
2537 : : }
2538 : 145100 : gel(Tp, 1) = t;
2539 [ + + ]: 249912 : for (i=2; i<=m; i++)
2540 : : {
2541 : 104812 : GEN u = gel(T, i-1), M = gel(T, i);
2542 : 104812 : GEN t = cgetg(lg(M), t_VEC);
2543 : 104812 : GEN v = gel(Tp, i-1);
2544 : 104812 : long n = lg(v)-1;
2545 [ + + ]: 245879 : for (j=1, k=1; k<n; j++, k+=2)
2546 : : {
2547 : 141067 : pari_sp av = avma;
2548 : 141067 : gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
2549 : 423201 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
2550 : : }
2551 [ + + ]: 104812 : if (k==n) gel(t, j) = gel(v, k);
2552 : 104812 : gel(Tp, i) = t;
2553 : : }
2554 : 145100 : return gmael(Tp,m,1);
2555 : : }
2556 : :
2557 : : static GEN
2558 : 0 : ncV_polint_tree(GEN T, GEN R, GEN xa, GEN Va)
2559 : : {
2560 : 0 : long i, j, l = lg(gel(Va,1)), n = lg(xa);
2561 : 0 : GEN V = cgetg(l, t_COL);
2562 [ # # ]: 0 : for(i=1; i < l; i++)
2563 : : {
2564 : 0 : pari_sp av = avma;
2565 : 0 : GEN ya = cgetg(n, t_VECSMALL);
2566 [ # # ]: 0 : for(j=1; j < n; j++)
2567 : 0 : ya[j] = mael(Va,j,i);
2568 : 0 : gel(V,i) = gerepilecopy(av, ZV_polint_tree(T, R, xa, ya));
2569 : : }
2570 : 0 : return V;
2571 : : }
2572 : :
2573 : : static GEN
2574 : 0 : nmV_polint_tree(GEN T, GEN R, GEN xa, GEN Ma)
2575 : : {
2576 : 0 : long i, j, l = lg(gel(Ma,1)), n = lg(xa);
2577 : 0 : GEN ya = cgetg(n, t_VEC);
2578 : 0 : GEN M = cgetg(l, t_MAT);
2579 [ # # ]: 0 : for(i=1; i < l; i++)
2580 : : {
2581 [ # # ]: 0 : for(j=1; j < n; j++)
2582 : 0 : gel(ya,j) = gmael(Ma,j,i);
2583 : 0 : gel(M,i) = ncV_polint_tree(T, R, xa, ya);
2584 : : }
2585 : 0 : return M;
2586 : : }
2587 : : GEN
2588 : 0 : Z_ZV_mod(GEN P, GEN xa)
2589 : : {
2590 : 0 : pari_sp av = avma;
2591 : 0 : GEN T = ZV_producttree(xa);
2592 : 0 : return gerepileuptoleaf(av, Z_ZV_mod_tree(P, xa, T));
2593 : : }
2594 : :
2595 : : GEN
2596 : 0 : Z_nv_mod(GEN P, GEN xa)
2597 : : {
2598 : 0 : return Z_ZV_mod(P, xa);
2599 : : }
2600 : :
2601 : : GEN
2602 : 57765 : ZX_nv_mod_tree(GEN P, GEN xa, GEN T)
2603 : : {
2604 : 57765 : long i, j, l = lg(P), n = lg(xa)-1;
2605 : 57765 : GEN V = cgetg(n+1, t_VEC);
2606 [ + + ]: 243006 : for (j=1; j <= n; j++)
2607 : : {
2608 : 185241 : gel(V, j) = cgetg(l, t_VECSMALL);
2609 : 185241 : mael(V, j, 1) = P[1]&VARNBITS;
2610 : : }
2611 [ + + ]: 678317 : for (i=2; i < l; i++)
2612 : : {
2613 : 620552 : GEN v = Z_ZV_mod_tree(gel(P, i), xa, T);
2614 [ + + ]: 2816846 : for (j=1; j <= n; j++)
2615 : 2196294 : mael(V, j, i) = v[j];
2616 : : }
2617 [ + + ]: 243006 : for (j=1; j <= n; j++)
2618 : 185241 : (void) Flx_renormalize(gel(V, j), l);
2619 : 57765 : return V;
2620 : : }
2621 : :
2622 : : static GEN
2623 : 145100 : ZV_sqr(GEN z)
2624 : : {
2625 : 145100 : long i,l = lg(z);
2626 : 145100 : GEN x = cgetg(l, t_VEC);
2627 [ + + ]: 145100 : if (typ(z)==t_VECSMALL)
2628 [ + + ]: 201951 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
2629 : : else
2630 [ + + ]: 468686 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
2631 : 145100 : return x;
2632 : : }
2633 : :
2634 : : static GEN
2635 : 844812 : ZT_sqr(GEN z)
2636 : : {
2637 [ + + ]: 844812 : if (typ(z) == t_INT)
2638 : 449800 : return sqri(z);
2639 : : else
2640 : : {
2641 : 395012 : long i,l = lg(z);
2642 : 395012 : GEN x = cgetg(l, t_VEC);
2643 [ + + ]: 1094724 : for (i=1; i<l; i++) gel(x,i) = ZT_sqr(gel(z,i));
2644 : 844812 : return x;
2645 : : }
2646 : : }
2647 : :
2648 : : static GEN
2649 : 145100 : ZV_invdivexact(GEN y, GEN x)
2650 : : {
2651 : 145100 : long i, l = lg(y);
2652 : 145100 : GEN z = cgetg(l,t_VEC);
2653 [ + + ]: 145100 : if (typ(x)==t_VECSMALL)
2654 [ + + ]: 201951 : for (i=1; i<l; i++)
2655 : : {
2656 : 156063 : pari_sp av = avma;
2657 : 156063 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
2658 : 156063 : avma = av;
2659 : 156063 : gel(z,i) = utoi(a);
2660 : : }
2661 : : else
2662 [ + + ]: 468686 : for (i=1; i<l; i++)
2663 : 369474 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
2664 : 145100 : return z;
2665 : : }
2666 : :
2667 : : static GEN
2668 : 145100 : ZV_chinesetree(GEN T, GEN xa)
2669 : : {
2670 : 145100 : GEN T2 = ZT_sqr(T), xa2 = ZV_sqr(xa);
2671 : 145100 : GEN mod = gmael(T,lg(T)-1,1);
2672 : 145100 : return ZV_invdivexact(Z_ZV_mod_tree(mod, xa2, T2), xa);
2673 : : }
2674 : :
2675 : : GEN
2676 : 45888 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN *pt_mod)
2677 : : {
2678 : 45888 : pari_sp av = avma;
2679 : 45888 : GEN R = ZV_chinesetree(T, P);
2680 : 45888 : GEN a = ZV_polint_tree(T, R, P, A);
2681 [ - + ]: 45888 : if (!pt_mod)
2682 : 0 : return gerepileuptoleaf(av, a);
2683 : : else
2684 : : {
2685 : 45888 : GEN mod = gmael(T, lg(T)-1, 1);
2686 : 45888 : gerepileall(av, 2, &a, &mod);
2687 : 45888 : *pt_mod = mod;
2688 : 45888 : return a;
2689 : : }
2690 : : }
2691 : :
2692 : : GEN
2693 : 99212 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
2694 : : {
2695 : 99212 : pari_sp av = avma;
2696 : 99212 : GEN T = ZV_producttree(P);
2697 : 99212 : GEN R = ZV_chinesetree(T, P);
2698 : 99212 : GEN a = ZV_polint_tree(T, R, P, A);
2699 [ - + ]: 99212 : if (!pt_mod)
2700 : 0 : return gerepileuptoint(av, a);
2701 : : else
2702 : : {
2703 : 99212 : GEN mod = gmael(T, lg(T)-1, 1);
2704 : 99212 : gerepileall(av, 2, &a, &mod);
2705 : 99212 : *pt_mod = mod;
2706 : 99212 : return a;
2707 : : }
2708 : : }
2709 : :
2710 : : GEN
2711 : 0 : nmV_chinese(GEN A, GEN P, GEN *pt_mod)
2712 : : {
2713 : 0 : pari_sp av = avma;
2714 : 0 : GEN T = ZV_producttree(P);
2715 : 0 : GEN R = ZV_chinesetree(T, P);
2716 : 0 : GEN a = nmV_polint_tree(T, R, P, A);
2717 [ # # ]: 0 : if (!pt_mod)
2718 : 0 : return gerepileupto(av, a);
2719 : : else
2720 : : {
2721 : 0 : GEN mod = gmael(T, lg(T)-1, 1);
2722 : 0 : gerepileall(av, 2, &a, &mod);
2723 : 0 : *pt_mod = mod;
2724 : 0 : return a;
2725 : : }
2726 : : }
2727 : :
2728 : : /**********************************************************************
2729 : : ** **
2730 : : ** Powering over (Z/NZ)^*, small N **
2731 : : ** **
2732 : : **********************************************************************/
2733 : :
2734 : : ulong
2735 : 15585001 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
2736 : : {
2737 : : ulong y, z, n;
2738 [ + + ]: 15585001 : if (n0 <= 1)
2739 : : { /* frequent special cases */
2740 [ + + ]: 251513 : if (n0 == 1) return x;
2741 [ + - ]: 96531 : if (n0 == 0) return 1;
2742 : : }
2743 [ + + ]: 15333488 : if (x <= 1) return x; /* 0 or 1 */
2744 : 15303929 : y = 1; z = x; n = n0;
2745 : : for(;;)
2746 : : {
2747 [ + + ]: 321480570 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
2748 [ + + ]: 321480570 : n>>=1; if (!n) return y;
2749 : 306176641 : z = Fl_sqr_pre(z,p,pi);
2750 : 321761642 : }
2751 : : }
2752 : :
2753 : : ulong
2754 : 27606784 : Fl_powu(ulong x, ulong n0, ulong p)
2755 : : {
2756 : : ulong y, z, n;
2757 [ + + ]: 27606784 : if (n0 <= 2)
2758 : : { /* frequent special cases */
2759 [ + + ]: 16261219 : if (n0 == 2) return Fl_sqr(x,p);
2760 [ + + ]: 1850796 : if (n0 == 1) return x;
2761 [ + - ]: 1126 : if (n0 == 0) return 1;
2762 : : }
2763 [ + + ]: 11345565 : if (x <= 1) return x; /* 0 or 1 */
2764 [ + + ]: 11318620 : if (!SMALL_ULONG(p))
2765 : 3925835 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
2766 : 7392785 : y = 1; z = x; n = n0;
2767 : : for(;;)
2768 : : {
2769 [ + + ]: 123291777 : if (n&1) y = Fl_mul(y,z,p);
2770 [ + + ]: 123304718 : n>>=1; if (!n) return y;
2771 : 115904744 : z = Fl_sqr(z,p);
2772 : 143512965 : }
2773 : : }
2774 : :
2775 : : /* Reduce data dependency to maximize internal parallelism */
2776 : : GEN
2777 : 959320 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
2778 : : {
2779 : : long i, k;
2780 : 959320 : GEN powers = cgetg(n + 2, t_VECSMALL);
2781 [ - + ]: 959320 : powers[1] = 1; if (n == 0) return powers;
2782 : 959320 : powers[2] = x;
2783 [ + + ]: 5382855 : for (i = 3, k=2; i <= n; i+=2, k++)
2784 : : {
2785 : 4423535 : powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
2786 : 4423535 : powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
2787 : : }
2788 [ + + ]: 959320 : if (i==n+1)
2789 : 847350 : powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
2790 : 959320 : return powers;
2791 : : }
2792 : :
2793 : : GEN
2794 : 104 : Fl_powers(ulong x, long n, ulong p)
2795 : : {
2796 : 104 : return Fl_powers_pre(x, n, p, get_Fl_red(p));
2797 : : }
2798 : :
2799 : : /**********************************************************************
2800 : : ** **
2801 : : ** Powering over (Z/NZ)^*, large N **
2802 : : ** **
2803 : : **********************************************************************/
2804 : :
2805 : : /* modified Barrett reduction with one fold */
2806 : : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
2807 : :
2808 : : static GEN
2809 : 6383 : Fp_invmBarrett(GEN p, long s)
2810 : : {
2811 : 6383 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
2812 : 6383 : return mkvec2(Q,R);
2813 : : }
2814 : :
2815 : : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
2816 : : * a = r (mod N) */
2817 : : static GEN
2818 : 324871 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
2819 : : {
2820 : 324871 : pari_sp av = avma;
2821 : 324871 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
2822 : 324871 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
2823 : 324871 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
2824 : 324871 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
2825 : 324871 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
2826 : 324871 : GEN r = subii(A, mulii(q, N));
2827 : 324871 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
2828 [ + + ]: 324871 : if (signe(sr)<0) return gerepileuptoint(av, r);
2829 : 311374 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
2830 [ + + ]: 311374 : if (signe(sr)<0) return gerepileuptoint(av, r);
2831 : 17878 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
2832 [ + + ]: 324871 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
2833 : : }
2834 : :
2835 : : /* Montgomery reduction */
2836 : :
2837 : : INLINE ulong
2838 : 227670 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
2839 : :
2840 : : typedef struct muldata {
2841 : : GEN N;
2842 : : GEN iM;
2843 : : ulong inv, s;
2844 : : GEN (*res)(struct muldata *,GEN);
2845 : : GEN (*mul2)(struct muldata *,GEN);
2846 : : } muldata;
2847 : :
2848 : : /* Montgomery reduction */
2849 : : static GEN
2850 : 12398187 : _montred(muldata *D, GEN x)
2851 : : {
2852 : 12398187 : return red_montgomery(x, D->N, D->inv);
2853 : : }
2854 : :
2855 : : static GEN
2856 : 4180845 : _remii(muldata *D, GEN x) { return remii(x, D->N); }
2857 : :
2858 : : static GEN
2859 : 324871 : _remiibar(muldata *D, GEN x) {
2860 : : #if DEBUG
2861 : : GEN r = Fp_rem_mBarrett(x, D->iM, D->s, D->N);
2862 : : if (cmpii(r, D->N) >= 0) pari_err_BUG("Rp_rem_mBarrett");
2863 : : return r;
2864 : : #else
2865 : 324871 : return Fp_rem_mBarrett(x, D->iM, D->s, D->N);
2866 : : #endif
2867 : : }
2868 : :
2869 : : /* 2x mod N */
2870 : : static GEN
2871 : 2484209 : _muli2red(muldata *D, GEN x)
2872 : : {
2873 : 2484209 : GEN z = shifti(x,1);
2874 [ + + ]: 2484212 : return (cmpii(z,D->N) >= 0)? subii(z,D->N): z;
2875 : : }
2876 : : static GEN
2877 : 1871757 : _muli2montred(muldata *D, GEN x)
2878 : : {
2879 : 1871757 : GEN z = _muli2red(D,x);
2880 : 1871760 : long l = lgefint(D->N);
2881 [ + + ]: 1873485 : while (lgefint(z) > l) z = subii(z,D->N);
2882 : 1871760 : return z;
2883 : : }
2884 : : static GEN
2885 : 849324 : _mul(void *data, GEN x, GEN y)
2886 : : {
2887 : 849324 : muldata *D = (muldata *)data;
2888 : 849324 : return D->res(D, mulii(x,y));
2889 : : }
2890 : : static GEN
2891 : 13338366 : _sqr(void *data, GEN x)
2892 : : {
2893 : 13338366 : muldata *D = (muldata *)data;
2894 : 13338366 : return D->res(D, sqri(x));
2895 : : }
2896 : : static GEN
2897 : 2484210 : _m2sqr(void *data, GEN x)
2898 : : {
2899 : 2484210 : muldata *D = (muldata *)data;
2900 : 2484210 : return D->mul2(D, D->res(D, sqri(x)));
2901 : : }
2902 : :
2903 : : static long
2904 : 303640 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D)
2905 : : {
2906 : 303640 : D->N = N;
2907 [ + + ][ + + ]: 303640 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
[ + + ]
2908 : : {
2909 : 6383 : D->mul2 = &_muli2red;
2910 : 6383 : D->res = &_remiibar;
2911 : 6383 : D->s = 1+(expi(N)>>1);
2912 : 6383 : D->iM = Fp_invmBarrett(N, D->s);
2913 : 6383 : return 0;
2914 : : }
2915 [ + + ][ + + ]: 297257 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
2916 : : {
2917 : 227671 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
2918 : 227670 : D->mul2 = &_muli2montred;
2919 : 227670 : D->res = &_montred;
2920 : 227670 : D->inv = init_montdata(N);
2921 : 227670 : return 1;
2922 : : }
2923 : : else
2924 : : {
2925 : 69587 : D->mul2 = &_muli2red;
2926 : 69587 : D->res = &_remii;
2927 : 303640 : return 0;
2928 : : }
2929 : : }
2930 : :
2931 : : GEN
2932 : 980436 : Fp_powu(GEN A, ulong k, GEN N)
2933 : : {
2934 : 980436 : long lN = lgefint(N), sA;
2935 : : int base_is_2, use_montgomery;
2936 : : muldata D;
2937 : : pari_sp av;
2938 : :
2939 [ + + ]: 980436 : if (lN == 3) {
2940 : 74223 : ulong n = uel(N,2);
2941 : 74223 : return utoi( Fl_powu(umodiu(A, n), k, n) );
2942 : : }
2943 [ + + ]: 906213 : if (k <= 2)
2944 : : { /* frequent special cases */
2945 [ + + ]: 704930 : if (k == 2) return Fp_sqr(A,N);
2946 [ + - ]: 112372 : if (k == 1) return A;
2947 [ # # ]: 0 : if (k == 0) return gen_1;
2948 : : }
2949 [ - + ][ # # ]: 201283 : sA = signe(A)==-1 && odd(k);
2950 : 201283 : base_is_2 = 0;
2951 [ + + ]: 201283 : if (lgefint(A) == 3) switch(A[2])
[ + + + ]
2952 : : {
2953 [ - + ]: 2316 : case 1: return sA ? gen_m1 : gen_1;
2954 : 46313 : case 2: base_is_2 = 1; break;
2955 : : }
2956 : :
2957 : : /* TODO: Move this out of here and use for general modular computations */
2958 : 198967 : av = avma;
2959 : 198967 : use_montgomery = Fp_select_red(&A, k, N, lN, &D);
2960 [ + + ]: 198967 : if (base_is_2)
2961 : 46313 : A = gen_powu_fold_i(A, k, (void*)&D, &_sqr, &_m2sqr);
2962 : : else
2963 : 152654 : A = gen_powu_i(A, k, (void*)&D, &_sqr, &_mul);
2964 [ + + ]: 198967 : if (use_montgomery)
2965 : : {
2966 : 171543 : A = _montred(&D, A);
2967 [ - + ]: 171543 : if (cmpii(A,N) >= 0) A = subii(A,N);
2968 [ - + ]: 171543 : if (sA) A = subii(N, A);
2969 : : }
2970 : 980436 : return gerepileuptoint(av, A);
2971 : : }
2972 : :
2973 : : GEN
2974 : 14336 : Fp_pows(GEN A, long k, GEN N)
2975 : : {
2976 [ + + ]: 14336 : if (lgefint(N) == 3) {
2977 : 309 : ulong n = N[2];
2978 : 309 : ulong a = umodiu(A, n);
2979 [ - + ]: 309 : if (k < 0) {
2980 : 0 : a = Fl_inv(a, n);
2981 : 0 : k = -k;
2982 : : }
2983 : 309 : return utoi( Fl_powu(a, (ulong)k, n) );
2984 : : }
2985 [ - + ]: 14027 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
2986 : 14336 : return Fp_powu(A, (ulong)k, N);
2987 : : }
2988 : :
2989 : : /* A^K mod N */
2990 : : GEN
2991 : 2351847 : Fp_pow(GEN A, GEN K, GEN N)
2992 : : {
2993 : 2351847 : pari_sp av = avma;
2994 : 2351847 : long t,s, lN = lgefint(N), sA;
2995 : : int base_is_2, use_montgomery;
2996 : : GEN y;
2997 : : muldata D;
2998 : :
2999 : 2351847 : s = signe(K);
3000 [ + + ]: 2351847 : if (!s)
3001 : : {
3002 : 5810 : t = signe(remii(A,N)); avma = av;
3003 [ + - ]: 5810 : return t? gen_1: gen_0;
3004 : : }
3005 [ + + ]: 2346037 : if (lN == 3)
3006 : : {
3007 : 2113832 : ulong k, n = N[2], a = umodiu(A, n);
3008 [ + + ]: 2113832 : if (s < 0) a = Fl_inv(a, n);
3009 [ + + ]: 2113776 : if (a <= 1) return utoi(a); /* 0 or 1 */
3010 [ + + ]: 2019959 : if (lgefint(K) > 3)
3011 : : { /* silly case : huge exponent, small modulus */
3012 : 21 : pari_warn(warner, "Mod(a,b)^n with n >> b : wasteful");
3013 [ + + ]: 21 : if (s > 0)
3014 : : {
3015 : 14 : ulong d = ugcd(a, n);
3016 [ + + ]: 14 : if (d != 1)
3017 : : { /* write n = n1 n2, with n2 maximal such that (n1,a) = 1 */
3018 : 7 : ulong n1 = ucoprime_part(n, d), n2 = n/n1;
3019 : :
3020 : 7 : k = umodiu(K, eulerphiu(n1));
3021 : : /* CRT: = a^K (mod n1), = 0 (mod n2)*/
3022 : 7 : return utoi( Fl_mul(Fl_powu(a, k, n1), n2 * Fl_inv(n2,n1), n) );
3023 : : }
3024 : : /* gcd(a,n) = 1 */
3025 : 7 : k = umodiu(K, eulerphiu(n));
3026 : : }
3027 : : else
3028 : 7 : k = umodiu(negi(K), eulerphiu(n));
3029 : : }
3030 : : else
3031 : 2019938 : k = uel(K,2);
3032 : 2019952 : return utoi(Fl_powu(a, k, n));
3033 : : }
3034 : :
3035 [ + + ]: 232205 : if (s < 0) y = Fp_inv(A,N);
3036 : : else
3037 : : {
3038 : 231854 : y = modii(A,N);
3039 [ + + ]: 231854 : if (!signe(y)) { avma = av; return gen_0; }
3040 : : }
3041 [ + + ]: 232198 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
3042 : :
3043 : 104714 : base_is_2 = 0;
3044 [ - + ][ # # ]: 104714 : sA = signe(y)==-1 && mod2(K);
3045 [ + + ]: 104714 : if (lgefint(y) == 3) switch(y[2])
[ + + + ]
3046 : : {
3047 [ - + ]: 41 : case 1: return sA ? gen_m1 : gen_1;
3048 : 82323 : case 2: base_is_2 = 1; break;
3049 : : }
3050 : :
3051 : : /* TODO: Move this out of here and use for general modular computations */
3052 : 104673 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D);
3053 [ + + ]: 104673 : if (base_is_2)
3054 : 82322 : y = gen_pow_fold_i(y, K, (void*)&D, &_sqr, &_m2sqr);
3055 : : else
3056 : 22351 : y = gen_pow_i(y, K, (void*)&D, &_sqr, &_mul);
3057 [ + + ]: 104674 : if (use_montgomery)
3058 : : {
3059 : 56128 : y = _montred(&D,y);
3060 [ - + ]: 56128 : if (cmpii(y,N) >= 0) y = subii(y,N);
3061 [ - + ]: 56128 : if (sA) y = subii(N, y);
3062 : : }
3063 : 2351792 : return gerepileuptoint(av,y);
3064 : : }
3065 : :
3066 : : static GEN
3067 : 142003 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
3068 : :
3069 : : static GEN
3070 : 60 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
3071 : :
3072 : : static GEN
3073 : 8 : _Fp_one(void *E) { (void) E; return gen_1; }
3074 : :
3075 : : GEN
3076 : 112 : Fp_powers(GEN x, long n, GEN p)
3077 : : {
3078 [ + + ]: 112 : if (lgefint(p) == 3)
3079 : 104 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
3080 : 112 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
3081 : : }
3082 : :
3083 : : static GEN
3084 : 485087 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
3085 : :
3086 : : static GEN
3087 : 406 : _Fp_rand(void *E) { return addis(randomi(subis((GEN)E,1)),1); }
3088 : :
3089 : : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
3090 : :
3091 : : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
3092 : : equalii,equali1,Fp_easylog};
3093 : :
3094 : : static GEN
3095 : 2639165 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
3096 : :
3097 : : static GEN
3098 : 25188318 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
3099 : :
3100 : : static GEN
3101 : 145974 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
3102 : :
3103 : : static GEN
3104 : 26186709 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
3105 : :
3106 : : static GEN
3107 : 60906 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
3108 : :
3109 : : static int
3110 : 1467321 : _Fp_equal0(GEN x) { return signe(x)==0; }
3111 : :
3112 : : static GEN
3113 : 125067 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
3114 : :
3115 : : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
3116 : : _Fp_inv,_Fp_equal0,_Fp_s};
3117 : :
3118 : 6750 : const struct bb_field *get_Fp_field(void **E, GEN p)
3119 : : {
3120 : 6750 : *E = (void*)p; return &Fp_field;
3121 : : }
3122 : :
3123 : : /*********************************************************************/
3124 : : /** **/
3125 : : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
3126 : : /** **/
3127 : : /*********************************************************************/
3128 : : ulong
3129 : 4872 : Fl_order(ulong a, ulong o, ulong p)
3130 : : {
3131 : 4872 : pari_sp av = avma;
3132 : : GEN m, P, E;
3133 : : long i;
3134 [ - + ]: 4872 : if (!o) o = p-1;
3135 : 4872 : m = factoru(o);
3136 : 4872 : P = gel(m,1);
3137 : 4872 : E = gel(m,2);
3138 [ + + ]: 12145 : for (i = lg(P)-1; i; i--)
3139 : : {
3140 : 7273 : ulong j, l=P[i], e=E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
3141 [ + + ]: 7273 : if (y == 1) o = t;
3142 : : else {
3143 [ + + ][ + + ]: 6727 : for (j = 1; j < e; j++) { y = Fl_powu(y, l, p); if (y == 1) break; }
3144 : 4739 : o = t * upowuu(l, j);
3145 : : }
3146 : : }
3147 : 4872 : avma = av; return o;
3148 : : }
3149 : :
3150 : : /*Find the exact order of a assuming a^o==1*/
3151 : : GEN
3152 : 11022 : Fp_order(GEN a, GEN o, GEN p) {
3153 [ + + ][ + + ]: 11022 : if (lgefint(p) == 3 && typ(o) == t_INT && lgefint(o)==3)
[ + - ]
3154 : : {
3155 : 21 : ulong pp = p[2], oo = o[2];
3156 : 21 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
3157 : : }
3158 : 11022 : return gen_order(a, o, (void*)p, &Fp_star);
3159 : : }
3160 : : GEN
3161 : 49 : Fp_factored_order(GEN a, GEN o, GEN p)
3162 : 49 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
3163 : :
3164 : : /* return order of a mod p^e, e > 0, pe = p^e */
3165 : : static GEN
3166 : 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
3167 : : {
3168 : : GEN ap, op;
3169 [ + + ]: 70 : if (equaliu(p, 2))
3170 : : {
3171 [ - + ]: 56 : if (e == 1) return gen_1;
3172 [ + + ][ - + ]: 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
3173 [ + + ]: 49 : if (mod4(a) == 1)
3174 : 14 : op = gen_1;
3175 : : else {
3176 : 35 : op = gen_2;
3177 : 35 : a = Fp_sqr(a, pe);
3178 : : }
3179 : : } else {
3180 [ - + ]: 14 : ap = (e == 1)? a: remii(a,p);
3181 : 14 : op = Fp_order(ap, subis(p,1), p);
3182 [ + - ]: 14 : if (e == 1) return op;
3183 : 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
3184 : : }
3185 [ + + ]: 49 : if (equali1(a)) return op;
3186 : 70 : return mulii(op, powiu(p, e - Z_pval(subis(a,1), p)));
3187 : : }
3188 : :
3189 : : GEN
3190 : 63 : znorder(GEN x, GEN o)
3191 : : {
3192 : 63 : pari_sp av = avma;
3193 : : GEN b, a;
3194 : :
3195 [ + + ]: 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
3196 : 56 : b = gel(x,1); a = gel(x,2);
3197 [ + + ]: 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
3198 [ + + ]: 49 : if (!o)
3199 : : {
3200 : 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
3201 : 35 : long i, l = lg(P);
3202 : 35 : o = gen_1;
3203 [ + + ]: 70 : for (i = 1; i < l; i++)
3204 : : {
3205 : 35 : GEN p = gel(P,i);
3206 : 35 : long e = itos(gel(E,i));
3207 : :
3208 [ + - ]: 35 : if (l == 2)
3209 : 35 : o = Zp_order(a, p, e, b);
3210 : : else {
3211 : 0 : GEN pe = powiu(p,e);
3212 : 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
3213 : : }
3214 : : }
3215 : 35 : return gerepileuptoint(av, o);
3216 : : }
3217 : 49 : return Fp_order(a, o, b);
3218 : : }
3219 : : GEN
3220 : 0 : order(GEN x) { return znorder(x, NULL); }
3221 : :
3222 : : /*********************************************************************/
3223 : : /** **/
3224 : : /** DISCRETE LOGARITHM in (Z/nZ)* **/
3225 : : /** **/
3226 : : /*********************************************************************/
3227 : : static GEN
3228 : 61585 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
3229 : : {
3230 : 61585 : pari_sp av = avma;
3231 : : GEN h1, h2, F, G;
3232 [ + + ]: 61585 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return NULL;
3233 [ + + ][ + + ]: 37103 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
3234 : : {
3235 : 214 : GEN M = cgetg(3, t_MAT);
3236 : 214 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
3237 : 214 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
3238 : 214 : return gerepileupto(av, M);
3239 : : }
3240 : 61585 : avma = av; return NULL;
3241 : : }
3242 : :
3243 : : static GEN
3244 : 61585 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
3245 : : {
3246 : : GEN rel;
3247 : : do
3248 : : {
3249 : 61585 : (*e)++; *g = Fp_mul(*g, b, p);
3250 : 61585 : rel = Fp_log_halfgcd(bnd, C, *g, p);
3251 [ + + ]: 61585 : } while (!rel);
3252 : 214 : return rel;
3253 : : }
3254 : :
3255 : : struct Fp_log_rel
3256 : : {
3257 : : GEN rel;
3258 : : long *sieve;
3259 : : ulong prmax;
3260 : : long nbrel, nbmax;
3261 : : };
3262 : :
3263 : : /* add u^e */
3264 : : static long
3265 : 36596 : addifsmooth1(struct Fp_log_rel *r, GEN h, long u, long e)
3266 : : {
3267 : 36596 : pari_sp av = avma;
3268 : : GEN z;
3269 [ + + ]: 36596 : if ((z = Z_issmooth_fact(h, r->prmax)))
3270 : : {
3271 : 2226 : long off = r->prmax+1;
3272 : 2226 : GEN F = cgetg(3, t_MAT);
3273 : 2226 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
3274 : 2226 : gel(F,2) = vecsmall_append(gel(z,2), e);
3275 : 2226 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3276 : : }
3277 : 36596 : return r->nbrel==r->nbmax;
3278 : : }
3279 : :
3280 : : /* add u^-1 v^-1 */
3281 : : static long
3282 : 248990 : addifsmooth2(struct Fp_log_rel *r, GEN h, long u, long v)
3283 : : {
3284 : 248990 : pari_sp av = avma;
3285 : : GEN z;
3286 [ + + ]: 248990 : if ((z = Z_issmooth_fact(h, r->prmax)))
3287 : : {
3288 : 93030 : long off = r->prmax+1;
3289 : 93030 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
3290 : 93030 : GEN F = cgetg(3, t_MAT);
3291 : 93030 : gel(F,1) = vecsmall_concat(gel(z,1), P);
3292 : 93030 : gel(F,2) = vecsmall_concat(gel(z,2), E);
3293 : 93030 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3294 : : }
3295 : 248990 : return r->nbrel==r->nbmax;
3296 : : }
3297 : :
3298 : : /*
3299 : : Let p=C^2+c
3300 : : Solve h = (C+x)*(C+a)-p = 0 [mod l]
3301 : : h= -c+x*(C+a)+C*a = 0 [mod l]
3302 : : x = (c-C*a)/(C+a) [mod l]
3303 : : h = -c+C*(x+a)+a*x
3304 : : */
3305 : :
3306 : : static void
3307 : 36127 : Fp_log_sieve_h(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, long a, GEN pr, GEN sz)
3308 : : {
3309 : 36127 : long th = expi(C), n = lg(pr)-1;
3310 : : long i,j;
3311 [ - + ]: 36127 : if (addifsmooth1(r, addis(C,a), a, -1)) return;
3312 [ + + ]: 25486818 : for(j=0; j<=a; j++)
3313 : 25450691 : r->sieve[j]=0;
3314 [ + + ]: 16596685 : for(i=1; i<=n; i++)
3315 : : {
3316 : 16560558 : ulong li = pr[i], s = sz[i], al = a % li;
3317 : 16560558 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
3318 [ + + ]: 16560558 : if (!iv) continue;
3319 : 16476712 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
3320 [ + + ]: 65275077 : for(j = u; j<=a; j+=li)
3321 : 48798365 : r->sieve[j] += s;
3322 : : }
3323 : 36127 : th = th - expu(th)-1;
3324 [ + + ]: 25427024 : for(j=0; j<a; j++)
3325 [ + + ]: 25390939 : if (r->sieve[j]>=th)
3326 : : {
3327 : 248990 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
3328 [ + + ]: 248990 : if (addifsmooth2(r, h, a, j)) return;
3329 : : }
3330 : : /* j = a */
3331 [ + + ]: 36085 : if (r->sieve[a]>=th)
3332 : : {
3333 : 469 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
3334 [ - + ]: 36127 : if (addifsmooth1(r, h, a, -2)) return;
3335 : : }
3336 : : }
3337 : :
3338 : : static GEN
3339 : 1234 : _psi(void*E, GEN y)
3340 : : {
3341 : 1234 : GEN lx = (GEN) E;
3342 : 1234 : long prec = lg(lx);
3343 : 1234 : GEN ly = glog(y, prec);
3344 : 1234 : GEN u = gdiv(lx, ly);
3345 : 1234 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
3346 : : }
3347 : :
3348 : : static GEN
3349 : 42 : opt_param(GEN x, long prec)
3350 : : {
3351 : 42 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
3352 : : }
3353 : :
3354 : : static GEN
3355 : 42 : check_kernel(long N, long prmax, GEN C, GEN M, GEN p, GEN m)
3356 : : {
3357 : 42 : pari_sp av = avma;
3358 : 42 : GEN K = FpMs_leftkernel_elt(M, N, m);
3359 : 42 : long i, f=0;
3360 : 42 : long l = lg(K), lm = lgefint(m);
3361 : 42 : GEN idx = diviiexact(subis(p,1),m), g;
3362 : : pari_timer ti;
3363 [ - + ]: 42 : if (DEBUGLEVEL) timer_start(&ti);
3364 [ + - ]: 84 : for(i=1; i<l; i++)
3365 [ + + ]: 84 : if (signe(gel(K,i)))
3366 : 42 : break;
3367 : 42 : g = utoi(i);
3368 : 42 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
3369 [ + + ]: 117467 : for(i=1; i<l; i++)
3370 : : {
3371 : 117425 : GEN k = gel(K,i);
3372 [ + + ]: 117425 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
3373 [ + + ][ - + ]: 117425 : if (signe(k)==0 || !equalii(Fp_pow(g, mulii(k,idx), p),
3374 : : Fp_pow(j, idx, p)))
3375 : 69608 : gel(K,i) = cgetineg(lm);
3376 : : else
3377 : 47817 : f++;
3378 : : }
3379 [ - + ]: 42 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld logs", f);
3380 : 42 : return gerepileupto(av, K);
3381 : : }
3382 : :
3383 : : static GEN
3384 : 84 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
3385 : : {
3386 : 84 : pari_sp av=avma;
3387 : 84 : GEN aa = gen_1;
3388 : 84 : long AV = 0;
3389 : : for(;;)
3390 : : {
3391 : 214 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
3392 : 214 : GEN F = gel(A,1), E = gel(A,2);
3393 : 214 : GEN Ao = gen_0;
3394 : 214 : long i, l = lg(F);
3395 [ + + ]: 1082 : for(i=1; i<l; i++)
3396 : : {
3397 : 998 : GEN Ki = gel(K,F[i]);
3398 [ + + ]: 998 : if (signe(Ki)<0) break;
3399 : 868 : Ao = addii(Ao, mulis(Ki, E[i]));
3400 : : }
3401 [ + + ]: 214 : if (i==l) return Fp_div(Ao, utoi(AV), m);
3402 : 130 : aa = gerepileuptoint(av, aa);
3403 : 130 : }
3404 : : }
3405 : :
3406 : : static GEN
3407 : 42 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
3408 : : {
3409 : 42 : pari_sp av = avma, av2;
3410 : : long i, nbi, nbrow;
3411 : : GEN C, c, Ci, ci, pr, sz, l, Ao, Bo, K, d, p_1;
3412 : : pari_timer ti;
3413 : : struct Fp_log_rel r;
3414 : 42 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
3415 : 42 : ulong bnd = 4*bnds;
3416 [ + - ][ - + ]: 42 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
3417 : :
3418 : 42 : p_1 = subiu(p,1);
3419 [ - + ]: 42 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
3420 : 0 : m = diviiexact(p_1, coprime_part(p_1, m));
3421 : 42 : pr = primes_upto_zv(bnd);
3422 : 42 : nbi = lg(pr)-1;
3423 [ - + ]: 42 : if (DEBUGLEVEL)
3424 : : {
3425 : 0 : err_printf("bnd=%lu Size FB=%ld\n", bnd, nbi);
3426 : 0 : timer_start(&ti);
3427 : : }
3428 : 42 : C = sqrtremi(p, &c);
3429 : 42 : av2 = avma;
3430 : 42 : Ci = cgetg(nbi+1,t_VECSMALL);
3431 : 42 : ci = cgetg(nbi+1,t_VECSMALL);
3432 : 42 : sz = cgetg(nbi+1,t_VECSMALL);
3433 [ + + ]: 11949 : for (i = 1; i <= nbi; ++i)
3434 : : {
3435 : 11907 : ulong lp = pr[i];
3436 : 11907 : Ci[i] = umodiu(C, lp);
3437 : 11907 : ci[i] = umodiu(c, lp);
3438 : 11907 : sz[i] = expu(lp);
3439 : : }
3440 : 42 : r.nbrel = 0;
3441 : 42 : r.nbmax = 8*nbi;
3442 : 42 : r.rel = cgetg(r.nbmax+1,t_VEC);
3443 : 42 : r.sieve = cgetg(r.nbmax+2,t_VECSMALL)+1;
3444 : 42 : r.prmax = pr[nbi];
3445 [ + + ]: 36169 : for(i=0; r.nbrel < r.nbmax; i++)
3446 : : {
3447 : 36127 : Fp_log_sieve_h(&r, C, c, Ci, ci, i, pr, sz);
3448 [ - + ][ # # ]: 36127 : if (DEBUGLEVEL && (i&127)==0)
3449 : 0 : err_printf("%ld%% ",100*r.nbrel/(r.nbmax));
3450 : : }
3451 : 42 : nbrow = r.prmax+i;
3452 [ - + ]: 42 : if (DEBUGLEVEL)
3453 : : {
3454 : 0 : err_printf("\n");
3455 : 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+i);
3456 : : }
3457 : 42 : setlg(r.rel,r.nbrel+1);
3458 : 42 : r.rel = gerepileupto(av2, r.rel);
3459 : 42 : K = check_kernel(nbrow, r.prmax, C, r.rel, p, m);
3460 [ - + ]: 42 : if (DEBUGLEVEL) timer_start(&ti);
3461 : 42 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
3462 [ - + ]: 42 : if (DEBUGLEVEL) timer_printf(&ti," log element");
3463 : 42 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
3464 [ - + ]: 42 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
3465 : 42 : d = gcdii(Ao,Bo);
3466 : 42 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
3467 [ - + ]: 42 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
3468 : 42 : return gerepileuptoint(av, l);
3469 : : }
3470 : :
3471 : : static int
3472 : 164731 : Fp_log_use_index(long e, long p)
3473 : : {
3474 [ + + ][ + - ]: 164731 : return (e >= 27 && 20*(p+6)<=e*e);
3475 : : }
3476 : :
3477 : : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
3478 : : static GEN
3479 : 181238 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
3480 : : {
3481 : 181238 : pari_sp av = avma;
3482 : 181238 : GEN p = (GEN)E;
3483 : : /* assume a reduced mod p, p not necessarily prime */
3484 [ + + ]: 181238 : if (equali1(a)) return gen_0;
3485 : : /* p > 2 */
3486 [ + + ]: 123033 : if (equalii(subis(p,1), a)) /* -1 */
3487 : : {
3488 : : pari_sp av2;
3489 : : GEN t;
3490 : 55904 : ord = dlog_get_ord(ord);
3491 [ + + ]: 55904 : if (mpodd(ord)) { avma = av; return cgetg(1, t_VEC); } /* no solution */
3492 : 55890 : t = shifti(ord,-1); /* only possible solution */
3493 : 55890 : av2 = avma;
3494 [ + + ]: 55890 : if (!equalii(Fp_pow(g, t, p), a)) { avma = av; return cgetg(1, t_VEC); }
3495 : 55862 : avma = av2; return gerepileuptoint(av, t);
3496 : : }
3497 [ + + ][ + + ]: 67129 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
[ + + ]
3498 : 42 : return Fp_log_index(a, g, ord, p);
3499 : 181238 : avma = av; return NULL; /* not easy */
3500 : : }
3501 : :
3502 : : GEN
3503 : 191706 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
3504 : : {
3505 : 191706 : GEN v = dlog_get_ordfa(ord);
3506 : 191678 : GEN F = gmael(v,2,1);
3507 : 191678 : long lF = lg(F)-1, lmax;
3508 [ + + ][ + + ]: 191678 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
3509 : 137407 : lmax = expi(gel(F,lF));
3510 [ + + ][ + + ]: 137407 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
3511 : 42 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
3512 : 191678 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
3513 : : }
3514 : :
3515 : : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
3516 : : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
3517 : : static GEN
3518 : 105 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
3519 : : {
3520 : 105 : long l = lg(P) - 1, e = E[l];
3521 [ + + ]: 105 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
3522 : : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
3523 : :
3524 [ + + ]: 105 : if (l == 1) {
3525 : 77 : hpe = h;
3526 : 77 : gpe = g;
3527 : : } else {
3528 : 28 : hpe = modii(h, pe);
3529 : 28 : gpe = modii(g, pe);
3530 : : }
3531 [ + + ]: 105 : if (e == 1) {
3532 : 21 : hp = hpe;
3533 : 21 : gp = gpe;
3534 : : } else {
3535 : 84 : hp = remii(hpe, p);
3536 : 84 : gp = remii(gpe, p);
3537 : : }
3538 [ + + ][ - + ]: 105 : if (hp == gen_0 || gp == gen_0) return NULL;
3539 [ + + ]: 84 : if (equaliu(p, 2))
3540 : : {
3541 : 35 : GEN n = int2n(e);
3542 : 35 : ogpe = Zp_order(gpe, gen_2, e, n);
3543 : 35 : a = Fp_log(hpe, gpe, ogpe, n);
3544 [ + + ]: 35 : if (typ(a) != t_INT) return NULL;
3545 : : }
3546 : : else
3547 : : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
3548 : : is trivial */
3549 : : /* [order(gp), factor(order(gp))] */
3550 : 49 : GEN v = Fp_factored_order(gp, subis(p,1), p);
3551 : 49 : GEN ogp = gel(v,1);
3552 [ - + ]: 49 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
3553 : 49 : a = Fp_log(hp, gp, v, p);
3554 [ - + ]: 49 : if (typ(a) != t_INT) return NULL;
3555 [ + + ]: 49 : if (e == 1) ogpe = ogp;
3556 : : else
3557 : : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
3558 : : /* use p-adic log: O(log p + e) mul*/
3559 : : long vpogpe, vpohpe;
3560 : :
3561 : 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
3562 : 28 : gpe = Fp_pow(gpe, ogp, pe);
3563 : : /* g,h = 1 mod p; compute b s.t. h = g^b */
3564 : :
3565 : : /* v_p(order g mod pe) */
3566 [ + - ]: 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subis(gpe,1), p);
3567 : : /* v_p(order h mod pe) */
3568 [ + - ]: 28 : vpohpe = equali1(gpe)? 0: e - Z_pval(subis(hpe,1), p);
3569 [ - + ]: 28 : if (vpohpe > vpogpe) return NULL;
3570 : :
3571 : 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
3572 [ - + ][ # # ]: 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
3573 : 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
3574 : 28 : a = addii(a, mulii(ogp, gtrunc(b)));
3575 : : }
3576 : : }
3577 : : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
3578 [ + + ]: 70 : if (l == 1) return a;
3579 : :
3580 : 28 : N = diviiexact(N, pe); /* make N coprime to p */
3581 : 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
3582 : 28 : g = Fp_pow(g, modii(ogpe, phi), N);
3583 : 28 : setlg(P, l); /* remove last element */
3584 : 28 : setlg(E, l);
3585 : 28 : b = znlog_rec(h, g, N, P, E, PHI);
3586 [ - + ]: 28 : if (!b) return NULL;
3587 : 105 : return addmulii(a, b, ogpe);
3588 : : }
3589 : :
3590 : : static GEN
3591 : 77 : get_PHI(GEN P, GEN E)
3592 : : {
3593 : 77 : long i, l = lg(P);
3594 : 77 : GEN PHI = cgetg(l, t_VEC);
3595 : 77 : gel(PHI,1) = gen_1;
3596 [ + + ]: 105 : for (i=1; i<l-1; i++)
3597 : : {
3598 : 28 : GEN t, p = gel(P,i);
3599 : 28 : long e = E[i];
3600 : 28 : t = mulii(powiu(p, e-1), subis(p,1));
3601 [ + + ]: 28 : if (i > 1) t = mulii(t, gel(PHI,i));
3602 : 28 : gel(PHI,i+1) = t;
3603 : : }
3604 : 77 : return PHI;
3605 : : }
3606 : :
3607 : : GEN
3608 : 203 : znlog(GEN h, GEN g, GEN o)
3609 : : {
3610 : 203 : pari_sp av = avma;
3611 : : GEN N, fa, P, E, x;
3612 [ + + - ]: 203 : switch (typ(g))
3613 : : {
3614 : : case t_PADIC:
3615 : : {
3616 : 28 : GEN p = gel(g,2);
3617 : 28 : long v = valp(g);
3618 [ - + ]: 28 : if (v < 0) pari_err_DIM("znlog");
3619 [ - + ]: 28 : if (v > 0) {
3620 : 0 : long k = gvaluation(h, p);
3621 [ # # ]: 0 : if (k % v) return cgetg(1,t_VEC);
3622 : 0 : k /= v;
3623 [ # # ]: 0 : if (!gequal(h, gpowgs(g,k))) { avma = av; return cgetg(1,t_VEC); }
3624 : 0 : avma = av; return stoi(k);
3625 : : }
3626 : 28 : N = gel(g,3);
3627 : 28 : g = Rg_to_Fp(g, N);
3628 : 28 : break;
3629 : : }
3630 : : case t_INTMOD:
3631 : 175 : N = gel(g,1);
3632 : 175 : g = gel(g,2); break;
3633 : 0 : default: pari_err_TYPE("znlog", g);
3634 : 0 : return NULL; /* not reached */
3635 : : }
3636 [ - + ]: 203 : if (equali1(N)) { avma = av; return gen_0; }
3637 : 203 : h = Rg_to_Fp(h, N);
3638 [ + + ]: 203 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
3639 : 77 : fa = Z_factor(N);
3640 : 77 : P = gel(fa,1);
3641 : 77 : E = vec_to_vecsmall(gel(fa,2));
3642 : 77 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
3643 [ + + ]: 77 : if (!x) { avma = av; return cgetg(1,t_VEC); }
3644 : 175 : return gerepileuptoint(av, x);
3645 : : }
3646 : :
3647 : : GEN
3648 : 1715 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
3649 : : {
3650 : 1715 : a = modii(a,p);
3651 [ + + ]: 1715 : if (!signe(a))
3652 : : {
3653 [ - + ]: 7 : if (zeta) *zeta = gen_1;
3654 [ + - ]: 7 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
3655 : 0 : return gen_0;
3656 : : }
3657 [ + + ]: 1708 : if (equaliu(n,2))
3658 : : {
3659 [ + + ]: 854 : if (zeta) *zeta = addis(p,-1);
3660 : 854 : return Fp_sqrt(a,p);
3661 : : }
3662 : 1708 : return gen_Shanks_sqrtn(a,n,addis(p,-1),zeta,(void*)p,&Fp_star);
3663 : : }
3664 : :
3665 : : /*********************************************************************/
3666 : : /** **/
3667 : : /** FUNDAMENTAL DISCRIMINANTS **/
3668 : : /** **/
3669 : : /*********************************************************************/
3670 : : long
3671 : 14021 : isfundamental(GEN x) {
3672 [ - + ]: 14021 : if (typ(x) != t_INT) pari_err_TYPE("isfundamental",x);
3673 : 14021 : return Z_isfundamental(x);
3674 : : }
3675 : :
3676 : : /* x fundamental ? */
3677 : : long
3678 : 6293 : uposisfundamental(ulong x)
3679 : : {
3680 : 6293 : ulong r = x & 15; /* x mod 16 */
3681 [ + + ]: 6293 : if (!r) return 0;
3682 [ + + + ]: 5915 : switch(r & 3)
3683 : : { /* x mod 4 */
3684 [ + + ]: 1184 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
3685 : 1689 : case 1: return uissquarefree(x);
3686 : 6293 : default: return 0;
3687 : : }
3688 : : }
3689 : : /* -x fundamental ? */
3690 : : long
3691 : 9807 : unegisfundamental(ulong x)
3692 : : {
3693 : 9807 : ulong r = x & 15; /* x mod 16 */
3694 [ + + ]: 9807 : if (!r) return 0;
3695 [ + + + ]: 9310 : switch(r & 3)
3696 : : { /* x mod 4 */
3697 [ + + ]: 1765 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
3698 : 4545 : case 3: return uissquarefree(x);
3699 : 9807 : default: return 0;
3700 : : }
3701 : : }
3702 : : long
3703 : 14574 : Z_isfundamental(GEN x)
3704 : : {
3705 : : long r;
3706 [ - + + ]: 14574 : switch(lgefint(x))
3707 : : {
3708 : 0 : case 2: return 0;
3709 : 12572 : case 3: return signe(x) < 0? unegisfundamental(x[2])
3710 [ + + ]: 12572 : : uposisfundamental(x[2]);
3711 : : }
3712 : 2002 : r = mod16(x);
3713 [ + + ]: 2002 : if (!r) return 0;
3714 [ + + ]: 1876 : if ((r & 3) == 0)
3715 : : {
3716 : : pari_sp av;
3717 : 376 : r >>= 2; /* |x|/4 mod 4 */
3718 [ + + ]: 376 : if (signe(x) < 0) r = 4-r;
3719 [ + + ]: 376 : if (r == 1) return 0;
3720 : 250 : av = avma;
3721 : 250 : r = Z_issquarefree( shifti(x,-2) );
3722 : 250 : avma = av; return r;
3723 : : }
3724 : 1500 : r &= 3; /* |x| mod 4 */
3725 [ + + ]: 1500 : if (signe(x) < 0) r = 4-r;
3726 [ + + ]: 14574 : return (r==1) ? Z_issquarefree(x) : 0;
3727 : : }
3728 : :
3729 : : GEN
3730 : 7 : quaddisc(GEN x)
3731 : : {
3732 : 7 : const pari_sp av = avma;
3733 : 7 : long i,r,tx=typ(x);
3734 : : GEN P,E,f,s;
3735 : :
3736 [ - + ]: 7 : if (!is_rational_t(tx)) pari_err_TYPE("quaddisc",x);
3737 : 7 : f = factor(x);
3738 : 7 : P = gel(f,1);
3739 : 7 : E = gel(f,2); s = gen_1;
3740 [ + + ]: 35 : for (i=1; i<lg(P); i++)
3741 [ + + ]: 28 : if (odd(mael(E,i,2))) s = mulii(s,gel(P,i));
3742 [ + - ]: 7 : r = mod4(s); if (gsigne(x) < 0) r = 4-r;
3743 [ - + ]: 7 : if (r>1) s = shifti(s,2);
3744 : 7 : return gerepileuptoint(av, s);
3745 : : }
3746 : :
3747 : : /*********************************************************************/
3748 : : /** **/
3749 : : /** FACTORIAL **/
3750 : : /** **/
3751 : : /*********************************************************************/
3752 : : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
3753 : : * first is slower ... ] */
3754 : : GEN
3755 : 158085 : mulu_interval(ulong a, ulong b)
3756 : : {
3757 : 158085 : pari_sp av = avma;
3758 : 158085 : ulong k, l, N, n = b - a + 1;
3759 : : long lx;
3760 : : GEN x;
3761 : :
3762 [ + + ]: 158085 : if (n < 61)
3763 : : {
3764 : 153361 : x = utoi(a);
3765 [ + + ]: 1349635 : for (k=a+1; k<=b; k++) x = mului(k,x);
3766 : 153361 : return gerepileuptoint(av, x);
3767 : : }
3768 : 4724 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3769 : 4724 : N = b + a;
3770 : 4724 : for (k = a;; k++)
3771 : : {
3772 [ + + ]: 628594 : l = N - k; if (l <= k) break;
3773 : 623870 : gel(x,lx++) = muluu(k,l);
3774 : 623870 : }
3775 [ + + ]: 4724 : if (l == k) gel(x,lx++) = utoipos(k);
3776 : 4724 : setlg(x, lx);
3777 : 158085 : return gerepileuptoint(av, divide_conquer_prod(x, mulii));
3778 : : }
3779 : :
3780 : : GEN
3781 : 75422 : mpfact(long n)
3782 : : {
3783 [ + + ]: 75422 : if (n < 2)
3784 : : {
3785 [ - + ]: 1904 : if (n < 0) pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
3786 : 1904 : return gen_1;
3787 : : }
3788 : 75422 : return mulu_interval(2UL, (ulong)n);
3789 : : }
3790 : :
3791 : : /*******************************************************************/
3792 : : /** **/
3793 : : /** LUCAS & FIBONACCI **/
3794 : : /** **/
3795 : : /*******************************************************************/
3796 : : static void
3797 : 56 : lucas(ulong n, GEN *a, GEN *b)
3798 : : {
3799 : : GEN z, t, zt;
3800 [ + + ]: 105 : if (!n) { *a = gen_2; *b = gen_1; return; }
3801 : 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
3802 [ + + + + : 49 : switch(n & 3) {
- ]
3803 : 14 : case 0: *a = addsi(-2,sqri(z)); *b = addsi(-1,zt); break;
3804 : 14 : case 1: *a = addsi(-1,zt); *b = addsi(2,sqri(t)); break;
3805 : 7 : case 2: *a = addsi(2,sqri(z)); *b = addsi(1,zt); break;
3806 : 49 : case 3: *a = addsi(1,zt); *b = addsi(-2,sqri(t));
3807 : : }
3808 : : }
3809 : :
3810 : : GEN
3811 : 7 : fibo(long n)
3812 : : {
3813 : 7 : pari_sp av = avma;
3814 : : GEN a, b;
3815 [ - + ]: 7 : if (!n) return gen_0;
3816 : 7 : lucas((ulong)(labs(n)-1), &a, &b);
3817 : 7 : a = diviuexact(addii(shifti(a,1),b), 5);
3818 [ - + ][ # # ]: 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
3819 : 7 : return gerepileuptoint(av, a);
3820 : : }
3821 : :
3822 : : /*******************************************************************/
3823 : : /* */
3824 : : /* CONTINUED FRACTIONS */
3825 : : /* */
3826 : : /*******************************************************************/
3827 : : static GEN
3828 : 2414400 : icopy_lg(GEN x, long l)
3829 : : {
3830 : 2414400 : long lx = lgefint(x);
3831 : : GEN y;
3832 : :
3833 [ + + ]: 2414400 : if (lx >= l) return icopy(x);
3834 : 2414400 : y = cgeti(l); affii(x, y); return y;
3835 : : }
3836 : :
3837 : : /* continued fraction of a/b. If y != NULL, stop when partial quotients
3838 : : * differ from y */
3839 : : static GEN
3840 : 2414689 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
3841 : : {
3842 : : GEN z, c;
3843 : 2414689 : ulong i, l, ly = lgefint(b);
3844 : :
3845 : : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
3846 : 2414689 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
3847 [ + + ][ + - ]: 2414689 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
[ + - ]
3848 [ - + ]: 2414689 : if (l > LGBITS) l = LGBITS;
3849 : :
3850 : 2414689 : z = cgetg(l,t_VEC);
3851 : 2414689 : l--;
3852 [ + + ]: 2414689 : if (y) {
3853 : 289 : pari_sp av = avma;
3854 [ + + ]: 289 : if (l >= (ulong)lg(y)) l = lg(y)-1;
3855 [ + + ]: 19578 : for (i = 1; i <= l; i++)
3856 : : {
3857 : 19289 : GEN q = gel(y,i);
3858 : 19289 : gel(z,i) = q;
3859 [ + + ]: 19289 : c = b; if (!gequal1(q)) c = mulii(q, b);
3860 : 19289 : c = subii(a, c);
3861 [ + + ]: 19289 : if (signe(c) < 0)
3862 : : { /* partial quotient too large */
3863 : 110 : c = addii(c, b);
3864 [ + + ]: 110 : if (signe(c) >= 0) i++; /* by 1 */
3865 : 110 : break;
3866 : : }
3867 [ + + ]: 19179 : if (cmpii(c, b) >= 0)
3868 : : { /* partial quotient too small */
3869 : 1 : c = subii(c, b);
3870 [ - + ]: 1 : if (cmpii(c, b) < 0) {
3871 : : /* by 1. If next quotient is 1 in y, add 1 */
3872 [ # # ][ # # ]: 0 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addis(q,1);
3873 : 0 : i++;
3874 : : }
3875 : 1 : break;
3876 : : }
3877 [ + + ]: 19178 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
3878 : 19178 : a = b; b = c;
3879 : : }
3880 : : } else {
3881 : 2414400 : a = icopy_lg(a, ly);
3882 : 2414400 : b = icopy(b);
3883 [ + + ]: 22421521 : for (i = 1; i <= l; i++)
3884 : : {
3885 : 22421251 : gel(z,i) = truedvmdii(a,b,&c);
3886 [ + + ]: 22421251 : if (c == gen_0) { i++; break; }
3887 : 20007121 : affii(c, a); cgiv(c); c = a;
3888 : 20007121 : a = b; b = c;
3889 : : }
3890 : : }
3891 : 2414689 : i--;
3892 [ + - ][ + + ]: 2414689 : if (i > 1 && gequal1(gel(z,i)))
3893 : : {
3894 : 79 : cgiv(gel(z,i)); --i;
3895 : 79 : gel(z,i) = addsi(1, gel(z,i)); /* unclean: leave old z[i] on stack */
3896 : : }
3897 : 2414689 : setlg(z,i+1); return z;
3898 : : }
3899 : :
3900 : : static GEN
3901 : 0 : sersfcont(GEN a, GEN b, long k)
3902 : : {
3903 [ # # ]: 0 : long i, l = typ(a) == t_POL? lg(a): 3;
3904 : : GEN y, c;
3905 [ # # ]: 0 : if (lg(b) > l) l = lg(b);
3906 [ # # ][ # # ]: 0 : if (k > 0 && l > k+1) l = k+1;
3907 : 0 : y = cgetg(l,t_VEC);
3908 [ # # ]: 0 : for (i=1; i<l; i++)
3909 : : {
3910 : 0 : gel(y,i) = poldivrem(a,b,&c);
3911 [ # # ]: 0 : if (gequal0(c)) { i++; break; }
3912 : 0 : a = b; b = c;
3913 : : }
3914 : 0 : setlg(y, i); return y;
3915 : : }
3916 : :
3917 : : GEN
3918 : 2414484 : gboundcf(GEN x, long k)
3919 : : {
3920 : : pari_sp av;
3921 : 2414484 : long tx = typ(x), e;
3922 : : GEN y, a, b, c;
3923 : :
3924 [ + + ]: 2414484 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
3925 [ + - ]: 2414477 : if (is_scalar_t(tx))
3926 : : {
3927 [ + + ]: 2414477 : if (gequal0(x)) return mkvec(gen_0);
3928 [ - + + - ]: 2414407 : switch(tx)
3929 : : {
3930 : 0 : case t_INT: return mkveccopy(x);
3931 : : case t_REAL:
3932 : 296 : av = avma;
3933 : 296 : c = mantissa_real(x,&e);
3934 [ + + ]: 296 : if (e < 0) pari_err_PREC("gboundcf");
3935 : 289 : y = int2n(e);
3936 : 289 : a = Qsfcont(c,y, NULL, k);
3937 : 289 : b = addsi(signe(x), c);
3938 : 289 : return gerepilecopy(av, Qsfcont(b,y, a, k));
3939 : :
3940 : : case t_FRAC:
3941 : 2414111 : av = avma;
3942 : 2414111 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
3943 : : }
3944 : 0 : pari_err_TYPE("gboundcf",x);
3945 : : }
3946 : :
3947 [ # # # # ]: 0 : switch(tx)
3948 : : {
3949 : 0 : case t_POL: return mkveccopy(x);
3950 : : case t_SER:
3951 : 0 : av = avma;
3952 : 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
3953 : : case t_RFRAC:
3954 : 0 : av = avma;
3955 : 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
3956 : : }
3957 : 0 : pari_err_TYPE("gboundcf",x);
3958 : 2414470 : return NULL; /* not reached */
3959 : : }
3960 : :
3961 : : static GEN
3962 : 14 : sfcont2(GEN b, GEN x, long k)
3963 : : {
3964 : 14 : pari_sp av = avma;
3965 : 14 : long lb = lg(b), tx = typ(x), i;
3966 : : GEN y,p1;
3967 : :
3968 [ + + ]: 14 : if (k)
3969 : : {
3970 [ + - ]: 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
3971 : 0 : lb = k+1;
3972 : : }
3973 : 7 : y = cgetg(lb,t_VEC);
3974 [ - + ]: 7 : if (lb==1) return y;
3975 [ + - ]: 7 : if (is_scalar_t(tx))
3976 : : {
3977 [ - + ][ # # ]: 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
3978 : : }
3979 [ # # ]: 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
3980 : :
3981 [ - + ]: 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
3982 : 7 : for (i = 1;;)
3983 : : {
3984 [ + - ]: 35 : if (tx == t_REAL)
3985 : : {
3986 : 35 : long e = expo(x);
3987 [ + + ][ - + ]: 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
3988 : 35 : gel(y,i) = floorr(x);
3989 : 35 : p1 = subri(x, gel(y,i));
3990 : : }
3991 : : else
3992 : : {
3993 : 0 : gel(y,i) = gfloor(x);
3994 : 0 : p1 = gsub(x, gel(y,i));
3995 : : }
3996 [ + + ]: 35 : if (++i >= lb) break;
3997 [ - + ]: 28 : if (gequal0(p1)) break;
3998 : 28 : x = gdiv(gel(b,i),p1);
3999 : 28 : }
4000 : 7 : setlg(y,i);
4001 : 7 : return gerepilecopy(av,y);
4002 : : }
4003 : :
4004 : :
4005 : : GEN
4006 : 0 : gcf(GEN x) { return gboundcf(x,0); }
4007 : : GEN
4008 : 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
4009 : : GEN
4010 : 49 : contfrac0(GEN x, GEN b, long nmax)
4011 : : {
4012 : : long tb;
4013 : :
4014 [ + + ]: 49 : if (!b) return gboundcf(x,nmax);
4015 : 28 : tb = typ(b);
4016 [ + + ]: 28 : if (tb == t_INT) return gboundcf(x,itos(b));
4017 [ - + ]: 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
4018 [ + + ]: 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
4019 : 28 : return sfcont2(b,x,nmax);
4020 : : }
4021 : :
4022 : : GEN
4023 : 126 : contfracpnqn(GEN x, long n)
4024 : : {
4025 : 126 : pari_sp av = avma;
4026 : 126 : long i, lx = lg(x);
4027 : : GEN M,A,B, p0,p1, q0,q1;
4028 : :
4029 [ + + ]: 126 : if (lx == 1)
4030 : : {
4031 [ + + ]: 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
4032 [ + + ]: 21 : if (n >= 0) return cgetg(1,t_MAT);
4033 : 7 : return matid(2);
4034 : : }
4035 [ + + - ]: 98 : switch(typ(x))
4036 : : {
4037 : 56 : case t_VEC: case t_COL: A = x; B = NULL; break;
4038 : : case t_MAT:
4039 [ - + + ]: 42 : switch(lgcols(x))
4040 : : {
4041 : 0 : case 2: A = row(x,1); B = NULL; break;
4042 : 35 : case 3: A = row(x,2); B = row(x,1); break;
4043 : 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
4044 : 0 : return NULL; /*not reached*/
4045 : : }
4046 : 35 : break;
4047 : 0 : default: pari_err_TYPE("pnqn",x);
4048 : 0 : return NULL; /*not reached*/
4049 : : }
4050 : 91 : p1 = gel(A,1);
4051 [ + + ]: 91 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
4052 [ + + ]: 91 : if (n >= 0)
4053 : : {
4054 : 56 : lx = minss(lx, n+2);
4055 [ + + ]: 56 : if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
4056 : : }
4057 [ + + ]: 35 : else if (lx == 2)
4058 : 7 : return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
4059 : : /* lx >= 3 */
4060 : 56 : p0 = gen_1;
4061 : 56 : q0 = gen_0; /* p[-1], q[-1] */
4062 : 56 : M = cgetg(lx, t_MAT);
4063 : 56 : gel(M,1) = mkcol2(p1,q1);
4064 [ + + ]: 217 : for (i=2; i<lx; i++)
4065 : : {
4066 : 161 : GEN a = gel(A,i), p2,q2;
4067 [ + + ]: 161 : if (B) {
4068 : 84 : GEN b = gel(B,i);
4069 : 84 : p0 = gmul(b,p0);
4070 : 84 : q0 = gmul(b,q0);
4071 : : }
4072 : 161 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
4073 : 161 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
4074 : 161 : gel(M,i) = mkcol2(p1,q1);
4075 : : }
4076 [ + + ]: 56 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
4077 : 112 : return gerepilecopy(av, M);
4078 : : }
4079 : : GEN
4080 : 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
4081 : : /* x = [a0, ..., an] from gboundcf, n >= 0;
4082 : : * return [[p0, ..., pn], [q0,...,qn]] */
4083 : : GEN
4084 : 423311 : ZV_allpnqn(GEN x)
4085 : : {
4086 : 423311 : long i, lx = lg(x);
4087 : 423311 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
4088 : :
4089 : 423311 : gel(v,1) = P = cgetg(lx, t_VEC);
4090 : 423311 : gel(v,2) = Q = cgetg(lx, t_VEC);
4091 : 423311 : p0 = gen_1; q0 = gen_0;
4092 : 423311 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
4093 [ + + ]: 1420678 : for (i=2; i<lx; i++)
4094 : : {
4095 : 997367 : GEN a = gel(x,i);
4096 : 997367 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
4097 : 997367 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
4098 : : }
4099 : 423311 : return v;
4100 : : }
4101 : :
4102 : : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
4103 : : static GEN
4104 : 35 : mod_to_frac(GEN x, GEN N, GEN B)
4105 : : {
4106 : : GEN a, b, A;
4107 [ + + ]: 35 : if (B)
4108 : : {
4109 : 21 : A = divii(shifti(N, -1), B);
4110 : : /* denominator bound useless, don't use it */
4111 [ - + ]: 21 : if (cmpii(A, B) < 0) B = NULL;
4112 : : }
4113 [ + + ]: 35 : if (!B)
4114 : : {
4115 : 14 : A = sqrti(shifti(N, -1));
4116 : 14 : B = A;
4117 : : }
4118 [ + + ][ - + ]: 35 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
4119 [ + - ]: 35 : return equali1(b)? a: mkfrac(a,b);
4120 : : }
4121 : :
4122 : : static GEN
4123 : 56 : mod_to_rfrac(GEN x, GEN N, long B)
4124 : : {
4125 : : GEN a, b;
4126 : 56 : long A, d = degpol(N);
4127 [ + + ]: 56 : if (B >= 0)
4128 : : {
4129 : 21 : A = d-1 - B;
4130 : : /* denominator bound useless, don't use it */
4131 [ + + ]: 21 : if (A < B) B = -1;
4132 : : }
4133 [ + + ]: 56 : if (B < 0)
4134 : : {
4135 : 42 : B = d >> 1;
4136 [ + + ]: 42 : A = odd(d)? B : B-1;
4137 : : }
4138 [ - + ]: 56 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
4139 [ - + ]: 56 : if (! RgXQ_ratlift(x, N, A, B, &a,&b)) return NULL;
4140 [ + + ]: 56 : if (degpol(RgX_gcd(a,b)) > 0) return NULL;
4141 : 56 : return gdiv(a,b);
4142 : : }
4143 : :
4144 : : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
4145 : : * of the continued fraction of x with b <= k maximal */
4146 : : static GEN
4147 : 0 : bestappr_frac(GEN x, GEN k)
4148 : : {
4149 : : pari_sp av;
4150 : : GEN p0, p1, p, q0, q1, q, a, y;
4151 : :
4152 [ # # ]: 0 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
4153 : 0 : av = avma; y = x;
4154 : 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
4155 : 0 : q1 = gen_0; q0 = gen_1;
4156 : 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
4157 : : for(;;)
4158 : : {
4159 : 0 : x = ginv(x); /* > 1 */
4160 [ # # ]: 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
4161 [ # # ]: 0 : if (cmpii(a,k) > 0)
4162 : : { /* next partial quotient will overflow limits */
4163 : : GEN n, d;
4164 : 0 : a = divii(subii(k, q1), q0);
4165 : 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4166 : 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4167 : : /* compare |y-p0/q0|, |y-p1/q1| */
4168 : 0 : n = gel(y,1);
4169 : 0 : d = gel(y,2);
4170 [ # # ]: 0 : if (absi_cmp(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
4171 : : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
4172 : 0 : { p1 = p0; q1 = q0; }
4173 : 0 : break;
4174 : : }
4175 : 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4176 : 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4177 : :
4178 [ # # ]: 0 : if (cmpii(q0,k) > 0) break;
4179 : 0 : x = gsub(x,a); /* 0 <= x < 1 */
4180 [ # # ]: 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
4181 : :
4182 : 0 : }
4183 : 0 : return gerepileupto(av, gdiv(p1,q1));
4184 : : }
4185 : : /* bestappr(t_REAL != 0), to maximal accuracy */
4186 : : static GEN
4187 : 28 : bestappr_real_max(GEN x)
4188 : : {
4189 : 28 : pari_sp av = avma;
4190 : : GEN p0, p1, p, q0, q1, q, a;
4191 : : long e;
4192 : 28 : p1 = gen_1; a = p0 = floorr(x); q1 = gen_0; q0 = gen_1;
4193 : 28 : x = subri(x,a); /* 0 <= x < 1 */
4194 : 28 : e = bit_prec(x) - expo(x);
4195 : : for(;;)
4196 : : {
4197 : : long d;
4198 [ + + ][ - + ]: 55 : if (!signe(x) || expi(q0) > e) { p1 = p0; q1 = q0; break; }
4199 : 34 : x = invr(x); /* > 1 */
4200 : 34 : d = nbits2prec(expo(x) + 1);
4201 [ + + ]: 34 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4202 : :
4203 : 27 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4204 : 27 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4205 : 27 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4206 : 27 : x = subri(x,a); /* 0 <= x < 1 */
4207 : 27 : }
4208 : 28 : return gerepileupto(av, gdiv(p1,q1));
4209 : : }
4210 : : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
4211 : : * of the continued fraction of x with b <= k maximal */
4212 : : static GEN
4213 : 195447 : bestappr_real(GEN x, GEN k)
4214 : : {
4215 : 195447 : pari_sp av = avma;
4216 : : GEN kr, p0, p1, p, q0, q1, q, a, y;
4217 : :
4218 : 195447 : y = x;
4219 : 195447 : p1 = gen_1; a = p0 = floorr(x);
4220 : 195447 : q1 = gen_0; q0 = gen_1;
4221 : 195447 : x = subri(x,a); /* 0 <= x < 1 */
4222 [ + + ]: 195447 : if (!signe(x)) { cgiv(x); return a; }
4223 : 191648 : kr = itor(k, realprec(x));
4224 : : for(;;)
4225 : : {
4226 : : long d;
4227 : 540810 : x = invr(x); /* > 1 */
4228 [ + + ]: 540810 : if (cmprr(x,kr) > 0)
4229 : : { /* next partial quotient will overflow limits */
4230 : 187418 : a = divii(subii(k, q1), q0);
4231 : 187418 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4232 : 187418 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4233 : : /* compare |y-p0/q0|, |y-p1/q1| */
4234 [ + + ]: 187418 : if (absr_cmp(mulir(q1, subri(mulir(q0,y), p0)),
4235 : : mulir(q0, subri(mulir(q1,y), p1))) < 0)
4236 : 3706 : { p1 = p0; q1 = q0; }
4237 : 187418 : break;
4238 : : }
4239 : 353392 : d = nbits2prec(expo(x) + 1);
4240 [ - + ]: 353392 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4241 : :
4242 : 353392 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4243 : 353392 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4244 : 353392 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4245 : :
4246 [ + + ]: 353392 : if (cmpii(q0,k) > 0) break;
4247 : 349345 : x = subri(x,a); /* 0 <= x < 1 */
4248 [ + + ]: 349345 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
4249 : 349162 : }
4250 : 195447 : return gerepileupto(av, gdiv(p1,q1));
4251 : : }
4252 : :
4253 : : /* k t_INT or NULL */
4254 : : static GEN
4255 : 299544 : bestappr_Q(GEN x, GEN k)
4256 : : {
4257 : 299544 : long lx, tx = typ(x), i;
4258 : : GEN a, y;
4259 : :
4260 [ - - + + : 299544 : switch(tx)
+ - + - ]
4261 : : {
4262 : 0 : case t_INT: return icopy(x);
4263 [ # # ]: 0 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
4264 : : case t_REAL:
4265 [ + + ]: 229503 : if (!signe(x)) return gen_0;
4266 [ + + ]: 195475 : return k? bestappr_real(x, k): bestappr_real_max(x);
4267 : :
4268 : : case t_INTMOD: {
4269 : 21 : pari_sp av = avma;
4270 [ + + ]: 21 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
4271 : 14 : return gerepilecopy(av, a);
4272 : : }
4273 : : case t_PADIC: {
4274 : 14 : pari_sp av = avma;
4275 : 14 : long v = valp(x);
4276 [ + + ]: 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
4277 [ - + ]: 7 : if (v) a = gmul(a, powis(gel(x,2), v));
4278 : 7 : return gerepilecopy(av, a);
4279 : : }
4280 : :
4281 : : case t_SER:
4282 [ # # ]: 0 : if (ser_isexactzero(x)) return gcopy(x);
4283 : : /* fall through */
4284 : : case t_COMPLEX: case t_POLMOD: case t_POL: case t_RFRAC:
4285 : : case t_VEC: case t_COL: case t_MAT:
4286 : 70006 : y = cgetg_copy(x, &lx);
4287 [ + - ]: 70006 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
4288 [ + + ]: 361713 : for (; i<lx; i++)
4289 : : {
4290 [ - + ]: 291707 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
4291 : 291707 : gel(y,i) = a;
4292 : : }
4293 [ - + ]: 70006 : if (tx == t_POL) return normalizepol(y);
4294 [ - + ]: 70006 : if (tx == t_SER) return normalize(y);
4295 : 70006 : return y;
4296 : : }
4297 : 0 : pari_err_TYPE("bestappr_Q",x);
4298 : 299544 : return NULL; /* not reached */
4299 : : }
4300 : :
4301 : : static GEN
4302 : 49 : bestappr_ser(GEN x, long B)
4303 : : {
4304 : 49 : long v = valp(x), lx = lg(x);
4305 : : GEN N, t;
4306 : 49 : x = normalizepol(ser2pol_i(x, lx));
4307 : 49 : N = monomial(gen_1, lx-2, varn(x));
4308 [ + + ]: 49 : t = mod_to_rfrac(x, N, B); if (!t) return NULL;
4309 [ + + ]: 42 : if (v)
4310 : : {
4311 : : GEN a, b;
4312 : : long vx;
4313 [ - + ]: 14 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
4314 : : /* t_RFRAC */
4315 : 14 : vx = varn(x);
4316 : 14 : a = gel(t,1);
4317 : 14 : b = gel(t,2);
4318 : 14 : v -= RgX_valrem(b, &b);
4319 [ - + ][ # # ]: 14 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
4320 [ + + ]: 14 : if (v < 0) b = RgX_shift(b, -v);
4321 [ + - ]: 7 : else if (v > 0) {
4322 [ - + ][ # # ]: 7 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
4323 : 7 : a = RgX_shift(a, v);
4324 : : }
4325 : 14 : t = mkrfraccopy(a, b);
4326 : : }
4327 : 49 : return t;
4328 : : }
4329 : : static GEN bestappr_RgX(GEN x, long B);
4330 : : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
4331 : : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
4332 : : static GEN
4333 : 63 : bestappr_RgX(GEN x, long B)
4334 : : {
4335 : 63 : long i, lx, tx = typ(x);
4336 : : GEN y, t;
4337 [ - + + + : 63 : switch(tx)
- - ]
4338 : : {
4339 : : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
4340 : : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
4341 : 0 : return gcopy(x);
4342 : :
4343 : : case t_RFRAC: {
4344 : 14 : pari_sp av = avma;
4345 [ + + ][ - + ]: 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
4346 : 7 : x = rfractoser(x, varn(gel(x,2)), 2*B+1);
4347 [ - + ]: 7 : t = bestappr_ser(x, B); if (!t) return NULL;
4348 : 7 : return gerepileupto(av, t);
4349 : : }
4350 : : case t_POLMOD: {
4351 : 7 : pari_sp av = avma;
4352 [ - + ]: 7 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
4353 : 7 : return gerepileupto(av, t);
4354 : : }
4355 : : case t_SER: {
4356 : 42 : pari_sp av = avma;
4357 [ + + ]: 42 : t = bestappr_ser(x, B); if (!t) return NULL;
4358 : 35 : return gerepileupto(av, t);
4359 : : }
4360 : :
4361 : : case t_VEC: case t_COL: case t_MAT:
4362 : 0 : y = cgetg_copy(x, &lx);
4363 [ # # ]: 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
4364 [ # # ]: 0 : for (; i<lx; i++)
4365 : : {
4366 [ # # ]: 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
4367 : 0 : gel(y,i) = t;
4368 : : }
4369 : 0 : return y;
4370 : : }
4371 : 0 : pari_err_TYPE("bestappr_RgX",x);
4372 : 63 : return NULL; /* not reached */
4373 : : }
4374 : :
4375 : : /* allow k = NULL: maximal accuracy */
4376 : : GEN
4377 : 7837 : bestappr(GEN x, GEN k)
4378 : : {
4379 : 7837 : pari_sp av = avma;
4380 [ + + ]: 7837 : if (k) { /* replace by floor(k) */
4381 [ + + - ]: 7795 : switch(typ(k))
4382 : : {
4383 : : case t_INT:
4384 : 42 : break;
4385 : : case t_REAL: case t_FRAC:
4386 : 7753 : k = floor_safe(k); /* left on stack for efficiency */
4387 [ - + ]: 7753 : if (!signe(k)) k = gen_1;
4388 : 7753 : break;
4389 : : default:
4390 : 0 : pari_err_TYPE("bestappr [bound type]", k);
4391 : 0 : break;
4392 : : }
4393 : : }
4394 : 7837 : x = bestappr_Q(x, k);
4395 [ + + ]: 7837 : if (!x) { avma = av; return cgetg(1,t_VEC); }
4396 : 7837 : return x;
4397 : : }
4398 : : GEN
4399 : 63 : bestapprPade(GEN x, long B)
4400 : : {
4401 : 63 : pari_sp av = avma;
4402 : 63 : GEN t = bestappr_RgX(x, B);
4403 [ + + ]: 63 : if (!t) { avma = av; return cgetg(1,t_VEC); }
4404 : 63 : return t;
4405 : : }
4406 : :
4407 : : /***********************************************************************/
4408 : : /** **/
4409 : : /** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
4410 : : /** **/
4411 : : /***********************************************************************/
4412 : :
4413 : : static GEN
4414 : 14 : get_quad(GEN f, GEN pol, long r)
4415 : : {
4416 : 14 : GEN p1 = gcoeff(f,1,2), q1 = gcoeff(f,2,2);
4417 [ + - ]: 14 : return mkquad(pol, r? subii(p1,q1): p1, q1);
4418 : : }
4419 : :
4420 : : /* replace f by f * [a,1; 1,0] */
4421 : : static void
4422 : 14 : update_f(GEN f, GEN a)
4423 : : {
4424 : : GEN p1;
4425 : 14 : p1 = gcoeff(f,1,1);
4426 : 14 : gcoeff(f,1,1) = addii(mulii(a,p1), gcoeff(f,1,2));
4427 : 14 : gcoeff(f,1,2) = p1;
4428 : :
4429 : 14 : p1 = gcoeff(f,2,1);
4430 : 14 : gcoeff(f,2,1) = addii(mulii(a,p1), gcoeff(f,2,2));
4431 : 14 : gcoeff(f,2,2) = p1;
4432 : 14 : }
4433 : :
4434 : : GEN
4435 : 7 : quadunit(GEN x)
4436 : : {
4437 : 7 : pari_sp av = avma, av2;
4438 : : GEN pol, y, a, u, v, sqd, f;
4439 : : long r;
4440 : :
4441 : 7 : check_quaddisc_real(x, &r, "quadunit");
4442 : 7 : pol = quadpoly(x);
4443 : 7 : sqd = sqrti(x); av2 = avma;
4444 : 7 : a = shifti(addsi(r,sqd),-1);
4445 : 7 : f = mkmat2(mkcol2(a, gen_1), mkcol2(gen_1, gen_0)); /* [a,0; 1,0] */
4446 : 7 : u = stoi(r); v = gen_2;
4447 : : for(;;)
4448 : : {
4449 : : GEN u1, v1;
4450 : 14 : u1 = subii(mulii(a,v),u);
4451 : 14 : v1 = divii(subii(x,sqri(u1)),v);
4452 [ + + ]: 14 : if ( equalii(v,v1) ) {
4453 : 7 : y = get_quad(f,pol,r);
4454 : 7 : update_f(f,a);
4455 : 7 : y = gdiv(get_quad(f,pol,r), gconj(y));
4456 : 7 : break;
4457 : : }
4458 : 7 : a = divii(addii(sqd,u1), v1);
4459 [ - + ]: 7 : if ( equalii(u,u1) ) {
4460 : 0 : y = get_quad(f,pol,r);
4461 : 0 : y = gdiv(y, gconj(y));
4462 : 0 : break;
4463 : : }
4464 : 7 : update_f(f,a);
4465 : 7 : u = u1; v = v1;
4466 [ - + ]: 7 : if (gc_needed(av2,2))
4467 : : {
4468 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit");
4469 : 0 : gerepileall(av2,4, &a,&f,&u,&v);
4470 : : }
4471 : 7 : }
4472 [ + - ]: 7 : if (signe(gel(y,3)) < 0) y = gneg(y);
4473 : 7 : return gerepileupto(av, y);
4474 : : }
4475 : :
4476 : : GEN
4477 : 21 : quadregulator(GEN x, long prec)
4478 : : {
4479 : 21 : pari_sp av = avma, av2;
4480 : : GEN R, rsqd, u, v, sqd;
4481 : : long r, Rexpo;
4482 : :
4483 : 21 : check_quaddisc_real(x, &r, "quadregulator");
4484 : 21 : sqd = sqrti(x);
4485 : 21 : rsqd = gsqrt(x,prec);
4486 : 21 : Rexpo = 0; R = real2n(1, prec); /* = 2 */
4487 : 21 : av2 = avma;
4488 : 21 : u = stoi(r); v = gen_2;
4489 : : for(;;)
4490 : : {
4491 : 70 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
4492 : 70 : GEN v1 = divii(subii(x,sqri(u1)),v);
4493 [ + + ]: 70 : if (equalii(v,v1))
4494 : : {
4495 : 7 : R = sqrr(R); shiftr_inplace(R, -1);
4496 : 7 : R = mulrr(R, divri(addir(u1,rsqd),v));
4497 : 7 : break;
4498 : : }
4499 [ + + ]: 63 : if (equalii(u,u1))
4500 : : {
4501 : 14 : R = sqrr(R); shiftr_inplace(R, -1);
4502 : 14 : break;
4503 : : }
4504 : 49 : R = mulrr(R, divri(addir(u1,rsqd),v));
4505 : 49 : Rexpo += expo(R); setexpo(R,0);
4506 : 49 : u = u1; v = v1;
4507 [ - + ]: 49 : if (Rexpo & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
4508 [ - + ]: 49 : if (gc_needed(av2,2))
4509 : : {
4510 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
4511 : 0 : gerepileall(av2,3, &R,&u,&v);
4512 : : }
4513 : 49 : }
4514 : 21 : R = logr_abs(divri(R,v));
4515 [ + - ]: 21 : if (Rexpo)
4516 : : {
4517 : 21 : GEN t = mulsr(Rexpo, mplog2(prec));
4518 : 21 : shiftr_inplace(t, 1);
4519 : 21 : R = addrr(R,t);
4520 : : }
4521 : 21 : return gerepileuptoleaf(av, R);
4522 : : }
4523 : :
4524 : : /*************************************************************************/
4525 : : /** **/
4526 : : /** CLASS NUMBER **/
4527 : : /** **/
4528 : : /*************************************************************************/
4529 : 1302 : static int qfb_is_1(GEN f) { return equali1(gel(f,1)); }
4530 : 2198 : static GEN qfb_pow(void *E, GEN f, GEN n) { (void)E; return powgi(f,n); }
4531 : 9779 : static GEN qfi_comp(void *E, GEN f, GEN g) { (void)E; return qficomp(f,g); }
4532 : 10437 : static ulong qfb_hash(GEN x)
4533 : : {
4534 : 10437 : GEN a = gel(x,1);
4535 : 10437 : GEN b = gel(x,2);
4536 : 10437 : ulong A = mod2BIL(a);
4537 [ + + ]: 10437 : ulong B = signe(b)? mod2BIL(b): 0;
4538 : 10437 : return (A << BITS_IN_HALFULONG) | B;
4539 : : }
4540 : : static const struct bb_group qfi_group={ qfi_comp,qfb_pow,NULL,qfb_hash,
4541 : : gidentical,qfb_is_1,NULL};
4542 : :
4543 : : GEN
4544 : 133 : qfbclassno0(GEN x,long flag)
4545 : : {
4546 [ + + - ]: 133 : switch(flag)
4547 : : {
4548 : 119 : case 0: return map_proto_G(classno,x);
4549 : 14 : case 1: return map_proto_G(classno2,x);
4550 : 0 : default: pari_err_FLAG("qfbclassno");
4551 : : }
4552 : 133 : return NULL; /* not reached */
4553 : : }
4554 : :
4555 : : /* f^h = 1, return order(f). Set *pfao to its factorization */
4556 : : static GEN
4557 : 441 : find_order(GEN f, GEN h, GEN *pfao)
4558 : : {
4559 : 441 : GEN v = gen_factored_order(f, h, NULL, &qfi_group);
4560 : 441 : *pfao = gel(v,2); return gel(v,1);
4561 : : }
4562 : :
4563 : : static int
4564 : 7 : ok_q(GEN q, GEN h, GEN d2, long r2)
4565 : : {
4566 [ + - ]: 7 : if (d2)
4567 : : {
4568 [ - + ][ # # ]: 7 : if (r2 <= 2 && !mpodd(q)) return 0;
4569 : 7 : return is_pm1(coprime_part(q,d2));
4570 : : }
4571 : : else
4572 : : {
4573 [ # # ][ # # ]: 0 : if (r2 <= 1 && !mpodd(q)) return 0;
4574 : 7 : return is_pm1(coprime_part(q,h));
4575 : : }
4576 : : }
4577 : :
4578 : : /* a,b given by their factorizations. Return factorization of lcm(a,b).
4579 : : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
4580 : : static GEN
4581 : 77 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
4582 : : {
4583 : 77 : GEN P = ZV_union_shallow(gel(Fa,1), gel(Fb,1));
4584 : 77 : GEN A = gen_1, B = gen_1;
4585 : 77 : long i, l = lg(P);
4586 : 77 : GEN E = cgetg(l, t_COL);
4587 [ + + ]: 287 : for (i=1; i<l; i++)
4588 : : {
4589 : 210 : GEN p = gel(P,i);
4590 : 210 : long va = Z_pval(a,p);
4591 : 210 : long vb = Z_pval(b,p);
4592 [ + + ]: 210 : if (va < vb)
4593 : : {
4594 : 77 : B = mulii(B,powiu(p,vb));
4595 : 77 : gel(E,i) = utoi(vb);
4596 : : }
4597 : : else
4598 : : {
4599 : 133 : A = mulii(A,powiu(p,va));
4600 : 133 : gel(E,i) = utoi(va);
4601 : : }
4602 : : }
4603 : 77 : *pA = A;
4604 : 77 : *pB = B; return mkmat2(P,E);
4605 : : }
4606 : :
4607 : : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
4608 : : static void
4609 : 77 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
4610 : : {
4611 : 77 : GEN A,B, g1 = *pg1, d1 = *pd1;
4612 : 77 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
4613 : 77 : *pg1 = gmul(powgi(g1, diviiexact(d1,A)), powgi(f, diviiexact(o,B)));
4614 : 77 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
4615 : 77 : }
4616 : :
4617 : : /* Write x = Df^2, where D = fundamental discriminant,
4618 : : * P^E = factorisation of conductor f, with E[i] >= 0 */
4619 : : static void
4620 : 147 : corediscfact(GEN x, long xmod4, GEN *ptD, GEN *ptP, GEN *ptE)
4621 : : {
4622 : 147 : long s = signe(x), l, i;
4623 : 147 : GEN fa = absi_factor(x);
4624 : 147 : GEN d, P = gel(fa,1), E = gtovecsmall(gel(fa,2));
4625 : :
4626 : 147 : l = lg(P); d = gen_1;
4627 [ + + ]: 518 : for (i=1; i<l; i++)
4628 : : {
4629 [ + + ]: 371 : if (E[i] & 1) d = mulii(d, gel(P,i));
4630 : 371 : E[i] >>= 1;
4631 : : }
4632 [ + + ][ + - ]: 147 : if (!xmod4 && mod4(d) != ((s < 0)? 3: 1)) { d = shifti(d,2); E[1]--; }
[ + - ]
4633 [ + + ]: 147 : *ptD = (s < 0)? negi(d): d;
4634 : 147 : *ptP = P;
4635 : 147 : *ptE = E;
4636 : 147 : }
4637 : :
4638 : : static GEN
4639 : 133 : conductor_part(GEN x, long xmod4, GEN *ptD, GEN *ptreg)
4640 : : {
4641 : 133 : long l, i, s = signe(x);
4642 : : GEN E, H, D, P, reg;
4643 : :
4644 : 133 : corediscfact(x, xmod4, &D, &P, &E);
4645 : 133 : H = gen_1; l = lg(P);
4646 : : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
4647 [ + + ]: 476 : for (i=1; i<l; i++)
4648 : : {
4649 : 343 : long e = E[i];
4650 [ - + ]: 343 : if (e)
4651 : : {
4652 : 0 : GEN p = gel(P,i);
4653 : 0 : H = mulii(H, subis(p, kronecker(D,p)));
4654 [ # # ]: 0 : if (e >= 2) H = mulii(H, powiu(p,e-1));
4655 : : }
4656 : : }
4657 : :
4658 : : /* divide by [ O_K^* : O^* ] */
4659 [ + + ]: 133 : if (s < 0)
4660 : : {
4661 : 119 : reg = NULL;
4662 [ - - + ]: 119 : switch(itou_or_0(D))
4663 : : {
4664 : 0 : case 4: H = shifti(H,-1); break;
4665 : 119 : case 3: H = divis(H,3); break;
4666 : : }
4667 : : } else {
4668 : 14 : reg = quadregulator(D,DEFAULTPREC);
4669 [ - + ]: 14 : if (!equalii(x,D))
4670 : 0 : H = divii(H, roundr(divrr(quadregulator(x,DEFAULTPREC), reg)));
4671 : : }
4672 [ + + ]: 133 : if (ptreg) *ptreg = reg;
4673 : 133 : *ptD = D; return H;
4674 : : }
4675 : :
4676 : : static long
4677 : 112 : two_rank(GEN x)
4678 : : {
4679 : 112 : GEN p = gel(absi_factor(x),1);
4680 : 112 : long l = lg(p)-1;
4681 : : #if 0 /* positive disc not needed */
4682 : : if (signe(x) > 0)
4683 : : {
4684 : : long i;
4685 : : for (i=1; i<=l; i++)
4686 : : if (mod4(gel(p,i)) == 3) { l--; break; }
4687 : : }
4688 : : #endif
4689 : 112 : return l-1;
4690 : : }
4691 : :
4692 : : static GEN
4693 : 2128 : sqr_primeform(GEN x, ulong p) { return redimag(qfisqr(primeform_u(x, p))); }
4694 : : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
4695 : : static GEN
4696 : 112 : get_forms(GEN D, GEN *pL)
4697 : : {
4698 : 112 : const long MAXFORM = 20;
4699 : 112 : GEN L, sqrtD = gsqrt(absi(D),DEFAULTPREC), forms = vectrunc_init(MAXFORM+1);
4700 : 112 : long s, nforms = 0;
4701 : : ulong p;
4702 : : forprime_t S;
4703 : 112 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
4704 : 112 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
4705 [ - + ]: 112 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
4706 [ - + ]: 112 : if (s < 10) s = 200;
4707 [ - + ]: 112 : else if (s < 20) s = 1000;
4708 [ + - ]: 112 : else if (s < 5000) s = 5000;
4709 : 112 : u_forprime_init(&S, 2, s);
4710 [ + + ]: 75040 : while ( (p = u_forprime_next(&S)) )
4711 : : {
4712 : 74928 : long d, k = kroiu(D,p);
4713 : : pari_sp av2;
4714 [ + + ]: 74928 : if (!k) continue;
4715 [ + + ]: 74725 : if (k > 0)
4716 : : {
4717 [ + + ]: 38157 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
4718 : 38157 : d = p-1;
4719 : : }
4720 : : else
4721 : 36568 : d = p+1;
4722 : 74725 : av2 = avma; affrr(divru(mulur(p,L),d), L); avma = av2;
4723 : : }
4724 : 112 : *pL = L; return forms;
4725 : : }
4726 : :
4727 : : /* h ~ #G, return o = order of f, set fao = its factorization */
4728 : : static GEN
4729 : 154 : Shanks_order(GEN f, GEN h, GEN *pfao)
4730 : : {
4731 : 154 : long s = minss(itos(sqrti(h)), 10000);
4732 : 154 : GEN T = gen_Shanks_init(f, s, NULL, &qfi_group);
4733 : 154 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, NULL, &qfi_group);
4734 : 154 : return find_order(f, addiu(v,1), pfao);
4735 : : }
4736 : :
4737 : : /* if g = 1 in G/<f> ? */
4738 : : static int
4739 : 518 : equal1(GEN T, ulong N, GEN g)
4740 : 518 : { return !!gen_Shanks(T, g, N, NULL, &qfi_group); }
4741 : :
4742 : : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
4743 : : * FIXME: should be gen_order, but equal1 has the wrong prototype */
4744 : : static GEN
4745 : 112 : relative_order(GEN a, GEN o, ulong N, GEN T)
4746 : : {
4747 : 112 : pari_sp av = avma;
4748 : : long i, l;
4749 : : GEN m;
4750 : :
4751 : 112 : m = dlog_get_ordfa(o);
4752 [ - + ]: 112 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
4753 : 112 : o = gel(m,1);
4754 : 112 : m = gel(m,2); l = lgcols(m);
4755 [ + + ]: 322 : for (i = l-1; i; i--)
4756 : : {
4757 : 210 : GEN t, y, p = gcoeff(m,i,1);
4758 : 210 : long j, e = itos(gcoeff(m,i,2));
4759 [ + + ]: 210 : if (l == 2) {
4760 : 35 : t = gen_1;
4761 : 35 : y = a;
4762 : : } else {
4763 : 175 : t = diviiexact(o, powiu(p,e));
4764 : 175 : y = powgi(a, t);
4765 : : }
4766 [ + + ]: 210 : if (equal1(T,N,y)) o = t;
4767 : : else {
4768 [ + + ]: 126 : for (j = 1; j < e; j++)
4769 : : {
4770 : 28 : y = powgi(y, p);
4771 [ + + ]: 28 : if (equal1(T,N,y)) break;
4772 : : }
4773 [ + + ]: 119 : if (j < e) {
4774 [ - + ]: 21 : if (j > 1) p = powiu(p, j);
4775 : 21 : o = mulii(t, p);
4776 : : }
4777 : : }
4778 : : }
4779 : 112 : return gerepilecopy(av, o);
4780 : : }
4781 : :
4782 : : /* h(x) for x<0 using Baby Step/Giant Step.
4783 : : * Assumes G is not too far from being cyclic.
4784 : : *
4785 : : * Compute G^2 instead of G so as to kill most of the non-cyclicity */
4786 : : GEN
4787 : 119 : classno(GEN x)
4788 : : {
4789 : 119 : pari_sp av = avma;
4790 : : long r2, k, s, i, l;
4791 : : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
4792 : :
4793 [ + + ]: 119 : if (signe(x) >= 0) return classno2(x);
4794 : :
4795 : 112 : check_quaddisc(x, &s, &k, "classno");
4796 [ - + ]: 112 : if (cmpiu(x,12) <= 0) return gen_1;
4797 : :
4798 : 112 : Hf = conductor_part(x, k, &D, NULL);
4799 [ - + ]: 112 : if (cmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
4800 : 112 : forms = get_forms(D, &L);
4801 : 112 : r2 = two_rank(D);
4802 : 112 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
4803 : :
4804 : 112 : l = lg(forms);
4805 : 112 : order_bound = const_vec(l-1, NULL);
4806 : 112 : g1 = gel(forms,1);
4807 : 112 : gel(order_bound,1) = d1 = Shanks_order(g1, hin, &fad1);
4808 : 112 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
4809 : 112 : d2 = NULL; /* not computed yet */
4810 [ + + ]: 112 : if (is_pm1(q)) goto END;
4811 [ + + ]: 1715 : for (i=2; i < l; i++)
4812 : : {
4813 : 1645 : GEN o, fao, a, F, fd, f = gel(forms,i);
4814 [ + + ]: 1645 : fd = powgi(f, d1); if (is_pm1(gel(fd,1))) continue;
4815 : 77 : F = powgi(fd, q);
4816 : 77 : a = gel(F,1);
4817 [ + + ]: 77 : o = is_pm1(a)? find_order(fd, q, &fao): Shanks_order(fd, q, &fao);
4818 : : /* f^(d1 q) = 1 */
4819 : 77 : fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
4820 : 77 : o = find_order(f, fao, &fao);
4821 : 77 : gel(order_bound,i) = o;
4822 : : /* o = order of f, fao = factor(o) */
4823 : 77 : update_g1(&g1,&d1,&fad1, f,o,fao);
4824 : 77 : q = diviiround(hin, d1);
4825 [ + + ]: 77 : if (is_pm1(q)) goto END;
4826 : : }
4827 : : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
4828 [ + - ]: 70 : if (expi(q) > 3)
4829 : : { /* q large: compute d2, 2nd elt divisor */
4830 : 70 : ulong N, n = 2*itou(sqrti(d1));
4831 : 70 : GEN D = d1, T = gen_Shanks_init(g1, n, NULL, &qfi_group);
4832 : 70 : d2 = gen_1;
4833 : 70 : N = itou( gceil(gdivgs(d1,n)) ); /* order(g1) <= n*N */
4834 [ + + ]: 287 : for (i = 1; i < l; i++)
4835 : : {
4836 : 280 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
4837 [ + + ]: 280 : if (!B) B = find_order(f, fad1, /*junk*/&d);
4838 : 280 : f = powgi(f,d2);
4839 [ + + ]: 280 : if (equal1(T,N,f)) continue;
4840 [ - + ]: 112 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
4841 : : /* f^B = 1 */
4842 : 112 : d = relative_order(f, B, N,T);
4843 : 112 : d2= mulii(d,d2);
4844 : 112 : D = mulii(d1,d2);
4845 : 112 : q = diviiround(hin,D);
4846 [ + + ]: 112 : if (is_pm1(q)) { d1 = D; goto END; }
4847 : : }
4848 : : /* very probably, d2 is the 2nd elementary divisor */
4849 : 7 : d1 = D; /* product of first two elt divisors */
4850 : : }
4851 : : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
4852 : : * 2-rank */
4853 [ - + ]: 7 : if (!ok_q(q,d1,d2,r2))
4854 : : {
4855 : 0 : GEN q0 = q;
4856 : : long d;
4857 [ # # ]: 0 : if (cmpii(mulii(q,d1), hin) < 0)
4858 : : { /* try q = q0+1,-1,+2,-2 */
4859 : 0 : d = 1;
4860 [ # # ][ # # ]: 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
4861 : : }
4862 : : else
4863 : : { /* q0-1,+1,-2,+2 */
4864 : 0 : d = -1;
4865 [ # # ][ # # ]: 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
4866 : : }
4867 : : }
4868 : 7 : d1 = mulii(d1,q);
4869 : :
4870 : : END:
4871 : 119 : return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
4872 : : }
4873 : :
4874 : : /* use Euler products */
4875 : : GEN
4876 : 21 : classno2(GEN x)
4877 : : {
4878 : 21 : pari_sp av = avma;
4879 : 21 : const long prec = DEFAULTPREC;
4880 : : long n, i, r, s;
4881 : : GEN p1, p2, S, p4, p5, p7, Hf, Pi, reg, logd, d, dr, D, half;
4882 : :
4883 : 21 : check_quaddisc(x, &s, &r, "classno2");
4884 [ + + ][ - + ]: 21 : if (s < 0 && cmpiu(x,12) <= 0) return gen_1;
4885 : :
4886 : 21 : Hf = conductor_part(x, r, &D, ®);
4887 [ + + ][ - + ]: 21 : if (s < 0 && cmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
4888 : :
4889 : 21 : Pi = mppi(prec);
4890 : 21 : d = absi(D); dr = itor(d, prec);
4891 : 21 : logd = logr_abs(dr);
4892 : 21 : p1 = sqrtr(divrr(mulir(d,logd), gmul2n(Pi,1)));
4893 [ + + ]: 21 : if (s > 0)
4894 : : {
4895 : 14 : GEN invlogd = invr(logd);
4896 : 14 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
4897 [ - + ]: 14 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
4898 : : }
4899 : 21 : n = itos_or_0( mptrunc(p1) );
4900 [ - + ]: 21 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
4901 : :
4902 : 21 : p4 = divri(Pi,d);
4903 : 21 : p7 = invr(sqrtr_abs(Pi));
4904 : 21 : half = real2n(-1, prec);
4905 [ + + ]: 21 : if (s > 0)
4906 : : { /* i = 1, shortcut */
4907 : 14 : p1 = sqrtr_abs(dr);
4908 : 14 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
4909 : 14 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
4910 [ + + ]: 546 : for (i=2; i<=n; i++)
4911 : : {
4912 [ + + ]: 532 : long k = kroiu(D,i); if (!k) continue;
4913 : 434 : p2 = mulir(sqru(i), p4);
4914 : 434 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
4915 : 434 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
4916 [ + + ]: 434 : S = (k>0)? addrr(S,p5): subrr(S,p5);
4917 : : }
4918 : 14 : S = shiftr(divrr(S,reg),-1);
4919 : : }
4920 : : else
4921 : : { /* i = 1, shortcut */
4922 : 7 : p1 = gdiv(sqrtr_abs(dr), Pi);
4923 : 7 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
4924 : 7 : S = addrr(p5, divrr(p1, mpexp(p4)));
4925 [ + + ]: 952 : for (i=2; i<=n; i++)
4926 : : {
4927 [ - + ]: 945 : long k = kroiu(D,i); if (!k) continue;
4928 : 945 : p2 = mulir(sqru(i), p4);
4929 : 945 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
4930 : 945 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
4931 [ + + ]: 945 : S = (k>0)? addrr(S,p5): subrr(S,p5);
4932 : : }
4933 : : }
4934 : 21 : return gerepileuptoint(av, mulii(Hf, roundr(S)));
4935 : : }
4936 : :
4937 : : static GEN
4938 : 14 : hclassno2(GEN x)
4939 : : {
4940 : : long i, l, s, xmod4;
4941 : : GEN Q, H, D, P, E;
4942 : :
4943 : 14 : x = negi(x);
4944 : 14 : check_quaddisc(x, &s, &xmod4, "hclassno");
4945 : 14 : corediscfact(x, xmod4, &D, &P, &E);
4946 : :
4947 : 14 : Q = quadclassunit0(D, 0, NULL, 0);
4948 : 14 : H = gel(Q,1); l = lg(P);
4949 : :
4950 : : /* H \prod_{p^e||f} (1 + (p^e-1)/(p-1))[ p - (D/p) ] */
4951 [ + + ]: 42 : for (i=1; i<l; i++)
4952 : : {
4953 : 28 : long e = E[i];
4954 [ + + ]: 28 : if (e)
4955 : : {
4956 : 14 : GEN p = gel(P,i), t = subis(p, kronecker(D,p));
4957 [ + + ]: 14 : if (e > 1) t = mulii(t, diviiexact(subis(powiu(p,e), 1), subis(p,1)));
4958 : 14 : H = mulii(H, addsi(1, t));
4959 : : }
4960 : : }
4961 [ - - + ]: 14 : switch( itou_or_0(D) )
4962 : : {
4963 : 0 : case 3: H = gdivgs(H, 3); break;
4964 : 0 : case 4: H = gdivgs(H, 2); break;
4965 : : }
4966 : 14 : return H;
4967 : : }
4968 : :
4969 : : GEN
4970 : 4627 : hclassno(GEN x)
4971 : : {
4972 : : ulong a, b, b2, d, h;
4973 : : long s;
4974 : : int f;
4975 : :
4976 [ - + ]: 4627 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
4977 : 4627 : s = signe(x);
4978 [ - + ]: 4627 : if (s < 0) return gen_0;
4979 [ - + ]: 4627 : if (!s) return gdivgs(gen_1, -12);
4980 : :
4981 [ + - ][ - + ]: 4627 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
4982 : :
4983 : 4627 : d = itou_or_0(x);
4984 [ + - ][ + + ]: 4627 : if (!d || d > 500000) return hclassno2(x);
4985 : :
4986 : 4613 : h = 0; b = d&1; b2 = (1+d)>>2; f=0;
4987 [ + + ]: 4613 : if (!b)
4988 : : {
4989 [ + + ]: 230013 : for (a=1; a*a<b2; a++)
4990 [ + + ]: 227381 : if (b2%a == 0) h++;
4991 : 2632 : f = (a*a==b2); b=2; b2=(4+d)>>2;
4992 : : }
4993 [ + + ]: 239008 : while (b2*3 < d)
4994 : : {
4995 [ + + ]: 234395 : if (b2%b == 0) h++;
4996 [ + + ]: 18017167 : for (a=b+1; a*a < b2; a++)
4997 [ + + ]: 17782772 : if (b2%a == 0) h += 2;
4998 [ + + ]: 234395 : if (a*a == b2) h++;
4999 : 234395 : b += 2; b2 = (b*b+d)>>2;
5000 : : }
5001 [ + + ]: 4613 : if (b2*3 == d)
5002 : : {
5003 : 105 : GEN y = cgetg(3,t_FRAC);
5004 : 105 : gel(y,1) = utoipos(3*h+1);
5005 : 105 : gel(y,2) = utoipos(3); return y;
5006 : : }
5007 [ + + ]: 4508 : if (f)
5008 : : {
5009 : 42 : GEN y = cgetg(3,t_FRAC);
5010 : 42 : gel(y,1) = utoipos(2*h+1);
5011 : 42 : gel(y,2) = gen_2; return y;
5012 : : }
5013 : 4627 : return utoipos(h);
5014 : : }
5015 : :
|
