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 43 : remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
30 : static ulong
31 42623 : u_remove2(ulong q) { return q >> vals(q); }
32 : GEN
33 43 : odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
34 : static GEN
35 42623 : 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 244 : is_gener_expo(GEN p, GEN L0)
40 : {
41 244 : GEN L, q = shifti(p,-1);
42 : long i, l;
43 244 : if (L0) {
44 224 : l = lg(L0);
45 224 : L = cgetg(l, t_VEC);
46 : } else {
47 20 : L0 = L = odd_prime_divisors(q);
48 20 : l = lg(L);
49 : }
50 244 : for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
51 244 : return L;
52 : }
53 : static GEN
54 41970 : u_is_gener_expo(ulong p, GEN L0)
55 : {
56 41970 : const ulong q = p >> 1;
57 : long i, l;
58 : GEN L;
59 41970 : if (L0) {
60 3498 : l = lg(L0);
61 3498 : L = cgetg(l, t_VECSMALL);
62 : } else {
63 38472 : L0 = L = u_odd_prime_divisors(q);
64 38475 : l = lg(L);
65 : }
66 41974 : for (i=1; i<l; i++) L[i] = q / uel(L0,i);
67 41974 : return L;
68 : }
69 :
70 : int
71 104998 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
72 : {
73 : long i;
74 104998 : if (krouu(x, p) >= 0) return 0;
75 119504 : for (i=lg(L)-1; i; i--)
76 : {
77 74902 : ulong t = Fl_powu(x, uel(L,i), p);
78 74918 : if (t == p_1 || t == 1) return 0;
79 : }
80 44602 : return 1;
81 : }
82 : /* assume p prime */
83 : ulong
84 136567 : pgener_Fl_local(ulong p, GEN L0)
85 : {
86 136567 : const pari_sp av = avma;
87 136567 : const ulong p_1 = p-1;
88 : long x;
89 : GEN L;
90 136567 : if (p <= 19) switch(p)
91 : { /* quick trivial cases */
92 21 : case 2: return 1;
93 : case 7:
94 17126 : case 17: return 3;
95 77479 : default: return 2;
96 : }
97 41941 : L = u_is_gener_expo(p,L0);
98 41944 : for (x=2;;x++) { if (is_gener_Fl(x,p,p_1,L)) break; }
99 41960 : avma = av; return x;
100 : }
101 : ulong
102 112495 : 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 524 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
108 : {
109 524 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
110 524 : if (t >= 0) return 0;
111 477 : for (i = lg(L)-1; i; i--)
112 : {
113 182 : GEN t = Fp_pow(x, gel(L,i), p);
114 182 : if (equalii(t, p_1) || equali1(t)) return 0;
115 : }
116 295 : return 1;
117 : }
118 :
119 : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
120 : GEN
121 24423 : pgener_Fp_local(GEN p, GEN L0)
122 : {
123 24423 : pari_sp av0 = avma;
124 : GEN x, p_1, L;
125 24423 : if (lgefint(p) == 3)
126 : {
127 : ulong z;
128 24184 : if (p[2] == 2) return gen_1;
129 18185 : if (L0) L0 = ZV_to_nv(L0);
130 18185 : z = pgener_Fl_local(uel(p,2), L0);
131 18185 : avma = av0; return utoipos(z);
132 : }
133 239 : p_1 = subiu(p,1); L = is_gener_expo(p, L0);
134 239 : x = utoipos(2);
135 452 : for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
136 239 : avma = av0; return utoipos(uel(x,2));
137 : }
138 :
139 : GEN
140 23590 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
141 :
142 : ulong
143 69204 : pgener_Zl(ulong p)
144 : {
145 69204 : 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 69204 : if (p == 40487) return 10;
148 : #ifndef LONG_IS_64BIT
149 9882 : return pgener_Fl(p);
150 : #else
151 59322 : if (p < (1UL<<32)) return pgener_Fl(p);
152 : else
153 : {
154 30 : const pari_sp av = avma;
155 30 : const ulong p_1 = p-1;
156 : long x ;
157 30 : GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
158 102 : for (x=2;;x++)
159 174 : if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2))) break;
160 30 : avma = av; return x;
161 : }
162 : #endif
163 : }
164 :
165 : /* p prime. Return a primitive root modulo p^e, e > 1 */
166 : GEN
167 69209 : pgener_Zp(GEN p)
168 : {
169 69209 : if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
170 : else
171 : {
172 5 : const pari_sp av = avma;
173 5 : GEN p_1 = subiu(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
174 5 : GEN x = utoipos(2);
175 12 : for (;; x[2]++)
176 29 : if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
177 5 : avma = av; return utoipos(uel(x,2));
178 : }
179 : }
180 :
181 : static GEN
182 231 : gener_Zp(GEN q, GEN F)
183 : {
184 231 : GEN p = NULL;
185 231 : long e = 0;
186 231 : 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 (absequaliu(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 217 : e = Z_isanypower(q, &p);
201 231 : return e > 1? pgener_Zp(p): pgener_Fp(q);
202 : }
203 :
204 : GEN
205 301 : znprimroot(GEN N)
206 : {
207 301 : pari_sp av = avma;
208 : GEN x, n, F;
209 :
210 301 : 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 294 : N = absi_shallow(N);
216 294 : if (abscmpiu(N, 4) <= 0) { avma = av; return mkintmodu(N[2]-1,N[2]); }
217 245 : 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 210 : x = gener_Zp(N,F);
228 210 : break;
229 : }
230 231 : 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(subiu(p,1), n), p); /* prim. n-th root of 1 */
241 0 : return gerepileuptoint(av, z);
242 : }
243 :
244 : GEN
245 217 : rootsof1u_Fp(ulong n, GEN p)
246 : {
247 217 : pari_sp av = avma;
248 217 : GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
249 217 : z = pgener_Fp_local(p, Flv_to_ZV(L));
250 217 : z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
251 217 : return gerepileuptoint(av, z);
252 : }
253 :
254 : ulong
255 3934 : rootsof1_Fl(ulong n, ulong p)
256 : {
257 3934 : pari_sp av = avma;
258 3934 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
259 3934 : ulong z = pgener_Fl_local(p, L);
260 3934 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
261 3934 : avma = av; return z;
262 : }
263 :
264 : /*********************************************************************/
265 : /** **/
266 : /** INVERSE TOTIENT FUNCTION **/
267 : /** **/
268 : /*********************************************************************/
269 : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
270 : * primes. Return factor(N) */
271 : GEN
272 350651 : Z_factor_listP(GEN N, GEN L)
273 : {
274 350651 : long i, k, l = lg(L);
275 350651 : GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
276 1346688 : for (i = k = 1; i < l; i++)
277 : {
278 996037 : GEN p = gel(L,i);
279 996037 : long v = Z_pvalrem(N, p, &N);
280 996037 : if (v)
281 : {
282 792176 : gel(P,k) = p;
283 792176 : gel(E,k) = utoipos(v);
284 792176 : k++;
285 : }
286 : }
287 350651 : setlg(P, k);
288 350651 : setlg(E, k); return mkmat2(P,E);
289 : }
290 :
291 : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
292 : * L is a list of primes containing all prime divisors of n. */
293 : static long
294 621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
295 : {
296 621565 : pari_sp av = avma, av2;
297 : GEN k, D;
298 : long i, v;
299 621565 : if (m && mod2(n))
300 : {
301 270914 : if (!equali1(n)) return 0;
302 69986 : if (px) *px = gen_1;
303 69986 : return 1;
304 : }
305 350651 : D = divisors(Z_factor_listP(shifti(n, -1), L));
306 : /* loop through primes p > m, d = p-1 | n */
307 350651 : av2 = avma;
308 350651 : if (!m)
309 : { /* special case p = 2, d = 1 */
310 69986 : k = n;
311 69986 : for (v = 1;; v++) {
312 69986 : if (istotient_i(k, gen_2, L, px)) {
313 69986 : if (px) *px = shifti(*px, v);
314 69986 : return 1;
315 : }
316 0 : if (mod2(k)) break;
317 0 : k = shifti(k,-1);
318 : }
319 0 : avma = av2;
320 : }
321 1099462 : for (i = 1; i < lg(D); ++i)
322 : {
323 1001588 : GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
324 1001588 : if (m && cmpii(d, m) < 0) continue;
325 677782 : p = addiu(d, 1);
326 677782 : if (!isprime(p)) continue;
327 442064 : k = diviiexact(n, d);
328 481593 : for (v = 1;; v++) {
329 : GEN r;
330 481593 : if (istotient_i(k, p, L, px)) {
331 182791 : if (px) *px = mulii(*px, powiu(p, v));
332 182791 : return 1;
333 : }
334 298802 : k = dvmdii(k, p, &r);
335 298802 : if (r != gen_0) break;
336 : }
337 259273 : avma = av2;
338 : }
339 97874 : avma = av; return 0;
340 : }
341 :
342 : /* find x such that phi(x) = n */
343 : long
344 70000 : istotient(GEN n, GEN *px)
345 : {
346 70000 : pari_sp av = avma;
347 70000 : if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
348 70000 : if (signe(n) < 1) return 0;
349 70000 : if (mod2(n))
350 : {
351 14 : if (!equali1(n)) return 0;
352 14 : if (px) *px = gen_1;
353 14 : return 1;
354 : }
355 69986 : if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
356 : {
357 69986 : if (!px) avma = av;
358 : else
359 69986 : *px = gerepileuptoint(av, *px);
360 69986 : return 1;
361 : }
362 0 : avma = av; return 0;
363 : }
364 :
365 : /*********************************************************************/
366 : /** **/
367 : /** INTEGRAL LOGARITHM **/
368 : /** **/
369 : /*********************************************************************/
370 :
371 : /* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
372 : * e = floor(log_y B). Set *ptq = y^e if non-NULL */
373 : long
374 253503 : ulogintall(ulong B, ulong y, ulong *ptq)
375 : {
376 : ulong r, r2;
377 : long e;
378 :
379 253503 : if (y == 2)
380 : {
381 4685 : long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
382 4685 : if (ptq) *ptq = 1UL << eB;
383 4685 : return eB;
384 : }
385 248818 : r = y, r2 = 1UL;
386 899962 : for (e=1;; e++)
387 : { /* here, r = y^e, r2 = y^(e-1) */
388 1551106 : if (r >= B)
389 : {
390 248707 : if (r != B) { e--; r = r2; }
391 248707 : if (ptq) *ptq = r;
392 248707 : return e;
393 : }
394 651255 : r2 = r;
395 651255 : r = umuluu_or_0(y, r);
396 651255 : if (!r)
397 : {
398 111 : if (ptq) *ptq = r2;
399 111 : return e;
400 : }
401 : }
402 : }
403 :
404 : /* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
405 : * e = floor(log_y B). Set *ptq = y^e if non-NULL */
406 : long
407 261183 : logintall(GEN B, GEN y, GEN *ptq)
408 : {
409 : pari_sp av;
410 261183 : long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
411 : GEN q, pow2;
412 :
413 261183 : if (lgefint(B) == 3)
414 : {
415 : ulong q;
416 253503 : if (lgefint(y) > 3)
417 : {
418 0 : if (ptq) *ptq = gen_1;
419 0 : return 0;
420 : }
421 253503 : if (!ptq) return ulogintall(B[2], y[2], NULL);
422 46929 : e = ulogintall(B[2], y[2], &q);
423 46929 : *ptq = utoi(q); return e;
424 : }
425 7680 : if (equaliu(y,2))
426 : {
427 159 : if (ptq) *ptq = int2n(eB);
428 159 : return eB;
429 : }
430 7521 : av = avma;
431 7521 : ey = expi(y);
432 : /* eB/(ey+1) - 1 < e <= eB/ey */
433 7521 : emax = eB/ey;
434 7521 : if (emax <= 13) /* e small, be naive */
435 : {
436 1492 : GEN r = y, r2 = gen_1;
437 15831 : for (e=1;; e++)
438 14339 : { /* here, r = y^e, r2 = y^(e-1) */
439 15831 : long fl = cmpii(r, B);
440 15831 : if (fl >= 0)
441 : {
442 1492 : if (fl) { e--; cgiv(r); r = r2; }
443 1492 : if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
444 1492 : return e;
445 : }
446 14339 : r2 = r; r = mulii(r,y);
447 : }
448 : }
449 : /* e >= 13 ey / (ey+1) >= 6.5 */
450 :
451 : /* binary splitting: compute bits of e one by one */
452 : /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
453 6029 : pow2 = new_chunk((long)log2(eB)+2);
454 6029 : gel(pow2,0) = y;
455 6029 : for (i=0, q=y;; )
456 28887 : {
457 34916 : GEN r = gel(pow2,i); /* r = y^2^i */
458 34916 : long fl = cmpii(r,B);
459 34916 : if (!fl)
460 : {
461 0 : e = 1L<<i;
462 0 : if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
463 0 : return e;
464 : }
465 34916 : if (fl > 0) { i--; break; }
466 32924 : q = r;
467 32924 : if (1L<<(i+1) > emax) break;
468 28887 : gel(pow2,++i) = sqri(q);
469 : }
470 :
471 6029 : for (e = 1L<<i;;)
472 26874 : { /* y^e = q < B < r = q * y^(2^i) */
473 32903 : pari_sp av2 = avma;
474 : long fl;
475 : GEN r;
476 32903 : if (--i < 0) break;
477 26881 : r = mulii(q, gel(pow2,i));
478 26881 : fl = cmpii(r, B);
479 26881 : if (fl > 0) avma = av2;
480 : else
481 : {
482 13036 : e += (1L<<i);
483 13036 : q = r;
484 13036 : if (!fl) break; /* B = r */
485 : }
486 : }
487 6029 : if (ptq) *ptq = gerepileuptoint(av, q); else avma = av;
488 6029 : return e;
489 : }
490 :
491 : long
492 56 : logint0(GEN B, GEN y, GEN *ptq)
493 : {
494 56 : if (typ(B) != t_INT) pari_err_TYPE("logint",B);
495 56 : if (signe(B) <= 0) pari_err_DOMAIN("logint", "x" ,"<=", gen_0, B);
496 56 : if (typ(y) != t_INT) pari_err_TYPE("logint",y);
497 56 : if (cmpis(y, 2) < 0) pari_err_DOMAIN("logint", "b" ,"<=", gen_1, y);
498 56 : return logintall(B,y,ptq);
499 : }
500 :
501 : /*********************************************************************/
502 : /** **/
503 : /** INTEGRAL SQUARE ROOT **/
504 : /** **/
505 : /*********************************************************************/
506 : GEN
507 30577 : sqrtint(GEN a)
508 : {
509 30577 : if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
510 30577 : switch (signe(a))
511 : {
512 30563 : case 1: return sqrti(a);
513 7 : case 0: return gen_0;
514 7 : default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
515 : }
516 : return NULL; /* LCOV_EXCL_LINE */
517 : }
518 :
519 : /*********************************************************************/
520 : /** **/
521 : /** PERFECT SQUARE **/
522 : /** **/
523 : /*********************************************************************/
524 : static int
525 13217157 : carremod(ulong A)
526 : {
527 13217157 : const int carresmod64[]={
528 : 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,
529 : 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};
530 13217157 : const int carresmod63[]={
531 : 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,
532 : 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};
533 13217157 : const int carresmod65[]={
534 : 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,
535 : 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};
536 13217157 : const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
537 13217157 : return (carresmod64[A & 0x3fUL]
538 5146443 : && carresmod63[A % 63UL]
539 3124032 : && carresmod65[A % 65UL]
540 15818194 : && carresmod11[A % 11UL]);
541 : }
542 :
543 : /* emulate Z_issquareall on single-word integers */
544 : long
545 11726515 : uissquareall(ulong A, ulong *sqrtA)
546 : {
547 11726515 : if (!A) { *sqrtA = 0; return 1; }
548 11726515 : if (carremod(A))
549 : {
550 1761212 : ulong a = usqrt(A);
551 1761206 : if (a * a == A) { *sqrtA = a; return 1; }
552 : }
553 10056335 : return 0;
554 : }
555 : long
556 120633 : uissquare(ulong A)
557 : {
558 120633 : if (!A) return 1;
559 120633 : if (carremod(A))
560 : {
561 3499 : ulong a = usqrt(A);
562 3499 : if (a * a == A) return 1;
563 : }
564 117158 : return 0;
565 : }
566 :
567 : long
568 6302343 : Z_issquareall(GEN x, GEN *pt)
569 : {
570 : pari_sp av;
571 : GEN y, r;
572 :
573 6302343 : switch(signe(x))
574 : {
575 2206566 : case -1: return 0;
576 1064 : case 0: if (pt) *pt=gen_0; return 1;
577 : }
578 4094713 : if (lgefint(x) == 3)
579 : {
580 2724727 : ulong u = uel(x,2), a;
581 2724727 : if (!pt) return uissquare(u);
582 2604094 : if (!uissquareall(u, &a)) return 0;
583 1392168 : *pt = utoipos(a); return 1;
584 : }
585 1369986 : if (!carremod(umodiu(x, 64*63*65*11))) return 0;
586 608705 : av = avma; y = sqrtremi(x, &r);
587 608705 : if (r != gen_0) { avma = av; return 0; }
588 17861 : if (pt) { *pt = y; avma = (pari_sp)y; } else avma = av;
589 17861 : return 1;
590 : }
591 :
592 : /* a t_INT, p prime */
593 : long
594 0 : Zp_issquare(GEN a, GEN p)
595 : {
596 : long v;
597 : GEN ap;
598 :
599 0 : if (!signe(a) || gequal1(a)) return 1;
600 0 : v = Z_pvalrem(a, p, &ap);
601 0 : if (v&1) return 0;
602 0 : return absequaliu(p, 2)? umodiu(ap, 8) == 1
603 0 : : kronecker(ap,p) == 1;
604 : }
605 :
606 : static long
607 3010 : polissquareall(GEN x, GEN *pt)
608 : {
609 : pari_sp av;
610 : long v;
611 : GEN y, a, b, p;
612 :
613 3010 : if (!signe(x))
614 : {
615 7 : if (pt) *pt = gcopy(x);
616 7 : return 1;
617 : }
618 3003 : if (odd(degpol(x))) return 0; /* odd degree */
619 2177 : av = avma;
620 2177 : v = RgX_valrem(x, &x);
621 2177 : if (v & 1) { avma = av; return 0; }
622 2170 : a = gel(x,2); /* test constant coeff */
623 2170 : if (!pt)
624 70 : { if (!issquare(a)) { avma = av; return 0; } }
625 : else
626 2100 : { if (!issquareall(a,&b)) { avma = av; return 0; } }
627 2170 : if (!degpol(x)) { /* constant polynomial */
628 77 : if (!pt) { avma = av; return 1; }
629 35 : y = scalarpol(b, varn(x)); goto END;
630 : }
631 2093 : p = characteristic(x);
632 2093 : if (signe(p) && !mod2(p))
633 : {
634 : long i, lx;
635 35 : if (!absequaliu(p,2)) pari_err_IMPL("issquare for even characteristic != 2");
636 28 : x = gmul(x, mkintmod(gen_1, gen_2));
637 28 : lx = lg(x);
638 28 : if ((lx-3) & 1) { avma = av; return 0; }
639 49 : for (i = 3; i < lx; i+=2)
640 28 : if (!gequal0(gel(x,i))) { avma = av; return 0; }
641 21 : if (pt) {
642 14 : y = cgetg((lx+3) / 2, t_POL);
643 49 : for (i = 2; i < lx; i+=2)
644 35 : if (!issquareall(gel(x,i), &gel(y, (i+2)>>1))) { avma = av; return 0; }
645 14 : y[1] = evalsigne(1) | evalvarn(varn(x));
646 14 : goto END;
647 : } else {
648 21 : for (i = 2; i < lx; i+=2)
649 14 : if (!issquare(gel(x,i))) { avma = av; return 0; }
650 7 : avma = av; return 1;
651 : }
652 : }
653 : else
654 : {
655 2058 : long m = 1;
656 2058 : x = RgX_Rg_div(x,a);
657 : /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
658 2058 : if (!signe(p)) x = RgX_deflate_max(x,&m);
659 2058 : y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
660 3325 : if (!RgX_equal(RgX_sqr(y), x)) { avma = av; return 0; }
661 798 : if (!pt) { avma = av; return 1; }
662 791 : if (!gequal1(a)) y = gmul(b, y);
663 791 : if (m != 1) y = RgX_inflate(y,m);
664 : }
665 : END:
666 840 : if (v) y = RgX_shift_shallow(y, v>>1);
667 840 : *pt = gerepilecopy(av, y); return 1;
668 : }
669 :
670 : /* b unit mod p */
671 : static int
672 287 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
673 : {
674 287 : if (d == 1)
675 : { /* mod p: faster */
676 203 : if (!Fp_ispower(b, K, p)) return 0;
677 203 : if (pt) *pt = Fp_sqrtn(b, K, p, NULL);
678 : }
679 : else
680 : { /* mod p^{2 +} */
681 84 : if (!ispower(cvtop(b, p, d), K, pt)) return 0;
682 63 : if (pt) *pt = gtrunc(*pt);
683 : }
684 266 : return 1;
685 : }
686 :
687 : /* We're studying whether a mod (q*p^e) is a K-th power, (q,p) = 1.
688 : * Decide mod p^e, then reduce a mod q unless q = NULL. */
689 : static int
690 427 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
691 : {
692 : GEN t, A;
693 427 : long v = Z_pvalrem(*pa, p, &A), d = e - v;
694 427 : if (d <= 0) t = gen_0;
695 : else
696 : {
697 : ulong r;
698 371 : v = uabsdivui_rem(v, K, &r);
699 371 : if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
700 266 : if (L && v) t = mulii(t, powiu(p, v));
701 : }
702 322 : if (q) *pa = modii(*pa, q);
703 322 : if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
704 322 : return 1;
705 : }
706 : long
707 329 : Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
708 : {
709 : GEN L, N;
710 : pari_sp av;
711 : long e, i, l;
712 : ulong pp;
713 : forprime_t S;
714 :
715 329 : if (!signe(a))
716 : {
717 21 : if (pt) {
718 21 : GEN t = cgetg(3, t_INTMOD);
719 21 : gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
720 : }
721 21 : return 1;
722 : }
723 : /* a != 0 */
724 308 : av = avma;
725 :
726 308 : if (typ(q) != t_INT) /* integer factorization */
727 : {
728 0 : GEN P = gel(q,1), E = gel(q,2);
729 0 : l = lg(P);
730 0 : L = pt? vectrunc_init(l): NULL;
731 0 : for (i = 1; i < l; i++)
732 : {
733 0 : GEN p = gel(P,i);
734 0 : long e = itos(gel(E,i));
735 0 : if (!handle_pe(&a, NULL, L, K, p, e)) { avma = av; return 0; }
736 : }
737 0 : goto END;
738 : }
739 308 : if (!mod2(K)
740 189 : && kronecker(a, shifti(q,-vali(q))) == -1) { avma = av; return 0; }
741 301 : L = pt? vectrunc_init(expi(q)+1): NULL;
742 301 : u_forprime_init(&S, 2, tridiv_bound(q));
743 301 : while ((pp = u_forprime_next(&S)))
744 : {
745 : int stop;
746 883407 : e = Z_lvalrem_stop(&q, pp, &stop);
747 883407 : if (!e) continue;
748 203 : if (!handle_pe(&a, q, L, K, utoipos(pp), e)) { avma = av; return 0; }
749 161 : if (stop)
750 : {
751 126 : if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
752 126 : goto END;
753 : }
754 : }
755 154 : l = lg(primetab);
756 154 : for (i = 1; i < l; i++)
757 : {
758 0 : GEN p = gel(primetab,i);
759 0 : e = Z_pvalrem(q, p, &q);
760 0 : if (!e) continue;
761 0 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
762 0 : if (is_pm1(q)) goto END;
763 : }
764 154 : N = gcdii(a,q);
765 154 : if (!is_pm1(N))
766 : {
767 112 : if (ifac_isprime(N))
768 : {
769 70 : e = Z_pvalrem(q, N, &q);
770 70 : if (!handle_pe(&a, q, L, K, N, e)) { avma = av; return 0; }
771 : }
772 : else
773 : {
774 42 : GEN part = ifac_start(N, 0);
775 : for(;;)
776 42 : {
777 : long e;
778 : GEN p;
779 84 : if (!ifac_next(&part, &p, &e)) break;
780 42 : e = Z_pvalrem(q, p, &q);
781 42 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
782 : }
783 : }
784 : }
785 84 : if (!is_pm1(q))
786 : {
787 84 : if (ifac_isprime(q))
788 : {
789 28 : if (!handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
790 : }
791 : else
792 : {
793 56 : GEN part = ifac_start(q, 0);
794 : for(;;)
795 84 : {
796 : long e;
797 : GEN p;
798 140 : if (!ifac_next(&part, &p, &e)) break;
799 98 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
800 : }
801 : }
802 : }
803 : END:
804 196 : if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
805 196 : return 1;
806 : }
807 :
808 : static long
809 49 : polmodispower(GEN x, GEN K, GEN *pt)
810 : {
811 49 : pari_sp av = avma;
812 49 : GEN p = NULL, T = NULL;
813 49 : if (Rg_is_FpXQ(x, &T,&p) && p)
814 : {
815 35 : x = liftall_shallow(x);
816 35 : if (!Fq_ispower(x, K, T, p)) { avma = av; return 0; }
817 28 : if (!pt) { avma = av; return 1; }
818 21 : x = Fq_sqrtn(x, K, T,p, NULL);
819 21 : if (typ(x) == t_INT)
820 7 : x = Fp_to_mod(x,p);
821 : else
822 14 : x = mkpolmod(FpX_to_mod(x,p), FpX_to_mod(T,p));
823 21 : *pt = gerepilecopy(av, x); return 1;
824 : }
825 14 : pari_err_IMPL("ispower for general t_POLMOD");
826 0 : return 0;
827 : }
828 :
829 : long
830 163142 : issquareall(GEN x, GEN *pt)
831 : {
832 163142 : long tx = typ(x);
833 : GEN F;
834 : pari_sp av;
835 :
836 163142 : if (!pt) return issquare(x);
837 20300 : switch(tx)
838 : {
839 2331 : case t_INT: return Z_issquareall(x, pt);
840 161 : case t_FRAC: av = avma;
841 161 : F = cgetg(3, t_FRAC);
842 161 : if ( !Z_issquareall(gel(x,1), &gel(F,1))
843 105 : || !Z_issquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
844 105 : *pt = F; return 1;
845 :
846 : case t_POLMOD:
847 21 : return polmodispower(x, gen_2, pt);
848 2919 : case t_POL: return polissquareall(x,pt);
849 14 : case t_RFRAC: av = avma;
850 14 : F = cgetg(3, t_RFRAC);
851 14 : if ( !issquareall(gel(x,1), &gel(F,1))
852 14 : || !polissquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
853 7 : *pt = F; return 1;
854 :
855 : case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
856 14756 : if (!issquare(x)) return 0;
857 14756 : *pt = gsqrt(x, DEFAULTPREC); return 1;
858 :
859 : case t_INTMOD:
860 63 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
861 :
862 35 : case t_FFELT: return FF_issquareall(x, pt);
863 :
864 : }
865 0 : pari_err_TYPE("issquareall",x);
866 : return 0; /* LCOV_EXCL_LINE */
867 : }
868 :
869 : long
870 157871 : issquare(GEN x)
871 : {
872 : pari_sp av;
873 : GEN a, p;
874 : long i, v;
875 :
876 157871 : switch(typ(x))
877 : {
878 : case t_INT:
879 142751 : return Z_issquare(x);
880 :
881 : case t_REAL:
882 14714 : return (signe(x)>=0);
883 :
884 : case t_INTMOD:
885 77 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, NULL);
886 :
887 : case t_FRAC:
888 21 : return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
889 :
890 7 : case t_FFELT: return FF_issquareall(x, NULL);
891 :
892 : case t_COMPLEX:
893 56 : return 1;
894 :
895 : case t_PADIC:
896 126 : a = gel(x,4); if (!signe(a)) return 1;
897 126 : if (valp(x)&1) return 0;
898 112 : p = gel(x,2);
899 112 : if (!absequaliu(p, 2)) return (kronecker(a,p) != -1);
900 :
901 42 : v = precp(x); /* here p=2, a is odd */
902 42 : if ((v>=3 && mod8(a) != 1 ) ||
903 0 : (v==2 && mod4(a) != 1)) return 0;
904 21 : return 1;
905 :
906 : case t_POLMOD:
907 14 : return polmodispower(x, gen_2, NULL);
908 :
909 : case t_POL:
910 77 : return polissquareall(x,NULL);
911 :
912 : case t_SER:
913 21 : if (!signe(x)) return 1;
914 14 : if (valp(x)&1) return 0;
915 7 : return issquare(gel(x,2));
916 :
917 : case t_RFRAC:
918 7 : av = avma; i = issquare(gmul(gel(x,1),gel(x,2)));
919 7 : avma = av; return i;
920 : }
921 0 : pari_err_TYPE("issquare",x);
922 : return 0; /* LCOV_EXCL_LINE */
923 : }
924 0 : GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
925 0 : GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
926 :
927 : long
928 1386 : ispolygonal(GEN x, GEN S, GEN *N)
929 : {
930 1386 : pari_sp av = avma;
931 : GEN D, d, n;
932 1386 : if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
933 1386 : if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
934 1386 : if (abscmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
935 1386 : if (signe(x) < 0) return 0;
936 1386 : if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
937 1260 : if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
938 : /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
939 1134 : if (abscmpiu(S, 1<<16) < 0) /* common case ! */
940 : {
941 441 : ulong s = S[2], r;
942 504 : if (s == 4) return Z_issquareall(x, N);
943 378 : if (s == 3)
944 0 : D = addiu(shifti(x, 3), 1);
945 : else
946 378 : D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
947 378 : if (!Z_issquareall(D, &d)) { avma = av; return 0; }
948 378 : if (s == 3)
949 0 : d = subiu(d, 1);
950 : else
951 378 : d = addiu(d, s - 4);
952 378 : n = absdiviu_rem(d, 2*s - 4, &r);
953 378 : if (r) { avma = av; return 0; }
954 : }
955 : else
956 : {
957 693 : GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
958 693 : D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
959 693 : if (!Z_issquareall(D, &d)) { avma = av; return 0; }
960 693 : d = addii(d, S_4);
961 693 : n = dvmdii(shifti(d,-1), S_2, &r);
962 693 : if (r != gen_0) { avma = av; return 0; }
963 : }
964 1071 : if (N) *N = gerepileuptoint(av, n); else avma = av;
965 1071 : return 1;
966 : }
967 :
968 : /*********************************************************************/
969 : /** **/
970 : /** PERFECT POWER **/
971 : /** **/
972 : /*********************************************************************/
973 : static long
974 721 : polispower(GEN x, GEN K, GEN *pt)
975 : {
976 : pari_sp av;
977 721 : long v, d, k = itos(K);
978 : GEN y, a, b;
979 721 : GEN T = NULL, p = NULL;
980 :
981 721 : if (!signe(x))
982 : {
983 7 : if (pt) *pt = gcopy(x);
984 7 : return 1;
985 : }
986 714 : d = degpol(x);
987 714 : if (d % k) return 0; /* degree not multiple of k */
988 707 : av = avma;
989 707 : if (RgX_is_FpXQX(x, &T, &p) && p)
990 : { /* over Fq */
991 336 : if (T && typ(T) == t_FFELT)
992 : {
993 126 : if (!FFX_ispower(x, k, T, pt)) { avma = av; return 0; }
994 105 : return 1;
995 : }
996 210 : x = RgX_to_FqX(x,T,p);
997 210 : if (!FqX_ispower(x, k, T,p, pt)) { avma = av; return 0; }
998 175 : if (pt)
999 175 : *pt = gerepileupto(av, FqX_to_mod(*pt, T, p));
1000 175 : return 1;
1001 : }
1002 371 : v = RgX_valrem(x, &x);
1003 371 : if (v % k) return 0;
1004 364 : v /= k;
1005 364 : a = gel(x,2); b = NULL;
1006 364 : if (!ispower(a, K, &b)) { avma = av; return 0; }
1007 350 : if (d)
1008 : {
1009 343 : GEN p = characteristic(x);
1010 343 : a = leading_coeff(x);
1011 343 : if (!ispower(a, K, &b)) { avma = av; return 0; }
1012 343 : x = RgX_normalize(x);
1013 343 : if (signe(p) && cmpii(p,K) <= 0)
1014 0 : pari_err_IMPL("ispower(general t_POL) in small characteristic");
1015 343 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
1016 343 : if (!RgX_equal(powgi(y, K), x)) { avma = av; return 0; }
1017 : }
1018 : else
1019 7 : y = pol_1(varn(x));
1020 350 : if (pt)
1021 : {
1022 350 : if (!gequal1(a))
1023 : {
1024 14 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
1025 14 : y = gmul(b,y);
1026 : }
1027 350 : if (v) y = RgX_shift_shallow(y, v);
1028 350 : *pt = gerepilecopy(av, y);
1029 : }
1030 0 : else avma = av;
1031 350 : return 1;
1032 : }
1033 :
1034 : long
1035 1771 : Z_ispowerall(GEN x, ulong k, GEN *pt)
1036 : {
1037 1771 : long s = signe(x);
1038 : ulong mask;
1039 1771 : if (!s) { if (pt) *pt = gen_0; return 1; }
1040 1771 : if (s > 0) {
1041 1645 : if (k == 2) return Z_issquareall(x, pt);
1042 973 : if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
1043 196 : if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
1044 168 : if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
1045 161 : return is_kth_power(x, k, pt);
1046 : }
1047 126 : if (!odd(k)) return 0;
1048 112 : if (Z_ispowerall(absi_shallow(x), k, pt))
1049 : {
1050 112 : if (pt) *pt = negi(*pt);
1051 112 : return 1;
1052 : };
1053 0 : return 0;
1054 : }
1055 :
1056 : /* is x a K-th power mod p ? Assume p prime. */
1057 : int
1058 224 : Fp_ispower(GEN x, GEN K, GEN p)
1059 : {
1060 224 : pari_sp av = avma;
1061 : GEN p_1;
1062 : long r;
1063 224 : x = modii(x, p);
1064 224 : if (!signe(x) || equali1(x)) { avma = av; return 1; }
1065 : /* implies p > 2 */
1066 126 : p_1 = subiu(p,1);
1067 126 : K = gcdii(K, p_1);
1068 126 : if (absequaliu(K, 2)) { r = kronecker(x,p); avma = av; return (r > 0); }
1069 49 : x = Fp_pow(x, diviiexact(p_1,K), p);
1070 49 : avma = av; return equali1(x);
1071 : }
1072 :
1073 : /* x unit defined modulo 2^e, e > 0, p prime */
1074 : static int
1075 2373 : U2_issquare(GEN x, long e)
1076 : {
1077 2373 : long r = signe(x)>=0?mod8(x):8-mod8(x);
1078 2373 : if (e==1) return 1;
1079 2373 : if (e==2) return (r&3L) == 1;
1080 2009 : return r == 1;
1081 : }
1082 : /* x unit defined modulo p^e, e > 0, p prime */
1083 : static int
1084 4690 : Up_issquare(GEN x, GEN p, long e)
1085 4690 : { return (absequaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
1086 :
1087 : long
1088 2548 : Zn_issquare(GEN d, GEN fn)
1089 : {
1090 : long j, np;
1091 2548 : if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
1092 2548 : if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
1093 : /* integer factorization */
1094 2548 : np = nbrows(fn);
1095 5320 : for (j = 1; j <= np; ++j)
1096 : {
1097 4970 : GEN r, p = gcoeff(fn, j, 1);
1098 4970 : long e = itos(gcoeff(fn, j, 2));
1099 4970 : long v = Z_pvalrem(d,p,&r);
1100 4970 : if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
1101 : }
1102 350 : return 1;
1103 : }
1104 :
1105 : static long
1106 1113 : Qp_ispower(GEN x, GEN K, GEN *pt)
1107 : {
1108 1113 : pari_sp av = avma;
1109 1113 : GEN z = Qp_sqrtn(x, K, NULL);
1110 1113 : if (!z) { avma = av; return 0; }
1111 819 : if (pt) *pt = z;
1112 819 : return 1;
1113 : }
1114 :
1115 : long
1116 7003164 : ispower(GEN x, GEN K, GEN *pt)
1117 : {
1118 : GEN z;
1119 :
1120 7003164 : if (!K) return gisanypower(x, pt);
1121 2982 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
1122 2982 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
1123 2982 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
1124 2933 : switch(typ(x)) {
1125 : case t_INT:
1126 826 : return Z_ispowerall(x, itou(K), pt);
1127 : case t_FRAC:
1128 : {
1129 21 : GEN a = gel(x,1), b = gel(x,2);
1130 21 : ulong k = itou(K);
1131 21 : if (pt) {
1132 14 : z = cgetg(3, t_FRAC);
1133 14 : if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
1134 14 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1135 : }
1136 0 : avma = (pari_sp)(z + 3); return 0;
1137 : }
1138 7 : return Z_ispower(a, k) && Z_ispower(b, k);
1139 : }
1140 : case t_INTMOD:
1141 189 : return Zn_ispower(gel(x,2), gel(x,1), K, pt);
1142 : case t_FFELT:
1143 28 : return FF_ispower(x, K, pt);
1144 :
1145 : case t_PADIC:
1146 1113 : return Qp_ispower(x, K, pt);
1147 : case t_POLMOD:
1148 14 : return polmodispower(x, K, pt);
1149 : case t_POL:
1150 714 : return polispower(x, K, pt);
1151 : case t_RFRAC: {
1152 7 : GEN a = gel(x,1), b = gel(x,2);
1153 7 : if (pt) {
1154 7 : z = cgetg(3, t_RFRAC);
1155 7 : if (ispower(a, K, &a) && polispower(b, K, &b)) {
1156 7 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1157 : }
1158 0 : avma = (pari_sp)(z + 3); return 0;
1159 : }
1160 0 : return (ispower(a, K, NULL) && polispower(b, K, NULL));
1161 : }
1162 : case t_REAL:
1163 7 : if (signe(x) < 0 && !mpodd(K)) return 0;
1164 : case t_COMPLEX:
1165 14 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1166 14 : return 1;
1167 :
1168 : case t_SER:
1169 7 : if (signe(x) && (!dvdsi(valp(x), K) || !ispower(gel(x,2), K, NULL)))
1170 0 : return 0;
1171 7 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1172 7 : return 1;
1173 : }
1174 0 : pari_err_TYPE("ispower",x);
1175 : return 0; /* LCOV_EXCL_LINE */
1176 : }
1177 :
1178 : long
1179 7000182 : gisanypower(GEN x, GEN *pty)
1180 : {
1181 7000182 : long tx = typ(x);
1182 : ulong k, h;
1183 7000182 : if (tx == t_INT) return Z_isanypower(x, pty);
1184 14 : if (tx == t_FRAC)
1185 : {
1186 14 : pari_sp av = avma;
1187 14 : GEN fa, P, E, a = gel(x,1), b = gel(x,2);
1188 : long i, j, p, e;
1189 14 : int sw = (abscmpii(a, b) > 0);
1190 :
1191 14 : if (sw) swap(a, b);
1192 14 : k = Z_isanypower(a, pty? &a: NULL);
1193 14 : if (!k)
1194 : { /* a = -1,1 or not a pure power */
1195 7 : if (!is_pm1(a)) { avma = av; return 0; }
1196 7 : if (signe(a) < 0) b = negi(b);
1197 7 : k = Z_isanypower(b, pty? &b: NULL);
1198 7 : if (!k || !pty) { avma = av; return k; }
1199 7 : *pty = gerepileupto(av, ginv(b));
1200 7 : return k;
1201 : }
1202 7 : fa = factoru(k);
1203 7 : P = gel(fa,1);
1204 7 : E = gel(fa,2); h = k;
1205 14 : for (i = lg(P) - 1; i > 0; i--)
1206 : {
1207 7 : p = P[i];
1208 7 : e = E[i];
1209 21 : for (j = 0; j < e; j++)
1210 14 : if (!is_kth_power(b, p, &b)) break;
1211 7 : if (j < e) k /= upowuu(p, e - j);
1212 : }
1213 7 : if (k == 1) { avma = av; return 0; }
1214 7 : if (!pty) { avma = av; return k; }
1215 0 : if (k != h) a = powiu(a, h/k);
1216 0 : *pty = gerepilecopy(av, mkfrac(a, b));
1217 0 : return k;
1218 : }
1219 0 : pari_err_TYPE("gisanypower", x);
1220 : return 0; /* LCOV_EXCL_LINE */
1221 : }
1222 :
1223 : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
1224 : * No need to optimize for 2,3,5,7 powers (done before) */
1225 : static long
1226 505715 : split_exponent(ulong e, GEN *x)
1227 : {
1228 : GEN fa, P, E;
1229 505715 : long i, j, l, k = 1;
1230 505715 : if (e == 1) return 1;
1231 14 : fa = factoru(e);
1232 14 : P = gel(fa,1);
1233 14 : E = gel(fa,2); l = lg(P);
1234 28 : for (i = 1; i < l; i++)
1235 : {
1236 14 : ulong p = P[i];
1237 28 : for (j = 0; j < E[i]; j++)
1238 : {
1239 : GEN y;
1240 14 : if (!is_kth_power(*x, p, &y)) break;
1241 14 : k *= p; *x = y;
1242 : }
1243 : }
1244 14 : return k;
1245 : }
1246 :
1247 : static long
1248 864647 : Z_isanypower_nosmalldiv(GEN *px)
1249 : { /* any prime divisor of x is > 102 */
1250 864647 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
1251 864647 : const double LOG103 = 4.6347; /* lower bound for log(103) */
1252 : forprime_t T;
1253 864647 : ulong mask = 7, e2;
1254 : long k, ex;
1255 864647 : GEN y, x = *px;
1256 :
1257 864647 : k = 1;
1258 864647 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
1259 864647 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
1260 864647 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
1261 864647 : if (u_forprime_init(&T, 11, e2))
1262 : {
1263 16961 : GEN logx = NULL;
1264 16961 : const ulong Q = 30011; /* prime */
1265 : ulong p, xmodQ;
1266 16961 : double dlogx = 0;
1267 : /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
1268 : * for large p the modular checks are no longer competitively fast */
1269 33964 : while ( (ex = is_pth_power(x, &y, &T, 30)) )
1270 : {
1271 42 : k *= ex; x = y;
1272 42 : e2 = (ulong)((expi(x) + 1) / LOG2_103);
1273 42 : u_forprime_restrict(&T, e2);
1274 : }
1275 16961 : if (DEBUGLEVEL>4) err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x));
1276 16961 : xmodQ = umodiu(x, Q);
1277 : /* test Q | x, just in case */
1278 16961 : if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
1279 : /* x^(1/p) < 2^31 */
1280 16947 : p = T.p;
1281 16947 : if (p <= e2)
1282 : {
1283 16933 : logx = logr_abs( itor(x, DEFAULTPREC) );
1284 16933 : dlogx = rtodbl(logx);
1285 16933 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1286 : }
1287 153965 : while (p && p <= e2)
1288 : { /* is x a p-th power ? By computing y = round(x^(1/p)).
1289 : * Check whether y^p = x, first mod Q, then exactly. */
1290 120071 : pari_sp av = avma;
1291 : long e;
1292 120071 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
1293 120071 : ulong ymodQ = umodiu(y,Q);
1294 120071 : if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
1295 21 : || !equalii(powiu(y, p), x)) avma = av;
1296 : else
1297 : {
1298 21 : k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
1299 21 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1300 21 : u_forprime_restrict(&T, e2);
1301 21 : continue; /* if success, retry same p */
1302 : }
1303 120050 : p = u_forprime_next(&T);
1304 : }
1305 : }
1306 864633 : *px = x; return k;
1307 : }
1308 :
1309 : static ulong tinyprimes[] = {
1310 : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
1311 : 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
1312 : 157, 163, 167, 173, 179, 181, 191, 193, 197, 199
1313 : };
1314 :
1315 : /* disregard the sign of x, caller will take care of x < 0 */
1316 : static long
1317 7000841 : Z_isanypower_aux(GEN x, GEN *pty)
1318 : {
1319 : long ex, v, i, l, k;
1320 : GEN y, P, E;
1321 7000841 : ulong mask, e = 0;
1322 :
1323 7000841 : if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
1324 :
1325 7000827 : if (signe(x) < 0) x = negi(x);
1326 7000827 : k = l = 1;
1327 7000827 : P = cgetg(26 + 1, t_VECSMALL);
1328 7000827 : E = cgetg(26 + 1, t_VECSMALL);
1329 : /* trial division */
1330 122942258 : for(i = 0; i < 26; i++)
1331 : {
1332 60136627 : ulong p = tinyprimes[i];
1333 : int stop;
1334 60136627 : v = Z_lvalrem_stop(&x, p, &stop);
1335 60136627 : if (v)
1336 : {
1337 7922348 : P[l] = p;
1338 7922348 : E[l] = v; l++;
1339 13588673 : e = ugcd(e, v); if (e == 1) goto END;
1340 : }
1341 54718340 : if (stop) {
1342 248038 : if (is_pm1(x)) k = e;
1343 248038 : goto END;
1344 : }
1345 : }
1346 :
1347 1334502 : if (e)
1348 : { /* Bingo. Result divides e */
1349 : long v3, v5, v7;
1350 505701 : ulong e2 = e;
1351 505701 : v = u_lvalrem(e2, 2, &e2);
1352 505701 : if (v)
1353 : {
1354 375249 : for (i = 0; i < v; i++)
1355 : {
1356 374171 : if (!Z_issquareall(x, &y)) break;
1357 1288 : k <<= 1; x = y;
1358 : }
1359 : }
1360 505701 : mask = 0;
1361 505701 : v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
1362 505701 : v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
1363 505701 : v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
1364 1011479 : while ( (ex = is_357_power(x, &y, &mask)) ) {
1365 77 : x = y;
1366 77 : switch(ex)
1367 : {
1368 28 : case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
1369 28 : case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
1370 21 : case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
1371 : }
1372 : }
1373 505701 : k *= split_exponent(e2, &x);
1374 : }
1375 : else
1376 828801 : k = Z_isanypower_nosmalldiv(&x);
1377 : END:
1378 7000827 : if (pty && k != 1)
1379 : {
1380 8029 : if (e)
1381 : { /* add missing small factors */
1382 6867 : y = powuu(P[1], E[1] / k);
1383 6867 : for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
1384 6867 : x = equali1(x)? y: mulii(x,y);
1385 : }
1386 8029 : *pty = x;
1387 : }
1388 7000827 : return k == 1? 0: k;
1389 : }
1390 :
1391 : long
1392 7000841 : Z_isanypower(GEN x, GEN *pty)
1393 : {
1394 7000841 : pari_sp av = avma;
1395 7000841 : long k = Z_isanypower_aux(x, pty);
1396 7000841 : if (!k) { avma = av; return 0; }
1397 8092 : if (signe(x) < 0)
1398 : {
1399 42 : long v = vals(k);
1400 42 : if (v)
1401 : {
1402 : GEN y;
1403 28 : k >>= v;
1404 28 : if (k == 1) { avma = av; return 0; }
1405 21 : if (!pty) { avma = av; return k; }
1406 14 : y = *pty;
1407 14 : y = powiu(y, 1<<v);
1408 14 : togglesign(y);
1409 14 : *pty = gerepileuptoint(av, y);
1410 14 : return k;
1411 : }
1412 14 : if (pty) togglesign_safe(pty);
1413 : }
1414 8064 : if (pty) *pty = gerepilecopy(av, *pty); else avma = av;
1415 8064 : return k;
1416 : }
1417 :
1418 : /* Faster than */
1419 : /* !cmpii(n, int2n(vali(n))) */
1420 : /* !cmpis(shifti(n, -vali(n)), 1) */
1421 : /* expi(n) == vali(n) */
1422 : /* hamming(n) == 1 */
1423 : /* even for single-word values, and much faster for multiword values. */
1424 : /* If all you have is a word, you can just use n & !(n & (n-1)). */
1425 : long
1426 80682 : Z_ispow2(GEN n)
1427 : {
1428 : GEN xp;
1429 : long i, lx;
1430 : ulong u;
1431 80682 : if (signe(n) != 1) return 0;
1432 80675 : xp = int_LSW(n);
1433 80675 : lx = lgefint(n);
1434 80675 : u = *xp;
1435 80959 : for (i = 3; i < lx; ++i)
1436 : {
1437 78140 : if (u) return 0;
1438 284 : xp = int_nextW(xp);
1439 284 : u = *xp;
1440 : }
1441 2819 : return !(u & (u-1));
1442 : }
1443 :
1444 : static long
1445 841788 : isprimepower_i(GEN n, GEN *pt, long flag)
1446 : {
1447 841788 : pari_sp av = avma;
1448 : long i, v;
1449 :
1450 841788 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1451 841788 : if (signe(n) <= 0) return 0;
1452 :
1453 841788 : if (lgefint(n) == 3)
1454 : {
1455 : ulong p;
1456 541121 : v = uisprimepower(n[2], &p);
1457 541121 : if (v)
1458 : {
1459 54909 : if (pt) *pt = utoipos(p);
1460 54909 : return v;
1461 : }
1462 486212 : return 0;
1463 : }
1464 1662340 : for (i = 0; i < 26; i++)
1465 : {
1466 1626494 : ulong p = tinyprimes[i];
1467 1626494 : v = Z_lvalrem(n, p, &n);
1468 1626491 : if (v)
1469 : {
1470 264818 : avma = av;
1471 264818 : if (!is_pm1(n)) return 0;
1472 424 : if (pt) *pt = utoipos(p);
1473 422 : return v;
1474 : }
1475 : }
1476 : /* p | n => p >= 103 */
1477 35846 : v = Z_isanypower_nosmalldiv(&n); /* expensive */
1478 35846 : if (!(flag? isprime(n): BPSW_psp(n))) { avma = av; return 0; }
1479 5534 : if (pt) *pt = gerepilecopy(av, n); else avma = av;
1480 5534 : return v;
1481 : }
1482 : long
1483 840098 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
1484 : long
1485 1690 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1486 :
1487 : long
1488 542220 : uisprimepower(ulong n, ulong *pp)
1489 : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1490 : * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1491 542220 : const ulong CUTOFF = 200UL;
1492 542220 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1493 542220 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1494 : #ifdef LONG_IS_64BIT
1495 : /* primes preceeding the appropriate root of ULONG_MAX. */
1496 481914 : const ulong ROOT9 = 137;
1497 481914 : const ulong ROOT8 = 251;
1498 481914 : const ulong ROOT7 = 563;
1499 481914 : const ulong ROOT5 = 7129;
1500 481914 : const ulong ROOT4 = 65521;
1501 : #else
1502 60306 : const ulong ROOT9 = 11;
1503 60306 : const ulong ROOT8 = 13;
1504 60306 : const ulong ROOT7 = 23;
1505 60306 : const ulong ROOT5 = 83;
1506 60306 : const ulong ROOT4 = 251;
1507 : #endif
1508 : ulong mask;
1509 : long v, i;
1510 : int e;
1511 542220 : if (n < 2) return 0;
1512 542206 : if (!odd(n)) {
1513 271187 : if (n & (n-1)) return 0;
1514 1383 : *pp = 2; return vals(n);
1515 : }
1516 271019 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1517 3653916 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1518 : {
1519 3594841 : ulong p = tinyprimes[i];
1520 3594841 : if (n % p == 0)
1521 : {
1522 211538 : v = u_lvalrem(n, p, &n);
1523 211538 : if (n == 1) { *pp = p; return v; }
1524 209438 : return 0;
1525 : }
1526 : }
1527 : /* p | n => p >= CUTOFF */
1528 :
1529 59075 : if (n < CUTOFF3)
1530 : {
1531 46354 : if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1532 0 : if (uissquareall(n, &n)) { *pp = n; return 2; }
1533 0 : return 0;
1534 : }
1535 :
1536 : /* Check for squares, fourth powers, and eighth powers as appropriate. */
1537 12721 : v = 1;
1538 12721 : if (uissquareall(n, &n)) {
1539 0 : v <<= 1;
1540 0 : if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1541 0 : v <<= 1;
1542 0 : if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1543 : }
1544 : }
1545 :
1546 12721 : if (CUTOFF > ROOT5) mask = 1;
1547 : else
1548 : {
1549 12720 : const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1550 12720 : if (n < CUTOFF5) mask = 1; else mask = 3;
1551 12720 : if (CUTOFF <= ROOT7)
1552 : {
1553 12720 : const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1554 12720 : if (n >= CUTOFF7) mask = 7;
1555 : }
1556 : }
1557 :
1558 12721 : if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1559 12721 : if ((e = uis_357_power(n, &n, &mask))) v *= e;
1560 :
1561 12721 : if (uisprime_101(n)) { *pp = n; return v; }
1562 6984 : return 0;
1563 : }
1564 :
1565 : /*********************************************************************/
1566 : /** **/
1567 : /** KRONECKER SYMBOL **/
1568 : /** **/
1569 : /*********************************************************************/
1570 : /* t = 3,5 mod 8 ? (= 2 not a square mod t) */
1571 : static int
1572 634128020 : ome(long t)
1573 : {
1574 634128020 : switch(t & 7)
1575 : {
1576 : case 3:
1577 360640202 : case 5: return 1;
1578 273487818 : default: return 0;
1579 : }
1580 : }
1581 : /* t a t_INT, is t = 3,5 mod 8 ? */
1582 : static int
1583 5371111 : gome(GEN t)
1584 5371111 : { return signe(t)? ome( mod2BIL(t) ): 0; }
1585 :
1586 : /* assume y odd, return kronecker(x,y) * s */
1587 : static long
1588 480718755 : krouu_s(ulong x, ulong y, long s)
1589 : {
1590 480718755 : ulong x1 = x, y1 = y, z;
1591 2606367622 : while (x1)
1592 : {
1593 1644998896 : long r = vals(x1);
1594 1645063303 : if (r)
1595 : {
1596 885758877 : if (odd(r) && ome(y1)) s = -s;
1597 885625686 : x1 >>= r;
1598 : }
1599 1644930112 : if (x1 & y1 & 2) s = -s;
1600 1644930112 : z = y1 % x1; y1 = x1; x1 = z;
1601 : }
1602 480649971 : return (y1 == 1)? s: 0;
1603 : }
1604 :
1605 : long
1606 6103807 : kronecker(GEN x, GEN y)
1607 : {
1608 6103807 : pari_sp av = avma;
1609 6103807 : long s = 1, r;
1610 : ulong xu, yu;
1611 :
1612 6103807 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
1613 6103807 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
1614 6103807 : switch (signe(y))
1615 : {
1616 0 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
1617 0 : case 0: return is_pm1(x);
1618 : }
1619 6103807 : r = vali(y);
1620 6103807 : if (r)
1621 : {
1622 12216 : if (!mpodd(x)) { avma = av; return 0; }
1623 10494 : if (odd(r) && gome(x)) s = -s;
1624 10494 : y = shifti(y,-r);
1625 : }
1626 6102085 : x = modii(x,y);
1627 13108465 : while (lgefint(x) > 3) /* x < y */
1628 : {
1629 : GEN z;
1630 904295 : r = vali(x);
1631 904295 : if (r)
1632 : {
1633 493412 : if (odd(r) && gome(y)) s = -s;
1634 493412 : x = shifti(x,-r);
1635 : }
1636 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1637 904295 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
1638 904295 : z = remii(y,x); y = x; x = z;
1639 904295 : if (gc_needed(av,2))
1640 : {
1641 0 : if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
1642 0 : gerepileall(av, 2, &x, &y);
1643 : }
1644 : }
1645 6102085 : xu = itou(x);
1646 6102085 : if (!xu) return is_pm1(y)? s: 0;
1647 6080904 : r = vals(xu);
1648 6080904 : if (r)
1649 : {
1650 4131440 : if (odd(r) && gome(y)) s = -s;
1651 4131440 : xu >>= r;
1652 : }
1653 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1654 6080904 : if (xu & mod2BIL(y) & 2) s = -s;
1655 6080904 : yu = umodiu(y, xu);
1656 6080904 : avma = av; return krouu_s(yu, xu, s);
1657 : }
1658 :
1659 : long
1660 30842 : krois(GEN x, long y)
1661 : {
1662 : ulong yu;
1663 30842 : long s = 1;
1664 :
1665 30842 : if (y <= 0)
1666 : {
1667 7 : if (y == 0) return is_pm1(x);
1668 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
1669 : }
1670 : else
1671 30835 : yu = (ulong)y;
1672 30835 : if (!odd(yu))
1673 : {
1674 : long r;
1675 14056 : if (!mpodd(x)) return 0;
1676 10472 : r = vals(yu); yu >>= r;
1677 10472 : if (odd(r) && gome(x)) s = -s;
1678 : }
1679 27251 : return krouu_s(umodiu(x, yu), yu, s);
1680 : }
1681 : /* assume y != 0 */
1682 : long
1683 342596563 : kroiu(GEN x, ulong y)
1684 : {
1685 : long r;
1686 342596563 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
1687 2132032 : if (!mpodd(x)) return 0;
1688 2109898 : r = vals(y); y >>= r;
1689 2109898 : return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
1690 : }
1691 :
1692 : /* assume y > 0, odd, return s * kronecker(x,y) */
1693 : static long
1694 120510 : krouodd(ulong x, GEN y, long s)
1695 : {
1696 : long r;
1697 120510 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
1698 41844 : if (!x) return 0; /* y != 1 */
1699 41844 : r = vals(x);
1700 41844 : if (r)
1701 : {
1702 7341 : if (odd(r) && gome(y)) s = -s;
1703 7341 : x >>= r;
1704 : }
1705 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1706 41844 : if (x & mod2BIL(y) & 2) s = -s;
1707 41844 : return krouu_s(umodiu(y,x), x, s);
1708 : }
1709 :
1710 : long
1711 119997 : krosi(long x, GEN y)
1712 : {
1713 119997 : const pari_sp av = avma;
1714 119997 : long s = 1, r;
1715 119997 : switch (signe(y))
1716 : {
1717 0 : case -1: y = negi(y); if (x < 0) s = -1; break;
1718 0 : case 0: return (x==1 || x==-1);
1719 : }
1720 119997 : r = vali(y);
1721 119997 : if (r)
1722 : {
1723 8421 : if (!odd(x)) { avma = av; return 0; }
1724 8421 : if (odd(r) && ome(x)) s = -s;
1725 8421 : y = shifti(y,-r);
1726 : }
1727 119997 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
1728 119997 : s = krouodd((ulong)x, y, s);
1729 119997 : avma = av; return s;
1730 : }
1731 :
1732 : long
1733 513 : kroui(ulong x, GEN y)
1734 : {
1735 513 : const pari_sp av = avma;
1736 513 : long s = 1, r;
1737 513 : switch (signe(y))
1738 : {
1739 0 : case -1: y = negi(y); break;
1740 0 : case 0: return x==1UL;
1741 : }
1742 513 : r = vali(y);
1743 513 : if (r)
1744 : {
1745 0 : if (!odd(x)) { avma = av; return 0; }
1746 0 : if (odd(r) && ome(x)) s = -s;
1747 0 : y = shifti(y,-r);
1748 : }
1749 513 : s = krouodd(x, y, s);
1750 513 : avma = av; return s;
1751 : }
1752 :
1753 : long
1754 73601538 : kross(long x, long y)
1755 : {
1756 : ulong yu;
1757 73601538 : long s = 1;
1758 :
1759 73601538 : if (y <= 0)
1760 : {
1761 441 : if (y == 0) return (labs(x)==1);
1762 413 : yu = (ulong)-y; if (x < 0) s = -1;
1763 : }
1764 : else
1765 73601097 : yu = (ulong)y;
1766 73601510 : if (!odd(yu))
1767 : {
1768 : long r;
1769 17156651 : if (!odd(x)) return 0;
1770 12332671 : r = vals(yu); yu >>= r;
1771 12332671 : if (odd(r) && ome(x)) s = -s;
1772 : }
1773 68777530 : x %= (long)yu; if (x < 0) x += yu;
1774 68777530 : return krouu_s((ulong)x, yu, s);
1775 : }
1776 :
1777 : long
1778 63022301 : krouu(ulong x, ulong y)
1779 : {
1780 : long r;
1781 63022301 : if (odd(y)) return krouu_s(x, y, 1);
1782 1674 : if (!odd(x)) return 0;
1783 1674 : r = vals(y); y >>= r;
1784 1674 : return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
1785 : }
1786 :
1787 : /*********************************************************************/
1788 : /** **/
1789 : /** HILBERT SYMBOL **/
1790 : /** **/
1791 : /*********************************************************************/
1792 : /* x,y are t_INT or t_REAL */
1793 : static long
1794 9982 : mphilbertoo(GEN x, GEN y)
1795 : {
1796 9982 : long sx = signe(x), sy = signe(y);
1797 9982 : if (!sx || !sy) return 0;
1798 9982 : return (sx < 0 && sy < 0)? -1: 1;
1799 : }
1800 :
1801 : long
1802 53144 : hilbertii(GEN x, GEN y, GEN p)
1803 : {
1804 : pari_sp av;
1805 : long oddvx, oddvy, z;
1806 :
1807 53144 : if (!p) return mphilbertoo(x,y);
1808 43183 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
1809 43183 : if (!signe(x) || !signe(y)) return 0;
1810 43162 : av = avma;
1811 43162 : oddvx = odd(Z_pvalrem(x,p,&x));
1812 43162 : oddvy = odd(Z_pvalrem(y,p,&y));
1813 : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
1814 43162 : if (absequaliu(p, 2))
1815 : {
1816 10689 : z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
1817 10689 : if (oddvx && gome(y)) z = -z;
1818 10689 : if (oddvy && gome(x)) z = -z;
1819 : }
1820 : else
1821 : {
1822 32473 : z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
1823 32473 : if (oddvx && kronecker(y,p) < 0) z = -z;
1824 32473 : if (oddvy && kronecker(x,p) < 0) z = -z;
1825 : }
1826 43162 : avma = av; return z;
1827 : }
1828 :
1829 : static void
1830 196 : err_prec(void) { pari_err_PREC("hilbert"); }
1831 : static void
1832 161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
1833 : static void
1834 56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
1835 :
1836 : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
1837 : * Return lift(x) provided it's p-adic accuracy is large enough to decide
1838 : * hilbert()'s value [ problem at p = 2 ] */
1839 : static GEN
1840 420 : lift_intmod(GEN x, GEN *pp)
1841 : {
1842 420 : GEN p = *pp, N = gel(x,1);
1843 420 : x = gel(x,2);
1844 420 : if (!p)
1845 : {
1846 266 : *pp = p = N;
1847 266 : switch(itos_or_0(p))
1848 : {
1849 : case 2:
1850 126 : case 4: err_prec();
1851 : }
1852 140 : return x;
1853 : }
1854 154 : if (!signe(p)) err_oo(N);
1855 112 : if (absequaliu(p,2))
1856 42 : { if (vali(N) <= 2) err_prec(); }
1857 : else
1858 70 : { if (!dvdii(N,p)) err_p(N,p); }
1859 28 : if (!signe(x)) err_prec();
1860 21 : return x;
1861 : }
1862 : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
1863 : * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
1864 : * to decide hilbert()'s value [ problem at p = 2 ]*/
1865 : static GEN
1866 210 : lift_padic(GEN x, GEN *pp)
1867 : {
1868 210 : GEN p = *pp, q = gel(x,2), y = gel(x,4);
1869 210 : if (!p) *pp = p = q;
1870 147 : else if (!equalii(p,q)) err_p(p, q);
1871 105 : if (absequaliu(p,2) && precp(x) <= 2) err_prec();
1872 70 : if (!signe(y)) err_prec();
1873 70 : return odd(valp(x))? mulii(p,y): y;
1874 : }
1875 :
1876 : long
1877 658 : hilbert(GEN x, GEN y, GEN p)
1878 : {
1879 658 : pari_sp av = avma;
1880 658 : long tx = typ(x), ty = typ(y), z;
1881 :
1882 658 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
1883 658 : if (tx == t_REAL)
1884 : {
1885 77 : if (p && signe(p)) err_oo(p);
1886 63 : switch (ty)
1887 : {
1888 : case t_INT:
1889 7 : case t_REAL: return mphilbertoo(x,y);
1890 0 : case t_FRAC: return mphilbertoo(x,gel(y,1));
1891 56 : default: pari_err_TYPE2("hilbert",x,y);
1892 : }
1893 : }
1894 581 : if (ty == t_REAL)
1895 : {
1896 14 : if (p && signe(p)) err_oo(p);
1897 14 : switch (tx)
1898 : {
1899 : case t_INT:
1900 14 : case t_REAL: return mphilbertoo(x,y);
1901 0 : case t_FRAC: return mphilbertoo(gel(x,1),y);
1902 0 : default: pari_err_TYPE2("hilbert",x,y);
1903 : }
1904 : }
1905 567 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
1906 364 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
1907 :
1908 308 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
1909 245 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
1910 :
1911 168 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
1912 168 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
1913 :
1914 168 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
1915 168 : if (p && !signe(p)) p = NULL;
1916 168 : z = hilbertii(x,y,p); avma = av; return z;
1917 : }
1918 :
1919 : /*******************************************************************/
1920 : /* */
1921 : /* SQUARE ROOT MODULO p */
1922 : /* */
1923 : /*******************************************************************/
1924 : static void
1925 4267163 : checkp(ulong q, ulong p)
1926 4267163 : { if (!q) pari_err_PRIME("Fl_nonsquare",utoipos(p)); }
1927 : /* p = 1 (mod 4) prime, return the first quadratic non-residue, a prime */
1928 : static ulong
1929 26768373 : nonsquare1_Fl(ulong p)
1930 : {
1931 : forprime_t S;
1932 : ulong q;
1933 26768373 : if ((p & 7UL) != 1) return 2UL;
1934 10951424 : q = p % 3; if (q == 2) return 3UL;
1935 3412517 : checkp(q, p);
1936 3412510 : q = p % 5; if (q == 2 || q == 3) return 5UL;
1937 518825 : checkp(q, p);
1938 518825 : q = p % 7; if (q != 4 && q >= 3) return 7UL;
1939 209946 : checkp(q, p);
1940 209946 : u_forprime_init(&S, 11, p);
1941 545767 : while ((q = u_forprime_next(&S)))
1942 : {
1943 335821 : long i = krouu(q, p);
1944 335821 : if (i < 0) return q;
1945 125875 : checkp(q, p);
1946 : }
1947 0 : checkp(0, p);
1948 : return 0; /*LCOV_EXCL_LINE*/
1949 : }
1950 : /* p > 2 a prime */
1951 : ulong
1952 7714 : nonsquare_Fl(ulong p)
1953 7714 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
1954 :
1955 : ulong
1956 150646 : Fl_2gener_pre(ulong p, ulong pi)
1957 : {
1958 150646 : ulong p1 = p-1;
1959 150646 : long e = vals(p1);
1960 150644 : if (e == 1) return p1;
1961 56689 : return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
1962 : }
1963 :
1964 : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
1965 : ulong
1966 64210127 : Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
1967 : {
1968 : long i, e, k;
1969 : ulong p1, q, v, w;
1970 :
1971 64210127 : if (!a) return 0;
1972 62880546 : p1 = p - 1; e = vals(p1);
1973 62876417 : if (e == 0) /* p = 2 */
1974 : {
1975 418894 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1976 418887 : return ((a & 1) == 0)? 0: 1;
1977 : }
1978 62457523 : if (e == 1)
1979 : {
1980 35744233 : v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
1981 35745609 : if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
1982 35713146 : p1 = p - v; if (v > p1) v = p1;
1983 35713146 : return v;
1984 : }
1985 26713290 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
1986 26713290 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
1987 26713290 : if (!p1) return 0;
1988 26713290 : v = Fl_mul_pre(a, p1, p, pi);
1989 26713290 : w = Fl_mul_pre(v, p1, p, pi);
1990 26713290 : if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
1991 75728625 : while (w != 1)
1992 : { /* a*w = v^2, y primitive 2^e-th root of 1
1993 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1994 22330313 : p1 = Fl_sqr_pre(w,p,pi);
1995 22330313 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
1996 22330313 : if (k == e) return ~0UL;
1997 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1998 22302059 : p1 = y;
1999 22302059 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
2000 22302059 : y = Fl_sqr_pre(p1, p, pi); e = k;
2001 22302059 : w = Fl_mul_pre(y, w, p, pi);
2002 22302059 : v = Fl_mul_pre(v, p1, p, pi);
2003 : }
2004 26685029 : p1 = p - v; if (v > p1) v = p1;
2005 26685029 : return v;
2006 : }
2007 :
2008 : ulong
2009 60759156 : Fl_sqrt(ulong a, ulong p)
2010 : {
2011 60759156 : ulong pi = get_Fl_red(p);
2012 60760814 : return Fl_sqrt_pre_i(a, 0, p, pi);
2013 : }
2014 :
2015 : ulong
2016 3416097 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
2017 : {
2018 3416097 : return Fl_sqrt_pre_i(a, 0, p, pi);
2019 : }
2020 :
2021 : static ulong
2022 46242 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
2023 : {
2024 : ulong x, y, m;
2025 46242 : ulong le1 = upowuu(l, e-1);
2026 72873 : for (x = 2; ; x++)
2027 : {
2028 99504 : y = Fl_powu_pre(x, r, p, pi);
2029 72873 : if (y==1) continue;
2030 56598 : m = Fl_powu_pre(y, le1, p, pi);
2031 56598 : if (m != 1) break;
2032 : }
2033 46242 : *pt_m = m;
2034 46242 : return y;
2035 : }
2036 :
2037 : /* solve x^l = a , l prime in G of order q.
2038 : *
2039 : * q = (l^e)*r, e >= 1, (r,l) = 1
2040 : * y generates the l-Sylow of G
2041 : * m = y^(l^(e-1)) != 1 */
2042 : static ulong
2043 110495 : Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
2044 : {
2045 : ulong p1, v, w, z, dl;
2046 : ulong u2;
2047 110495 : if (a==0) return a;
2048 110495 : u2 = Fl_inv(l%r, r);
2049 110496 : v = Fl_powu_pre(a, u2, p, pi);
2050 110495 : w = Fl_powu_pre(v, l, p, pi);
2051 110496 : w = Fl_mul_pre(w, Fl_inv(a, p), p, pi);
2052 110482 : if (w==1) return v;
2053 45309 : if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
2054 109623 : while (w!=1)
2055 : {
2056 49554 : ulong k = 0;
2057 49554 : p1 = w;
2058 : do
2059 : {
2060 73427 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
2061 73427 : k++;
2062 73427 : } while (p1!=1);
2063 49554 : if (k==e) return ULONG_MAX;
2064 19005 : dl = Fl_log_pre(z, m, l, p, pi);
2065 19005 : dl = Fl_neg(dl, l);
2066 19005 : p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
2067 19005 : m = Fl_powu_pre(m, dl, p, pi);
2068 19005 : e = k;
2069 19005 : v = Fl_mul_pre(p1,v,p,pi);
2070 19005 : y = Fl_powu_pre(p1,l,p,pi);
2071 19005 : w = Fl_mul_pre(y,w,p,pi);
2072 : }
2073 14760 : return v;
2074 : }
2075 :
2076 : static ulong
2077 109789 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
2078 : {
2079 109789 : ulong r, e = u_lvalrem(p-1, l, &r);
2080 109788 : return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
2081 : }
2082 :
2083 : ulong
2084 109788 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
2085 : {
2086 109788 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2087 : }
2088 :
2089 : ulong
2090 0 : Fl_sqrtl(ulong a, ulong l, ulong p)
2091 : {
2092 0 : ulong pi = get_Fl_red(p);
2093 0 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2094 : }
2095 :
2096 : ulong
2097 64969 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
2098 : {
2099 64969 : ulong m, q = p-1, z;
2100 64969 : ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
2101 64969 : if (a==0)
2102 : {
2103 48118 : if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
2104 48111 : if (zetan) *zetan = 1UL;
2105 48111 : return 0;
2106 : }
2107 16851 : if (n==1)
2108 : {
2109 0 : if (zetan) *zetan = 1;
2110 0 : return n < 0? Fl_inv(a,p): a;
2111 : }
2112 16851 : if (n==2)
2113 : {
2114 3444 : if (zetan) *zetan = p-1;
2115 3444 : return Fl_sqrt_pre_i(a, 0, p, pi);
2116 : }
2117 13407 : if (a == 1 && !zetan) return a;
2118 7478 : m = ugcd(nn, q);
2119 7478 : z = 1;
2120 7478 : if (m!=1)
2121 : {
2122 912 : GEN F = factoru(m);
2123 : long i, j, e;
2124 : ulong r, zeta, y, l;
2125 1915 : for (i = nbrows(F); i; i--)
2126 : {
2127 1017 : l = ucoeff(F,i,1);
2128 1017 : j = ucoeff(F,i,2);
2129 1017 : e = u_lvalrem(q,l, &r);
2130 1017 : y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
2131 1017 : if (zetan)
2132 464 : z = Fl_mul_pre(z, Fl_powu_pre(y, upowuu(l,e-j), p, pi), p, pi);
2133 1017 : if (a!=1)
2134 : do
2135 : {
2136 707 : a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
2137 693 : if (!a) return ULONG_MAX;
2138 693 : } while (--j);
2139 : }
2140 : }
2141 7464 : if (m != nn)
2142 : {
2143 6594 : ulong qm = q/m, nm = nn/m;
2144 6594 : a = Fl_powu_pre(a, Fl_inv(nm%qm, qm), p, pi);
2145 : }
2146 7464 : if (n < 0) a = Fl_inv(a, p);
2147 7464 : if (zetan) *zetan = z;
2148 7464 : return a;
2149 : }
2150 :
2151 : ulong
2152 64969 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
2153 : {
2154 64969 : ulong pi = get_Fl_red(p);
2155 64969 : return Fl_sqrtn_pre(a, n, p, pi, zetan);
2156 : }
2157 :
2158 : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
2159 : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
2160 : * and in average is about 2 or 2.5 times worse. But try both algorithms for
2161 : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
2162 : *
2163 : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
2164 : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
2165 : * If (a|p)=1, then sqrt(a) is in F_p.
2166 : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
2167 :
2168 : /* compute y^2, y = y[1] + y[2] X */
2169 : static GEN
2170 449 : sqrt_Cipolla_sqr(void *data, GEN y)
2171 : {
2172 449 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
2173 449 : GEN u2 = sqri(u), v2 = sqri(v);
2174 449 : v = subii(sqri(addii(v,u)), addii(u2,v2));
2175 449 : u = addii(u2, mulii(v2,n));
2176 : /* NOT mkvec2: must be gerepileupto-able */
2177 449 : retmkvec2(modii(u,p), modii(v,p));
2178 : }
2179 : /* compute (t+X) y^2 */
2180 : static GEN
2181 23 : sqrt_Cipolla_msqr(void *data, GEN y)
2182 : {
2183 23 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
2184 23 : ulong t = gt[2];
2185 23 : GEN d = addii(u, mului(t,v)), d2= sqri(d);
2186 23 : GEN b = remii(mulii(a,v), p);
2187 23 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
2188 23 : v = subii(d2, mulii(b,v));
2189 : /* NOT mkvec2: must be gerepileupto-able */
2190 23 : retmkvec2(modii(u,p), modii(v,p));
2191 : }
2192 : /* assume a reduced mod p [ otherwise correct but inefficient ] */
2193 : static GEN
2194 8 : sqrt_Cipolla(GEN a, GEN p)
2195 : {
2196 : pari_sp av1;
2197 : GEN u, v, n, y, pov2;
2198 : ulong t;
2199 :
2200 8 : if (kronecker(a, p) < 0) return NULL;
2201 8 : pov2 = shifti(p,-1);
2202 8 : if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
2203 :
2204 8 : av1 = avma;
2205 41 : for(t=1; ; t++)
2206 : {
2207 74 : n = subsi((long)(t*t), a);
2208 41 : if (kronecker(n, p) < 0) break;
2209 33 : avma = av1;
2210 : }
2211 :
2212 : /* compute (t+X)^((p-1)/2) =: u+vX */
2213 8 : u = utoipos(t);
2214 8 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
2215 : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
2216 : /* Now u+vX = (t+X)^((p-1)/2); thus
2217 : * (u+vX)(t+X) = sqrt(a) + 0 X
2218 : * Whence,
2219 : * sqrt(a) = (u+vt)t - v*a
2220 : * 0 = (u+vt)
2221 : * Thus a square root is v*a */
2222 :
2223 8 : v = Fp_mul(gel(y, 2), a, p);
2224 8 : if (cmpii(v,pov2) > 0) v = subii(p,v);
2225 8 : return v;
2226 : }
2227 :
2228 : /* Return NULL if p is found to be composite */
2229 : static GEN
2230 2800 : Fp_2gener_all(long e, GEN p)
2231 : {
2232 : GEN y, m;
2233 : long k;
2234 2800 : GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
2235 2800 : if (e==0 && !equaliu(p,2)) return NULL;
2236 9265 : for (k=2; ; k++)
2237 6465 : {
2238 9265 : long i = krosi(k, p);
2239 9265 : if (i >= 0)
2240 : {
2241 6465 : if (i) continue;
2242 0 : return NULL;
2243 : }
2244 2800 : y = m = Fp_pow(utoi(k), q, p);
2245 9402 : for (i=1; i<e; i++)
2246 6602 : if (equali1(m = Fp_sqr(m, p))) break;
2247 2800 : if (i == e) break; /* success */
2248 : }
2249 2800 : return y;
2250 : }
2251 :
2252 : /* Return NULL if p is found to be composite */
2253 : GEN
2254 980 : Fp_2gener(GEN p)
2255 980 : { return Fp_2gener_all(vali(subis(p,1)),p); }
2256 :
2257 : /* smallest square root */
2258 : static GEN
2259 27543 : choose_sqrt(GEN v, GEN p)
2260 : {
2261 27543 : pari_sp av = avma;
2262 27543 : GEN q = subii(p,v);
2263 27543 : if (cmpii(v,q) > 0) v = q; else avma = av;
2264 27543 : return v;
2265 : }
2266 : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
2267 : GEN
2268 2361838 : Fp_sqrt_i(GEN a, GEN y, GEN p)
2269 : {
2270 2361838 : pari_sp av = avma;
2271 : long i, k, e;
2272 : GEN p1, q, v, w;
2273 :
2274 2361838 : if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
2275 2361838 : if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
2276 2361838 : if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
2277 2361838 : if (lgefint(p) == 3)
2278 : {
2279 2334188 : ulong pp = uel(p,2), u = Fl_sqrt(umodiu(a, pp), pp);
2280 2334174 : if (u == ~0UL) return NULL;
2281 2334132 : return utoi(u);
2282 : }
2283 :
2284 27650 : a = modii(a, p); if (!signe(a)) { avma = av; return gen_0; }
2285 27559 : p1 = subiu(p,1); e = vali(p1);
2286 27559 : if (e <= 2)
2287 : { /* direct formulas more efficient */
2288 : pari_sp av2;
2289 22553 : if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
2290 22553 : if (e == 1)
2291 : {
2292 13346 : q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
2293 13346 : v = Fp_pow(a, q, p);
2294 : }
2295 : else
2296 : { /* Atkin's formula */
2297 9207 : GEN i, a2 = shifti(a,1);
2298 9207 : if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
2299 9207 : q = shifti(p1, -3); /* (p-5)/8 */
2300 9207 : v = Fp_pow(a2, q, p);
2301 9207 : i = Fp_mul(a2, Fp_sqr(v,p), p); /* i^2 = -1 */
2302 9207 : v = Fp_mul(a, Fp_mul(v, subiu(i,1), p), p);
2303 : }
2304 22553 : av2 = avma;
2305 : /* must check equality in case (a/p) = -1 or p not prime */
2306 22553 : e = equalii(Fp_sqr(v,p), a); avma = av2;
2307 22553 : return e? gerepileuptoint(av,choose_sqrt(v,p)): NULL;
2308 : }
2309 : /* On average, Cipolla is better than Tonelli/Shanks if and only if
2310 : * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
2311 5006 : if (e*(e-1) > 20 + 8 * expi(p))
2312 : {
2313 8 : v = sqrt_Cipolla(a,p); if (!v) { avma = av; return NULL; }
2314 8 : return gerepileuptoint(av,v);
2315 : }
2316 4998 : if (!y)
2317 : {
2318 1820 : y = Fp_2gener_all(e, p);
2319 1820 : if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
2320 : }
2321 4998 : q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
2322 4998 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
2323 4998 : v = Fp_mul(a, p1, p);
2324 4998 : w = Fp_mul(v, p1, p);
2325 17035 : while (!equali1(w))
2326 : { /* a*w = v^2, y primitive 2^e-th root of 1
2327 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
2328 7039 : p1 = Fp_sqr(w,p);
2329 7039 : for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
2330 7039 : if (k == e) { avma=av; return NULL; } /* p composite or (a/p) != 1 */
2331 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
2332 7039 : p1 = y;
2333 7039 : for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
2334 7039 : y = Fp_sqr(p1, p); e = k;
2335 7039 : w = Fp_mul(y, w, p);
2336 7039 : v = Fp_mul(v, p1, p);
2337 7039 : if (gc_needed(av,1))
2338 : {
2339 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
2340 0 : gerepileall(av,3, &y,&w,&v);
2341 : }
2342 : }
2343 4998 : return gerepileuptoint(av, choose_sqrt(v,p));
2344 : }
2345 :
2346 : GEN
2347 2347411 : Fp_sqrt(GEN a, GEN p)
2348 : {
2349 2347411 : return Fp_sqrt_i(a, NULL, p);
2350 : }
2351 :
2352 : /*********************************************************************/
2353 : /** **/
2354 : /** GCD & BEZOUT **/
2355 : /** **/
2356 : /*********************************************************************/
2357 :
2358 : GEN
2359 19872396 : lcmii(GEN x, GEN y)
2360 : {
2361 : pari_sp av;
2362 : GEN a, b;
2363 19872396 : if (!signe(x) || !signe(y)) return gen_0;
2364 19872396 : av = avma;
2365 19872396 : a = gcdii(x,y); if (!equali1(a)) y = diviiexact(y,a);
2366 19872396 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
2367 : }
2368 :
2369 : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
2370 : * set *pd = gcd(x,N) */
2371 : GEN
2372 3894053 : Fp_invgen(GEN x, GEN N, GEN *pd)
2373 : {
2374 : GEN d, d0, e, v;
2375 3894053 : if (lgefint(N) == 3)
2376 : {
2377 3358966 : ulong dd, NN = N[2], xx = umodiu(x,NN);
2378 3358966 : if (!xx) { *pd = N; return gen_0; }
2379 3358966 : xx = Fl_invgen(xx, NN, &dd);
2380 3358966 : *pd = utoi(dd); return utoi(xx);
2381 : }
2382 535087 : *pd = d = bezout(x, N, &v, NULL);
2383 535087 : if (equali1(d)) return v;
2384 : /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
2385 442936 : e = diviiexact(N,d);
2386 442936 : d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
2387 442936 : if (equali1(d0)) return v;
2388 311849 : if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
2389 311849 : return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
2390 : }
2391 :
2392 : /*********************************************************************/
2393 : /** **/
2394 : /** CHINESE REMAINDERS **/
2395 : /** **/
2396 : /*********************************************************************/
2397 :
2398 : /* Chinese Remainder Theorem. x and y must have the same type (integermod,
2399 : * polymod, or polynomial/vector/matrix recursively constructed with these
2400 : * as coefficients). Creates (with the same type) a z in the same residue
2401 : * class as x and the same residue class as y, if it is possible.
2402 : *
2403 : * We also allow (during recursion) two identical objects even if they are
2404 : * not integermod or polymod. For example:
2405 : *
2406 : * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
2407 : * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
2408 : * ? chinese(x, y)
2409 : * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
2410 :
2411 : static GEN
2412 327170 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
2413 : {
2414 327170 : GEN z = gassoc_proto(f,x,NULL);
2415 327163 : if (z == gen_1) retmkintmod(gen_0,gen_1);
2416 327128 : return z;
2417 : }
2418 :
2419 : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
2420 : * call chinese: makes Mod(0,1) a better "neutral" element */
2421 : static GEN
2422 21 : chinese_intpol(GEN x,GEN y)
2423 : {
2424 21 : pari_sp av = avma;
2425 21 : GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
2426 21 : return gerepileupto(av, chinese(z, y));
2427 : }
2428 :
2429 : GEN
2430 49 : chinese1(GEN x) { return gen_chinese(x,chinese); }
2431 :
2432 : GEN
2433 16520 : chinese(GEN x, GEN y)
2434 : {
2435 : pari_sp av;
2436 16520 : long tx = typ(x), ty;
2437 : GEN z,p1,p2,d,u,v;
2438 :
2439 16520 : if (!y) return chinese1(x);
2440 16471 : if (gidentical(x,y)) return gcopy(x);
2441 16464 : ty = typ(y);
2442 16464 : if (tx == ty) switch(tx)
2443 : {
2444 : case t_POLMOD:
2445 : {
2446 28 : GEN A = gel(x,1), B = gel(y,1);
2447 28 : GEN a = gel(x,2), b = gel(y,2);
2448 28 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
2449 28 : if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
2450 28 : av = avma;
2451 28 : d = RgX_extgcd(A,B,&u,&v);
2452 28 : p2 = gsub(b, a);
2453 28 : if (!gequal0(gmod(p2, d))) break;
2454 28 : p1 = gdiv(A,d);
2455 28 : p2 = gadd(a, gmul(gmul(u,p1), p2));
2456 :
2457 28 : z = cgetg(3, t_POLMOD);
2458 28 : gel(z,1) = gmul(p1,B);
2459 28 : gel(z,2) = gmod(p2,gel(z,1));
2460 28 : return gerepileupto(av, z);
2461 : }
2462 : case t_INTMOD:
2463 : {
2464 16401 : GEN A = gel(x,1), B = gel(y,1);
2465 16401 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
2466 16401 : z = cgetg(3,t_INTMOD);
2467 16401 : Z_chinese_pre(A, B, &C, &U, &d);
2468 16401 : c = Z_chinese_post(a, b, C, U, d);
2469 16401 : if (!c) pari_err_OP("chinese", x,y);
2470 16401 : gel(z,1) = icopy_avma(C, (pari_sp)z);
2471 16401 : gel(z,2) = icopy_avma(c, (pari_sp)gel(z,1));
2472 16401 : avma = (pari_sp)gel(z,2); return z;
2473 : }
2474 : case t_POL:
2475 : {
2476 7 : long i, lx = lg(x), ly = lg(y);
2477 7 : if (varn(x) != varn(y)) break;
2478 7 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
2479 7 : z = cgetg(lx, t_POL); z[1] = x[1];
2480 7 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2481 7 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
2482 7 : return z;
2483 : }
2484 :
2485 : case t_VEC: case t_COL: case t_MAT:
2486 : {
2487 : long i, lx;
2488 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
2489 7 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2490 7 : return z;
2491 : }
2492 : }
2493 21 : if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
2494 7 : if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
2495 0 : pari_err_OP("chinese",x,y);
2496 : return NULL; /* LCOV_EXCL_LINE */
2497 : }
2498 :
2499 : /* init chinese(Mod(.,A), Mod(.,B)) */
2500 : void
2501 237890 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
2502 : {
2503 237890 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
2504 237890 : GEN t = diviiexact(A,d);
2505 237890 : *pU = mulii(u, t);
2506 237890 : *pC = mulii(t, B);
2507 237890 : if (pd) *pd = d;
2508 237890 : }
2509 : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
2510 : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
2511 : * If d not NULL, check wether a = b mod d. */
2512 : GEN
2513 847109 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
2514 : {
2515 : GEN b_a;
2516 847109 : if (!signe(a))
2517 : {
2518 314985 : if (d && !dvdii(b, d)) return NULL;
2519 314985 : return Fp_mul(b, U, C);
2520 : }
2521 532124 : b_a = subii(b,a);
2522 532124 : if (d && !dvdii(b_a, d)) return NULL;
2523 532124 : return modii(addii(a, mulii(U, b_a)), C);
2524 : }
2525 : static ulong
2526 2159133 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
2527 : {
2528 2159133 : if (!a) return Fl_mul(b, U, C);
2529 2159133 : return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
2530 : }
2531 :
2532 : GEN
2533 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
2534 : {
2535 2142 : pari_sp av = avma;
2536 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
2537 2142 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
2538 : }
2539 : GEN
2540 219291 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
2541 : {
2542 219291 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
2543 219291 : return Z_chinese_post(a,b, *pC, U, NULL);
2544 : }
2545 :
2546 : /* return lift(chinese(a mod A, b mod B))
2547 : * assume(A,B)=1, a,b,A,B integers and C = A*B */
2548 : GEN
2549 312899 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
2550 : {
2551 312899 : pari_sp av = avma;
2552 312899 : GEN U = mulii(Fp_inv(A,B), A);
2553 312899 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2554 : }
2555 : ulong
2556 2159133 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
2557 2159133 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
2558 :
2559 : /* chinese1 for coprime moduli in Z */
2560 : static GEN
2561 296054 : chinese1_coprime_Z_aux(GEN x, GEN y)
2562 : {
2563 296054 : GEN z = cgetg(3, t_INTMOD);
2564 296054 : GEN A = gel(x,1), a = gel(x, 2);
2565 296054 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
2566 296054 : pari_sp av = avma;
2567 296054 : GEN U = mulii(Fp_inv(A,B), A);
2568 296054 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2569 296054 : gel(z,1) = C; return z;
2570 : }
2571 : GEN
2572 327121 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
2573 :
2574 : /*********************************************************************/
2575 : /** **/
2576 : /** MODULAR EXPONENTIATION **/
2577 : /** **/
2578 : /*********************************************************************/
2579 :
2580 : /* xa, ya = t_VECSMALL */
2581 : GEN
2582 484963 : ZV_producttree(GEN xa)
2583 : {
2584 484963 : long n = lg(xa)-1;
2585 484963 : long m = n==1 ? 1: expu(n-1)+1;
2586 484962 : GEN T = cgetg(m+1, t_VEC), t;
2587 : long i, j, k;
2588 484955 : t = cgetg(((n+1)>>1)+1, t_VEC);
2589 484957 : if (typ(xa)==t_VECSMALL)
2590 : {
2591 972846 : for (j=1, k=1; k<n; j++, k+=2)
2592 565110 : gel(t, j) = muluu(xa[k], xa[k+1]);
2593 407736 : if (k==n) gel(t, j) = utoi(xa[k]);
2594 : } else {
2595 377629 : for (j=1, k=1; k<n; j++, k+=2)
2596 300404 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
2597 77225 : if (k==n) gel(t, j) = icopy(gel(xa,k));
2598 : }
2599 484959 : gel(T,1) = t;
2600 917532 : for (i=2; i<=m; i++)
2601 : {
2602 432571 : GEN u = gel(T, i-1);
2603 432571 : long n = lg(u)-1;
2604 432571 : t = cgetg(((n+1)>>1)+1, t_VEC);
2605 1058516 : for (j=1, k=1; k<n; j++, k+=2)
2606 625943 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
2607 432573 : if (k==n) gel(t, j) = gel(u, k);
2608 432573 : gel(T, i) = t;
2609 : }
2610 484961 : return T;
2611 : }
2612 :
2613 : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
2614 : GEN
2615 10295675 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
2616 : {
2617 : long i,j,k;
2618 10295675 : long m = lg(T)-1, n = lg(P)-1;
2619 : GEN t;
2620 10295675 : GEN Tp = cgetg(m+1, t_VEC);
2621 10294828 : gel(Tp, m) = mkvec(A);
2622 17572713 : for (i=m-1; i>=1; i--)
2623 : {
2624 7277907 : GEN u = gel(T, i);
2625 7277907 : GEN v = gel(Tp, i+1);
2626 7277907 : long n = lg(u)-1;
2627 7277907 : t = cgetg(n+1, t_VEC);
2628 14834696 : for (j=1, k=1; k<n; j++, k+=2)
2629 : {
2630 7556041 : gel(t, k) = modii(gel(v, j), gel(u, k));
2631 7556633 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
2632 : }
2633 7278655 : if (k==n) gel(t, k) = gel(v, j);
2634 7278655 : gel(Tp, i) = t;
2635 : }
2636 : {
2637 10294806 : GEN u = gel(T, i+1);
2638 10294806 : GEN v = gel(Tp, i+1);
2639 10294806 : long l = lg(u)-1;
2640 10294806 : if (typ(P)==t_VECSMALL)
2641 : {
2642 9809980 : GEN R = cgetg(n+1, t_VECSMALL);
2643 26550034 : for (j=1, k=1; j<=l; j++, k+=2)
2644 : {
2645 16739505 : uel(R,k) = umodiu(gel(v, j), P[k]);
2646 16740248 : if (k < n)
2647 11764821 : uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
2648 : }
2649 9810529 : return R;
2650 : }
2651 : else
2652 : {
2653 484826 : GEN R = cgetg(n+1, t_VEC);
2654 1595444 : for (j=1, k=1; j<=l; j++, k+=2)
2655 : {
2656 1110619 : gel(R,k) = modii(gel(v, j), gel(P,k));
2657 1110667 : if (k < n)
2658 865402 : gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
2659 : }
2660 484825 : return R;
2661 : }
2662 : }
2663 : }
2664 :
2665 : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
2666 : GEN
2667 5960748 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
2668 : {
2669 5960748 : long m = lg(T)-1, n = lg(A)-1;
2670 : long i,j,k;
2671 5960748 : GEN Tp = cgetg(m+1, t_VEC);
2672 5933367 : GEN M = gel(T, 1);
2673 5933367 : GEN t = cgetg(lg(M), t_VEC);
2674 5925385 : if (typ(P)==t_VECSMALL)
2675 : {
2676 19340919 : for (j=1, k=1; k<n; j++, k+=2)
2677 : {
2678 16653758 : pari_sp av = avma;
2679 16653758 : GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
2680 16522071 : GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
2681 16561942 : gel(t, j) = gerepileuptoint(av, tj);
2682 : }
2683 2687161 : if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
2684 : } else
2685 : {
2686 7546216 : for (j=1, k=1; k<n; j++, k+=2)
2687 : {
2688 4312142 : pari_sp av = avma;
2689 4312142 : GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
2690 4305370 : GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
2691 4333613 : gel(t, j) = gerepileuptoint(av, tj);
2692 : }
2693 3234074 : if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
2694 : }
2695 5914788 : gel(Tp, 1) = t;
2696 13685426 : for (i=2; i<=m; i++)
2697 : {
2698 7751518 : GEN u = gel(T, i-1), M = gel(T, i);
2699 7751518 : GEN t = cgetg(lg(M), t_VEC);
2700 7815718 : GEN v = gel(Tp, i-1);
2701 7815718 : long n = lg(v)-1;
2702 24346022 : for (j=1, k=1; k<n; j++, k+=2)
2703 : {
2704 16575384 : pari_sp av = avma;
2705 66301536 : gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
2706 49726152 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
2707 : }
2708 7770638 : if (k==n) gel(t, j) = gel(v, k);
2709 7770638 : gel(Tp, i) = t;
2710 : }
2711 5933908 : return gmael(Tp,m,1);
2712 : }
2713 :
2714 : static GEN
2715 221026 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2716 : {
2717 221026 : long i, l = lg(gel(vA,1)), n = lg(P);
2718 221026 : GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
2719 5418736 : for (i=1; i < l; i++)
2720 : {
2721 5198111 : pari_sp av = avma;
2722 5198111 : GEN c, A = cgetg(n, typ(P));
2723 : long j;
2724 5222465 : for (j=1; j < n; j++) A[j] = mael(vA,j,i);
2725 5222465 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
2726 5195291 : gel(V,i) = gerepileuptoint(av, c);
2727 : }
2728 220625 : return V;
2729 : }
2730 :
2731 : static GEN
2732 90043 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2733 : {
2734 90043 : long i, j, l, n = lg(P);
2735 90043 : GEN mod = gmael(T, lg(T)-1, 1), V, w;
2736 90043 : w = cgetg(n, t_VECSMALL);
2737 90029 : for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
2738 90029 : l = vecsmall_max(w);
2739 90059 : V = cgetg(l, t_POL);
2740 90117 : V[1] = mael(vA,1,1);
2741 391328 : for (i=2; i < l; i++)
2742 : {
2743 301375 : pari_sp av = avma;
2744 301375 : GEN c, A = cgetg(n, typ(P));
2745 300501 : if (typ(P)==t_VECSMALL)
2746 5298 : for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
2747 : else
2748 295203 : for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
2749 300501 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
2750 301167 : gel(V,i) = gerepileuptoint(av, c);
2751 : }
2752 89953 : return ZX_renormalize(V, l);
2753 : }
2754 :
2755 : static GEN
2756 4493 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2757 : {
2758 4493 : long i, j, l = lg(gel(vA,1)), n = lg(P);
2759 4493 : GEN A = cgetg(n, t_VEC);
2760 4493 : GEN V = cgetg(l, t_COL);
2761 92009 : for (i=1; i < l; i++)
2762 : {
2763 87519 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2764 87519 : gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
2765 : }
2766 4490 : return V;
2767 : }
2768 :
2769 : static GEN
2770 18651 : polint_chinese(GEN worker, GEN mA, GEN P)
2771 : {
2772 18651 : long i, j, l = lg(gel(mA,1)), n = lg(P);
2773 18651 : long pending = 0, workid, cnt = 0;
2774 : struct pari_mt pt;
2775 : GEN done, va, M;
2776 18651 : GEN A = cgetg(n, t_VEC);
2777 18651 : va = mkvec(gen_0);
2778 18651 : M = cgetg(l, t_MAT);
2779 18651 : if (DEBUGLEVEL>2) err_printf("Start parallel Chinese remainder: ");
2780 18651 : mt_queue_start_lim(&pt, worker, l-1);
2781 252863 : for (i=1; i<l || pending; i++)
2782 : {
2783 234212 : for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
2784 234212 : gel(va, 1) = A;
2785 234212 : mt_queue_submit(&pt, i, i<l? va: NULL);
2786 234212 : done = mt_queue_get(&pt, &workid, &pending);
2787 234212 : if (done)
2788 : {
2789 213085 : gel(M,workid) = done;
2790 213085 : if (DEBUGLEVEL>2) err_printf("%ld%% ",(++cnt)*100/(l-1));
2791 : }
2792 : }
2793 18651 : if (DEBUGLEVEL>2) err_printf("\n");
2794 18651 : mt_queue_end(&pt);
2795 18651 : return M;
2796 : }
2797 :
2798 : GEN
2799 4178 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
2800 : {
2801 4178 : return nxCV_polint_center_tree(vA, P, T, R, m2);
2802 : }
2803 :
2804 : static GEN
2805 66 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2806 : {
2807 66 : long i, j, l = lg(gel(vA,1)), n = lg(P);
2808 66 : GEN A = cgetg(n, t_VEC);
2809 66 : GEN V = cgetg(l, t_MAT);
2810 381 : for (i=1; i < l; i++)
2811 : {
2812 315 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2813 315 : gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
2814 : }
2815 66 : return V;
2816 : }
2817 :
2818 : static GEN
2819 291 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
2820 : {
2821 291 : GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
2822 291 : return polint_chinese(worker, mA, P);
2823 : }
2824 :
2825 : static GEN
2826 757 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2827 : {
2828 757 : long i, j, l = lg(gel(vA,1)), n = lg(P);
2829 757 : GEN A = cgetg(n, t_VEC);
2830 757 : GEN V = cgetg(l, t_MAT);
2831 9193 : for (i=1; i < l; i++)
2832 : {
2833 8436 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2834 8436 : gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
2835 : }
2836 757 : return V;
2837 : }
2838 :
2839 : GEN
2840 208834 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
2841 : {
2842 208834 : return ncV_polint_center_tree(vA, P, T, R, m2);
2843 : }
2844 :
2845 : static GEN
2846 18360 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
2847 : {
2848 18360 : GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
2849 18360 : return polint_chinese(worker, mA, P);
2850 : }
2851 :
2852 : /* return [A mod P[i], i=1..#P] */
2853 : GEN
2854 0 : Z_ZV_mod(GEN A, GEN P)
2855 : {
2856 0 : pari_sp av = avma;
2857 0 : return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
2858 : }
2859 : /* P a t_VECSMALL */
2860 : GEN
2861 0 : Z_nv_mod(GEN A, GEN P)
2862 : {
2863 0 : pari_sp av = avma;
2864 0 : return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
2865 : }
2866 : /* B a ZX, T = ZV_producttree(P) */
2867 : GEN
2868 408813 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
2869 : {
2870 : pari_sp av;
2871 408813 : long i, j, l = lg(B), n = lg(A)-1;
2872 408813 : GEN V = cgetg(n+1, t_VEC);
2873 1627550 : for (j=1; j <= n; j++)
2874 : {
2875 1218627 : gel(V, j) = cgetg(l, t_VECSMALL);
2876 1218624 : mael(V, j, 1) = B[1]&VARNBITS;
2877 : }
2878 408923 : av = avma;
2879 9413109 : for (i=2; i < l; i++)
2880 : {
2881 9004297 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
2882 35413610 : for (j=1; j <= n; j++)
2883 26409424 : mael(V, j, i) = v[j];
2884 9004186 : avma = av;
2885 : }
2886 1627459 : for (j=1; j <= n; j++)
2887 1218643 : (void) Flx_renormalize(gel(V, j), l);
2888 408816 : return V;
2889 : }
2890 :
2891 : static GEN
2892 3159 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
2893 :
2894 : GEN
2895 405 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
2896 : {
2897 405 : pari_sp av = avma;
2898 405 : long i, j, l = lg(P), n = lg(xa)-1, vP = varn(P);
2899 405 : GEN V = cgetg(n+1, t_VEC);
2900 1435 : for (j=1; j <= n; j++)
2901 : {
2902 1030 : gel(V, j) = cgetg(l, t_POL);
2903 1030 : mael(V, j, 1) = vP;
2904 : }
2905 1881 : for (i=2; i < l; i++)
2906 : {
2907 1476 : GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
2908 5091 : for (j=1; j <= n; j++)
2909 3615 : gmael(V, j, i) = gel(v,j);
2910 : }
2911 1435 : for (j=1; j <= n; j++)
2912 1030 : (void) FlxX_renormalize(gel(V, j), l);
2913 405 : return gerepilecopy(av, V);
2914 : }
2915 :
2916 : GEN
2917 315 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
2918 : {
2919 315 : pari_sp av = avma;
2920 315 : long i, j, l = lg(C), n = lg(xa)-1;
2921 315 : GEN V = cgetg(n+1, t_VEC);
2922 1073 : for (j = 1; j <= n; j++)
2923 758 : gel(V, j) = cgetg(l, t_COL);
2924 1998 : for (i = 1; i < l; i++)
2925 : {
2926 1683 : GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
2927 5773 : for (j = 1; j <= n; j++)
2928 4090 : gmael(V, j, i) = gel(v,j);
2929 : }
2930 315 : return gerepilecopy(av, V);
2931 : }
2932 :
2933 : GEN
2934 66 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
2935 : {
2936 66 : pari_sp av = avma;
2937 66 : long i, j, l = lg(M), n = lg(xa)-1;
2938 66 : GEN V = cgetg(n+1, t_VEC);
2939 222 : for (j=1; j <= n; j++)
2940 156 : gel(V, j) = cgetg(l, t_MAT);
2941 381 : for (i=1; i < l; i++)
2942 : {
2943 315 : GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
2944 1073 : for (j=1; j <= n; j++)
2945 758 : gmael(V, j, i) = gel(v,j);
2946 : }
2947 66 : return gerepilecopy(av, V);
2948 : }
2949 :
2950 : /* B a ZX, T = ZV_producttree(P) */
2951 : GEN
2952 96236 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
2953 : {
2954 : pari_sp av;
2955 96236 : long i, j, l = lg(B), n = lg(A)-1;
2956 96236 : GEN V = cgetg(n+1, t_VEC);
2957 343977 : for (j=1; j <= n; j++)
2958 247756 : gel(V, j) = cgetg(l, t_VECSMALL);
2959 96221 : av = avma;
2960 902867 : for (i=1; i < l; i++)
2961 : {
2962 806629 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
2963 2904085 : for (j=1; j <= n; j++)
2964 2097439 : mael(V, j, i) = v[j];
2965 806646 : avma = av;
2966 : }
2967 96238 : return V;
2968 : }
2969 :
2970 : GEN
2971 12017 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
2972 : {
2973 12017 : pari_sp av = avma;
2974 12017 : long i, j, l = lg(M), n = lg(xa)-1;
2975 12017 : GEN V = cgetg(n+1, t_VEC);
2976 42606 : for (j=1; j <= n; j++)
2977 30591 : gel(V, j) = cgetg(l, t_MAT);
2978 108253 : for (i=1; i < l; i++)
2979 : {
2980 96236 : GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
2981 344023 : for (j=1; j <= n; j++)
2982 247785 : gmael(V, j, i) = gel(v,j);
2983 : }
2984 12017 : return gerepilecopy(av, V);
2985 : }
2986 :
2987 : static GEN
2988 481476 : ZV_sqr(GEN z)
2989 : {
2990 481476 : long i,l = lg(z);
2991 481476 : GEN x = cgetg(l, t_VEC);
2992 481421 : if (typ(z)==t_VECSMALL)
2993 407605 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
2994 : else
2995 73816 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
2996 481398 : return x;
2997 : }
2998 :
2999 : static GEN
3000 3198192 : ZT_sqr(GEN z)
3001 : {
3002 3198192 : if (typ(z) == t_INT)
3003 1809218 : return sqri(z);
3004 : else
3005 : {
3006 1388974 : long i,l = lg(z);
3007 1388974 : GEN x = cgetg(l, t_VEC);
3008 1388892 : for (i=1; i<l; i++) gel(x,i) = ZT_sqr(gel(z,i));
3009 1390014 : return x;
3010 : }
3011 : }
3012 :
3013 : static GEN
3014 481404 : ZV_invdivexact(GEN y, GEN x)
3015 : {
3016 481404 : long i, l = lg(y);
3017 481404 : GEN z = cgetg(l,t_VEC);
3018 481473 : if (typ(x)==t_VECSMALL)
3019 1747984 : for (i=1; i<l; i++)
3020 : {
3021 1340422 : pari_sp av = avma;
3022 1340422 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
3023 1340568 : avma = av;
3024 1340568 : gel(z,i) = utoipos(a);
3025 : }
3026 : else
3027 684908 : for (i=1; i<l; i++)
3028 611092 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
3029 481378 : return z;
3030 : }
3031 :
3032 : /* P t_VECSMALL or t_VEC of t_INT */
3033 : GEN
3034 481026 : ZV_chinesetree(GEN P, GEN T)
3035 : {
3036 481026 : GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
3037 481403 : GEN mod = gmael(T,lg(T)-1,1);
3038 481403 : return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
3039 : }
3040 :
3041 : static GEN
3042 80135 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
3043 : {
3044 80135 : if (!pt_mod)
3045 2559 : return gerepileupto(av, a);
3046 : else
3047 : {
3048 77576 : GEN mod = gmael(T, lg(T)-1, 1);
3049 77576 : gerepileall(av, 2, &a, &mod);
3050 77576 : *pt_mod = mod;
3051 77576 : return a;
3052 : }
3053 : }
3054 :
3055 : GEN
3056 42630 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3057 : {
3058 42630 : pari_sp av = avma;
3059 42630 : GEN T = ZV_producttree(P);
3060 42630 : GEN R = ZV_chinesetree(P, T);
3061 42630 : GEN a = ZV_chinese_tree(A, P, T, R);
3062 42630 : GEN mod = gmael(T, lg(T)-1, 1);
3063 42630 : GEN ca = Fp_center(a, mod, shifti(mod,-1));
3064 42630 : return gc_chinese(av, T, ca, pt_mod);
3065 : }
3066 :
3067 : GEN
3068 13073 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
3069 : {
3070 13073 : pari_sp av = avma;
3071 13073 : GEN T = ZV_producttree(P);
3072 13073 : GEN R = ZV_chinesetree(P, T);
3073 13073 : GEN a = ZV_chinese_tree(A, P, T, R);
3074 13073 : return gc_chinese(av, T, a, pt_mod);
3075 : }
3076 :
3077 : GEN
3078 405 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
3079 : {
3080 405 : pari_sp av = avma;
3081 405 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3082 405 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
3083 405 : return gerepileupto(av, a);
3084 : }
3085 :
3086 : GEN
3087 2021 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3088 : {
3089 2021 : pari_sp av = avma;
3090 2021 : GEN T = ZV_producttree(P);
3091 2021 : GEN R = ZV_chinesetree(P, T);
3092 2021 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3093 2021 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
3094 2021 : return gc_chinese(av, T, a, pt_mod);
3095 : }
3096 :
3097 : GEN
3098 3760 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3099 : {
3100 3760 : pari_sp av = avma;
3101 3760 : GEN T = ZV_producttree(P);
3102 3760 : GEN R = ZV_chinesetree(P, T);
3103 3760 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3104 3760 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
3105 3760 : return gc_chinese(av, T, a, pt_mod);
3106 : }
3107 :
3108 : GEN
3109 0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
3110 : {
3111 0 : pari_sp av = avma;
3112 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3113 0 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
3114 0 : return gerepileupto(av, a);
3115 : }
3116 :
3117 : GEN
3118 757 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
3119 : {
3120 757 : pari_sp av = avma;
3121 757 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3122 757 : GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
3123 757 : return gerepileupto(av, a);
3124 : }
3125 :
3126 : GEN
3127 18360 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3128 : {
3129 18360 : pari_sp av = avma;
3130 18360 : GEN T = ZV_producttree(P);
3131 18360 : GEN R = ZV_chinesetree(P, T);
3132 18360 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3133 18360 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
3134 18360 : return gc_chinese(av, T, a, pt_mod);
3135 : }
3136 :
3137 : GEN
3138 66 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
3139 : {
3140 66 : pari_sp av = avma;
3141 66 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3142 66 : GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
3143 66 : return gerepileupto(av, a);
3144 : }
3145 :
3146 : GEN
3147 291 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3148 : {
3149 291 : pari_sp av = avma;
3150 291 : GEN T = ZV_producttree(P);
3151 291 : GEN R = ZV_chinesetree(P, T);
3152 291 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3153 291 : GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
3154 291 : return gc_chinese(av, T, a, pt_mod);
3155 : }
3156 :
3157 : /**********************************************************************
3158 : ** **
3159 : ** Powering over (Z/NZ)^*, small N **
3160 : ** **
3161 : **********************************************************************/
3162 :
3163 : /* 2^n mod p; assume n > 1 */
3164 : static ulong
3165 20781911 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
3166 : {
3167 20781911 : ulong y = 2;
3168 20781911 : int j = 1+bfffo(n);
3169 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
3170 20781911 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
3171 421759804 : for (; j; n<<=1,j--)
3172 : {
3173 400977871 : y = Fl_sqr_pre(y,p,pi);
3174 400977895 : if (n & HIGHBIT) y = Fl_double(y, p);
3175 : }
3176 20781933 : return y;
3177 : }
3178 :
3179 : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
3180 : static ulong
3181 2596993 : Fl_2powu(ulong n, ulong p)
3182 : {
3183 2596993 : ulong y = 2;
3184 2596993 : int j = 1+bfffo(n);
3185 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
3186 2596993 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
3187 30090880 : for (; j; n<<=1,j--)
3188 : {
3189 27475269 : y = (y*y) % p;
3190 27475269 : if (n & HIGHBIT) y = Fl_double(y, p);
3191 : }
3192 2615611 : return y;
3193 : }
3194 :
3195 : ulong
3196 96725544 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
3197 : {
3198 : ulong y, z, n;
3199 96725544 : if (n0 <= 1)
3200 : { /* frequent special cases */
3201 10891664 : if (n0 == 1) return x;
3202 2984013 : if (n0 == 0) return 1;
3203 : }
3204 85833880 : if (x <= 2)
3205 : {
3206 22574265 : if (x == 2) return Fl_2powu_pre(n0, p, pi);
3207 1792341 : return x; /* 0 or 1 */
3208 : }
3209 63259615 : y = 1; z = x; n = n0;
3210 : for(;;)
3211 : {
3212 1061916445 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
3213 562662814 : n>>=1; if (!n) return y;
3214 499403998 : z = Fl_sqr_pre(z,p,pi);
3215 : }
3216 : }
3217 :
3218 : ulong
3219 43949131 : Fl_powu(ulong x, ulong n0, ulong p)
3220 : {
3221 : ulong y, z, n;
3222 43949131 : if (n0 <= 2)
3223 : { /* frequent special cases */
3224 32238929 : if (n0 == 2) return Fl_sqr(x,p);
3225 3630529 : if (n0 == 1) return x;
3226 61963 : if (n0 == 0) return 1;
3227 : }
3228 11710202 : if (x <= 1) return x; /* 0 or 1 */
3229 11657890 : if (p & HIGHMASK)
3230 6299654 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
3231 5358236 : if (x == 2) return Fl_2powu(n0, p);
3232 2762157 : y = 1; z = x; n = n0;
3233 : for(;;)
3234 : {
3235 49285999 : if (n&1) y = (y*z) % p;
3236 26024078 : n>>=1; if (!n) return y;
3237 23261921 : z = (z*z) % p;
3238 : }
3239 : }
3240 :
3241 : /* Reduce data dependency to maximize internal parallelism */
3242 : GEN
3243 11052688 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
3244 : {
3245 : long i, k;
3246 11052688 : GEN powers = cgetg(n + 2, t_VECSMALL);
3247 11047329 : powers[1] = 1; if (n == 0) return powers;
3248 11047329 : powers[2] = x;
3249 46570770 : for (i = 3, k=2; i <= n; i+=2, k++)
3250 : {
3251 35511298 : powers[i] = Fl_sqr_pre(powers[k], p, pi);
3252 35517499 : powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
3253 : }
3254 11059472 : if (i==n+1)
3255 9662217 : powers[i] = Fl_sqr_pre(powers[k], p, pi);
3256 11059394 : return powers;
3257 : }
3258 :
3259 : GEN
3260 2961 : Fl_powers(ulong x, long n, ulong p)
3261 : {
3262 2961 : return Fl_powers_pre(x, n, p, get_Fl_red(p));
3263 : }
3264 :
3265 : /**********************************************************************
3266 : ** **
3267 : ** Powering over (Z/NZ)^*, large N **
3268 : ** **
3269 : **********************************************************************/
3270 :
3271 : static GEN
3272 4259221 : Fp_dblsqr(GEN x, GEN N)
3273 : {
3274 4259221 : GEN z = shifti(Fp_sqr(x, N), 1);
3275 4259221 : return cmpii(z, N) >= 0? subii(z, N): z;
3276 : }
3277 :
3278 : typedef struct muldata {
3279 : GEN (*sqr)(void * E, GEN x);
3280 : GEN (*mul)(void * E, GEN x, GEN y);
3281 : GEN (*mul2)(void * E, GEN x);
3282 : } muldata;
3283 :
3284 : /* modified Barrett reduction with one fold */
3285 : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
3286 :
3287 : static GEN
3288 10092 : Fp_invmBarrett(GEN p, long s)
3289 : {
3290 10092 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
3291 10092 : return mkvec2(Q,R);
3292 : }
3293 :
3294 : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
3295 : * a = r (mod N) */
3296 : static GEN
3297 4235005 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
3298 : {
3299 4235005 : pari_sp av = avma;
3300 4235005 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
3301 4235005 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
3302 4235005 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
3303 4235005 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
3304 4235005 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
3305 4235005 : GEN r = subii(A, mulii(q, N));
3306 4235005 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
3307 4235005 : if (signe(sr)<0) return gerepileuptoint(av, r);
3308 2548401 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
3309 2548401 : if (signe(sr)<0) return gerepileuptoint(av, r);
3310 98169 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
3311 98169 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
3312 : }
3313 :
3314 : /* Montgomery reduction */
3315 :
3316 : INLINE ulong
3317 271268 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
3318 :
3319 : struct montred
3320 : {
3321 : GEN N;
3322 : ulong inv;
3323 : };
3324 :
3325 : /* Montgomery reduction */
3326 : static GEN
3327 29575561 : _sqr_montred(void * E, GEN x)
3328 : {
3329 29575561 : struct montred * D = (struct montred *) E;
3330 29575561 : return red_montgomery(sqri(x), D->N, D->inv);
3331 : }
3332 :
3333 : /* Montgomery reduction */
3334 : static GEN
3335 2460301 : _mul_montred(void * E, GEN x, GEN y)
3336 : {
3337 2460301 : struct montred * D = (struct montred *) E;
3338 2460301 : return red_montgomery(mulii(x, y), D->N, D->inv);
3339 : }
3340 :
3341 : static GEN
3342 5816539 : _mul2_montred(void * E, GEN x)
3343 : {
3344 5816539 : struct montred * D = (struct montred *) E;
3345 5816539 : GEN z = shifti(_sqr_montred(E, x), 1);
3346 5816521 : long l = lgefint(D->N);
3347 5816521 : while (lgefint(z) > l) z = subii(z, D->N);
3348 5816522 : return z;
3349 : }
3350 :
3351 : static GEN
3352 10533844 : _sqr_remii(void* N, GEN x)
3353 10533844 : { return remii(sqri(x), (GEN) N); }
3354 :
3355 : static GEN
3356 938683 : _mul_remii(void* N, GEN x, GEN y)
3357 938683 : { return remii(mulii(x, y), (GEN) N); }
3358 :
3359 : static GEN
3360 2997472 : _mul2_remii(void* N, GEN x)
3361 2997472 : { return Fp_dblsqr(x, (GEN) N); }
3362 :
3363 : struct redbarrett
3364 : {
3365 : GEN iM, N;
3366 : long s;
3367 : };
3368 :
3369 : static GEN
3370 3823998 : _sqr_remiibar(void *E, GEN x)
3371 : {
3372 3823998 : struct redbarrett * D = (struct redbarrett *) E;
3373 3823998 : return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
3374 : }
3375 :
3376 : static GEN
3377 411007 : _mul_remiibar(void *E, GEN x, GEN y)
3378 : {
3379 411007 : struct redbarrett * D = (struct redbarrett *) E;
3380 411007 : return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
3381 : }
3382 :
3383 : static GEN
3384 1261749 : _mul2_remiibar(void *E, GEN x)
3385 : {
3386 1261749 : struct redbarrett * D = (struct redbarrett *) E;
3387 1261749 : return Fp_dblsqr(x, D->N);
3388 : }
3389 :
3390 : static long
3391 358012 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
3392 : {
3393 358012 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
3394 : {
3395 10092 : struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
3396 10092 : D->sqr = &_sqr_remiibar;
3397 10092 : D->mul = &_mul_remiibar;
3398 10092 : D->mul2 = &_mul2_remiibar;
3399 10092 : E->N = N;
3400 10092 : E->s = 1+(expi(N)>>1);
3401 10092 : E->iM = Fp_invmBarrett(N, E->s);
3402 10092 : *pt_E = (void*) E;
3403 10092 : return 0;
3404 : }
3405 347920 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
3406 : {
3407 271264 : struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
3408 271264 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
3409 271267 : D->sqr = &_sqr_montred;
3410 271267 : D->mul = &_mul_montred;
3411 271267 : D->mul2 = &_mul2_montred;
3412 271267 : E->N = N;
3413 271267 : E->inv = init_montdata(N);
3414 271267 : *pt_E = (void*) E;
3415 271267 : return 1;
3416 : }
3417 : else
3418 : {
3419 76655 : D->sqr = &_sqr_remii;
3420 76655 : D->mul = &_mul_remii;
3421 76655 : D->mul2 = &_mul2_remii;
3422 76655 : *pt_E = (void*) N;
3423 76655 : return 0;
3424 : }
3425 : }
3426 :
3427 : GEN
3428 882897 : Fp_powu(GEN A, ulong k, GEN N)
3429 : {
3430 882897 : long lN = lgefint(N);
3431 : int base_is_2, use_montgomery;
3432 : muldata D;
3433 : void *E;
3434 : pari_sp av;
3435 :
3436 882897 : if (lN == 3) {
3437 88481 : ulong n = uel(N,2);
3438 88481 : return utoi( Fl_powu(umodiu(A, n), k, n) );
3439 : }
3440 794416 : if (k <= 2)
3441 : { /* frequent special cases */
3442 565241 : if (k == 2) return Fp_sqr(A,N);
3443 131277 : if (k == 1) return A;
3444 0 : if (k == 0) return gen_1;
3445 : }
3446 229175 : av = avma; A = modii(A,N);
3447 229176 : base_is_2 = 0;
3448 229176 : if (lgefint(A) == 3) switch(A[2])
3449 : {
3450 523 : case 1: avma = av; return gen_1;
3451 34202 : case 2: base_is_2 = 1; break;
3452 : }
3453 :
3454 : /* TODO: Move this out of here and use for general modular computations */
3455 228653 : use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
3456 228653 : if (base_is_2)
3457 34202 : A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
3458 : else
3459 194451 : A = gen_powu_i(A, k, E, D.sqr, D.mul);
3460 228652 : if (use_montgomery)
3461 : {
3462 198180 : A = red_montgomery(A, N, ((struct montred *) E)->inv);
3463 198178 : if (cmpii(A, N) >= 0) A = subii(A,N);
3464 : }
3465 228650 : return gerepileuptoint(av, A);
3466 : }
3467 :
3468 : GEN
3469 21623 : Fp_pows(GEN A, long k, GEN N)
3470 : {
3471 21623 : if (lgefint(N) == 3) {
3472 7372 : ulong n = N[2];
3473 7372 : ulong a = umodiu(A, n);
3474 7372 : if (k < 0) {
3475 126 : a = Fl_inv(a, n);
3476 126 : k = -k;
3477 : }
3478 7372 : return utoi( Fl_powu(a, (ulong)k, n) );
3479 : }
3480 14251 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
3481 14251 : return Fp_powu(A, (ulong)k, N);
3482 : }
3483 :
3484 : /* A^K mod N */
3485 : GEN
3486 5528198 : Fp_pow(GEN A, GEN K, GEN N)
3487 : {
3488 : pari_sp av;
3489 5528198 : long s, lN = lgefint(N), sA;
3490 : int base_is_2, use_montgomery;
3491 : GEN y;
3492 : muldata D;
3493 : void *E;
3494 :
3495 5528198 : s = signe(K);
3496 5528198 : if (!s) return dvdii(A,N)? gen_0: gen_1;
3497 5451086 : if (lN == 3 && lgefint(K) == 3)
3498 : {
3499 5200168 : ulong n = N[2], a = umodiu(A, n);
3500 5200168 : if (s < 0) a = Fl_inv(a, n);
3501 5200168 : if (a <= 1) return utoi(a); /* 0 or 1 */
3502 4867932 : return utoi(Fl_powu(a, uel(K,2), n));
3503 : }
3504 :
3505 250918 : av = avma;
3506 250918 : if (s < 0) y = Fp_inv(A,N);
3507 : else
3508 : {
3509 250450 : y = modii(A,N);
3510 250452 : if (!signe(y)) { avma = av; return gen_0; }
3511 : }
3512 250918 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
3513 :
3514 129441 : base_is_2 = 0;
3515 129441 : sA = signe(y)==-1 && mod2(K);
3516 129441 : if (lgefint(y) == 3) switch(y[2])
3517 : {
3518 82 : case 1: return sA ? gen_m1 : gen_1;
3519 87415 : case 2: base_is_2 = 1; break;
3520 : }
3521 :
3522 : /* TODO: Move this out of here and use for general modular computations */
3523 129359 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
3524 129361 : if (base_is_2)
3525 87415 : y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
3526 : else
3527 41946 : y = gen_pow_i(y, K, E, D.sqr, D.mul);
3528 129359 : if (use_montgomery)
3529 : {
3530 73084 : y = red_montgomery(y, N, ((struct montred *) E)->inv);
3531 73084 : if (cmpii(y,N) >= 0) y = subii(y,N);
3532 73082 : if (sA) y = subii(N, y);
3533 : }
3534 129357 : return gerepileuptoint(av,y);
3535 : }
3536 :
3537 : static GEN
3538 1247536 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
3539 :
3540 : static GEN
3541 23500 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
3542 :
3543 : static GEN
3544 55230 : _Fp_one(void *E) { (void) E; return gen_1; }
3545 :
3546 : GEN
3547 84 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
3548 : {
3549 84 : return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul);
3550 : }
3551 :
3552 : GEN
3553 55090 : Fp_pow_table(GEN R, GEN n, GEN p)
3554 : {
3555 55090 : return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul);
3556 : }
3557 :
3558 : GEN
3559 2016 : Fp_powers(GEN x, long n, GEN p)
3560 : {
3561 2016 : if (lgefint(p) == 3)
3562 1876 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
3563 140 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
3564 : }
3565 :
3566 : static GEN
3567 3179631 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
3568 :
3569 : static GEN
3570 112 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
3571 :
3572 : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
3573 :
3574 : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
3575 : equalii,equali1,Fp_easylog};
3576 :
3577 : static GEN
3578 1052565 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
3579 :
3580 : static GEN
3581 1504160 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
3582 :
3583 : static GEN
3584 268331 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
3585 :
3586 : static GEN
3587 1665214 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
3588 :
3589 : static GEN
3590 27616 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
3591 :
3592 : static int
3593 385860 : _Fp_equal0(GEN x) { return signe(x)==0; }
3594 :
3595 : static GEN
3596 150543 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
3597 :
3598 : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
3599 : _Fp_inv,_Fp_equal0,_Fp_s};
3600 :
3601 7371 : const struct bb_field *get_Fp_field(void **E, GEN p)
3602 : {
3603 7371 : *E = (void*)p; return &Fp_field;
3604 : }
3605 :
3606 : /*********************************************************************/
3607 : /** **/
3608 : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
3609 : /** **/
3610 : /*********************************************************************/
3611 : ulong
3612 12096 : Fl_order(ulong a, ulong o, ulong p)
3613 : {
3614 12096 : pari_sp av = avma;
3615 : GEN m, P, E;
3616 : long i;
3617 12096 : if (!o) o = p-1;
3618 12096 : m = factoru(o);
3619 12096 : P = gel(m,1);
3620 12096 : E = gel(m,2);
3621 28952 : for (i = lg(P)-1; i; i--)
3622 : {
3623 16856 : ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
3624 16856 : if (y == 1) o = t;
3625 14861 : else for (j = 1; j < e; j++)
3626 : {
3627 4284 : y = Fl_powu(y, l, p);
3628 4284 : if (y == 1) { o = t * upowuu(l, j); break; }
3629 : }
3630 : }
3631 12096 : avma = av; return o;
3632 : }
3633 :
3634 : /*Find the exact order of a assuming a^o==1*/
3635 : GEN
3636 10676 : Fp_order(GEN a, GEN o, GEN p) {
3637 10676 : if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
3638 : {
3639 21 : ulong pp = p[2], oo = (o && lgefint(o)==3)? o[2]: pp-1;
3640 21 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
3641 : }
3642 10655 : return gen_order(a, o, (void*)p, &Fp_star);
3643 : }
3644 : GEN
3645 56 : Fp_factored_order(GEN a, GEN o, GEN p)
3646 56 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
3647 :
3648 : /* return order of a mod p^e, e > 0, pe = p^e */
3649 : static GEN
3650 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
3651 : {
3652 : GEN ap, op;
3653 70 : if (absequaliu(p, 2))
3654 : {
3655 56 : if (e == 1) return gen_1;
3656 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
3657 49 : if (mod4(a) == 1)
3658 14 : op = gen_1;
3659 : else {
3660 35 : op = gen_2;
3661 35 : a = Fp_sqr(a, pe);
3662 : }
3663 : } else {
3664 14 : ap = (e == 1)? a: remii(a,p);
3665 14 : op = Fp_order(ap, subiu(p,1), p);
3666 14 : if (e == 1) return op;
3667 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
3668 : }
3669 49 : if (equali1(a)) return op;
3670 7 : return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
3671 : }
3672 :
3673 : GEN
3674 63 : znorder(GEN x, GEN o)
3675 : {
3676 63 : pari_sp av = avma;
3677 : GEN b, a;
3678 :
3679 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
3680 56 : b = gel(x,1); a = gel(x,2);
3681 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
3682 49 : if (!o)
3683 : {
3684 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
3685 35 : long i, l = lg(P);
3686 35 : o = gen_1;
3687 70 : for (i = 1; i < l; i++)
3688 : {
3689 35 : GEN p = gel(P,i);
3690 35 : long e = itos(gel(E,i));
3691 :
3692 35 : if (l == 2)
3693 35 : o = Zp_order(a, p, e, b);
3694 : else {
3695 0 : GEN pe = powiu(p,e);
3696 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
3697 : }
3698 : }
3699 35 : return gerepileuptoint(av, o);
3700 : }
3701 14 : return Fp_order(a, o, b);
3702 : }
3703 : GEN
3704 0 : order(GEN x) { return znorder(x, NULL); }
3705 :
3706 : /*********************************************************************/
3707 : /** **/
3708 : /** DISCRETE LOGARITHM in (Z/nZ)* **/
3709 : /** **/
3710 : /*********************************************************************/
3711 : static GEN
3712 59520 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
3713 : {
3714 59520 : pari_sp av = avma;
3715 : GEN h1, h2, F, G;
3716 59520 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return NULL;
3717 35723 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
3718 : {
3719 226 : GEN M = cgetg(3, t_MAT);
3720 226 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
3721 226 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
3722 226 : return gerepileupto(av, M);
3723 : }
3724 35497 : avma = av; return NULL;
3725 : }
3726 :
3727 : static GEN
3728 59520 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
3729 : {
3730 : GEN rel;
3731 : do
3732 : {
3733 59520 : (*e)++; *g = Fp_mul(*g, b, p);
3734 59520 : rel = Fp_log_halfgcd(bnd, C, *g, p);
3735 59520 : } while (!rel);
3736 226 : return rel;
3737 : }
3738 :
3739 : struct Fp_log_rel
3740 : {
3741 : GEN rel;
3742 : ulong prmax;
3743 : long nbrel, nbmax, nbgen;
3744 : };
3745 :
3746 : /* add u^e */
3747 : static void
3748 2583 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
3749 : {
3750 2583 : pari_sp av = avma;
3751 2583 : long off = r->prmax+1;
3752 2583 : GEN F = cgetg(3, t_MAT);
3753 2583 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
3754 2583 : gel(F,2) = vecsmall_append(gel(z,2), e);
3755 2583 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3756 2583 : }
3757 :
3758 : /* add u^-1 v^-1 */
3759 : static void
3760 99869 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
3761 : {
3762 99869 : pari_sp av = avma;
3763 99869 : long off = r->prmax+1;
3764 99869 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
3765 99869 : GEN F = cgetg(3, t_MAT);
3766 99869 : gel(F,1) = vecsmall_concat(gel(z,1), P);
3767 99869 : gel(F,2) = vecsmall_concat(gel(z,2), E);
3768 99869 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3769 99869 : }
3770 :
3771 : /*
3772 : Let p=C^2+c
3773 : Solve h = (C+x)*(C+a)-p = 0 [mod l]
3774 : h= -c+x*(C+a)+C*a = 0 [mod l]
3775 : x = (c-C*a)/(C+a) [mod l]
3776 : h = -c+C*(x+a)+a*x
3777 : */
3778 :
3779 : GEN
3780 39020 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
3781 : {
3782 39020 : pari_sp ltop = avma;
3783 39020 : long th = expi(mulis(C,a)), n = lg(pi)-1;
3784 : long i, j;
3785 39017 : GEN sieve = zero_zv(a+2)+1;
3786 39036 : GEN L = cgetg(1+a+2, t_VEC);
3787 39022 : pari_sp av = avma;
3788 39022 : long rel = 1;
3789 : GEN z, h;
3790 39022 : h = addis(C,a);
3791 39004 : if ((z = Z_issmooth_fact(h, prmax)))
3792 : {
3793 2415 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
3794 2415 : av = avma;
3795 : }
3796 16761861 : for(i=1; i<=n; i++)
3797 : {
3798 16722896 : ulong li = pi[i], s = sz[i], al = a % li;
3799 16722896 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
3800 17246358 : if (!iv) continue;
3801 16783209 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
3802 76677994 : for(j = u; j<=a; j+=li)
3803 60411034 : sieve[j] += s;
3804 : }
3805 38965 : th = th - expu(th)-1;
3806 27880928 : for(j=0; j<a; j++)
3807 27841894 : if (sieve[j]>=th)
3808 : {
3809 116710 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
3810 111523 : if ((z = Z_issmooth_fact(h, prmax)))
3811 : {
3812 104741 : gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
3813 104908 : av = avma;
3814 6644 : } else avma = av;
3815 : }
3816 : /* j = a */
3817 39034 : if (sieve[a]>=th)
3818 : {
3819 420 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
3820 420 : if ((z = Z_issmooth_fact(h, prmax)))
3821 : {
3822 343 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
3823 343 : av = avma;
3824 : }
3825 : }
3826 39034 : setlg(L, rel);
3827 39036 : return gerepilecopy(ltop, L);
3828 : }
3829 :
3830 : static long
3831 49 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
3832 : {
3833 : struct pari_mt pt;
3834 49 : long i, j, nb = 0;
3835 49 : GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
3836 : mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
3837 49 : long running, pending = 0;
3838 49 : GEN W = zerovec(r->nbgen);
3839 49 : mt_queue_start_lim(&pt, worker, r->nbgen);
3840 39305 : for (i = 0; (running = (i < r->nbgen)) || pending; i++)
3841 : {
3842 : GEN done;
3843 : long idx;
3844 39256 : mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
3845 39256 : done = mt_queue_get(&pt, &idx, &pending);
3846 39256 : if (!done || lg(done)==1) continue;
3847 34146 : gel(W, idx+1) = done;
3848 34146 : nb += lg(done)-1;
3849 34146 : if (DEBUGLEVEL && (i&127)==0)
3850 0 : err_printf("%ld%% ",100*nb/r->nbmax);
3851 : }
3852 49 : mt_queue_end(&pt);
3853 37870 : for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
3854 : {
3855 : long ll, m;
3856 37821 : GEN L = gel(W,j);
3857 37821 : if (isintzero(L)) continue;
3858 32942 : ll = lg(L);
3859 135394 : for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
3860 : {
3861 102452 : GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
3862 102452 : if (v[1] == 1)
3863 2583 : addifsmooth1(r, h, v[2], v[3]);
3864 : else
3865 99869 : addifsmooth2(r, h, v[2], v[3]);
3866 : }
3867 : }
3868 49 : return j;
3869 : }
3870 :
3871 : static GEN
3872 525 : ECP_psi(GEN x, GEN y)
3873 : {
3874 525 : long prec = realprec(x);
3875 525 : GEN lx = glog(x, prec), ly = glog(y, prec);
3876 525 : GEN u = gdiv(lx, ly);
3877 525 : return gpow(u, gneg(u),prec);
3878 : }
3879 :
3880 : struct computeG
3881 : {
3882 : GEN C;
3883 : long bnd, nbi;
3884 : };
3885 :
3886 : static GEN
3887 525 : _computeG(void *E, GEN gen)
3888 : {
3889 525 : struct computeG * d = (struct computeG *) E;
3890 525 : GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
3891 525 : return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
3892 : }
3893 :
3894 : static long
3895 49 : compute_nbgen(GEN C, long bnd, long nbi)
3896 : {
3897 : struct computeG d;
3898 49 : d.C = shifti(C, 1);
3899 49 : d.bnd = bnd;
3900 49 : d.nbi = nbi;
3901 49 : return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
3902 : }
3903 :
3904 : static GEN
3905 1367 : _psi(void*E, GEN y)
3906 : {
3907 1367 : GEN lx = (GEN) E;
3908 1367 : long prec = realprec(lx);
3909 1367 : GEN ly = glog(y, prec);
3910 1367 : GEN u = gdiv(lx, ly);
3911 1367 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
3912 : }
3913 :
3914 : static GEN
3915 49 : opt_param(GEN x, long prec)
3916 : {
3917 49 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
3918 : }
3919 :
3920 : static GEN
3921 49 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
3922 : {
3923 49 : pari_sp av = avma;
3924 49 : long lM = lg(M)-1, nbcol = lM;
3925 49 : long tbs = maxss(1, expu(nbg/expi(m)));
3926 : for (;;)
3927 35 : {
3928 84 : GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
3929 : GEN tab;
3930 84 : long i, f=0;
3931 84 : long l = lg(K), lm = lgefint(m);
3932 84 : GEN idx = diviiexact(subiu(p,1),m), g;
3933 : pari_timer ti;
3934 84 : if (DEBUGLEVEL) timer_start(&ti);
3935 140 : for(i=1; i<l; i++)
3936 140 : if (signe(gel(K,i)))
3937 84 : break;
3938 84 : g = Fp_pow(utoi(i), idx, p);
3939 84 : tab = Fp_pow_init(g, p, tbs, p);
3940 84 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
3941 130438 : for(i=1; i<l; i++)
3942 : {
3943 130354 : GEN k = gel(K,i);
3944 130354 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
3945 130354 : if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
3946 80542 : gel(K,i) = cgetineg(lm);
3947 : else
3948 49812 : f++;
3949 : }
3950 84 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
3951 133 : if(f > (nbg>>1)) return gerepileupto(av, K);
3952 9877 : for(i=1; i<=nbcol; i++)
3953 : {
3954 9842 : long a = 1+random_Fl(lM);
3955 9842 : swap(gel(M,a),gel(M,i));
3956 : }
3957 35 : if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
3958 : }
3959 : }
3960 :
3961 : static GEN
3962 98 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
3963 : {
3964 98 : pari_sp av=avma;
3965 98 : GEN aa = gen_1;
3966 98 : long AV = 0;
3967 : for(;;)
3968 128 : {
3969 226 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
3970 226 : GEN F = gel(A,1), E = gel(A,2);
3971 226 : GEN Ao = gen_0;
3972 226 : long i, l = lg(F);
3973 1145 : for(i=1; i<l; i++)
3974 : {
3975 1047 : GEN Ki = gel(K,F[i]);
3976 1047 : if (signe(Ki)<0) break;
3977 919 : Ao = addii(Ao, mulis(Ki, E[i]));
3978 : }
3979 324 : if (i==l) return Fp_div(Ao, utoi(AV), m);
3980 128 : aa = gerepileuptoint(av, aa);
3981 : }
3982 : }
3983 :
3984 : static GEN
3985 49 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
3986 : {
3987 49 : pari_sp av = avma, av2;
3988 49 : long i, j, nbi, nbr = 0, nbrow, nbg;
3989 : GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
3990 : pari_timer ti;
3991 : struct Fp_log_rel r;
3992 49 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
3993 49 : ulong bnd = 4*bnds;
3994 49 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
3995 :
3996 49 : p_1 = subiu(p,1);
3997 49 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
3998 0 : m = diviiexact(p_1, Z_ppo(p_1, m));
3999 49 : pr = primes_upto_zv(bnd);
4000 49 : nbi = lg(pr)-1;
4001 49 : C = sqrtremi(p, &c);
4002 49 : av2 = avma;
4003 12236 : for (i = 1; i <= nbi; ++i)
4004 : {
4005 12187 : ulong lp = pr[i];
4006 37793 : while (lp <= bnd)
4007 : {
4008 13419 : nbr++;
4009 13419 : lp *= pr[i];
4010 : }
4011 : }
4012 49 : pi = cgetg(nbr+1,t_VECSMALL);
4013 49 : Ci = cgetg(nbr+1,t_VECSMALL);
4014 49 : ci = cgetg(nbr+1,t_VECSMALL);
4015 49 : sz = cgetg(nbr+1,t_VECSMALL);
4016 12236 : for (i = 1, j = 1; i <= nbi; ++i)
4017 : {
4018 12187 : ulong lp = pr[i], sp = expu(2*lp-1);
4019 37793 : while (lp <= bnd)
4020 : {
4021 13419 : pi[j] = lp;
4022 13419 : Ci[j] = umodiu(C, lp);
4023 13419 : ci[j] = umodiu(c, lp);
4024 13419 : sz[j] = sp;
4025 13419 : lp *= pr[i];
4026 13419 : j++;
4027 : }
4028 : }
4029 49 : r.nbrel = 0;
4030 49 : r.nbgen = compute_nbgen(C, bnd, nbi);
4031 49 : r.nbmax = 2*(nbi+r.nbgen);
4032 49 : r.rel = cgetg(r.nbmax+1,t_VEC);
4033 49 : r.prmax = pr[nbi];
4034 49 : if (DEBUGLEVEL)
4035 : {
4036 0 : err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
4037 0 : timer_start(&ti);
4038 : }
4039 49 : nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
4040 49 : nbrow = r.prmax + nbg;
4041 49 : if (DEBUGLEVEL)
4042 : {
4043 0 : err_printf("\n");
4044 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
4045 : }
4046 49 : setlg(r.rel,r.nbrel+1);
4047 49 : r.rel = gerepilecopy(av2, r.rel);
4048 49 : K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
4049 49 : if (DEBUGLEVEL) timer_start(&ti);
4050 49 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
4051 49 : if (DEBUGLEVEL) timer_printf(&ti," log element");
4052 49 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
4053 49 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
4054 49 : d = gcdii(Ao,Bo);
4055 49 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
4056 49 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
4057 49 : return gerepileuptoint(av, l);
4058 : }
4059 :
4060 : static int
4061 713406 : Fp_log_use_index(long e, long p)
4062 : {
4063 713406 : return (e >= 27 && 20*(p+6)<=e*e);
4064 : }
4065 :
4066 : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
4067 : static GEN
4068 1076092 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
4069 : {
4070 1076092 : pari_sp av = avma;
4071 1076092 : GEN p = (GEN)E;
4072 : /* assume a reduced mod p, p not necessarily prime */
4073 1076092 : if (equali1(a)) return gen_0;
4074 : /* p > 2 */
4075 662590 : if (equalii(subiu(p,1), a)) /* -1 */
4076 : {
4077 : pari_sp av2;
4078 : GEN t;
4079 218013 : ord = get_arith_Z(ord);
4080 218013 : if (mpodd(ord)) { avma = av; return cgetg(1, t_VEC); } /* no solution */
4081 217999 : t = shifti(ord,-1); /* only possible solution */
4082 217999 : av2 = avma;
4083 217999 : if (!equalii(Fp_pow(g, t, p), a)) { avma = av; return cgetg(1, t_VEC); }
4084 217803 : avma = av2; return gerepileuptoint(av, t);
4085 : }
4086 444577 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
4087 49 : return Fp_log_index(a, g, ord, p);
4088 444528 : avma = av; return NULL; /* not easy */
4089 : }
4090 :
4091 : GEN
4092 747600 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
4093 : {
4094 747600 : GEN v = get_arith_ZZM(ord);
4095 747572 : GEN F = gmael(v,2,1);
4096 747572 : long lF = lg(F)-1, lmax;
4097 747572 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
4098 631276 : lmax = expi(gel(F,lF));
4099 631276 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
4100 49 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
4101 631276 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
4102 : }
4103 :
4104 : static ulong
4105 0 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
4106 : {
4107 0 : ulong i, h=1;
4108 0 : for(i=0; i<ord; i++, h = Fl_mul(h, g, p))
4109 0 : if(a==h) return i;
4110 0 : return ~0UL;
4111 : }
4112 :
4113 : static ulong
4114 19005 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
4115 : {
4116 19005 : ulong i, h=1;
4117 47866 : for(i=0; i<ord; i++, h = Fl_mul_pre(h, g, p, pi))
4118 47866 : if(a==h) return i;
4119 0 : return ~0UL;
4120 : }
4121 :
4122 : static ulong
4123 0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
4124 : {
4125 0 : pari_sp av = avma;
4126 0 : GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
4127 0 : ulong z = typ(r)==t_INT ? itou(r): ~0L;
4128 0 : avma = av; return z;
4129 : }
4130 :
4131 : ulong
4132 19005 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
4133 : {
4134 19005 : if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
4135 0 : return Fl_log_Fp(a, g, ord, p);
4136 : }
4137 :
4138 : ulong
4139 0 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
4140 : {
4141 0 : if (ord <= 200)
4142 0 : return (p&HIGHMASK) ? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
4143 0 : : Fl_log_naive(a, g, ord, p);
4144 0 : return Fl_log_Fp(a, g, ord, p);
4145 : }
4146 :
4147 : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
4148 : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
4149 : static GEN
4150 112 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
4151 : {
4152 112 : long l = lg(P) - 1, e = E[l];
4153 112 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
4154 : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
4155 :
4156 112 : if (l == 1) {
4157 84 : hpe = h;
4158 84 : gpe = g;
4159 : } else {
4160 28 : hpe = modii(h, pe);
4161 28 : gpe = modii(g, pe);
4162 : }
4163 112 : if (e == 1) {
4164 28 : hp = hpe;
4165 28 : gp = gpe;
4166 : } else {
4167 84 : hp = remii(hpe, p);
4168 84 : gp = remii(gpe, p);
4169 : }
4170 112 : if (hp == gen_0 || gp == gen_0) return NULL;
4171 91 : if (absequaliu(p, 2))
4172 : {
4173 35 : GEN n = int2n(e);
4174 35 : ogpe = Zp_order(gpe, gen_2, e, n);
4175 35 : a = Fp_log(hpe, gpe, ogpe, n);
4176 35 : if (typ(a) != t_INT) return NULL;
4177 : }
4178 : else
4179 : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
4180 : is trivial */
4181 : /* [order(gp), factor(order(gp))] */
4182 56 : GEN v = Fp_factored_order(gp, subiu(p,1), p);
4183 56 : GEN ogp = gel(v,1);
4184 56 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
4185 56 : a = Fp_log(hp, gp, v, p);
4186 56 : if (typ(a) != t_INT) return NULL;
4187 56 : if (e == 1) ogpe = ogp;
4188 : else
4189 : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
4190 : /* use p-adic log: O(log p + e) mul*/
4191 : long vpogpe, vpohpe;
4192 :
4193 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
4194 28 : gpe = Fp_pow(gpe, ogp, pe);
4195 : /* g,h = 1 mod p; compute b s.t. h = g^b */
4196 :
4197 : /* v_p(order g mod pe) */
4198 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
4199 : /* v_p(order h mod pe) */
4200 28 : vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
4201 28 : if (vpohpe > vpogpe) return NULL;
4202 :
4203 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
4204 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
4205 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
4206 28 : a = addii(a, mulii(ogp, padic_to_Q(b)));
4207 : }
4208 : }
4209 : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
4210 77 : if (l == 1) return a;
4211 :
4212 28 : N = diviiexact(N, pe); /* make N coprime to p */
4213 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
4214 28 : g = Fp_pow(g, modii(ogpe, phi), N);
4215 28 : setlg(P, l); /* remove last element */
4216 28 : setlg(E, l);
4217 28 : b = znlog_rec(h, g, N, P, E, PHI);
4218 28 : if (!b) return NULL;
4219 28 : return addmulii(a, b, ogpe);
4220 : }
4221 :
4222 : static GEN
4223 84 : get_PHI(GEN P, GEN E)
4224 : {
4225 84 : long i, l = lg(P);
4226 84 : GEN PHI = cgetg(l, t_VEC);
4227 84 : gel(PHI,1) = gen_1;
4228 112 : for (i=1; i<l-1; i++)
4229 : {
4230 28 : GEN t, p = gel(P,i);
4231 28 : long e = E[i];
4232 28 : t = mulii(powiu(p, e-1), subiu(p,1));
4233 28 : if (i > 1) t = mulii(t, gel(PHI,i));
4234 28 : gel(PHI,i+1) = t;
4235 : }
4236 84 : return PHI;
4237 : }
4238 :
4239 : GEN
4240 224 : znlog(GEN h, GEN g, GEN o)
4241 : {
4242 224 : pari_sp av = avma;
4243 : GEN N, fa, P, E, x;
4244 224 : switch (typ(g))
4245 : {
4246 : case t_PADIC:
4247 : {
4248 28 : GEN p = gel(g,2);
4249 28 : long v = valp(g);
4250 28 : if (v < 0) pari_err_DIM("znlog");
4251 28 : if (v > 0) {
4252 0 : long k = gvaluation(h, p);
4253 0 : if (k % v) return cgetg(1,t_VEC);
4254 0 : k /= v;
4255 0 : if (!gequal(h, gpowgs(g,k))) { avma = av; return cgetg(1,t_VEC); }
4256 0 : avma = av; return stoi(k);
4257 : }
4258 28 : N = gel(g,3);
4259 28 : g = Rg_to_Fp(g, N);
4260 28 : break;
4261 : }
4262 : case t_INTMOD:
4263 189 : N = gel(g,1);
4264 189 : g = gel(g,2); break;
4265 7 : default: pari_err_TYPE("znlog", g);
4266 : return NULL; /* LCOV_EXCL_LINE */
4267 : }
4268 217 : if (equali1(N)) { avma = av; return gen_0; }
4269 217 : h = Rg_to_Fp(h, N);
4270 210 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
4271 84 : fa = Z_factor(N);
4272 84 : P = gel(fa,1);
4273 84 : E = vec_to_vecsmall(gel(fa,2));
4274 84 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
4275 84 : if (!x) { avma = av; return cgetg(1,t_VEC); }
4276 49 : return gerepileuptoint(av, x);
4277 : }
4278 :
4279 : GEN
4280 61301 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
4281 : {
4282 61301 : if (lgefint(p)==3)
4283 : {
4284 60909 : long nn = itos_or_0(n);
4285 60909 : if (nn)
4286 : {
4287 60909 : ulong pp = p[2];
4288 : ulong uz;
4289 60909 : ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
4290 60888 : if (r==ULONG_MAX) return NULL;
4291 60853 : if (zeta) *zeta = utoi(uz);
4292 60853 : return utoi(r);
4293 : }
4294 : }
4295 392 : a = modii(a,p);
4296 392 : if (!signe(a))
4297 : {
4298 0 : if (zeta) *zeta = gen_1;
4299 0 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
4300 0 : return gen_0;
4301 : }
4302 392 : if (absequaliu(n,2))
4303 : {
4304 224 : if (zeta) *zeta = subiu(p,1);
4305 224 : return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
4306 : }
4307 168 : return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
4308 : }
4309 :
4310 : /*********************************************************************/
4311 : /** **/
4312 : /** FUNDAMENTAL DISCRIMINANTS **/
4313 : /** **/
4314 : /*********************************************************************/
4315 : static int
4316 1407 : fa_isfundamental(GEN F)
4317 : {
4318 1407 : GEN P = gel(F,1), E = gel(F,2);
4319 1407 : long i, s, l = lg(P);
4320 :
4321 1407 : if (l == 1) return 1;
4322 1400 : s = signe(gel(P,1)); /* = signe(x) */
4323 1400 : if (!s) return 0;
4324 1393 : if (s < 0) { l--; P = vecslice(P,2,l); E = vecslice(E,2,l); }
4325 1393 : if (l == 1) return 0;
4326 1386 : if (!absequaliu(gel(P,1), 2))
4327 686 : i = 1; /* need x = 1 mod 4 */
4328 : else
4329 : {
4330 700 : i = 2;
4331 700 : switch(itou(gel(E,1)))
4332 : {
4333 182 : case 2: s = -s; break; /* need x/4 = 3 mod 4 */
4334 84 : case 3: s = 0; break; /* no condition mod 4 */
4335 434 : default: return 0;
4336 : }
4337 : }
4338 1974 : for(; i < l; i++)
4339 : {
4340 1190 : if (!equali1(gel(E,i))) return 0;
4341 1022 : if (s && Mod4(gel(P,i)) == 3) s = -s;
4342 : }
4343 784 : return s >= 0;
4344 : }
4345 : long
4346 17563 : isfundamental(GEN x)
4347 : {
4348 17563 : if (typ(x) != t_INT)
4349 : {
4350 1407 : pari_sp av = avma;
4351 1407 : int v = fa_isfundamental(check_arith_all(x,"isfundamental"));
4352 1407 : avma = av; return v;
4353 : }
4354 16156 : return Z_isfundamental(x);
4355 : }
4356 :
4357 : /* x fundamental ? */
4358 : long
4359 10583 : uposisfundamental(ulong x)
4360 : {
4361 10583 : ulong r = x & 15; /* x mod 16 */
4362 10583 : if (!r) return 0;
4363 10002 : switch(r & 3)
4364 : { /* x mod 4 */
4365 1947 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
4366 3109 : case 1: return uissquarefree(x);
4367 4946 : default: return 0;
4368 : }
4369 : }
4370 : /* -x fundamental ? */
4371 : long
4372 21211 : unegisfundamental(ulong x)
4373 : {
4374 21211 : ulong r = x & 15; /* x mod 16 */
4375 21211 : if (!r) return 0;
4376 20210 : switch(r & 3)
4377 : { /* x mod 4 */
4378 3676 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
4379 9698 : case 3: return uissquarefree(x);
4380 6836 : default: return 0;
4381 : }
4382 : }
4383 : long
4384 10353 : sisfundamental(long x)
4385 10353 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
4386 :
4387 : long
4388 16723 : Z_isfundamental(GEN x)
4389 : {
4390 : long r;
4391 16723 : switch(lgefint(x))
4392 : {
4393 7 : case 2: return 0;
4394 22175 : case 3: return signe(x) < 0? unegisfundamental(x[2])
4395 22175 : : uposisfundamental(x[2]);
4396 : }
4397 2010 : r = mod16(x);
4398 2010 : if (!r) return 0;
4399 1884 : if ((r & 3) == 0)
4400 : {
4401 : pari_sp av;
4402 376 : r >>= 2; /* |x|/4 mod 4 */
4403 376 : if (signe(x) < 0) r = 4-r;
4404 376 : if (r == 1) return 0;
4405 250 : av = avma;
4406 250 : r = Z_issquarefree( shifti(x,-2) );
4407 250 : avma = av; return r;
4408 : }
4409 1508 : r &= 3; /* |x| mod 4 */
4410 1508 : if (signe(x) < 0) r = 4-r;
4411 1508 : return (r==1) ? Z_issquarefree(x) : 0;
4412 : }
4413 :
4414 : static GEN
4415 2821 : fa_quaddisc(GEN f)
4416 : {
4417 2821 : GEN P = gel(f,1), E = gel(f,2), s = gen_1;
4418 2821 : long i, l = lg(P);
4419 9051 : for (i = 1; i < l; i++) /* possibly including -1 */
4420 6230 : if (mpodd(gel(E,i))) s = mulii(s, gel(P,i));
4421 2821 : if (Mod4(s) > 1) s = shifti(s,2);
4422 2821 : return s;
4423 : }
4424 :
4425 : GEN
4426 2821 : quaddisc(GEN x)
4427 : {
4428 2821 : const pari_sp av = avma;
4429 2821 : if (is_rational_t(typ(x))) x = factor(x);
4430 1407 : else x = check_arith_all(x,"quaddisc");
4431 2821 : return gerepileuptoint(av, fa_quaddisc(x));
4432 : }
4433 :
4434 : /*********************************************************************/
4435 : /** **/
4436 : /** FACTORIAL **/
4437 : /** **/
4438 : /*********************************************************************/
4439 : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
4440 : * first is slower ... ] */
4441 : GEN
4442 1208921 : mulu_interval(ulong a, ulong b)
4443 : {
4444 1208921 : pari_sp av = avma;
4445 : ulong k, l, N, n;
4446 : long lx;
4447 : GEN x;
4448 :
4449 1208921 : if (!a) return gen_0;
4450 1208921 : n = b - a + 1;
4451 1208921 : if (n < 61)
4452 : {
4453 1200794 : if (n == 1) return utoi(a);
4454 574631 : x = muluu(a,a+1); if (n == 2) return x;
4455 284064 : for (k=a+2; k<=b; k++) x = mului(k,x);
4456 284064 : return gerepileuptoint(av, x);
4457 : }
4458 8127 : lx = 1; x = cgetg(2 + n/2, t_VEC);
4459 8127 : N = b + a;
4460 980334 : for (k = a;; k++)
4461 : {
4462 1952541 : l = N - k; if (l <= k) break;
4463 972207 : gel(x,lx++) = muluu(k,l);
4464 : }
4465 8127 : if (l == k) gel(x,lx++) = utoipos(k);
4466 8127 : setlg(x, lx);
4467 8127 : return gerepileuptoint(av, ZV_prod(x));
4468 : }
4469 : GEN
4470 448 : muls_interval(long a, long b)
4471 : {
4472 448 : pari_sp av = avma;
4473 448 : long lx, k, l, N, n = b - a + 1;
4474 : GEN x;
4475 :
4476 448 : if (a <= 0 && b >= 0) return gen_0;
4477 287 : if (n < 61)
4478 : {
4479 287 : x = stoi(a);
4480 287 : for (k=a+1; k<=b; k++) x = mulsi(k,x);
4481 287 : return gerepileuptoint(av, x);
4482 : }
4483 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
4484 0 : N = b + a;
4485 0 : for (k = a;; k++)
4486 : {
4487 0 : l = N - k; if (l <= k) break;
4488 0 : gel(x,lx++) = mulss(k,l);
4489 : }
4490 0 : if (l == k) gel(x,lx++) = stoi(k);
4491 0 : setlg(x, lx);
4492 0 : return gerepileuptoint(av, ZV_prod(x));
4493 : }
4494 :
4495 : GEN
4496 2934676 : mpfact(long n)
4497 : {
4498 2934676 : if (n < 2)
4499 : {
4500 1848295 : if (n < 0) pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
4501 1848295 : return gen_1;
4502 : }
4503 1086381 : return mulu_interval(2UL, (ulong)n);
4504 : }
4505 :
4506 : /*******************************************************************/
4507 : /** **/
4508 : /** LUCAS & FIBONACCI **/
4509 : /** **/
4510 : /*******************************************************************/
4511 : static void
4512 56 : lucas(ulong n, GEN *a, GEN *b)
4513 : {
4514 : GEN z, t, zt;
4515 56 : if (!n) { *a = gen_2; *b = gen_1; return; }
4516 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
4517 49 : switch(n & 3) {
4518 14 : case 0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
4519 14 : case 1: *a = subiu(zt,1); *b = addiu(sqri(t),2); break;
4520 7 : case 2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
4521 14 : case 3: *a = addiu(zt,1); *b = subiu(sqri(t),2);
4522 : }
4523 : }
4524 :
4525 : GEN
4526 7 : fibo(long n)
4527 : {
4528 7 : pari_sp av = avma;
4529 : GEN a, b;
4530 7 : if (!n) return gen_0;
4531 7 : lucas((ulong)(labs(n)-1), &a, &b);
4532 7 : a = diviuexact(addii(shifti(a,1),b), 5);
4533 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
4534 7 : return gerepileuptoint(av, a);
4535 : }
4536 :
4537 : /*******************************************************************/
4538 : /* */
4539 : /* CONTINUED FRACTIONS */
4540 : /* */
4541 : /*******************************************************************/
4542 : static GEN
4543 488413 : icopy_lg(GEN x, long l)
4544 : {
4545 488413 : long lx = lgefint(x);
4546 : GEN y;
4547 :
4548 488413 : if (lx >= l) return icopy(x);
4549 35 : y = cgeti(l); affii(x, y); return y;
4550 : }
4551 :
4552 : /* continued fraction of a/b. If y != NULL, stop when partial quotients
4553 : * differ from y */
4554 : static GEN
4555 488709 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
4556 : {
4557 : GEN z, c;
4558 488709 : ulong i, l, ly = lgefint(b);
4559 :
4560 : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
4561 488709 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
4562 488709 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
4563 488709 : if (l > LGBITS) l = LGBITS;
4564 :
4565 488709 : z = cgetg(l,t_VEC);
4566 488709 : l--;
4567 488709 : if (y) {
4568 296 : pari_sp av = avma;
4569 296 : if (l >= (ulong)lg(y)) l = lg(y)-1;
4570 19488 : for (i = 1; i <= l; i++)
4571 : {
4572 19303 : GEN q = gel(y,i);
4573 19303 : gel(z,i) = q;
4574 19303 : c = b; if (!gequal1(q)) c = mulii(q, b);
4575 19303 : c = subii(a, c);
4576 19303 : if (signe(c) < 0)
4577 : { /* partial quotient too large */
4578 110 : c = addii(c, b);
4579 110 : if (signe(c) >= 0) i++; /* by 1 */
4580 110 : break;
4581 : }
4582 19193 : if (cmpii(c, b) >= 0)
4583 : { /* partial quotient too small */
4584 1 : c = subii(c, b);
4585 1 : if (cmpii(c, b) < 0) {
4586 : /* by 1. If next quotient is 1 in y, add 1 */
4587 0 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
4588 0 : i++;
4589 : }
4590 1 : break;
4591 : }
4592 19192 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
4593 19192 : a = b; b = c;
4594 : }
4595 : } else {
4596 488413 : a = icopy_lg(a, ly);
4597 488413 : b = icopy(b);
4598 1738365 : for (i = 1; i <= l; i++)
4599 : {
4600 1738101 : gel(z,i) = truedvmdii(a,b,&c);
4601 1738101 : if (c == gen_0) { i++; break; }
4602 1249952 : affii(c, a); cgiv(c); c = a;
4603 1249952 : a = b; b = c;
4604 : }
4605 : }
4606 488709 : i--;
4607 488709 : if (i > 1 && gequal1(gel(z,i)))
4608 : {
4609 85 : cgiv(gel(z,i)); --i;
4610 85 : gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
4611 : }
4612 488709 : setlg(z,i+1); return z;
4613 : }
4614 :
4615 : static GEN
4616 0 : sersfcont(GEN a, GEN b, long k)
4617 : {
4618 0 : long i, l = typ(a) == t_POL? lg(a): 3;
4619 : GEN y, c;
4620 0 : if (lg(b) > l) l = lg(b);
4621 0 : if (k > 0 && l > k+1) l = k+1;
4622 0 : y = cgetg(l,t_VEC);
4623 0 : for (i=1; i<l; i++)
4624 : {
4625 0 : gel(y,i) = poldivrem(a,b,&c);
4626 0 : if (gequal0(c)) { i++; break; }
4627 0 : a = b; b = c;
4628 : }
4629 0 : setlg(y, i); return y;
4630 : }
4631 :
4632 : GEN
4633 488973 : gboundcf(GEN x, long k)
4634 : {
4635 : pari_sp av;
4636 488973 : long tx = typ(x), e;
4637 : GEN y, a, b, c;
4638 :
4639 488973 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
4640 488966 : if (is_scalar_t(tx))
4641 : {
4642 488966 : if (gequal0(x)) return mkvec(gen_0);
4643 488917 : switch(tx)
4644 : {
4645 497 : case t_INT: return mkveccopy(x);
4646 : case t_REAL:
4647 303 : av = avma;
4648 303 : c = mantissa_real(x,&e);
4649 303 : if (e < 0) pari_err_PREC("gboundcf");
4650 296 : y = int2n(e);
4651 296 : a = Qsfcont(c,y, NULL, k);
4652 296 : b = addsi(signe(x), c);
4653 296 : return gerepilecopy(av, Qsfcont(b,y, a, k));
4654 :
4655 : case t_FRAC:
4656 488117 : av = avma;
4657 488117 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
4658 : }
4659 0 : pari_err_TYPE("gboundcf",x);
4660 : }
4661 :
4662 0 : switch(tx)
4663 : {
4664 0 : case t_POL: return mkveccopy(x);
4665 : case t_SER:
4666 0 : av = avma;
4667 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
4668 : case t_RFRAC:
4669 0 : av = avma;
4670 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
4671 : }
4672 0 : pari_err_TYPE("gboundcf",x);
4673 : return NULL; /* LCOV_EXCL_LINE */
4674 : }
4675 :
4676 : static GEN
4677 14 : sfcont2(GEN b, GEN x, long k)
4678 : {
4679 14 : pari_sp av = avma;
4680 14 : long lb = lg(b), tx = typ(x), i;
4681 : GEN y,p1;
4682 :
4683 14 : if (k)
4684 : {
4685 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
4686 0 : lb = k+1;
4687 : }
4688 7 : y = cgetg(lb,t_VEC);
4689 7 : if (lb==1) return y;
4690 7 : if (is_scalar_t(tx))
4691 : {
4692 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
4693 : }
4694 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
4695 :
4696 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
4697 7 : for (i = 1;;)
4698 : {
4699 63 : if (tx == t_REAL)
4700 : {
4701 35 : long e = expo(x);
4702 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
4703 35 : gel(y,i) = floorr(x);
4704 35 : p1 = subri(x, gel(y,i));
4705 : }
4706 : else
4707 : {
4708 0 : gel(y,i) = gfloor(x);
4709 0 : p1 = gsub(x, gel(y,i));
4710 : }
4711 35 : if (++i >= lb) break;
4712 28 : if (gequal0(p1)) break;
4713 28 : x = gdiv(gel(b,i),p1);
4714 : }
4715 7 : setlg(y,i);
4716 7 : return gerepilecopy(av,y);
4717 : }
4718 :
4719 :
4720 : GEN
4721 98 : gcf(GEN x) { return gboundcf(x,0); }
4722 : GEN
4723 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
4724 : GEN
4725 49 : contfrac0(GEN x, GEN b, long nmax)
4726 : {
4727 : long tb;
4728 :
4729 49 : if (!b) return gboundcf(x,nmax);
4730 28 : tb = typ(b);
4731 28 : if (tb == t_INT) return gboundcf(x,itos(b));
4732 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
4733 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
4734 14 : return sfcont2(b,x,nmax);
4735 : }
4736 :
4737 : GEN
4738 238 : contfracpnqn(GEN x, long n)
4739 : {
4740 238 : pari_sp av = avma;
4741 238 : long i, lx = lg(x);
4742 : GEN M,A,B, p0,p1, q0,q1;
4743 :
4744 238 : if (lx == 1)
4745 : {
4746 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
4747 21 : if (n >= 0) return cgetg(1,t_MAT);
4748 7 : return matid(2);
4749 : }
4750 210 : switch(typ(x))
4751 : {
4752 168 : case t_VEC: case t_COL: A = x; B = NULL; break;
4753 : case t_MAT:
4754 42 : switch(lgcols(x))
4755 : {
4756 0 : case 2: A = row(x,1); B = NULL; break;
4757 35 : case 3: A = row(x,2); B = row(x,1); break;
4758 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
4759 : return NULL; /*LCOV_EXCL_LINE*/
4760 : }
4761 35 : break;
4762 0 : default: pari_err_TYPE("pnqn",x);
4763 : return NULL; /*LCOV_EXCL_LINE*/
4764 : }
4765 203 : p1 = gel(A,1);
4766 203 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
4767 203 : if (n >= 0)
4768 : {
4769 168 : lx = minss(lx, n+2);
4770 168 : if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
4771 : }
4772 35 : else if (lx == 2)
4773 7 : return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
4774 : /* lx >= 3 */
4775 105 : p0 = gen_1;
4776 105 : q0 = gen_0; /* p[-1], q[-1] */
4777 105 : M = cgetg(lx, t_MAT);
4778 105 : gel(M,1) = mkcol2(p1,q1);
4779 350 : for (i=2; i<lx; i++)
4780 : {
4781 245 : GEN a = gel(A,i), p2,q2;
4782 245 : if (B) {
4783 84 : GEN b = gel(B,i);
4784 84 : p0 = gmul(b,p0);
4785 84 : q0 = gmul(b,q0);
4786 : }
4787 245 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
4788 245 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
4789 245 : gel(M,i) = mkcol2(p1,q1);
4790 : }
4791 105 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
4792 105 : return gerepilecopy(av, M);
4793 : }
4794 : GEN
4795 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
4796 : /* x = [a0, ..., an] from gboundcf, n >= 0;
4797 : * return [[p0, ..., pn], [q0,...,qn]] */
4798 : GEN
4799 480802 : ZV_allpnqn(GEN x)
4800 : {
4801 480802 : long i, lx = lg(x);
4802 480802 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
4803 :
4804 480802 : gel(v,1) = P = cgetg(lx, t_VEC);
4805 480802 : gel(v,2) = Q = cgetg(lx, t_VEC);
4806 480802 : p0 = gen_1; q0 = gen_0;
4807 480802 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
4808 1685516 : for (i=2; i<lx; i++)
4809 : {
4810 1204714 : GEN a = gel(x,i);
4811 1204714 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
4812 1204714 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
4813 : }
4814 480802 : return v;
4815 : }
4816 :
4817 : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
4818 : static GEN
4819 42 : mod_to_frac(GEN x, GEN N, GEN B)
4820 : {
4821 : GEN a, b, A;
4822 42 : if (B) A = divii(shifti(N, -1), B);
4823 : else
4824 : {
4825 14 : A = sqrti(shifti(N, -1));
4826 14 : B = A;
4827 : }
4828 42 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
4829 28 : return equali1(b)? a: mkfrac(a,b);
4830 : }
4831 :
4832 : static GEN
4833 70 : mod_to_rfrac(GEN x, GEN N, long B)
4834 : {
4835 : GEN a, b;
4836 70 : long A, d = degpol(N);
4837 70 : if (B >= 0) A = d-1 - B;
4838 : else
4839 : {
4840 42 : B = d >> 1;
4841 42 : A = odd(d)? B : B-1;
4842 : }
4843 70 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
4844 70 : if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
4845 56 : return gdiv(a,b);
4846 : }
4847 :
4848 : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
4849 : * of the continued fraction of x with b <= k maximal */
4850 : static GEN
4851 0 : bestappr_frac(GEN x, GEN k)
4852 : {
4853 : pari_sp av;
4854 : GEN p0, p1, p, q0, q1, q, a, y;
4855 :
4856 0 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
4857 0 : av = avma; y = x;
4858 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
4859 0 : q1 = gen_0; q0 = gen_1;
4860 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
4861 : for(;;)
4862 : {
4863 0 : x = ginv(x); /* > 1 */
4864 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
4865 0 : if (cmpii(a,k) > 0)
4866 : { /* next partial quotient will overflow limits */
4867 : GEN n, d;
4868 0 : a = divii(subii(k, q1), q0);
4869 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4870 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4871 : /* compare |y-p0/q0|, |y-p1/q1| */
4872 0 : n = gel(y,1);
4873 0 : d = gel(y,2);
4874 0 : if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
4875 : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
4876 0 : { p1 = p0; q1 = q0; }
4877 0 : break;
4878 : }
4879 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4880 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4881 :
4882 0 : if (cmpii(q0,k) > 0) break;
4883 0 : x = gsub(x,a); /* 0 <= x < 1 */
4884 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
4885 :
4886 : }
4887 0 : return gerepileupto(av, gdiv(p1,q1));
4888 : }
4889 : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
4890 : * of the continued fraction of x with b <= k maximal */
4891 : static GEN
4892 277862 : bestappr_real(GEN x, GEN k)
4893 : {
4894 277862 : pari_sp av = avma;
4895 277862 : GEN kr, p0, p1, p, q0, q1, q, a, y = x;
4896 :
4897 277862 : p1 = gen_1; a = p0 = floorr(x);
4898 277862 : q1 = gen_0; q0 = gen_1;
4899 277862 : x = subri(x,a); /* 0 <= x < 1 */
4900 277862 : if (!signe(x)) { cgiv(x); return a; }
4901 270724 : kr = itor(k, realprec(x));
4902 : for(;;)
4903 468363 : {
4904 : long d;
4905 739087 : x = invr(x); /* > 1 */
4906 739087 : if (cmprr(x,kr) > 0)
4907 : { /* next partial quotient will overflow limits */
4908 265909 : a = divii(subii(k, q1), q0);
4909 265909 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4910 265909 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4911 : /* compare |y-p0/q0|, |y-p1/q1| */
4912 265909 : if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
4913 : mulir(q0, subri(mulir(q1,y), p1))) < 0)
4914 5055 : { p1 = p0; q1 = q0; }
4915 265909 : break;
4916 : }
4917 473178 : d = nbits2prec(expo(x) + 1);
4918 473178 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4919 :
4920 472993 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4921 472993 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4922 472993 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4923 :
4924 472993 : if (cmpii(q0,k) > 0) break;
4925 468820 : x = subri(x,a); /* 0 <= x < 1 */
4926 468820 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
4927 : }
4928 270724 : if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
4929 270724 : return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
4930 : }
4931 :
4932 : /* k t_INT or NULL */
4933 : static GEN
4934 426445 : bestappr_Q(GEN x, GEN k)
4935 : {
4936 426445 : long lx, tx = typ(x), i;
4937 : GEN a, y;
4938 :
4939 426445 : switch(tx)
4940 : {
4941 77 : case t_INT: return icopy(x);
4942 0 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
4943 : case t_REAL:
4944 317529 : if (!signe(x)) return gen_0;
4945 : /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
4946 277863 : i = bit_prec(x); if (i <= expo(x)) return NULL;
4947 277862 : return bestappr_real(x, k? k: int2n(i));
4948 :
4949 : case t_INTMOD: {
4950 28 : pari_sp av = avma;
4951 28 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
4952 21 : return gerepilecopy(av, a);
4953 : }
4954 : case t_PADIC: {
4955 14 : pari_sp av = avma;
4956 14 : long v = valp(x);
4957 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
4958 7 : if (v) a = gmul(a, powis(gel(x,2), v));
4959 7 : return gerepilecopy(av, a);
4960 : }
4961 :
4962 : case t_COMPLEX: {
4963 126 : pari_sp av = avma;
4964 126 : y = cgetg(3, t_COMPLEX);
4965 126 : gel(y,2) = bestappr(gel(x,2), k);
4966 126 : gel(y,1) = bestappr(gel(x,1), k);
4967 126 : if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
4968 0 : return y;
4969 : }
4970 : case t_SER:
4971 0 : if (ser_isexactzero(x)) return gcopy(x);
4972 : /* fall through */
4973 : case t_POLMOD: case t_POL: case t_RFRAC:
4974 : case t_VEC: case t_COL: case t_MAT:
4975 108671 : y = cgetg_copy(x, &lx);
4976 108671 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
4977 523941 : for (; i<lx; i++)
4978 : {
4979 415272 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
4980 415270 : gel(y,i) = a;
4981 : }
4982 108669 : if (tx == t_POL) return normalizepol(y);
4983 108655 : if (tx == t_SER) return normalize(y);
4984 108655 : return y;
4985 : }
4986 0 : pari_err_TYPE("bestappr_Q",x);
4987 : return NULL; /* LCOV_EXCL_LINE */
4988 : }
4989 :
4990 : static GEN
4991 56 : bestappr_ser(GEN x, long B)
4992 : {
4993 56 : long dN, v = valp(x), lx = lg(x);
4994 : GEN t;
4995 56 : x = normalizepol(ser2pol_i(x, lx));
4996 56 : dN = lx-2;
4997 56 : if (v > 0)
4998 : {
4999 14 : x = RgX_shift_shallow(x, v);
5000 14 : dN += v;
5001 : }
5002 42 : else if (v < 0)
5003 : {
5004 7 : if (B >= 0) B = maxss(B+v, 0);
5005 : }
5006 56 : t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
5007 56 : if (!t) return NULL;
5008 42 : if (v < 0)
5009 : {
5010 : GEN a, b;
5011 : long vx;
5012 7 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
5013 : /* t_RFRAC */
5014 7 : vx = varn(x);
5015 7 : a = gel(t,1);
5016 7 : b = gel(t,2);
5017 7 : v -= RgX_valrem(b, &b);
5018 7 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
5019 7 : if (v < 0) b = RgX_shift(b, -v);
5020 0 : else if (v > 0) {
5021 0 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
5022 0 : a = RgX_shift(a, v);
5023 : }
5024 7 : t = mkrfraccopy(a, b);
5025 : }
5026 42 : return t;
5027 : }
5028 : static GEN bestappr_RgX(GEN x, long B);
5029 : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
5030 : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
5031 : static GEN
5032 77 : bestappr_RgX(GEN x, long B)
5033 : {
5034 77 : long i, lx, tx = typ(x);
5035 : GEN y, t;
5036 77 : switch(tx)
5037 : {
5038 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
5039 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
5040 0 : return gcopy(x);
5041 :
5042 : case t_RFRAC: {
5043 14 : pari_sp av = avma;
5044 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
5045 7 : x = rfrac_to_ser(x, 2*B+1);
5046 7 : t = bestappr_ser(x, B); if (!t) return NULL;
5047 0 : return gerepileupto(av, t);
5048 : }
5049 : case t_POLMOD: {
5050 14 : pari_sp av = avma;
5051 14 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
5052 14 : return gerepileupto(av, t);
5053 : }
5054 : case t_SER: {
5055 49 : pari_sp av = avma;
5056 49 : t = bestappr_ser(x, B); if (!t) return NULL;
5057 42 : return gerepileupto(av, t);
5058 : }
5059 :
5060 : case t_VEC: case t_COL: case t_MAT:
5061 0 : y = cgetg_copy(x, &lx);
5062 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
5063 0 : for (; i<lx; i++)
5064 : {
5065 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
5066 0 : gel(y,i) = t;
5067 : }
5068 0 : return y;
5069 : }
5070 0 : pari_err_TYPE("bestappr_RgX",x);
5071 : return NULL; /* LCOV_EXCL_LINE */
5072 : }
5073 :
5074 : /* allow k = NULL: maximal accuracy */
5075 : GEN
5076 11173 : bestappr(GEN x, GEN k)
5077 : {
5078 11173 : pari_sp av = avma;
5079 11173 : if (k) { /* replace by floor(k) */
5080 10900 : switch(typ(k))
5081 : {
5082 : case t_INT:
5083 1225 : break;
5084 : case t_REAL: case t_FRAC:
5085 9675 : k = floor_safe(k); /* left on stack for efficiency */
5086 9675 : if (!signe(k)) k = gen_1;
5087 9675 : break;
5088 : default:
5089 0 : pari_err_TYPE("bestappr [bound type]", k);
5090 0 : break;
5091 : }
5092 : }
5093 11173 : x = bestappr_Q(x, k);
5094 11173 : if (!x) { avma = av; return cgetg(1,t_VEC); }
5095 11158 : return x;
5096 : }
5097 : GEN
5098 77 : bestapprPade(GEN x, long B)
5099 : {
5100 77 : pari_sp av = avma;
5101 77 : GEN t = bestappr_RgX(x, B);
5102 77 : if (!t) { avma = av; return cgetg(1,t_VEC); }
5103 63 : return t;
5104 : }
5105 :
5106 : /***********************************************************************/
5107 : /** **/
5108 : /** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
5109 : /** **/
5110 : /***********************************************************************/
5111 :
5112 : static GEN
5113 14 : get_quad(GEN f, GEN pol, long r)
5114 : {
5115 14 : GEN p1 = gcoeff(f,1,2), q1 = gcoeff(f,2,2);
5116 14 : return mkquad(pol, r? subii(p1,q1): p1, q1);
5117 : }
5118 :
5119 : /* replace f by f * [a,1; 1,0] */
5120 : static void
5121 14 : update_f(GEN f, GEN a)
5122 : {
5123 : GEN p1;
5124 14 : p1 = gcoeff(f,1,1);
5125 14 : gcoeff(f,1,1) = addii(mulii(a,p1), gcoeff(f,1,2));
5126 14 : gcoeff(f,1,2) = p1;
5127 :
5128 14 : p1 = gcoeff(f,2,1);
5129 14 : gcoeff(f,2,1) = addii(mulii(a,p1), gcoeff(f,2,2));
5130 14 : gcoeff(f,2,2) = p1;
5131 14 : }
5132 :
5133 : GEN
5134 7 : quadunit(GEN x)
5135 : {
5136 7 : pari_sp av = avma, av2;
5137 : GEN pol, y, a, u, v, sqd, f;
5138 : long r;
5139 :
5140 7 : check_quaddisc_real(x, &r, "quadunit");
5141 7 : pol = quadpoly(x);
5142 7 : sqd = sqrti(x); av2 = avma;
5143 7 : a = shifti(addui(r,sqd),-1);
5144 7 : f = mkmat2(mkcol2(a, gen_1), mkcol2(gen_1, gen_0)); /* [a,0; 1,0] */
5145 7 : u = stoi(r); v = gen_2;
5146 : for(;;)
5147 7 : {
5148 : GEN u1, v1;
5149 14 : u1 = subii(mulii(a,v),u);
5150 14 : v1 = divii(subii(x,sqri(u1)),v);
5151 14 : if ( equalii(v,v1) ) {
5152 7 : y = get_quad(f,pol,r);
5153 7 : update_f(f,a);
5154 7 : y = gdiv(get_quad(f,pol,r), conj_i(y));
5155 7 : break;
5156 : }
5157 7 : a = divii(addii(sqd,u1), v1);
5158 7 : if ( equalii(u,u1) ) {
5159 0 : y = get_quad(f,pol,r);
5160 0 : y = gdiv(y, conj_i(y));
5161 0 : break;
5162 : }
5163 7 : update_f(f,a);
5164 7 : u = u1; v = v1;
5165 7 : if (gc_needed(av2,2))
5166 : {
5167 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit");
5168 0 : gerepileall(av2,4, &a,&f,&u,&v);
5169 : }
5170 : }
5171 7 : if (signe(gel(y,3)) < 0) y = gneg(y);
5172 7 : return gerepileupto(av, y);
5173 : }
5174 :
5175 : GEN
5176 7 : quadunit0(GEN x, long v)
5177 : {
5178 7 : GEN y = quadunit(x);
5179 7 : if (v==-1) v = fetch_user_var("w");
5180 7 : setvarn(gel(y,1), v);
5181 7 : return y;
5182 : }
5183 :
5184 : GEN
5185 21 : quadregulator(GEN x, long prec)
5186 : {
5187 21 : pari_sp av = avma, av2;
5188 : GEN R, rsqd, u, v, sqd;
5189 : long r, Rexpo;
5190 :
5191 21 : check_quaddisc_real(x, &r, "quadregulator");
5192 21 : sqd = sqrti(x);
5193 21 : rsqd = gsqrt(x,prec);
5194 21 : Rexpo = 0; R = real2n(1, prec); /* = 2 */
5195 21 : av2 = avma;
5196 21 : u = stoi(r); v = gen_2;
5197 : for(;;)
5198 49 : {
5199 70 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
5200 70 : GEN v1 = divii(subii(x,sqri(u1)),v);
5201 70 : if (equalii(v,v1))
5202 : {
5203 7 : R = sqrr(R); shiftr_inplace(R, -1);
5204 7 : R = mulrr(R, divri(addir(u1,rsqd),v));
5205 7 : break;
5206 : }
5207 63 : if (equalii(u,u1))
5208 : {
5209 14 : R = sqrr(R); shiftr_inplace(R, -1);
5210 14 : break;
5211 : }
5212 49 : R = mulrr(R, divri(addir(u1,rsqd),v));
5213 49 : Rexpo += expo(R); setexpo(R,0);
5214 49 : u = u1; v = v1;
5215 49 : if (Rexpo & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
5216 49 : if (gc_needed(av2,2))
5217 : {
5218 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
5219 0 : gerepileall(av2,3, &R,&u,&v);
5220 : }
5221 : }
5222 21 : R = logr_abs(divri(R,v));
5223 21 : if (Rexpo)
5224 : {
5225 21 : GEN t = mulsr(Rexpo, mplog2(prec));
5226 21 : shiftr_inplace(t, 1);
5227 21 : R = addrr(R,t);
5228 : }
5229 21 : return gerepileuptoleaf(av, R);
5230 : }
5231 :
5232 : /*************************************************************************/
5233 : /** **/
5234 : /** CLASS NUMBER **/
5235 : /** **/
5236 : /*************************************************************************/
5237 :
5238 : int
5239 12943797 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
5240 :
5241 18316834 : static GEN qfi_pow(void *E, GEN f, GEN n)
5242 18316834 : { return E? nupow(f,n,(GEN)E): powgi(f,n); }
5243 23025703 : static GEN qfi_comp(void *E, GEN f, GEN g)
5244 23025703 : { return E? nucomp(f,g,(GEN)E): qficomp(f,g); }
5245 : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
5246 : gidentical,qfb_equal1,NULL};
5247 :
5248 : GEN
5249 2941241 : qfi_order(GEN q, GEN o)
5250 2941241 : { return gen_order(q, o, NULL, &qfi_group); }
5251 :
5252 : GEN
5253 0 : qfi_log(GEN a, GEN g, GEN o)
5254 0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
5255 :
5256 : GEN
5257 626556 : qfi_Shanks(GEN a, GEN g, long n)
5258 : {
5259 626556 : pari_sp av = avma;
5260 : GEN T, X;
5261 : long rt_n, c;
5262 :
5263 626556 : a = redimag(a);
5264 626556 : g = redimag(g);
5265 :
5266 626556 : rt_n = sqrt((double)n);
5267 626556 : c = n / rt_n;
5268 626556 : c = (c * rt_n < n + 1) ? c + 1 : c;
5269 :
5270 626556 : T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
5271 626556 : X = gen_Shanks(T, a, c, NULL, &qfi_group);
5272 :
5273 626556 : if (!X) { avma = av; return X; }
5274 332549 : return gerepileuptoint(av, X);
5275 : }
5276 :
5277 : GEN
5278 140 : qfbclassno0(GEN x,long flag)
5279 : {
5280 140 : switch(flag)
5281 : {
5282 126 : case 0: return map_proto_G(classno,x);
5283 14 : case 1: return map_proto_G(classno2,x);
5284 0 : default: pari_err_FLAG("qfbclassno");
5285 : }
5286 : return NULL; /* LCOV_EXCL_LINE */
5287 : }
5288 :
5289 : /* f^h = 1, return order(f). Set *pfao to its factorization */
5290 : static GEN
5291 2776995 : find_order(void *E, GEN f, GEN h, GEN *pfao)
5292 : {
5293 2776995 : GEN v = gen_factored_order(f, h, E, &qfi_group);
5294 2776995 : *pfao = gel(v,2); return gel(v,1);
5295 : }
5296 :
5297 : static int
5298 6705 : ok_q(GEN q, GEN h, GEN d2, long r2)
5299 : {
5300 6705 : if (d2)
5301 : {
5302 7 : if (r2 <= 2 && !mpodd(q)) return 0;
5303 7 : return is_pm1(Z_ppo(q,d2));
5304 : }
5305 : else
5306 : {
5307 6698 : if (r2 <= 1 && !mpodd(q)) return 0;
5308 6698 : return is_pm1(Z_ppo(q,h));
5309 : }
5310 : }
5311 :
5312 : /* a,b given by their factorizations. Return factorization of lcm(a,b).
5313 : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
5314 : static GEN
5315 362162 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
5316 : {
5317 362162 : GEN P = ZC_union_shallow(gel(Fa,1), gel(Fb,1));
5318 362162 : GEN A = gen_1, B = gen_1;
5319 362162 : long i, l = lg(P);
5320 362162 : GEN E = cgetg(l, t_COL);
5321 1071040 : for (i=1; i<l; i++)
5322 : {
5323 708878 : GEN p = gel(P,i);
5324 708878 : long va = Z_pval(a,p);
5325 708878 : long vb = Z_pval(b,p);
5326 708878 : if (va < vb)
5327 : {
5328 364289 : B = mulii(B,powiu(p,vb));
5329 364289 : gel(E,i) = utoi(vb);
5330 : }
5331 : else
5332 : {
5333 344589 : A = mulii(A,powiu(p,va));
5334 344589 : gel(E,i) = utoi(va);
5335 : }
5336 : }
5337 362162 : *pA = A;
5338 362162 : *pB = B; return mkmat2(P,E);
5339 : }
5340 :
5341 : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
5342 : static void
5343 362162 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
5344 : {
5345 362162 : GEN A,B, g1 = *pg1, d1 = *pd1;
5346 362162 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
5347 362162 : *pg1 = gmul(powgi(g1, diviiexact(d1,A)), powgi(f, diviiexact(o,B)));
5348 362162 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
5349 362162 : }
5350 :
5351 : /* Write x = Df^2, where D = fundamental discriminant,
5352 : * P^E = factorisation of conductor f, with E[i] >= 0 */
5353 : static void
5354 2060327 : corediscfact(GEN x, long xmod4, GEN *ptD, GEN *ptP, GEN *ptE)
5355 : {
5356 2060327 : long s = signe(x), l, i;
5357 2060327 : GEN fa = absZ_factor(x);
5358 2060327 : GEN d, P = gel(fa,1), E = gtovecsmall(gel(fa,2));
5359 :
5360 2060327 : l = lg(P); d = gen_1;
5361 5362951 : for (i=1; i<l; i++)
5362 : {
5363 3302624 : if (E[i] & 1) d = mulii(d, gel(P,i));
5364 3302624 : E[i] >>= 1;
5365 : }
5366 2060327 : if (!xmod4 && mod4(d) != ((s < 0)? 3: 1)) { d = shifti(d,2); E[1]--; }
5367 2060327 : *ptD = (s < 0)? negi(d): d;
5368 2060327 : *ptP = P;
5369 2060327 : *ptE = E;
5370 2060327 : }
5371 :
5372 : static GEN
5373 2052517 : conductor_part(GEN x, long xmod4, GEN *ptD, GEN *ptreg)
5374 : {
5375 2052517 : long l, i, s = signe(x);
5376 : GEN E, H, D, P, reg;
5377 :
5378 2052517 : corediscfact(x, xmod4, &D, &P, &E);
5379 2052517 : H = gen_1; l = lg(P);
5380 : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
5381 5329269 : for (i=1; i<l; i++)
5382 : {
5383 3276752 : long e = E[i];
5384 3276752 : if (e)
5385 : {
5386 7 : GEN p = gel(P,i);
5387 7 : H = mulii(H, subis(p, kronecker(D,p)));
5388 7 : if (e >= 2) H = mulii(H, powiu(p,e-1));
5389 : }
5390 : }
5391 :
5392 : /* divide by [ O_K^* : O^* ] */
5393 2052517 : if (s < 0)
5394 : {
5395 2052503 : reg = NULL;
5396 2052503 : switch(itou_or_0(D))
5397 : {
5398 0 : case 4: H = shifti(H,-1); break;
5399 0 : case 3: H = divis(H,3); break;
5400 : }
5401 : } else {
5402 14 : reg = quadregulator(D,DEFAULTPREC);
5403 14 : if (!equalii(x,D))
5404 0 : H = divii(H, roundr(divrr(quadregulator(x,DEFAULTPREC), reg)));
5405 : }
5406 2052517 : if (ptreg) *ptreg = reg;
5407 2052517 : *ptD = D; return H;
5408 : }
5409 :
5410 : static long
5411 2052496 : two_rank(GEN x)
5412 : {
5413 2052496 : GEN p = gel(absZ_factor(x),1);
5414 2052496 : long l = lg(p)-1;
5415 : #if 0 /* positive disc not needed */
5416 : if (signe(x) > 0)
5417 : {
5418 : long i;
5419 : for (i=1; i<=l; i++)
5420 : if (mod4(gel(p,i)) == 3) { l--; break; }
5421 : }
5422 : #endif
5423 2052496 : return l-1;
5424 : }
5425 :
5426 : static GEN
5427 38995259 : sqr_primeform(GEN x, ulong p) { return redimag(qfisqr(primeform_u(x, p))); }
5428 : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
5429 : static GEN
5430 2052496 : get_forms(GEN D, GEN *pL)
5431 : {
5432 2052496 : const long MAXFORM = 20;
5433 2052496 : GEN L, sqrtD = gsqrt(absi_shallow(D),DEFAULTPREC);
5434 2052496 : GEN forms = vectrunc_init(MAXFORM+1);
5435 2052496 : long s, nforms = 0;
5436 : ulong p;
5437 : forprime_t S;
5438 2052496 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
5439 2052496 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
5440 2052496 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
5441 2052496 : if (s < 10) s = 200;
5442 1901369 : else if (s < 20) s = 1000;
5443 1477 : else if (s < 5000) s = 5000;
5444 2052496 : u_forprime_init(&S, 2, s);
5445 344232110 : while ( (p = u_forprime_next(&S)) )
5446 : {
5447 340127118 : long d, k = kroiu(D,p);
5448 : pari_sp av2;
5449 340127118 : if (!k) continue;
5450 337897177 : if (k > 0)
5451 : {
5452 169485505 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
5453 169485505 : d = p-1;
5454 : }
5455 : else
5456 168411672 : d = p+1;
5457 337897177 : av2 = avma; affrr(divru(mulur(p,L),d), L); avma = av2;
5458 : }
5459 2052496 : *pL = L; return forms;
5460 : }
5461 :
5462 : /* h ~ #G, return o = order of f, set fao = its factorization */
5463 : static GEN
5464 2052545 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
5465 : {
5466 2052545 : long s = minss(itos(sqrti(h)), 10000);
5467 2052545 : GEN T = gen_Shanks_init(f, s, E, &qfi_group);
5468 2052545 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
5469 2052545 : return find_order(E, f, addiu(v,1), pfao);
5470 : }
5471 :
5472 : /* if g = 1 in G/<f> ? */
5473 : static int
5474 518 : equal1(void *E, GEN T, ulong N, GEN g)
5475 518 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
5476 :
5477 : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
5478 : * FIXME: should be gen_order, but equal1 has the wrong prototype */
5479 : static GEN
5480 112 : relative_order(void *E, GEN a, GEN o, ulong N, GEN T)
5481 : {
5482 112 : pari_sp av = avma;
5483 : long i, l;
5484 : GEN m;
5485 :
5486 112 : m = get_arith_ZZM(o);
5487 112 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
5488 112 : o = gel(m,1);
5489 112 : m = gel(m,2); l = lgcols(m);
5490 322 : for (i = l-1; i; i--)
5491 : {
5492 210 : GEN t, y, p = gcoeff(m,i,1);
5493 210 : long j, e = itos(gcoeff(m,i,2));
5494 210 : if (l == 2) {
5495 35 : t = gen_1;
5496 35 : y = a;
5497 : } else {
5498 175 : t = diviiexact(o, powiu(p,e));
5499 175 : y = powgi(a, t);
5500 : }
5501 210 : if (equal1(E, T,N,y)) o = t;
5502 : else {
5503 126 : for (j = 1; j < e; j++)
5504 : {
5505 28 : y = powgi(y, p);
5506 28 : if (equal1(E, T,N,y)) break;
5507 : }
5508 119 : if (j < e) {
5509 21 : if (j > 1) p = powiu(p, j);
5510 21 : o = mulii(t, p);
5511 : }
5512 : }
5513 : }
5514 112 : return gerepilecopy(av, o);
5515 : }
5516 :
5517 : /* h(x) for x<0 using Baby Step/Giant Step.
5518 : * Assumes G is not too far from being cyclic.
5519 : *
5520 : * Compute G^2 instead of G so as to kill most of the non-cyclicity */
5521 : GEN
5522 2054893 : classno(GEN x)
5523 : {
5524 2054893 : pari_sp av = avma;
5525 : long r2, k, s, i, l;
5526 : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
5527 : void *E;
5528 :
5529 2054893 : if (signe(x) >= 0) return classno2(x);
5530 :
5531 2054886 : check_quaddisc(x, &s, &k, "classno");
5532 2054886 : if (abscmpiu(x,12) <= 0) return gen_1;
5533 :
5534 2052496 : Hf = conductor_part(x, k, &D, NULL);
5535 2052496 : if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
5536 2052496 : forms = get_forms(D, &L);
5537 2052496 : r2 = two_rank(D);
5538 2052496 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
5539 :
5540 2052496 : l = lg(forms);
5541 2052496 : order_bound = const_vec(l-1, NULL);
5542 2052496 : E = expi(D) > 60? (void*)sqrtnint(shifti(absi_shallow(D),-2),4): NULL;
5543 2052496 : g1 = gel(forms,1);
5544 2052496 : gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
5545 2052496 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
5546 2052496 : d2 = NULL; /* not computed yet */
5547 2052496 : if (is_pm1(q)) goto END;
5548 509036 : for (i=2; i < l; i++)
5549 : {
5550 502268 : GEN o, fao, a, F, fd, f = gel(forms,i);
5551 502268 : fd = powgi(f, d1); if (is_pm1(gel(fd,1))) continue;
5552 362162 : F = powgi(fd, q);
5553 362162 : a = gel(F,1);
5554 362162 : o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
5555 : /* f^(d1 q) = 1 */
5556 362162 : fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
5557 362162 : o = find_order(E, f, fao, &fao);
5558 362162 : gel(order_bound,i) = o;
5559 : /* o = order of f, fao = factor(o) */
5560 362162 : update_g1(&g1,&d1,&fad1, f,o,fao);
5561 362162 : q = diviiround(hin, d1);
5562 362162 : if (is_pm1(q)) goto END;
5563 : }
5564 : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
5565 6768 : if (expi(q) > 3)
5566 : { /* q large: compute d2, 2nd elt divisor */
5567 70 : ulong N, n = 2*itou(sqrti(d1));
5568 70 : GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
5569 70 : d2 = gen_1;
5570 70 : N = itou( gceil(gdivgs(d1,n)) ); /* order(g1) <= n*N */
5571 287 : for (i = 1; i < l; i++)
5572 : {
5573 280 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
5574 280 : if (!B) B = find_order(E, f, fad1, /*junk*/&d);
5575 280 : f = powgi(f,d2);
5576 280 : if (equal1(E, T,N,f)) continue;
5577 112 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
5578 : /* f^B = 1 */
5579 112 : d = relative_order(E, f, B, N,T);
5580 112 : d2= mulii(d,d2);
5581 112 : D = mulii(d1,d2);
5582 112 : q = diviiround(hin,D);
5583 112 : if (is_pm1(q)) { d1 = D; goto END; }
5584 : }
5585 : /* very probably, d2 is the 2nd elementary divisor */
5586 7 : d1 = D; /* product of first two elt divisors */
5587 : }
5588 : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
5589 : * 2-rank */
5590 6705 : if (!ok_q(q,d1,d2,r2))
5591 : {
5592 0 : GEN q0 = q;
5593 : long d;
5594 0 : if (cmpii(mulii(q,d1), hin) < 0)
5595 : { /* try q = q0+1,-1,+2,-2 */
5596 0 : d = 1;
5597 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
5598 : }
5599 : else
5600 : { /* q0-1,+1,-2,+2 */
5601 0 : d = -1;
5602 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
5603 : }
5604 : }
5605 6705 : d1 = mulii(d1,q);
5606 :
5607 : END:
5608 2052496 : return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
5609 : }
5610 :
5611 : GEN
5612 0 : quadclassno(GEN x)
5613 : {
5614 0 : pari_sp av = avma;
5615 : GEN Hf, D;
5616 : long s, r;
5617 0 : check_quaddisc(x, &s, &r, "quadclassno");
5618 0 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
5619 0 : Hf = conductor_part(x, r, &D, NULL);
5620 0 : return gerepileuptoint(av, mulii(Hf, gel(quadclassunit0(D,0,NULL,0),1)));
5621 : }
5622 :
5623 : /* use Euler products */
5624 : GEN
5625 21 : classno2(GEN x)
5626 : {
5627 21 : pari_sp av = avma;
5628 21 : const long prec = DEFAULTPREC;
5629 : long n, i, r, s;
5630 : GEN p1, p2, S, p4, p5, p7, Hf, Pi, reg, logd, d, dr, D, half;
5631 :
5632 21 : check_quaddisc(x, &s, &r, "classno2");
5633 21 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
5634 :
5635 21 : Hf = conductor_part(x, r, &D, ®);
5636 21 : if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
5637 :
5638 21 : Pi = mppi(prec);
5639 21 : d = absi_shallow(D); dr = itor(d, prec);
5640 21 : logd = logr_abs(dr);
5641 21 : p1 = sqrtr(divrr(mulir(d,logd), gmul2n(Pi,1)));
5642 21 : if (s > 0)
5643 : {
5644 14 : GEN invlogd = invr(logd);
5645 14 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
5646 14 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
5647 : }
5648 21 : n = itos_or_0( mptrunc(p1) );
5649 21 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
5650 :
5651 21 : p4 = divri(Pi,d);
5652 21 : p7 = invr(sqrtr_abs(Pi));
5653 21 : half = real2n(-1, prec);
5654 21 : if (s > 0)
5655 : { /* i = 1, shortcut */
5656 14 : p1 = sqrtr_abs(dr);
5657 14 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5658 14 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
5659 546 : for (i=2; i<=n; i++)
5660 : {
5661 532 : long k = kroiu(D,i); if (!k) continue;
5662 434 : p2 = mulir(sqru(i), p4);
5663 434 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5664 434 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
5665 434 : S = (k>0)? addrr(S,p5): subrr(S,p5);
5666 : }
5667 14 : S = shiftr(divrr(S,reg),-1);
5668 : }
5669 : else
5670 : { /* i = 1, shortcut */
5671 7 : p1 = gdiv(sqrtr_abs(dr), Pi);
5672 7 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5673 7 : S = addrr(p5, divrr(p1, mpexp(p4)));
5674 952 : for (i=2; i<=n; i++)
5675 : {
5676 945 : long k = kroiu(D,i); if (!k) continue;
5677 945 : p2 = mulir(sqru(i), p4);
5678 945 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5679 945 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
5680 945 : S = (k>0)? addrr(S,p5): subrr(S,p5);
5681 : }
5682 : }
5683 21 : return gerepileuptoint(av, mulii(Hf, roundr(S)));
5684 : }
5685 :
5686 : /* 1 + q + ... + q^v, v > 0 */
5687 : static GEN
5688 120 : geomsumu(ulong q, long v)
5689 : {
5690 120 : GEN u = utoipos(1+q);
5691 120 : for (; v > 1; v--) u = addui(1, mului(q, u));
5692 120 : return u;
5693 : }
5694 : static GEN
5695 120 : geomsum(GEN q, long v)
5696 : {
5697 : GEN u;
5698 120 : if (lgefint(q) == 3) return geomsumu(q[2], v);
5699 0 : u = addiu(q,1);
5700 0 : for (; v > 1; v--) u = addui(1, mulii(q, u));
5701 0 : return u;
5702 : }
5703 :
5704 : static GEN
5705 7810 : hclassno6_large(GEN x)
5706 : {
5707 : long i, l, s, xmod4;
5708 : GEN Q, H, D, P, E;
5709 :
5710 7810 : x = negi(x);
5711 7810 : check_quaddisc(x, &s, &xmod4, "hclassno");
5712 7810 : corediscfact(x, xmod4, &D, &P, &E);
5713 :
5714 7810 : Q = quadclassunit0(D, 0, NULL, 0);
5715 7810 : H = gel(Q,1); l = lg(P);
5716 :
5717 : /* H \prod_{p^e||f} (1 + (p^e-1)/(p-1))[ p - (D/p) ] */
5718 33682 : for (i=1; i<l; i++)
5719 : {
5720 25872 : long e = E[i], s;
5721 : GEN p, t;
5722 25872 : if (!e) continue;
5723 5003 : p = gel(P,i); s = kronecker(D,p);
5724 5003 : if (e == 1) t = addiu(p, 1-s);
5725 1000 : else if (s == 1) t = powiu(p,e);
5726 120 : else t = addui(1, mulii(subis(p, s), geomsum(p,e-1)));
5727 5003 : H = mulii(H,t);
5728 : }
5729 7810 : switch( itou_or_0(D) )
5730 : {
5731 0 : case 3: H = shifti(H,1);break;
5732 0 : case 4: H = muliu(H,3); break;
5733 7810 : default:H = muliu(H,6); break;
5734 : }
5735 7810 : return H;
5736 : }
5737 :
5738 : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
5739 : GEN
5740 121870 : hclassno6(GEN x)
5741 : {
5742 121870 : ulong d = itou_or_0(x);
5743 121870 : if (!d || d > 500000) return hclassno6_large(x);
5744 114060 : return utoipos(hclassno6u(d));
5745 : }
5746 :
5747 : GEN
5748 46088 : hclassno(GEN x)
5749 : {
5750 : long a, s;
5751 46088 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
5752 46088 : s = signe(x);
5753 46088 : if (s < 0) return gen_0;
5754 46088 : if (!s) return gdivgs(gen_1, -12);
5755 46088 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
5756 46088 : return gdivgs(hclassno6(x), 6);
5757 : }
5758 : /******************************************************************/
5759 : /* */
5760 : /* RAMANUJAN's TAU FUNCTION */
5761 : /* */
5762 : /******************************************************************/
5763 : /* 4|N > 0, not fundamental at 2; 6 * Hurwitz class number in level 2,
5764 : * equal to 6*(H(N)+2H(N/4)), H=qfbhclassno */
5765 : static GEN
5766 36750 : Hspec(GEN N)
5767 : {
5768 36750 : long v2 = Z_lvalrem(N, 2, &N), v2f = v2 >> 1;
5769 : GEN t;
5770 36750 : if (odd(v2)) { v2f--; N = shifti(N,3); }
5771 32557 : else if (mod4(N)!=3) { v2f--; N = shifti(N,2); }
5772 : /* N fundamental at 2, v2f = v2(f) s.t. N = f^2 D, D fundamental */
5773 36750 : t = addui(3, muliu(subiu(int2n(v2f+1), 3), 2 - kroiu(N,2)));
5774 36750 : return mulii(t, hclassno6(N));
5775 : }
5776 :
5777 : /* Ramanujan tau function for p prime */
5778 : static GEN
5779 14903 : tauprime(GEN p)
5780 : {
5781 14903 : pari_sp av = avma, av2;
5782 : GEN s, p2, p2_7, p_9, T;
5783 : ulong lim, t, tin;
5784 :
5785 14903 : if (absequaliu(p, 2)) return utoineg(24);
5786 : /* p > 2 */
5787 11396 : p2 = sqri(p);
5788 11396 : p2_7 = mului(7, p2);
5789 11396 : p_9 = mului(9, p);
5790 11396 : av2 = avma;
5791 11396 : lim = itou(sqrtint(p));
5792 11396 : tin = mod4(p) == 3? 1: 0;
5793 11396 : s = gen_0;
5794 87178 : for (t = 1; t <= lim; ++t)
5795 : {
5796 75782 : GEN h, a, t2 = sqru(t), D = shifti(subii(p, t2), 2); /* 4(p-t^2) */
5797 : /* t mod 2 != tin <=> D not fundamental at 2 */
5798 75782 : h = ((t&1UL) == tin)? hclassno6(D): Hspec(D);
5799 75782 : a = mulii(powiu(t2,3), addii(p2_7, mulii(t2, subii(shifti(t2,2), p_9))));
5800 75782 : s = addii(s, mulii(a,h));
5801 75782 : if (!(t & 255)) s = gerepileuptoint(av2, s);
5802 : }
5803 : /* 28p^3 - 28p^2 - 90p - 35 */
5804 11396 : T = subii(shifti(mulii(p2_7, subiu(p,1)), 2), addiu(mului(90,p), 35));
5805 11396 : s = shifti(diviuexact(s, 3), 6);
5806 11396 : return gerepileuptoint(av, subii(mulii(mulii(p2,p),T), addui(1, s)));
5807 : }
5808 :
5809 : /* Ramanujan tau function, return 0 for <= 0 */
5810 : GEN
5811 7035 : ramanujantau(GEN n)
5812 : {
5813 7035 : pari_sp ltop = avma;
5814 : GEN T, F, P, E;
5815 : long j, lP;
5816 :
5817 7035 : if (!(F = check_arith_all(n,"ramanujantau")))
5818 : {
5819 7014 : if (signe(n) <= 0) return gen_0;
5820 7007 : F = Z_factor(n);
5821 : }
5822 : else
5823 : {
5824 21 : P = gel(F,1);
5825 21 : if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
5826 : }
5827 :
5828 7014 : P = gel(F,1);
5829 7014 : E = gel(F,2); lP = lg(P);
5830 7014 : T = gen_1;
5831 21917 : for (j = 1; j < lP; j++)
5832 : {
5833 14903 : GEN p = gel(P,j), tp = tauprime(p), t1 = tp, t0 = gen_1;
5834 14903 : long k, e = itou(gel(E,j));
5835 20160 : for (k = 1; k < e; k++)
5836 : {
5837 5257 : GEN t2 = subii(mulii(tp, t1), mulii(powiu(p, 11), t0));
5838 5257 : t0 = t1; t1 = t2;
5839 : }
5840 14903 : T = mulii(T, t1);
5841 : }
5842 7014 : return gerepileuptoint(ltop, T);
5843 : }
|
