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 44140 : 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 44140 : 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 43491 : u_is_gener_expo(ulong p, GEN L0)
55 : {
56 43491 : const ulong q = p >> 1;
57 : long i, l;
58 : GEN L;
59 43491 : if (L0) {
60 3466 : l = lg(L0);
61 3466 : L = cgetg(l, t_VECSMALL);
62 : } else {
63 40025 : L0 = L = u_odd_prime_divisors(q);
64 40020 : l = lg(L);
65 : }
66 43487 : for (i=1; i<l; i++) L[i] = q / uel(L0,i);
67 43487 : return L;
68 : }
69 :
70 : int
71 114099 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
72 : {
73 : long i;
74 114099 : if (krouu(x, p) >= 0) return 0;
75 125091 : for (i=lg(L)-1; i; i--)
76 : {
77 78980 : ulong t = Fl_powu(x, uel(L,i), p);
78 78984 : if (t == p_1 || t == 1) return 0;
79 : }
80 46111 : return 1;
81 : }
82 : /* assume p prime */
83 : ulong
84 137138 : pgener_Fl_local(ulong p, GEN L0)
85 : {
86 137138 : const pari_sp av = avma;
87 137138 : const ulong p_1 = p-1;
88 : long x;
89 : GEN L;
90 137138 : if (p <= 19) switch(p)
91 : { /* quick trivial cases */
92 21 : case 2: return 1;
93 : case 7:
94 17114 : case 17: return 3;
95 76542 : default: return 2;
96 : }
97 43461 : L = u_is_gener_expo(p,L0);
98 43457 : for (x=2;;x++) { if (is_gener_Fl(x,p,p_1,L)) break; }
99 43464 : avma = av; return x;
100 : }
101 : ulong
102 113245 : 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 514 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
108 : {
109 514 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
110 514 : if (t >= 0) return 0;
111 450 : for (i = lg(L)-1; i; i--)
112 : {
113 161 : GEN t = Fp_pow(x, gel(L,i), p);
114 161 : if (equalii(t, p_1) || equali1(t)) return 0;
115 : }
116 289 : return 1;
117 : }
118 :
119 : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
120 : GEN
121 24199 : pgener_Fp_local(GEN p, GEN L0)
122 : {
123 24199 : pari_sp av0 = avma;
124 : GEN x, p_1, L;
125 24199 : if (lgefint(p) == 3)
126 : {
127 : ulong z;
128 23960 : if (p[2] == 2) return gen_1;
129 18038 : if (L0) L0 = ZV_to_nv(L0);
130 18038 : z = pgener_Fl_local(uel(p,2), L0);
131 18038 : 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 23436 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
141 :
142 : ulong
143 68077 : pgener_Zl(ulong p)
144 : {
145 68077 : 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 68077 : if (p == 40487) return 10;
148 : #ifndef LONG_IS_64BIT
149 9721 : return pgener_Fl(p);
150 : #else
151 58356 : 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 102 : if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2))) break;
160 72 : avma = av; return x;
161 : }
162 : #endif
163 : }
164 :
165 : /* p prime. Return a primitive root modulo p^e, e > 1 */
166 : GEN
167 68082 : pgener_Zp(GEN p)
168 : {
169 68082 : 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 17 : if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
177 12 : 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 : if (signe(N) < 0) N = absi(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 3899 : rootsof1_Fl(ulong n, ulong p)
256 : {
257 3899 : pari_sp av = avma;
258 3899 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
259 3899 : ulong z = pgener_Fl_local(p, L);
260 3899 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
261 3899 : 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 0 : }
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 39529 : }
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 254326 : ulogintall(ulong B, ulong y, ulong *ptq)
375 : {
376 : ulong r, r2;
377 : long e;
378 :
379 254326 : if (y == 2)
380 : {
381 4699 : long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
382 4699 : if (ptq) *ptq = 1UL << eB;
383 4699 : return eB;
384 : }
385 249627 : r = y, r2 = 1UL;
386 901080 : for (e=1;; e++)
387 : { /* here, r = y^e, r2 = y^(e-1) */
388 901080 : if (r >= B)
389 : {
390 249516 : if (r != B) { e--; r = r2; }
391 249516 : if (ptq) *ptq = r;
392 249516 : return e;
393 : }
394 651564 : r2 = r;
395 651564 : r = umuluu_or_0(y, r);
396 651564 : if (!r)
397 : {
398 111 : if (ptq) *ptq = r2;
399 111 : return e;
400 : }
401 651453 : }
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 262056 : logintall(GEN B, GEN y, GEN *ptq)
408 : {
409 : pari_sp av;
410 262056 : long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
411 : GEN q, pow2;
412 :
413 262056 : if (lgefint(B) == 3)
414 : {
415 : ulong q;
416 254326 : if (lgefint(y) > 3)
417 : {
418 0 : if (ptq) *ptq = gen_1;
419 0 : return 0;
420 : }
421 254326 : if (!ptq) return ulogintall(B[2], y[2], NULL);
422 46440 : e = ulogintall(B[2], y[2], &q);
423 46440 : *ptq = utoi(q); return e;
424 : }
425 7730 : if (equaliu(y,2))
426 : {
427 159 : if (ptq) *ptq = int2n(eB);
428 159 : return eB;
429 : }
430 7571 : av = avma;
431 7571 : ey = expi(y);
432 : /* eB/(ey+1) - 1 < e <= eB/ey */
433 7571 : emax = eB/ey;
434 7571 : if (emax <= 13) /* e small, be naive */
435 : {
436 1518 : GEN r = y, r2 = gen_1;
437 16124 : for (e=1;; e++)
438 : { /* here, r = y^e, r2 = y^(e-1) */
439 16124 : long fl = cmpii(r, B);
440 16124 : if (fl >= 0)
441 : {
442 1518 : if (fl) { e--; cgiv(r); r = r2; }
443 1518 : if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
444 1518 : return e;
445 : }
446 14606 : r2 = r; r = mulii(r,y);
447 14606 : }
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 6053 : pow2 = new_chunk((long)log2(eB)+2);
454 6053 : gel(pow2,0) = y;
455 6053 : for (i=0, q=y;; )
456 : {
457 35059 : GEN r = gel(pow2,i); /* r = y^2^i */
458 35059 : long fl = cmpii(r,B);
459 35059 : 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 35059 : if (fl > 0) { i--; break; }
466 33048 : q = r;
467 33048 : if (1L<<(i+1) > emax) break;
468 29006 : gel(pow2,++i) = sqri(q);
469 29006 : }
470 :
471 6053 : for (e = 1L<<i;;)
472 : { /* y^e = q < B < r = q * y^(2^i) */
473 33027 : pari_sp av2 = avma;
474 : long fl;
475 : GEN r;
476 33027 : if (--i < 0) break;
477 26981 : r = mulii(q, gel(pow2,i));
478 26981 : fl = cmpii(r, B);
479 26981 : if (fl > 0) avma = av2;
480 : else
481 : {
482 13108 : e += (1L<<i);
483 13108 : q = r;
484 13108 : if (!fl) break; /* B = r */
485 : }
486 26974 : }
487 6053 : if (ptq) *ptq = gerepileuptoint(av, q); else avma = av;
488 6053 : 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 30563 : sqrtint(GEN a)
508 : {
509 30563 : if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
510 30563 : switch (signe(a))
511 : {
512 30549 : 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 13185619 : carremod(ulong A)
526 : {
527 13185619 : 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 13185619 : 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 13185619 : 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 13185619 : const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
537 26371238 : return (carresmod64[A & 0x3fUL]
538 5128332 : && carresmod63[A % 63UL]
539 3107832 : && carresmod65[A % 65UL]
540 15770577 : && carresmod11[A % 11UL]);
541 : }
542 :
543 : /* emulate Z_issquareall on single-word integers */
544 : long
545 11695097 : uissquareall(ulong A, ulong *sqrtA)
546 : {
547 11695097 : if (!A) { *sqrtA = 0; return 1; }
548 11695097 : if (carremod(A))
549 : {
550 1745286 : ulong a = usqrt(A);
551 1745274 : if (a * a == A) { *sqrtA = a; return 1; }
552 : }
553 10040838 : return 0;
554 : }
555 : long
556 120605 : uissquare(ulong A)
557 : {
558 120605 : if (!A) return 1;
559 120605 : if (carremod(A))
560 : {
561 3485 : ulong a = usqrt(A);
562 3485 : if (a * a == A) return 1;
563 : }
564 117144 : return 0;
565 : }
566 :
567 : long
568 6272042 : Z_issquareall(GEN x, GEN *pt)
569 : {
570 : pari_sp av;
571 : GEN y, r;
572 :
573 6272042 : switch(signe(x))
574 : {
575 2190935 : case -1: return 0;
576 1064 : case 0: if (pt) *pt=gen_0; return 1;
577 : }
578 4080043 : if (lgefint(x) == 3)
579 : {
580 2710154 : ulong u = uel(x,2), a;
581 2710154 : if (!pt) return uissquare(u);
582 2589549 : if (!uissquareall(u, &a)) return 0;
583 1377785 : *pt = utoipos(a); return 1;
584 : }
585 1369889 : if (!carremod(umodiu(x, 64*63*65*11))) return 0;
586 608580 : av = avma; y = sqrtremi(x, &r);
587 608580 : if (r != gen_0) { avma = av; return 0; }
588 17732 : if (pt) { *pt = y; avma = (pari_sp)y; } else avma = av;
589 17732 : 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 2975 : polissquareall(GEN x, GEN *pt)
608 : {
609 : pari_sp av;
610 : long v;
611 : GEN y, a, b, p;
612 :
613 2975 : if (!signe(x))
614 : {
615 7 : if (pt) *pt = gcopy(x);
616 7 : return 1;
617 : }
618 2968 : 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 273 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
673 : {
674 273 : if (d == 1)
675 : { /* mod p: faster */
676 189 : if (!Fp_ispower(b, K, p)) return 0;
677 189 : 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 252 : 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 413 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
691 : {
692 : GEN t, A;
693 413 : long v = Z_pvalrem(*pa, p, &A), d = e - v;
694 413 : if (d <= 0) t = gen_0;
695 : else
696 : {
697 : ulong r;
698 357 : v = udivui_rem(v, K, &r);
699 357 : if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
700 252 : if (L && v) t = mulii(t, powiu(p, v));
701 : }
702 308 : if (q) *pa = modii(*pa, q);
703 308 : if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
704 308 : return 1;
705 : }
706 : long
707 308 : 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 308 : if (!signe(a))
716 : {
717 14 : if (pt) {
718 14 : GEN t = cgetg(3, t_INTMOD);
719 14 : gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
720 : }
721 14 : return 1;
722 : }
723 : /* a != 0 */
724 294 : av = avma;
725 :
726 294 : 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 294 : if (!mod2(K)
740 189 : && kronecker(a, shifti(q,-vali(q))) == -1) { avma = av; return 0; }
741 287 : L = pt? vectrunc_init(expi(q)+1): NULL;
742 287 : u_forprime_init(&S, 2, tridiv_bound(q));
743 287 : while ((pp = u_forprime_next(&S)))
744 : {
745 : int stop;
746 883351 : e = Z_lvalrem_stop(&q, pp, &stop);
747 883351 : if (!e) continue;
748 189 : if (!handle_pe(&a, q, L, K, utoipos(pp), e)) { avma = av; return 0; }
749 147 : if (stop)
750 : {
751 112 : if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
752 112 : 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 : {
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 42 : }
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 : {
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 84 : }
801 : }
802 : }
803 : END:
804 182 : if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
805 182 : 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 163100 : issquareall(GEN x, GEN *pt)
831 : {
832 163100 : long tx = typ(x);
833 : GEN F;
834 : pari_sp av;
835 :
836 163100 : if (!pt) return issquare(x);
837 20258 : switch(tx)
838 : {
839 2324 : 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 2891 : case t_POL: return polissquareall(x,pt);
849 7 : case t_RFRAC: av = avma;
850 7 : F = cgetg(3, t_RFRAC);
851 7 : if ( !issquareall(gel(x,1), &gel(F,1))
852 7 : || !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 157857 : issquare(GEN x)
871 : {
872 : pari_sp av;
873 : GEN a, p;
874 : long i, v;
875 :
876 157857 : switch(typ(x))
877 : {
878 : case t_INT:
879 142737 : 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 = diviu_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 658 : polispower(GEN x, GEN K, GEN *pt)
975 : {
976 : pari_sp av;
977 658 : long v, d, k = itos(K);
978 : GEN y, a, b;
979 658 : GEN T = NULL, p = NULL;
980 :
981 658 : if (!signe(x))
982 : {
983 7 : if (pt) *pt = gcopy(x);
984 7 : return 1;
985 : }
986 651 : d = degpol(x);
987 651 : if (d % k) return 0; /* degree not multiple of k */
988 644 : av = avma;
989 644 : 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 308 : v = RgX_valrem(x, &x);
1003 308 : if (v % k) return 0;
1004 301 : v /= k;
1005 301 : a = gel(x,2); b = NULL;
1006 301 : if (!ispower(a, K, &b)) { avma = av; return 0; }
1007 287 : if (d)
1008 : {
1009 280 : GEN p = characteristic(x);
1010 280 : a = leading_coeff(x);
1011 280 : if (!ispower(a, K, &b)) { avma = av; return 0; }
1012 280 : x = RgX_normalize(x);
1013 280 : if (signe(p) && cmpii(p,K) <= 0)
1014 0 : pari_err_IMPL("ispower(general t_POL) in small characteristic");
1015 280 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
1016 280 : if (!RgX_equal(powgi(y, K), x)) { avma = av; return 0; }
1017 : }
1018 : else
1019 7 : y = pol_1(varn(x));
1020 287 : if (pt)
1021 : {
1022 287 : if (!gequal1(a))
1023 : {
1024 14 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
1025 14 : y = gmul(b,y);
1026 : }
1027 287 : if (v) y = RgX_shift_shallow(y, v);
1028 287 : *pt = gerepilecopy(av, y);
1029 : }
1030 0 : else avma = av;
1031 287 : return 1;
1032 : }
1033 :
1034 : long
1035 1547 : Z_ispowerall(GEN x, ulong k, GEN *pt)
1036 : {
1037 1547 : long s = signe(x);
1038 : ulong mask;
1039 1547 : if (!s) { if (pt) *pt = gen_0; return 1; }
1040 1547 : if (s > 0) {
1041 1421 : 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(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 210 : Fp_ispower(GEN x, GEN K, GEN p)
1059 : {
1060 210 : pari_sp av = avma;
1061 : GEN p_1;
1062 : long r;
1063 210 : x = modii(x, p);
1064 210 : if (!signe(x) || equali1(x)) { avma = av; return 1; }
1065 : /* implies p > 2 */
1066 119 : p_1 = subiu(p,1);
1067 119 : K = gcdii(K, p_1);
1068 119 : if (absequaliu(K, 2)) { r = kronecker(x,p); avma = av; return (r > 0); }
1069 42 : x = Fp_pow(x, diviiexact(p_1,K), p);
1070 42 : 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 7002947 : ispower(GEN x, GEN K, GEN *pt)
1117 : {
1118 : GEN z;
1119 :
1120 7002947 : if (!K) return gisanypower(x, pt);
1121 2765 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
1122 2765 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
1123 2765 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
1124 2716 : switch(typ(x)) {
1125 : case t_INT:
1126 693 : 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 168 : 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 651 : 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 864619 : Z_isanypower_nosmalldiv(GEN *px)
1249 : { /* any prime divisor of x is > 102 */
1250 864619 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
1251 864619 : const double LOG103 = 4.6347; /* lower bound for log(103) */
1252 : forprime_t T;
1253 864619 : ulong mask = 7, e2;
1254 : long k, ex;
1255 864619 : GEN y, x = *px;
1256 :
1257 864619 : k = 1;
1258 864619 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
1259 864619 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
1260 864619 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
1261 864619 : if (u_forprime_init(&T, 11, e2))
1262 : {
1263 16954 : GEN logx = NULL;
1264 16954 : const ulong Q = 30011; /* prime */
1265 : ulong p, xmodQ;
1266 16954 : 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 33950 : 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 16954 : if (DEBUGLEVEL>4) err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x));
1276 16954 : xmodQ = umodiu(x, Q);
1277 : /* test Q | x, just in case */
1278 16954 : if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
1279 : /* x^(1/p) < 2^31 */
1280 16940 : p = T.p;
1281 16940 : if (p <= e2)
1282 : {
1283 16926 : logx = logr_abs( itor(x, DEFAULTPREC) );
1284 16926 : dlogx = rtodbl(logx);
1285 16926 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1286 : }
1287 153902 : 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 120022 : pari_sp av = avma;
1291 : long e;
1292 120022 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
1293 120022 : ulong ymodQ = umodiu(y,Q);
1294 120022 : 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 120001 : p = u_forprime_next(&T);
1304 : }
1305 : }
1306 864605 : *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 7000813 : Z_isanypower_aux(GEN x, GEN *pty)
1318 : {
1319 : long ex, v, i, l, k;
1320 : GEN y, P, E;
1321 7000813 : ulong mask, e = 0;
1322 :
1323 7000813 : if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
1324 :
1325 7000799 : if (signe(x) < 0) x = negi(x);
1326 7000799 : k = l = 1;
1327 7000799 : P = cgetg(26 + 1, t_VECSMALL);
1328 7000799 : E = cgetg(26 + 1, t_VECSMALL);
1329 : /* trial division */
1330 122940746 : for(i = 0; i < 26; i++)
1331 : {
1332 60135899 : ulong p = tinyprimes[i];
1333 : int stop;
1334 60135899 : v = Z_lvalrem_stop(&x, p, &stop);
1335 60135899 : 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 54717612 : if (stop) {
1342 248038 : if (is_pm1(x)) k = e;
1343 248038 : goto END;
1344 : }
1345 : }
1346 :
1347 1334474 : 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 828773 : k = Z_isanypower_nosmalldiv(&x);
1377 : END:
1378 7000799 : if (pty && k != 1)
1379 : {
1380 8008 : 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 8008 : *pty = x;
1387 : }
1388 7000799 : return k == 1? 0: k;
1389 : }
1390 :
1391 : long
1392 7000813 : Z_isanypower(GEN x, GEN *pty)
1393 : {
1394 7000813 : pari_sp av = avma;
1395 7000813 : long k = Z_isanypower_aux(x, pty);
1396 7000813 : if (!k) { avma = av; return 0; }
1397 8071 : 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 8043 : if (pty) *pty = gerepilecopy(av, *pty); else avma = av;
1415 8043 : 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 80670 : Z_ispow2(GEN n)
1427 : {
1428 : GEN xp;
1429 : long i, lx;
1430 : ulong u;
1431 80670 : if (signe(n) != 1) return 0;
1432 80663 : xp = int_LSW(n);
1433 80663 : lx = lgefint(n);
1434 80663 : u = *xp;
1435 80947 : for (i = 3; i < lx; ++i)
1436 : {
1437 78128 : 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 842008 : isprimepower_i(GEN n, GEN *pt, long flag)
1446 : {
1447 842008 : pari_sp av = avma;
1448 : long i, v;
1449 :
1450 842008 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1451 842007 : if (signe(n) <= 0) return 0;
1452 :
1453 842007 : if (lgefint(n) == 3)
1454 : {
1455 : ulong p;
1456 541115 : v = uisprimepower(n[2], &p);
1457 541115 : if (v)
1458 : {
1459 54903 : if (pt) *pt = utoipos(p);
1460 54903 : return v;
1461 : }
1462 486212 : return 0;
1463 : }
1464 1662564 : for (i = 0; i < 26; i++)
1465 : {
1466 1626718 : ulong p = tinyprimes[i];
1467 1626718 : v = Z_lvalrem(n, p, &n);
1468 1626719 : if (v)
1469 : {
1470 265047 : avma = av;
1471 265047 : if (!is_pm1(n)) return 0;
1472 652 : if (pt) *pt = utoipos(p);
1473 651 : 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 1910 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1486 :
1487 : long
1488 541752 : 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 541752 : const ulong CUTOFF = 200UL;
1492 541752 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1493 541752 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1494 : #ifdef LONG_IS_64BIT
1495 : /* primes preceeding the appropriate root of ULONG_MAX. */
1496 481512 : const ulong ROOT9 = 137;
1497 481512 : const ulong ROOT8 = 251;
1498 481512 : const ulong ROOT7 = 563;
1499 481512 : const ulong ROOT5 = 7129;
1500 481512 : const ulong ROOT4 = 65521;
1501 : #else
1502 60240 : const ulong ROOT9 = 11;
1503 60240 : const ulong ROOT8 = 13;
1504 60240 : const ulong ROOT7 = 23;
1505 60240 : const ulong ROOT5 = 83;
1506 60240 : const ulong ROOT4 = 251;
1507 : #endif
1508 : ulong mask;
1509 : long v, i;
1510 : int e;
1511 541752 : if (n < 2) return 0;
1512 541738 : if (!odd(n)) {
1513 270725 : if (n & (n-1)) return 0;
1514 921 : *pp = 2; return vals(n);
1515 : }
1516 271013 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1517 3653910 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1518 : {
1519 3594835 : ulong p = tinyprimes[i];
1520 3594835 : if (n % p == 0)
1521 : {
1522 211532 : v = u_lvalrem(n, p, &n);
1523 211532 : 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 669708731 : ome(long t)
1573 : {
1574 669708731 : switch(t & 7)
1575 : {
1576 : case 3:
1577 385236455 : case 5: return 1;
1578 284472276 : default: return 0;
1579 : }
1580 : }
1581 : /* t a t_INT, is t = 3,5 mod 8 ? */
1582 : static int
1583 5355910 : gome(GEN t)
1584 5355910 : { return signe(t)? ome( mod2BIL(t) ): 0; }
1585 :
1586 : /* assume y odd, return kronecker(x,y) * s */
1587 : static long
1588 524201049 : krouu_s(ulong x, ulong y, long s)
1589 : {
1590 524201049 : ulong x1 = x, y1 = y, z;
1591 2751025014 : while (x1)
1592 : {
1593 1702745768 : long r = vals(x1);
1594 1702809557 : if (r)
1595 : {
1596 924754336 : if (odd(r) && ome(y1)) s = -s;
1597 924567695 : x1 >>= r;
1598 : }
1599 1702622916 : if (x1 & y1 & 2) s = -s;
1600 1702622916 : z = y1 % x1; y1 = x1; x1 = z;
1601 : }
1602 524078197 : return (y1 == 1)? s: 0;
1603 : }
1604 :
1605 : long
1606 6103539 : kronecker(GEN x, GEN y)
1607 : {
1608 6103539 : pari_sp av = avma;
1609 6103539 : long s = 1, r;
1610 : ulong xu, yu;
1611 :
1612 6103539 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
1613 6103539 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
1614 6103539 : 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 6103539 : r = vali(y);
1620 6103539 : if (r)
1621 : {
1622 12202 : if (!mpodd(x)) { avma = av; return 0; }
1623 10487 : if (odd(r) && gome(x)) s = -s;
1624 10487 : y = shifti(y,-r);
1625 : }
1626 6101824 : x = modii(x,y);
1627 13108316 : while (lgefint(x) > 3) /* x < y */
1628 : {
1629 : GEN z;
1630 904668 : r = vali(x);
1631 904668 : if (r)
1632 : {
1633 493700 : if (odd(r) && gome(y)) s = -s;
1634 493700 : x = shifti(x,-r);
1635 : }
1636 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1637 904668 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
1638 904668 : z = remii(y,x); y = x; x = z;
1639 904668 : 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 6101824 : xu = itou(x);
1646 6101824 : if (!xu) return is_pm1(y)? s: 0;
1647 6080650 : r = vals(xu);
1648 6080650 : if (r)
1649 : {
1650 4131487 : if (odd(r) && gome(y)) s = -s;
1651 4131487 : xu >>= r;
1652 : }
1653 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1654 6080650 : if (xu & mod2BIL(y) & 2) s = -s;
1655 6080650 : yu = umodiu(y, xu);
1656 6080650 : avma = av; return krouu_s(yu, xu, s);
1657 : }
1658 :
1659 : long
1660 30737 : krois(GEN x, long y)
1661 : {
1662 : ulong yu;
1663 30737 : long s = 1;
1664 :
1665 30737 : 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 30730 : yu = (ulong)y;
1672 30730 : if (!odd(yu))
1673 : {
1674 : long r;
1675 14021 : if (!mpodd(x)) return 0;
1676 10437 : r = vals(yu); yu >>= r;
1677 10437 : if (odd(r) && gome(x)) s = -s;
1678 : }
1679 27146 : return krouu_s(umodiu(x, yu), yu, s);
1680 : }
1681 : /* assume y != 0 */
1682 : long
1683 340149923 : kroiu(GEN x, ulong y)
1684 : {
1685 : long r;
1686 340149923 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
1687 2116590 : if (!mpodd(x)) return 0;
1688 2094484 : r = vals(y); y >>= r;
1689 2094484 : 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 120548 : krouodd(ulong x, GEN y, long s)
1695 : {
1696 : long r;
1697 120548 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
1698 41882 : if (!x) return 0; /* y != 1 */
1699 41882 : r = vals(x);
1700 41882 : if (r)
1701 : {
1702 7368 : if (odd(r) && gome(y)) s = -s;
1703 7368 : x >>= r;
1704 : }
1705 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1706 41882 : if (x & mod2BIL(y) & 2) s = -s;
1707 41882 : return krouu_s(umodiu(y,x), x, s);
1708 : }
1709 :
1710 : long
1711 120045 : krosi(long x, GEN y)
1712 : {
1713 120045 : const pari_sp av = avma;
1714 120045 : long s = 1, r;
1715 120045 : 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 120045 : r = vali(y);
1721 120045 : 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 120045 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
1728 120045 : s = krouodd((ulong)x, y, s);
1729 120045 : avma = av; return s;
1730 : }
1731 :
1732 : long
1733 503 : kroui(ulong x, GEN y)
1734 : {
1735 503 : const pari_sp av = avma;
1736 503 : long s = 1, r;
1737 503 : switch (signe(y))
1738 : {
1739 0 : case -1: y = negi(y); break;
1740 0 : case 0: return x==1UL;
1741 : }
1742 503 : r = vali(y);
1743 503 : 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 503 : s = krouodd(x, y, s);
1750 503 : avma = av; return s;
1751 : }
1752 :
1753 : long
1754 73218203 : kross(long x, long y)
1755 : {
1756 : ulong yu;
1757 73218203 : long s = 1;
1758 :
1759 73218203 : if (y <= 0)
1760 : {
1761 427 : if (y == 0) return (labs(x)==1);
1762 399 : yu = (ulong)-y; if (x < 0) s = -1;
1763 : }
1764 : else
1765 73217776 : yu = (ulong)y;
1766 73218175 : if (!odd(yu))
1767 : {
1768 : long r;
1769 17118753 : if (!odd(x)) return 0;
1770 12295984 : r = vals(yu); yu >>= r;
1771 12295984 : if (odd(r) && ome(x)) s = -s;
1772 : }
1773 68395406 : x %= (long)yu; if (x < 0) x += yu;
1774 68395406 : return krouu_s((ulong)x, yu, s);
1775 : }
1776 :
1777 : long
1778 109298287 : krouu(ulong x, ulong y)
1779 : {
1780 : long r;
1781 109298287 : if (odd(y)) return krouu_s(x, y, 1);
1782 1746 : if (!odd(x)) return 0;
1783 1746 : r = vals(y); y >>= r;
1784 1746 : 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 :
1925 : static ulong
1926 26681469 : Fl_2gener_pre_all(long e, ulong p, ulong pi)
1927 : {
1928 : ulong y, m;
1929 : long k, i;
1930 26681469 : ulong q = (p-1) >> e; /* q = (p-1)/2^oo is odd */
1931 46538440 : for (k=2; ; k++)
1932 : { /* loop terminates for k < p (even if p composite) */
1933 46538440 : i = krouu(k, p);
1934 46538447 : if (i >= 0)
1935 : {
1936 19856973 : if (i) continue;
1937 7 : pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1938 : }
1939 26681474 : y = m = Fl_powu_pre(k, q, p, pi);
1940 72248261 : for (i=1; i<e; i++)
1941 45566703 : if ((m = Fl_sqr_pre(m, p, pi)) == 1) break;
1942 26681558 : if (i == e) break; /* success */
1943 19856971 : }
1944 26681553 : return y;
1945 : }
1946 :
1947 : ulong
1948 149698 : Fl_2gener_pre(ulong p, ulong pi)
1949 149698 : { return Fl_2gener_pre_all(vals(p-1), p, pi); }
1950 :
1951 : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
1952 : ulong
1953 63922839 : Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
1954 : {
1955 : long i, e, k;
1956 : ulong p1, q, v, w;
1957 :
1958 63922839 : if (!a) return 0;
1959 62597040 : p1 = p - 1; e = vals(p1);
1960 62602321 : if (e == 0) /* p = 2 */
1961 : {
1962 418894 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1963 418887 : return ((a & 1) == 0)? 0: 1;
1964 : }
1965 62183427 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
1966 62183427 : if (e == 1) y = p1;
1967 26535563 : else if (y==0) y = Fl_2gener_pre_all(e, p, pi);
1968 62183420 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
1969 62161333 : if (!p1) return 0;
1970 62161333 : v = Fl_mul_pre(a, p1, p, pi);
1971 62167246 : w = Fl_mul_pre(v, p1, p, pi);
1972 146466282 : while (w != 1)
1973 : { /* a*w = v^2, y primitive 2^e-th root of 1
1974 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1975 22199261 : p1 = Fl_sqr_pre(w,p,pi);
1976 22199262 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
1977 22199262 : if (k == e) return ~0UL;
1978 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1979 22135077 : p1 = y;
1980 22135077 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
1981 22135077 : y = Fl_sqr_pre(p1, p, pi); e = k;
1982 22135077 : w = Fl_mul_pre(y, w, p, pi);
1983 22135077 : v = Fl_mul_pre(v, p1, p, pi);
1984 : }
1985 62100911 : p1 = p - v; if (v > p1) v = p1;
1986 62100911 : return v;
1987 : }
1988 :
1989 : ulong
1990 60506084 : Fl_sqrt(ulong a, ulong p)
1991 : {
1992 60506084 : ulong pi = get_Fl_red(p);
1993 60503834 : return Fl_sqrt_pre_i(a, 0, p, pi);
1994 : }
1995 :
1996 : ulong
1997 3387692 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
1998 : {
1999 3387692 : return Fl_sqrt_pre_i(a, 0, p, pi);
2000 : }
2001 :
2002 : static ulong
2003 46039 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
2004 : {
2005 : ulong x, y, m;
2006 46039 : ulong le1 = upowuu(l, e-1);
2007 72798 : for (x = 2; ; x++)
2008 : {
2009 72798 : y = Fl_powu_pre(x, r, p, pi);
2010 72798 : if (y==1) continue;
2011 56442 : m = Fl_powu_pre(y, le1, p, pi);
2012 56442 : if (m != 1) break;
2013 26759 : }
2014 46039 : *pt_m = m;
2015 46039 : return y;
2016 : }
2017 :
2018 : /* solve x^l = a , l prime in G of order q.
2019 : *
2020 : * q = (l^e)*r, e >= 1, (r,l) = 1
2021 : * y generates the l-Sylow of G
2022 : * m = y^(l^(e-1)) != 1 */
2023 : static ulong
2024 110079 : Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
2025 : {
2026 : ulong p1, v, w, z, dl;
2027 : ulong u2;
2028 110079 : if (a==0) return a;
2029 110079 : u2 = Fl_inv(l%r, r);
2030 110080 : v = Fl_powu_pre(a, u2, p, pi);
2031 110080 : w = Fl_powu_pre(v, l, p, pi);
2032 110082 : w = Fl_mul_pre(w, Fl_inv(a, p), p, pi);
2033 110066 : if (w==1) return v;
2034 45092 : if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
2035 109272 : while (w!=1)
2036 : {
2037 49440 : ulong k = 0;
2038 49440 : p1 = w;
2039 : do
2040 : {
2041 73204 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
2042 73204 : k++;
2043 73204 : } while (p1!=1);
2044 49440 : if (k==e) return ULONG_MAX;
2045 19088 : dl = Fl_log_pre(z, m, l, p, pi);
2046 19088 : dl = Fl_neg(dl, l);
2047 19088 : p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
2048 19088 : m = Fl_powu_pre(m, dl, p, pi);
2049 19088 : e = k;
2050 19088 : v = Fl_mul_pre(p1,v,p,pi);
2051 19088 : y = Fl_powu_pre(p1,l,p,pi);
2052 19088 : w = Fl_mul_pre(y,w,p,pi);
2053 : }
2054 14740 : return v;
2055 : }
2056 :
2057 : static ulong
2058 109402 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
2059 : {
2060 109402 : ulong r, e = u_lvalrem(p-1, l, &r);
2061 109401 : return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
2062 : }
2063 :
2064 : ulong
2065 109401 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
2066 : {
2067 109401 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2068 : }
2069 :
2070 : ulong
2071 0 : Fl_sqrtl(ulong a, ulong l, ulong p)
2072 : {
2073 0 : ulong pi = get_Fl_red(p);
2074 0 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2075 : }
2076 :
2077 : ulong
2078 64934 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
2079 : {
2080 64934 : ulong m, q = p-1, z;
2081 64934 : ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
2082 64934 : if (a==0)
2083 : {
2084 48118 : if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
2085 48111 : if (zetan) *zetan = 1UL;
2086 48111 : return 0;
2087 : }
2088 16816 : if (n==1)
2089 : {
2090 0 : if (zetan) *zetan = 1;
2091 0 : return n < 0? Fl_inv(a,p): a;
2092 : }
2093 16816 : if (n==2)
2094 : {
2095 3430 : if (zetan) *zetan = p-1;
2096 3430 : return Fl_sqrt_pre_i(a, 0, p, pi);
2097 : }
2098 13386 : if (a == 1 && !zetan) return a;
2099 7471 : m = ugcd(nn, q);
2100 7471 : z = 1;
2101 7471 : if (m!=1)
2102 : {
2103 912 : GEN F = factoru(m);
2104 : long i, j, e;
2105 : ulong r, zeta, y, l;
2106 1915 : for (i = nbrows(F); i; i--)
2107 : {
2108 1017 : l = ucoeff(F,i,1);
2109 1017 : j = ucoeff(F,i,2);
2110 1017 : e = u_lvalrem(q,l, &r);
2111 1017 : y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
2112 1017 : if (zetan)
2113 478 : z = Fl_mul_pre(z, Fl_powu_pre(y, upowuu(l,e-j), p, pi), p, pi);
2114 1017 : if (a!=1)
2115 : do
2116 : {
2117 679 : a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
2118 665 : if (!a) return ULONG_MAX;
2119 665 : } while (--j);
2120 : }
2121 : }
2122 7457 : if (m != nn)
2123 : {
2124 6587 : ulong qm = q/m, nm = nn/m;
2125 6587 : a = Fl_powu_pre(a, Fl_inv(nm%qm, qm), p, pi);
2126 : }
2127 7457 : if (n < 0) a = Fl_inv(a, p);
2128 7457 : if (zetan) *zetan = z;
2129 7457 : return a;
2130 : }
2131 :
2132 : ulong
2133 64934 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
2134 : {
2135 64934 : ulong pi = get_Fl_red(p);
2136 64934 : return Fl_sqrtn_pre(a, n, p, pi, zetan);
2137 : }
2138 :
2139 : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
2140 : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
2141 : * and in average is about 2 or 2.5 times worse. But try both algorithms for
2142 : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
2143 : *
2144 : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
2145 : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
2146 : * If (a|p)=1, then sqrt(a) is in F_p.
2147 : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
2148 :
2149 : /* compute y^2, y = y[1] + y[2] X */
2150 : static GEN
2151 449 : sqrt_Cipolla_sqr(void *data, GEN y)
2152 : {
2153 449 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
2154 449 : GEN u2 = sqri(u), v2 = sqri(v);
2155 449 : v = subii(sqri(addii(v,u)), addii(u2,v2));
2156 449 : u = addii(u2, mulii(v2,n));
2157 : /* NOT mkvec2: must be gerepileupto-able */
2158 449 : retmkvec2(modii(u,p), modii(v,p));
2159 : }
2160 : /* compute (t+X) y^2 */
2161 : static GEN
2162 23 : sqrt_Cipolla_msqr(void *data, GEN y)
2163 : {
2164 23 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
2165 23 : ulong t = gt[2];
2166 23 : GEN d = addii(u, mului(t,v)), d2= sqri(d);
2167 23 : GEN b = remii(mulii(a,v), p);
2168 23 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
2169 23 : v = subii(d2, mulii(b,v));
2170 : /* NOT mkvec2: must be gerepileupto-able */
2171 23 : retmkvec2(modii(u,p), modii(v,p));
2172 : }
2173 : /* assume a reduced mod p [ otherwise correct but inefficient ] */
2174 : static GEN
2175 8 : sqrt_Cipolla(GEN a, GEN p)
2176 : {
2177 : pari_sp av1;
2178 : GEN u, v, n, y, pov2;
2179 : ulong t;
2180 :
2181 8 : if (kronecker(a, p) < 0) return NULL;
2182 8 : pov2 = shifti(p,-1);
2183 8 : if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
2184 :
2185 8 : av1 = avma;
2186 41 : for(t=1; ; t++)
2187 : {
2188 41 : n = subsi((long)(t*t), a);
2189 41 : if (kronecker(n, p) < 0) break;
2190 33 : avma = av1;
2191 33 : }
2192 :
2193 : /* compute (t+X)^((p-1)/2) =: u+vX */
2194 8 : u = utoipos(t);
2195 8 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
2196 : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
2197 : /* Now u+vX = (t+X)^((p-1)/2); thus
2198 : * (u+vX)(t+X) = sqrt(a) + 0 X
2199 : * Whence,
2200 : * sqrt(a) = (u+vt)t - v*a
2201 : * 0 = (u+vt)
2202 : * Thus a square root is v*a */
2203 :
2204 8 : v = Fp_mul(gel(y, 2), a, p);
2205 8 : if (cmpii(v,pov2) > 0) v = subii(p,v);
2206 8 : return v;
2207 : }
2208 :
2209 : /* Return NULL if p is found to be composite */
2210 : static GEN
2211 2800 : Fp_2gener_all(long e, GEN p)
2212 : {
2213 : GEN y, m;
2214 : long k;
2215 2800 : GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
2216 2800 : if (e==0 && !equaliu(p,2)) return NULL;
2217 9265 : for (k=2; ; k++)
2218 : {
2219 9265 : long i = krosi(k, p);
2220 9265 : if (i >= 0)
2221 : {
2222 6465 : if (i) continue;
2223 0 : return NULL;
2224 : }
2225 2800 : y = m = Fp_pow(utoi(k), q, p);
2226 9402 : for (i=1; i<e; i++)
2227 6602 : if (equali1(m = Fp_sqr(m, p))) break;
2228 2800 : if (i == e) break; /* success */
2229 6465 : }
2230 2800 : return y;
2231 : }
2232 :
2233 : /* Return NULL if p is found to be composite */
2234 : GEN
2235 980 : Fp_2gener(GEN p)
2236 980 : { return Fp_2gener_all(vali(subis(p,1)),p); }
2237 :
2238 : /* smallest square root */
2239 : static GEN
2240 27537 : choose_sqrt(GEN v, GEN p)
2241 : {
2242 27537 : pari_sp av = avma;
2243 27537 : GEN q = subii(p,v);
2244 27537 : if (cmpii(v,q) > 0) v = q; else avma = av;
2245 27537 : return v;
2246 : }
2247 : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
2248 : GEN
2249 2347867 : Fp_sqrt_i(GEN a, GEN y, GEN p)
2250 : {
2251 2347867 : pari_sp av = avma;
2252 : long i, k, e;
2253 : GEN p1, q, v, w;
2254 :
2255 2347867 : if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
2256 2347867 : if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
2257 2347867 : if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
2258 2347867 : if (lgefint(p) == 3)
2259 : {
2260 2320223 : ulong pp = uel(p,2), u = Fl_sqrt(umodiu(a, pp), pp);
2261 2320209 : if (u == ~0UL) return NULL;
2262 2320167 : return utoi(u);
2263 : }
2264 :
2265 27644 : a = modii(a, p); if (!signe(a)) { avma = av; return gen_0; }
2266 27553 : p1 = subiu(p,1); e = vali(p1);
2267 27553 : if (e <= 2)
2268 : { /* direct formulas more efficient */
2269 : pari_sp av2;
2270 22547 : if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
2271 22547 : if (e == 1)
2272 : {
2273 13340 : q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
2274 13340 : v = Fp_pow(a, q, p);
2275 : }
2276 : else
2277 : { /* Atkin's formula */
2278 9207 : GEN i, a2 = shifti(a,1);
2279 9207 : if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
2280 9207 : q = shifti(p1, -3); /* (p-5)/8 */
2281 9207 : v = Fp_pow(a2, q, p);
2282 9207 : i = Fp_mul(a2, Fp_sqr(v,p), p); /* i^2 = -1 */
2283 9207 : v = Fp_mul(a, Fp_mul(v, subiu(i,1), p), p);
2284 : }
2285 22547 : av2 = avma;
2286 : /* must check equality in case (a/p) = -1 or p not prime */
2287 22547 : e = equalii(Fp_sqr(v,p), a); avma = av2;
2288 22547 : return e? gerepileuptoint(av,choose_sqrt(v,p)): NULL;
2289 : }
2290 : /* On average, Cipolla is better than Tonelli/Shanks if and only if
2291 : * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
2292 5006 : if (e*(e-1) > 20 + 8 * expi(p))
2293 : {
2294 8 : v = sqrt_Cipolla(a,p); if (!v) { avma = av; return NULL; }
2295 8 : return gerepileuptoint(av,v);
2296 : }
2297 4998 : if (!y)
2298 : {
2299 1820 : y = Fp_2gener_all(e, p);
2300 1820 : if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
2301 : }
2302 4998 : q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
2303 4998 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
2304 4998 : v = Fp_mul(a, p1, p);
2305 4998 : w = Fp_mul(v, p1, p);
2306 17035 : while (!equali1(w))
2307 : { /* a*w = v^2, y primitive 2^e-th root of 1
2308 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
2309 7039 : p1 = Fp_sqr(w,p);
2310 7039 : for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
2311 7039 : if (k == e) { avma=av; return NULL; } /* p composite or (a/p) != 1 */
2312 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
2313 7039 : p1 = y;
2314 7039 : for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
2315 7039 : y = Fp_sqr(p1, p); e = k;
2316 7039 : w = Fp_mul(y, w, p);
2317 7039 : v = Fp_mul(v, p1, p);
2318 7039 : if (gc_needed(av,1))
2319 : {
2320 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
2321 0 : gerepileall(av,3, &y,&w,&v);
2322 : }
2323 : }
2324 4998 : return gerepileuptoint(av, choose_sqrt(v,p));
2325 : }
2326 :
2327 : GEN
2328 2333440 : Fp_sqrt(GEN a, GEN p)
2329 : {
2330 2333440 : return Fp_sqrt_i(a, NULL, p);
2331 : }
2332 :
2333 : /*********************************************************************/
2334 : /** **/
2335 : /** GCD & BEZOUT **/
2336 : /** **/
2337 : /*********************************************************************/
2338 :
2339 : GEN
2340 19119610 : lcmii(GEN x, GEN y)
2341 : {
2342 : pari_sp av;
2343 : GEN a, b;
2344 19119610 : if (!signe(x) || !signe(y)) return gen_0;
2345 19119610 : av = avma;
2346 19119610 : a = gcdii(x,y); if (!equali1(a)) y = diviiexact(y,a);
2347 19119610 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
2348 : }
2349 :
2350 : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
2351 : * set *pd = gcd(x,N) */
2352 : GEN
2353 2671985 : Fp_invgen(GEN x, GEN N, GEN *pd)
2354 : {
2355 : GEN d, d0, e, v;
2356 2671985 : if (lgefint(N) == 3)
2357 : {
2358 2378130 : ulong dd, NN = N[2], xx = umodiu(x,NN);
2359 2378130 : if (!xx) { *pd = N; return gen_0; }
2360 2378130 : xx = Fl_invgen(xx, NN, &dd);
2361 2378130 : *pd = utoi(dd); return utoi(xx);
2362 : }
2363 293855 : *pd = d = bezout(x, N, &v, NULL);
2364 293855 : if (equali1(d)) return v;
2365 : /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
2366 231061 : e = diviiexact(N,d);
2367 231061 : d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
2368 231061 : if (equali1(d0)) return v;
2369 109541 : if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
2370 109541 : return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
2371 : }
2372 :
2373 : /*********************************************************************/
2374 : /** **/
2375 : /** CHINESE REMAINDERS **/
2376 : /** **/
2377 : /*********************************************************************/
2378 :
2379 : /* Chinese Remainder Theorem. x and y must have the same type (integermod,
2380 : * polymod, or polynomial/vector/matrix recursively constructed with these
2381 : * as coefficients). Creates (with the same type) a z in the same residue
2382 : * class as x and the same residue class as y, if it is possible.
2383 : *
2384 : * We also allow (during recursion) two identical objects even if they are
2385 : * not integermod or polymod. For example:
2386 : *
2387 : * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
2388 : * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
2389 : * ? chinese(x, y)
2390 : * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
2391 :
2392 : static GEN
2393 309314 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
2394 : {
2395 309314 : GEN z = gassoc_proto(f,x,NULL);
2396 309307 : if (z == gen_1) retmkintmod(gen_0,gen_1);
2397 309272 : return z;
2398 : }
2399 :
2400 : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
2401 : * call chinese: makes Mod(0,1) a better "neutral" element */
2402 : static GEN
2403 21 : chinese_intpol(GEN x,GEN y)
2404 : {
2405 21 : pari_sp av = avma;
2406 21 : GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
2407 21 : return gerepileupto(av, chinese(z, y));
2408 : }
2409 :
2410 : GEN
2411 49 : chinese1(GEN x) { return gen_chinese(x,chinese); }
2412 :
2413 : GEN
2414 16520 : chinese(GEN x, GEN y)
2415 : {
2416 : pari_sp av;
2417 16520 : long tx = typ(x), ty;
2418 : GEN z,p1,p2,d,u,v;
2419 :
2420 16520 : if (!y) return chinese1(x);
2421 16471 : if (gidentical(x,y)) return gcopy(x);
2422 16464 : ty = typ(y);
2423 16464 : if (tx == ty) switch(tx)
2424 : {
2425 : case t_POLMOD:
2426 : {
2427 28 : GEN A = gel(x,1), B = gel(y,1);
2428 28 : GEN a = gel(x,2), b = gel(y,2);
2429 28 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
2430 28 : if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
2431 28 : av = avma;
2432 28 : d = RgX_extgcd(A,B,&u,&v);
2433 28 : p2 = gsub(b, a);
2434 28 : if (!gequal0(gmod(p2, d))) break;
2435 28 : p1 = gdiv(A,d);
2436 28 : p2 = gadd(a, gmul(gmul(u,p1), p2));
2437 :
2438 28 : z = cgetg(3, t_POLMOD);
2439 28 : gel(z,1) = gmul(p1,B);
2440 28 : gel(z,2) = gmod(p2,gel(z,1));
2441 28 : return gerepileupto(av, z);
2442 : }
2443 : case t_INTMOD:
2444 : {
2445 16401 : GEN A = gel(x,1), B = gel(y,1);
2446 16401 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
2447 16401 : z = cgetg(3,t_INTMOD);
2448 16401 : Z_chinese_pre(A, B, &C, &U, &d);
2449 16401 : c = Z_chinese_post(a, b, C, U, d);
2450 16401 : if (!c) pari_err_OP("chinese", x,y);
2451 16401 : gel(z,1) = icopy_avma(C, (pari_sp)z);
2452 16401 : gel(z,2) = icopy_avma(c, (pari_sp)gel(z,1));
2453 16401 : avma = (pari_sp)gel(z,2); return z;
2454 : }
2455 : case t_POL:
2456 : {
2457 7 : long i, lx = lg(x), ly = lg(y);
2458 7 : if (varn(x) != varn(y)) break;
2459 7 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
2460 7 : z = cgetg(lx, t_POL); z[1] = x[1];
2461 7 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2462 7 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
2463 7 : return z;
2464 : }
2465 :
2466 : case t_VEC: case t_COL: case t_MAT:
2467 : {
2468 : long i, lx;
2469 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
2470 7 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2471 7 : return z;
2472 : }
2473 : }
2474 21 : if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
2475 7 : if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
2476 0 : pari_err_OP("chinese",x,y);
2477 : return NULL; /* LCOV_EXCL_LINE */
2478 : }
2479 :
2480 : /* init chinese(Mod(.,A), Mod(.,B)) */
2481 : void
2482 237218 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
2483 : {
2484 237218 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
2485 237218 : GEN t = diviiexact(A,d);
2486 237218 : *pU = mulii(u, t);
2487 237218 : *pC = mulii(t, B);
2488 237218 : if (pd) *pd = d;
2489 237218 : }
2490 : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
2491 : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
2492 : * If d not NULL, check wether a = b mod d. */
2493 : GEN
2494 637304 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
2495 : {
2496 : GEN b_a;
2497 637304 : if (!signe(a))
2498 : {
2499 312591 : if (d && !dvdii(b, d)) return NULL;
2500 312591 : return Fp_mul(b, U, C);
2501 : }
2502 324713 : b_a = subii(b,a);
2503 324713 : if (d && !dvdii(b_a, d)) return NULL;
2504 324713 : return modii(addii(a, mulii(U, b_a)), C);
2505 : }
2506 : static ulong
2507 1104325 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
2508 : {
2509 1104325 : if (!a) return Fl_mul(b, U, C);
2510 1104325 : return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
2511 : }
2512 :
2513 : GEN
2514 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
2515 : {
2516 2142 : pari_sp av = avma;
2517 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
2518 2142 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
2519 : }
2520 : GEN
2521 218619 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
2522 : {
2523 218619 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
2524 218619 : return Z_chinese_post(a,b, *pC, U, NULL);
2525 : }
2526 :
2527 : /* return lift(chinese(a mod A, b mod B))
2528 : * assume(A,B)=1, a,b,A,B integers and C = A*B */
2529 : GEN
2530 110591 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
2531 : {
2532 110591 : pari_sp av = avma;
2533 110591 : GEN U = mulii(Fp_inv(A,B), A);
2534 110591 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2535 : }
2536 : ulong
2537 1104325 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
2538 1104325 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
2539 :
2540 : /* chinese1 for coprime moduli in Z */
2541 : static GEN
2542 289229 : chinese1_coprime_Z_aux(GEN x, GEN y)
2543 : {
2544 289229 : GEN z = cgetg(3, t_INTMOD);
2545 289229 : GEN A = gel(x,1), a = gel(x, 2);
2546 289229 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
2547 289229 : pari_sp av = avma;
2548 289229 : GEN U = mulii(Fp_inv(A,B), A);
2549 289229 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2550 289229 : gel(z,1) = C; return z;
2551 : }
2552 : GEN
2553 309265 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
2554 :
2555 : /*********************************************************************/
2556 : /** **/
2557 : /** MODULAR EXPONENTIATION **/
2558 : /** **/
2559 : /*********************************************************************/
2560 :
2561 : /* xa, ya = t_VECSMALL */
2562 : GEN
2563 471589 : ZV_producttree(GEN xa)
2564 : {
2565 471589 : long n = lg(xa)-1;
2566 471589 : long m = n==1 ? 1: expu(n-1)+1;
2567 471584 : GEN T = cgetg(m+1, t_VEC), t;
2568 : long i, j, k;
2569 471559 : t = cgetg(((n+1)>>1)+1, t_VEC);
2570 471582 : if (typ(xa)==t_VECSMALL)
2571 : {
2572 952120 : for (j=1, k=1; k<n; j++, k+=2)
2573 553704 : gel(t, j) = muluu(xa[k], xa[k+1]);
2574 398416 : if (k==n) gel(t, j) = utoi(xa[k]);
2575 : } else {
2576 369128 : for (j=1, k=1; k<n; j++, k+=2)
2577 295963 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
2578 73165 : if (k==n) gel(t, j) = icopy(gel(xa,k));
2579 : }
2580 471583 : gel(T,1) = t;
2581 894836 : for (i=2; i<=m; i++)
2582 : {
2583 423252 : GEN u = gel(T, i-1);
2584 423252 : long n = lg(u)-1;
2585 423252 : t = cgetg(((n+1)>>1)+1, t_VEC);
2586 1039428 : for (j=1, k=1; k<n; j++, k+=2)
2587 616175 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
2588 423253 : if (k==n) gel(t, j) = gel(u, k);
2589 423253 : gel(T, i) = t;
2590 : }
2591 471584 : return T;
2592 : }
2593 :
2594 : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
2595 : GEN
2596 10237240 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
2597 : {
2598 : long i,j,k;
2599 10237240 : long m = lg(T)-1, n = lg(P)-1;
2600 : GEN t;
2601 10237240 : GEN Tp = cgetg(m+1, t_VEC);
2602 10236609 : gel(Tp, m) = mkvec(A);
2603 17481256 : for (i=m-1; i>=1; i--)
2604 : {
2605 7245049 : GEN u = gel(T, i);
2606 7245049 : GEN v = gel(Tp, i+1);
2607 7245049 : long n = lg(u)-1;
2608 7245049 : t = cgetg(n+1, t_VEC);
2609 14767967 : for (j=1, k=1; k<n; j++, k+=2)
2610 : {
2611 7522494 : gel(t, k) = modii(gel(v, j), gel(u, k));
2612 7522879 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
2613 : }
2614 7245473 : if (k==n) gel(t, k) = gel(v, j);
2615 7245473 : gel(Tp, i) = t;
2616 : }
2617 : {
2618 10236207 : GEN u = gel(T, i+1);
2619 10236207 : GEN v = gel(Tp, i+1);
2620 10236207 : long l = lg(u)-1;
2621 10236207 : if (typ(P)==t_VECSMALL)
2622 : {
2623 9764622 : GEN R = cgetg(n+1, t_VECSMALL);
2624 26436107 : for (j=1, k=1; j<=l; j++, k+=2)
2625 : {
2626 16670124 : uel(R,k) = umodiu(gel(v, j), P[k]);
2627 16671651 : if (k < n)
2628 11711946 : uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
2629 : }
2630 9765983 : return R;
2631 : }
2632 : else
2633 : {
2634 471585 : GEN R = cgetg(n+1, t_VEC);
2635 1559357 : for (j=1, k=1; j<=l; j++, k+=2)
2636 : {
2637 1087750 : gel(R,k) = modii(gel(v, j), gel(P,k));
2638 1087775 : if (k < n)
2639 849690 : gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
2640 : }
2641 471607 : return R;
2642 : }
2643 : }
2644 : }
2645 :
2646 : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
2647 : GEN
2648 5976828 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
2649 : {
2650 5976828 : long m = lg(T)-1, n = lg(A)-1;
2651 : long i,j,k;
2652 5976828 : GEN Tp = cgetg(m+1, t_VEC);
2653 5954178 : GEN M = gel(T, 1);
2654 5954178 : GEN t = cgetg(lg(M), t_VEC);
2655 5994240 : if (typ(P)==t_VECSMALL)
2656 : {
2657 19895987 : for (j=1, k=1; k<n; j++, k+=2)
2658 : {
2659 17153262 : pari_sp av = avma;
2660 17153262 : GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
2661 16976783 : GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
2662 16957121 : gel(t, j) = gerepileuptoint(av, tj);
2663 : }
2664 2742725 : if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
2665 : } else
2666 : {
2667 7526767 : for (j=1, k=1; k<n; j++, k+=2)
2668 : {
2669 4315270 : pari_sp av = avma;
2670 4315270 : GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
2671 4298136 : GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
2672 4325424 : gel(t, j) = gerepileuptoint(av, tj);
2673 : }
2674 3211497 : if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
2675 : }
2676 5950521 : gel(Tp, 1) = t;
2677 13882894 : for (i=2; i<=m; i++)
2678 : {
2679 7937549 : GEN u = gel(T, i-1), M = gel(T, i);
2680 7937549 : GEN t = cgetg(lg(M), t_VEC);
2681 8021618 : GEN v = gel(Tp, i-1);
2682 8021618 : long n = lg(v)-1;
2683 24870667 : for (j=1, k=1; k<n; j++, k+=2)
2684 : {
2685 16938294 : pari_sp av = avma;
2686 67753176 : gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
2687 50814882 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
2688 : }
2689 7932373 : if (k==n) gel(t, j) = gel(v, k);
2690 7932373 : gel(Tp, i) = t;
2691 : }
2692 5945345 : return gmael(Tp,m,1);
2693 : }
2694 :
2695 : static GEN
2696 219108 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2697 : {
2698 219108 : long i, l = lg(gel(vA,1)), n = lg(P);
2699 219108 : GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
2700 5441036 : for (i=1; i < l; i++)
2701 : {
2702 5221799 : pari_sp av = avma;
2703 5221799 : GEN c, A = cgetg(n, typ(P));
2704 : long j;
2705 5236752 : for (j=1; j < n; j++) A[j] = mael(vA,j,i);
2706 5236752 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
2707 5221500 : gel(V,i) = gerepileuptoint(av, c);
2708 : }
2709 219237 : return V;
2710 : }
2711 :
2712 : static GEN
2713 89763 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2714 : {
2715 89763 : long i, j, l, n = lg(P);
2716 89763 : GEN mod = gmael(T, lg(T)-1, 1), V, w;
2717 89763 : w = cgetg(n, t_VECSMALL);
2718 89762 : for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
2719 89762 : l = vecsmall_max(w);
2720 89800 : V = cgetg(l, t_POL);
2721 89973 : V[1] = mael(vA,1,1);
2722 389601 : for (i=2; i < l; i++)
2723 : {
2724 299870 : pari_sp av = avma;
2725 299870 : GEN c, A = cgetg(n, typ(P));
2726 299301 : if (typ(P)==t_VECSMALL)
2727 5327 : for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
2728 : else
2729 293974 : for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
2730 299301 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
2731 300026 : gel(V,i) = gerepileuptoint(av, c);
2732 : }
2733 89731 : return ZX_renormalize(V, l);
2734 : }
2735 :
2736 : static GEN
2737 4428 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2738 : {
2739 4428 : long i, j, l = lg(gel(vA,1)), n = lg(P);
2740 4428 : GEN A = cgetg(n, t_VEC);
2741 4426 : GEN V = cgetg(l, t_COL);
2742 91726 : for (i=1; i < l; i++)
2743 : {
2744 87293 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2745 87293 : gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
2746 : }
2747 4433 : return V;
2748 : }
2749 :
2750 : static GEN
2751 18563 : polint_chinese(GEN worker, GEN mA, GEN P)
2752 : {
2753 18563 : long i, j, l = lg(gel(mA,1)), n = lg(P);
2754 18563 : long pending = 0, workid, cnt = 0;
2755 : struct pari_mt pt;
2756 : GEN done, va, M;
2757 18563 : GEN A = cgetg(n, t_VEC);
2758 18563 : va = mkvec(gen_0);
2759 18563 : M = cgetg(l, t_MAT);
2760 18563 : if (DEBUGLEVEL>2) err_printf("Start parallel Chinese remainder: ");
2761 18563 : mt_queue_start_lim(&pt, worker, l-1);
2762 248426 : for (i=1; i<l || pending; i++)
2763 : {
2764 229863 : for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
2765 229863 : gel(va, 1) = A;
2766 229863 : mt_queue_submit(&pt, i, i<l? va: NULL);
2767 229863 : done = mt_queue_get(&pt, &workid, &pending);
2768 229863 : if (done)
2769 : {
2770 211409 : gel(M,workid) = done;
2771 211409 : if (DEBUGLEVEL>2) err_printf("%ld%% ",(++cnt)*100/(l-1));
2772 : }
2773 : }
2774 18563 : if (DEBUGLEVEL>2) err_printf("\n");
2775 18563 : mt_queue_end(&pt);
2776 18563 : return M;
2777 : }
2778 :
2779 : GEN
2780 4113 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
2781 : {
2782 4113 : return nxCV_polint_center_tree(vA, P, T, R, m2);
2783 : }
2784 :
2785 : static GEN
2786 66 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2787 : {
2788 66 : long i, j, l = lg(gel(vA,1)), n = lg(P);
2789 66 : GEN A = cgetg(n, t_VEC);
2790 66 : GEN V = cgetg(l, t_MAT);
2791 381 : for (i=1; i < l; i++)
2792 : {
2793 315 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2794 315 : gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
2795 : }
2796 66 : return V;
2797 : }
2798 :
2799 : static GEN
2800 284 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
2801 : {
2802 284 : GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
2803 284 : return polint_chinese(worker, mA, P);
2804 : }
2805 :
2806 : static GEN
2807 757 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2808 : {
2809 757 : long i, j, l = lg(gel(vA,1)), n = lg(P);
2810 757 : GEN A = cgetg(n, t_VEC);
2811 757 : GEN V = cgetg(l, t_MAT);
2812 9197 : for (i=1; i < l; i++)
2813 : {
2814 8440 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
2815 8440 : gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
2816 : }
2817 757 : return V;
2818 : }
2819 :
2820 : GEN
2821 206945 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
2822 : {
2823 206945 : return ncV_polint_center_tree(vA, P, T, R, m2);
2824 : }
2825 :
2826 : static GEN
2827 18279 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
2828 : {
2829 18279 : GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
2830 18279 : return polint_chinese(worker, mA, P);
2831 : }
2832 :
2833 : /* return [A mod P[i], i=1..#P] */
2834 : GEN
2835 0 : Z_ZV_mod(GEN A, GEN P)
2836 : {
2837 0 : pari_sp av = avma;
2838 0 : return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
2839 : }
2840 : /* P a t_VECSMALL */
2841 : GEN
2842 0 : Z_nv_mod(GEN A, GEN P)
2843 : {
2844 0 : pari_sp av = avma;
2845 0 : return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
2846 : }
2847 : /* B a ZX, T = ZV_producttree(P) */
2848 : GEN
2849 390345 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
2850 : {
2851 : pari_sp av;
2852 390345 : long i, j, l = lg(B), n = lg(A)-1;
2853 390345 : GEN V = cgetg(n+1, t_VEC);
2854 1558283 : for (j=1; j <= n; j++)
2855 : {
2856 1167860 : gel(V, j) = cgetg(l, t_VECSMALL);
2857 1167843 : mael(V, j, 1) = B[1]&VARNBITS;
2858 : }
2859 390423 : av = avma;
2860 9349753 : for (i=2; i < l; i++)
2861 : {
2862 8959339 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
2863 35246610 : for (j=1; j <= n; j++)
2864 26287280 : mael(V, j, i) = v[j];
2865 8959330 : avma = av;
2866 : }
2867 1558182 : for (j=1; j <= n; j++)
2868 1167863 : (void) Flx_renormalize(gel(V, j), l);
2869 390319 : return V;
2870 : }
2871 :
2872 : static GEN
2873 3159 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
2874 :
2875 : GEN
2876 405 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
2877 : {
2878 405 : pari_sp av = avma;
2879 405 : long i, j, l = lg(P), n = lg(xa)-1, vP = varn(P);
2880 405 : GEN V = cgetg(n+1, t_VEC);
2881 1435 : for (j=1; j <= n; j++)
2882 : {
2883 1030 : gel(V, j) = cgetg(l, t_POL);
2884 1030 : mael(V, j, 1) = vP;
2885 : }
2886 1881 : for (i=2; i < l; i++)
2887 : {
2888 1476 : GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
2889 5091 : for (j=1; j <= n; j++)
2890 3615 : gmael(V, j, i) = gel(v,j);
2891 : }
2892 1435 : for (j=1; j <= n; j++)
2893 1030 : (void) FlxX_renormalize(gel(V, j), l);
2894 405 : return gerepilecopy(av, V);
2895 : }
2896 :
2897 : GEN
2898 315 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
2899 : {
2900 315 : pari_sp av = avma;
2901 315 : long i, j, l = lg(C), n = lg(xa)-1;
2902 315 : GEN V = cgetg(n+1, t_VEC);
2903 1073 : for (j = 1; j <= n; j++)
2904 758 : gel(V, j) = cgetg(l, t_COL);
2905 1998 : for (i = 1; i < l; i++)
2906 : {
2907 1683 : GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
2908 5773 : for (j = 1; j <= n; j++)
2909 4090 : gmael(V, j, i) = gel(v,j);
2910 : }
2911 315 : return gerepilecopy(av, V);
2912 : }
2913 :
2914 : GEN
2915 66 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
2916 : {
2917 66 : pari_sp av = avma;
2918 66 : long i, j, l = lg(M), n = lg(xa)-1;
2919 66 : GEN V = cgetg(n+1, t_VEC);
2920 222 : for (j=1; j <= n; j++)
2921 156 : gel(V, j) = cgetg(l, t_MAT);
2922 381 : for (i=1; i < l; i++)
2923 : {
2924 315 : GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
2925 1073 : for (j=1; j <= n; j++)
2926 758 : gmael(V, j, i) = gel(v,j);
2927 : }
2928 66 : return gerepilecopy(av, V);
2929 : }
2930 :
2931 : /* B a ZX, T = ZV_producttree(P) */
2932 : GEN
2933 96244 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
2934 : {
2935 : pari_sp av;
2936 96244 : long i, j, l = lg(B), n = lg(A)-1;
2937 96244 : GEN V = cgetg(n+1, t_VEC);
2938 344019 : for (j=1; j <= n; j++)
2939 247723 : gel(V, j) = cgetg(l, t_VECSMALL);
2940 96296 : av = avma;
2941 902679 : for (i=1; i < l; i++)
2942 : {
2943 806438 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
2944 2903044 : for (j=1; j <= n; j++)
2945 2096661 : mael(V, j, i) = v[j];
2946 806383 : avma = av;
2947 : }
2948 96241 : return V;
2949 : }
2950 :
2951 : GEN
2952 12017 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
2953 : {
2954 12017 : pari_sp av = avma;
2955 12017 : long i, j, l = lg(M), n = lg(xa)-1;
2956 12017 : GEN V = cgetg(n+1, t_VEC);
2957 42605 : for (j=1; j <= n; j++)
2958 30591 : gel(V, j) = cgetg(l, t_MAT);
2959 108256 : for (i=1; i < l; i++)
2960 : {
2961 96239 : GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
2962 344036 : for (j=1; j <= n; j++)
2963 247794 : gmael(V, j, i) = gel(v,j);
2964 : }
2965 12017 : return gerepilecopy(av, V);
2966 : }
2967 :
2968 : static GEN
2969 468190 : ZV_sqr(GEN z)
2970 : {
2971 468190 : long i,l = lg(z);
2972 468190 : GEN x = cgetg(l, t_VEC);
2973 468208 : if (typ(z)==t_VECSMALL)
2974 398452 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
2975 : else
2976 69756 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
2977 468179 : return x;
2978 : }
2979 :
2980 : static GEN
2981 3130408 : ZT_sqr(GEN z)
2982 : {
2983 3130408 : if (typ(z) == t_INT)
2984 1776292 : return sqri(z);
2985 : else
2986 : {
2987 1354116 : long i,l = lg(z);
2988 1354116 : GEN x = cgetg(l, t_VEC);
2989 1354148 : for (i=1; i<l; i++) gel(x,i) = ZT_sqr(gel(z,i));
2990 1354264 : return x;
2991 : }
2992 : }
2993 :
2994 : static GEN
2995 468184 : ZV_invdivexact(GEN y, GEN x)
2996 : {
2997 468184 : long i, l = lg(y);
2998 468184 : GEN z = cgetg(l,t_VEC);
2999 468205 : if (typ(x)==t_VECSMALL)
3000 1713428 : for (i=1; i<l; i++)
3001 : {
3002 1314987 : pari_sp av = avma;
3003 1314987 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
3004 1315004 : avma = av;
3005 1315004 : gel(z,i) = utoipos(a);
3006 : }
3007 : else
3008 668119 : for (i=1; i<l; i++)
3009 598363 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
3010 468197 : return z;
3011 : }
3012 :
3013 : /* P t_VECSMALL or t_VEC of t_INT */
3014 : GEN
3015 468118 : ZV_chinesetree(GEN P, GEN T)
3016 : {
3017 468118 : GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
3018 468174 : GEN mod = gmael(T,lg(T)-1,1);
3019 468174 : return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
3020 : }
3021 :
3022 : static GEN
3023 76033 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
3024 : {
3025 76033 : if (!pt_mod)
3026 2545 : return gerepileupto(av, a);
3027 : else
3028 : {
3029 73488 : GEN mod = gmael(T, lg(T)-1, 1);
3030 73488 : gerepileall(av, 2, &a, &mod);
3031 73488 : *pt_mod = mod;
3032 73488 : return a;
3033 : }
3034 : }
3035 :
3036 : GEN
3037 38813 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3038 : {
3039 38813 : pari_sp av = avma;
3040 38813 : GEN T = ZV_producttree(P);
3041 38813 : GEN R = ZV_chinesetree(P, T);
3042 38813 : GEN a = ZV_chinese_tree(A, P, T, R);
3043 38813 : GEN mod = gmael(T, lg(T)-1, 1);
3044 38813 : GEN ca = Fp_center(a, mod, shifti(mod,-1));
3045 38813 : return gc_chinese(av, T, ca, pt_mod);
3046 : }
3047 :
3048 : GEN
3049 12931 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
3050 : {
3051 12931 : pari_sp av = avma;
3052 12931 : GEN T = ZV_producttree(P);
3053 12931 : GEN R = ZV_chinesetree(P, T);
3054 12931 : GEN a = ZV_chinese_tree(A, P, T, R);
3055 12931 : return gc_chinese(av, T, a, pt_mod);
3056 : }
3057 :
3058 : GEN
3059 405 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
3060 : {
3061 405 : pari_sp av = avma;
3062 405 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3063 405 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
3064 405 : return gerepileupto(av, a);
3065 : }
3066 :
3067 : GEN
3068 1994 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3069 : {
3070 1994 : pari_sp av = avma;
3071 1994 : GEN T = ZV_producttree(P);
3072 1994 : GEN R = ZV_chinesetree(P, T);
3073 1994 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3074 1994 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
3075 1994 : return gc_chinese(av, T, a, pt_mod);
3076 : }
3077 :
3078 : GEN
3079 3732 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3080 : {
3081 3732 : pari_sp av = avma;
3082 3732 : GEN T = ZV_producttree(P);
3083 3732 : GEN R = ZV_chinesetree(P, T);
3084 3732 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3085 3732 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
3086 3732 : return gc_chinese(av, T, a, pt_mod);
3087 : }
3088 :
3089 : GEN
3090 0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
3091 : {
3092 0 : pari_sp av = avma;
3093 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3094 0 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
3095 0 : return gerepileupto(av, a);
3096 : }
3097 :
3098 : GEN
3099 757 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
3100 : {
3101 757 : pari_sp av = avma;
3102 757 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3103 757 : GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
3104 757 : return gerepileupto(av, a);
3105 : }
3106 :
3107 : GEN
3108 18279 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3109 : {
3110 18279 : pari_sp av = avma;
3111 18279 : GEN T = ZV_producttree(P);
3112 18279 : GEN R = ZV_chinesetree(P, T);
3113 18279 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3114 18279 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
3115 18279 : return gc_chinese(av, T, a, pt_mod);
3116 : }
3117 :
3118 : GEN
3119 66 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
3120 : {
3121 66 : pari_sp av = avma;
3122 66 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3123 66 : GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
3124 66 : return gerepileupto(av, a);
3125 : }
3126 :
3127 : GEN
3128 284 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
3129 : {
3130 284 : pari_sp av = avma;
3131 284 : GEN T = ZV_producttree(P);
3132 284 : GEN R = ZV_chinesetree(P, T);
3133 284 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
3134 284 : GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
3135 284 : return gc_chinese(av, T, a, pt_mod);
3136 : }
3137 :
3138 : /**********************************************************************
3139 : ** **
3140 : ** Powering over (Z/NZ)^*, small N **
3141 : ** **
3142 : **********************************************************************/
3143 :
3144 : /* 2^n mod p; assume n > 1 */
3145 : static ulong
3146 20078323 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
3147 : {
3148 20078323 : ulong y = 2;
3149 20078323 : int j = 1+bfffo(n);
3150 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
3151 20078323 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
3152 415887059 : for (; j; n<<=1,j--)
3153 : {
3154 395808820 : y = Fl_sqr_pre(y,p,pi);
3155 395809002 : if (n & HIGHBIT) y = Fl_double(y, p);
3156 : }
3157 20078239 : return y;
3158 : }
3159 :
3160 : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
3161 : static ulong
3162 4661649 : Fl_2powu(ulong n, ulong p)
3163 : {
3164 4661649 : ulong y = 2;
3165 4661649 : int j = 1+bfffo(n);
3166 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
3167 4661649 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
3168 59990178 : for (; j; n<<=1,j--)
3169 : {
3170 55343632 : y = (y*y) % p;
3171 55343632 : if (n & HIGHBIT) y = Fl_double(y, p);
3172 : }
3173 4646546 : return y;
3174 : }
3175 :
3176 : ulong
3177 96228926 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
3178 : {
3179 : ulong y, z, n;
3180 96228926 : if (n0 <= 1)
3181 : { /* frequent special cases */
3182 12483201 : if (n0 == 1) return x;
3183 4714960 : if (n0 == 0) return 1;
3184 : }
3185 83745725 : if (x <= 2)
3186 : {
3187 21294015 : if (x == 2) return Fl_2powu_pre(n0, p, pi);
3188 1215808 : return x; /* 0 or 1 */
3189 : }
3190 62451710 : y = 1; z = x; n = n0;
3191 : for(;;)
3192 : {
3193 562364265 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
3194 562786443 : n>>=1; if (!n) return y;
3195 500349695 : z = Fl_sqr_pre(z,p,pi);
3196 499912555 : }
3197 : }
3198 :
3199 : ulong
3200 45730366 : Fl_powu(ulong x, ulong n0, ulong p)
3201 : {
3202 : ulong y, z, n;
3203 45730366 : if (n0 <= 2)
3204 : { /* frequent special cases */
3205 32112247 : if (n0 == 2) return Fl_sqr(x,p);
3206 3471277 : if (n0 == 1) return x;
3207 43735 : if (n0 == 0) return 1;
3208 : }
3209 13618119 : if (x <= 1) return x; /* 0 or 1 */
3210 13568503 : if (p & HIGHMASK)
3211 6216497 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
3212 7352006 : if (x == 2) return Fl_2powu(n0, p);
3213 2715266 : y = 1; z = x; n = n0;
3214 : for(;;)
3215 : {
3216 25870429 : if (n&1) y = (y*z) % p;
3217 25870429 : n>>=1; if (!n) return y;
3218 23155163 : z = (z*z) % p;
3219 23155163 : }
3220 : }
3221 :
3222 : /* Reduce data dependency to maximize internal parallelism */
3223 : GEN
3224 11040246 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
3225 : {
3226 : long i, k;
3227 11040246 : GEN powers = cgetg(n + 2, t_VECSMALL);
3228 11045031 : powers[1] = 1; if (n == 0) return powers;
3229 11045031 : powers[2] = x;
3230 46642056 : for (i = 3, k=2; i <= n; i+=2, k++)
3231 : {
3232 35602627 : powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
3233 35600458 : powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
3234 : }
3235 11039429 : if (i==n+1)
3236 9645339 : powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
3237 11039759 : return powers;
3238 : }
3239 :
3240 : GEN
3241 2821 : Fl_powers(ulong x, long n, ulong p)
3242 : {
3243 2821 : return Fl_powers_pre(x, n, p, get_Fl_red(p));
3244 : }
3245 :
3246 : /**********************************************************************
3247 : ** **
3248 : ** Powering over (Z/NZ)^*, large N **
3249 : ** **
3250 : **********************************************************************/
3251 :
3252 : static GEN
3253 4258758 : Fp_dblsqr(GEN x, GEN N)
3254 : {
3255 4258758 : GEN z = shifti(Fp_sqr(x, N), 1);
3256 4258758 : return cmpii(z, N) >= 0? subii(z, N): z;
3257 : }
3258 :
3259 : typedef struct muldata {
3260 : GEN (*sqr)(void * E, GEN x);
3261 : GEN (*mul)(void * E, GEN x, GEN y);
3262 : GEN (*mul2)(void * E, GEN x);
3263 : } muldata;
3264 :
3265 : /* modified Barrett reduction with one fold */
3266 : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
3267 :
3268 : static GEN
3269 10090 : Fp_invmBarrett(GEN p, long s)
3270 : {
3271 10090 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
3272 10090 : return mkvec2(Q,R);
3273 : }
3274 :
3275 : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
3276 : * a = r (mod N) */
3277 : static GEN
3278 4234830 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
3279 : {
3280 4234830 : pari_sp av = avma;
3281 4234830 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
3282 4234830 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
3283 4234830 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
3284 4234830 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
3285 4234830 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
3286 4234830 : GEN r = subii(A, mulii(q, N));
3287 4234830 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
3288 4234830 : if (signe(sr)<0) return gerepileuptoint(av, r);
3289 2548273 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
3290 2548273 : if (signe(sr)<0) return gerepileuptoint(av, r);
3291 98163 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
3292 98163 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
3293 : }
3294 :
3295 : /* Montgomery reduction */
3296 :
3297 : INLINE ulong
3298 328971 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
3299 :
3300 : struct montred
3301 : {
3302 : GEN N;
3303 : ulong inv;
3304 : };
3305 :
3306 : /* Montgomery reduction */
3307 : static GEN
3308 32522218 : _sqr_montred(void * E, GEN x)
3309 : {
3310 32522218 : struct montred * D = (struct montred *) E;
3311 32522218 : return red_montgomery(sqri(x), D->N, D->inv);
3312 : }
3313 :
3314 : /* Montgomery reduction */
3315 : static GEN
3316 3045670 : _mul_montred(void * E, GEN x, GEN y)
3317 : {
3318 3045670 : struct montred * D = (struct montred *) E;
3319 3045670 : return red_montgomery(mulii(x, y), D->N, D->inv);
3320 : }
3321 :
3322 : static GEN
3323 5837186 : _mul2_montred(void * E, GEN x)
3324 : {
3325 5837186 : struct montred * D = (struct montred *) E;
3326 5837186 : GEN z = shifti(_sqr_montred(E, x), 1);
3327 5837196 : long l = lgefint(D->N);
3328 5837196 : while (lgefint(z) > l) z = subii(z, D->N);
3329 5837196 : return z;
3330 : }
3331 :
3332 : static GEN
3333 10891387 : _sqr_remii(void* N, GEN x)
3334 10891387 : { return remii(sqri(x), (GEN) N); }
3335 :
3336 : static GEN
3337 1046021 : _mul_remii(void* N, GEN x, GEN y)
3338 1046021 : { return remii(mulii(x, y), (GEN) N); }
3339 :
3340 : static GEN
3341 2997200 : _mul2_remii(void* N, GEN x)
3342 2997200 : { return Fp_dblsqr(x, (GEN) N); }
3343 :
3344 : struct redbarrett
3345 : {
3346 : GEN iM, N;
3347 : long s;
3348 : };
3349 :
3350 : static GEN
3351 3823824 : _sqr_remiibar(void *E, GEN x)
3352 : {
3353 3823824 : struct redbarrett * D = (struct redbarrett *) E;
3354 3823824 : return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
3355 : }
3356 :
3357 : static GEN
3358 411006 : _mul_remiibar(void *E, GEN x, GEN y)
3359 : {
3360 411006 : struct redbarrett * D = (struct redbarrett *) E;
3361 411006 : return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
3362 : }
3363 :
3364 : static GEN
3365 1261558 : _mul2_remiibar(void *E, GEN x)
3366 : {
3367 1261558 : struct redbarrett * D = (struct redbarrett *) E;
3368 1261558 : return Fp_dblsqr(x, D->N);
3369 : }
3370 :
3371 : static long
3372 422257 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
3373 : {
3374 422257 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
3375 : {
3376 10090 : struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
3377 10090 : D->sqr = &_sqr_remiibar;
3378 10090 : D->mul = &_mul_remiibar;
3379 10090 : D->mul2 = &_mul2_remiibar;
3380 10090 : E->N = N;
3381 10090 : E->s = 1+(expi(N)>>1);
3382 10090 : E->iM = Fp_invmBarrett(N, E->s);
3383 10090 : *pt_E = (void*) E;
3384 10090 : return 0;
3385 : }
3386 412167 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
3387 : {
3388 328973 : struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
3389 328972 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
3390 328973 : D->sqr = &_sqr_montred;
3391 328973 : D->mul = &_mul_montred;
3392 328973 : D->mul2 = &_mul2_montred;
3393 328973 : E->N = N;
3394 328973 : E->inv = init_montdata(N);
3395 328970 : *pt_E = (void*) E;
3396 328970 : return 1;
3397 : }
3398 : else
3399 : {
3400 83195 : D->sqr = &_sqr_remii;
3401 83195 : D->mul = &_mul_remii;
3402 83195 : D->mul2 = &_mul2_remii;
3403 83195 : *pt_E = (void*) N;
3404 83195 : return 0;
3405 : }
3406 : }
3407 :
3408 : GEN
3409 1002948 : Fp_powu(GEN A, ulong k, GEN N)
3410 : {
3411 1002948 : long lN = lgefint(N), sA;
3412 : int base_is_2, use_montgomery;
3413 : muldata D;
3414 : void *E;
3415 : pari_sp av;
3416 :
3417 1002948 : if (lN == 3) {
3418 87683 : ulong n = uel(N,2);
3419 87683 : return utoi( Fl_powu(umodiu(A, n), k, n) );
3420 : }
3421 915265 : if (k <= 2)
3422 : { /* frequent special cases */
3423 630342 : if (k == 2) return Fp_sqr(A,N);
3424 175912 : if (k == 1) return A;
3425 0 : if (k == 0) return gen_1;
3426 : }
3427 284923 : sA = signe(A)==-1 && odd(k);
3428 284923 : base_is_2 = 0;
3429 284923 : if (lgefint(A) == 3) switch(A[2])
3430 : {
3431 523 : case 1: return sA ? gen_m1 : gen_1;
3432 34408 : case 2: base_is_2 = 1; break;
3433 : }
3434 :
3435 : /* TODO: Move this out of here and use for general modular computations */
3436 284400 : av = avma;
3437 284400 : use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
3438 284399 : if (base_is_2)
3439 34408 : A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
3440 : else
3441 249991 : A = gen_powu_i(A, k, E, D.sqr, D.mul);
3442 284398 : if (use_montgomery)
3443 : {
3444 253965 : A = red_montgomery(A, N, ((struct montred *) E)->inv);
3445 253966 : if (cmpii(A, N) >= 0) A = subii(A,N);
3446 253966 : if (sA) A = subii(N, A);
3447 : }
3448 284399 : return gerepileuptoint(av, A);
3449 : }
3450 :
3451 : GEN
3452 21567 : Fp_pows(GEN A, long k, GEN N)
3453 : {
3454 21567 : if (lgefint(N) == 3) {
3455 7316 : ulong n = N[2];
3456 7316 : ulong a = umodiu(A, n);
3457 7316 : if (k < 0) {
3458 126 : a = Fl_inv(a, n);
3459 126 : k = -k;
3460 : }
3461 7316 : return utoi( Fl_powu(a, (ulong)k, n) );
3462 : }
3463 14251 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
3464 14251 : return Fp_powu(A, (ulong)k, N);
3465 : }
3466 :
3467 : /* A^K mod N */
3468 : GEN
3469 5333665 : Fp_pow(GEN A, GEN K, GEN N)
3470 : {
3471 : pari_sp av;
3472 5333665 : long s, lN = lgefint(N), sA;
3473 : int base_is_2, use_montgomery;
3474 : GEN y;
3475 : muldata D;
3476 : void *E;
3477 :
3478 5333665 : s = signe(K);
3479 5333665 : if (!s) return dvdii(A,N)? gen_0: gen_1;
3480 5261495 : if (lN == 3 && lgefint(K) == 3)
3481 : {
3482 4968761 : ulong n = N[2], a = umodiu(A, n);
3483 4968761 : if (s < 0) a = Fl_inv(a, n);
3484 4968761 : if (a <= 1) return utoi(a); /* 0 or 1 */
3485 4643707 : return utoi(Fl_powu(a, uel(K,2), n));
3486 : }
3487 :
3488 292734 : av = avma;
3489 292734 : if (s < 0) y = Fp_inv(A,N);
3490 : else
3491 : {
3492 292270 : y = modii(A,N);
3493 292268 : if (!signe(y)) { avma = av; return gen_0; }
3494 : }
3495 292733 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
3496 :
3497 137939 : base_is_2 = 0;
3498 137939 : sA = signe(y)==-1 && mod2(K);
3499 137939 : if (lgefint(y) == 3) switch(y[2])
3500 : {
3501 82 : case 1: return sA ? gen_m1 : gen_1;
3502 88525 : case 2: base_is_2 = 1; break;
3503 : }
3504 :
3505 : /* TODO: Move this out of here and use for general modular computations */
3506 137857 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
3507 137856 : if (base_is_2)
3508 88525 : y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
3509 : else
3510 49331 : y = gen_pow_i(y, K, E, D.sqr, D.mul);
3511 137859 : if (use_montgomery)
3512 : {
3513 75007 : y = red_montgomery(y, N, ((struct montred *) E)->inv);
3514 75009 : if (cmpii(y,N) >= 0) y = subii(y,N);
3515 75006 : if (sA) y = subii(N, y);
3516 : }
3517 137858 : return gerepileuptoint(av,y);
3518 : }
3519 :
3520 : static GEN
3521 639712 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
3522 :
3523 : static GEN
3524 1128 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
3525 :
3526 : static GEN
3527 140 : _Fp_one(void *E) { (void) E; return gen_1; }
3528 :
3529 : GEN
3530 2016 : Fp_powers(GEN x, long n, GEN p)
3531 : {
3532 2016 : if (lgefint(p) == 3)
3533 1876 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
3534 140 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
3535 : }
3536 :
3537 : static GEN
3538 2981073 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
3539 :
3540 : static GEN
3541 112 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
3542 :
3543 : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
3544 :
3545 : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
3546 : equalii,equali1,Fp_easylog};
3547 :
3548 : static GEN
3549 1106416 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
3550 :
3551 : static GEN
3552 1572954 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
3553 :
3554 : static GEN
3555 318458 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
3556 :
3557 : static GEN
3558 1746877 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
3559 :
3560 : static GEN
3561 29472 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
3562 :
3563 : static int
3564 395272 : _Fp_equal0(GEN x) { return signe(x)==0; }
3565 :
3566 : static GEN
3567 151001 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
3568 :
3569 : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
3570 : _Fp_inv,_Fp_equal0,_Fp_s};
3571 :
3572 7619 : const struct bb_field *get_Fp_field(void **E, GEN p)
3573 : {
3574 7619 : *E = (void*)p; return &Fp_field;
3575 : }
3576 :
3577 : /*********************************************************************/
3578 : /** **/
3579 : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
3580 : /** **/
3581 : /*********************************************************************/
3582 : ulong
3583 12082 : Fl_order(ulong a, ulong o, ulong p)
3584 : {
3585 12082 : pari_sp av = avma;
3586 : GEN m, P, E;
3587 : long i;
3588 12082 : if (!o) o = p-1;
3589 12082 : m = factoru(o);
3590 12082 : P = gel(m,1);
3591 12082 : E = gel(m,2);
3592 28910 : for (i = lg(P)-1; i; i--)
3593 : {
3594 16828 : ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
3595 16828 : if (y == 1) o = t;
3596 14840 : else for (j = 1; j < e; j++)
3597 : {
3598 4284 : y = Fl_powu(y, l, p);
3599 4284 : if (y == 1) { o = t * upowuu(l, j); break; }
3600 : }
3601 : }
3602 12082 : avma = av; return o;
3603 : }
3604 :
3605 : /*Find the exact order of a assuming a^o==1*/
3606 : GEN
3607 10676 : Fp_order(GEN a, GEN o, GEN p) {
3608 10676 : if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
3609 : {
3610 21 : ulong pp = p[2], oo = (o && lgefint(o)==3)? o[2]: pp-1;
3611 21 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
3612 : }
3613 10655 : return gen_order(a, o, (void*)p, &Fp_star);
3614 : }
3615 : GEN
3616 56 : Fp_factored_order(GEN a, GEN o, GEN p)
3617 56 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
3618 :
3619 : /* return order of a mod p^e, e > 0, pe = p^e */
3620 : static GEN
3621 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
3622 : {
3623 : GEN ap, op;
3624 70 : if (absequaliu(p, 2))
3625 : {
3626 56 : if (e == 1) return gen_1;
3627 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
3628 49 : if (mod4(a) == 1)
3629 14 : op = gen_1;
3630 : else {
3631 35 : op = gen_2;
3632 35 : a = Fp_sqr(a, pe);
3633 : }
3634 : } else {
3635 14 : ap = (e == 1)? a: remii(a,p);
3636 14 : op = Fp_order(ap, subiu(p,1), p);
3637 14 : if (e == 1) return op;
3638 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
3639 : }
3640 49 : if (equali1(a)) return op;
3641 7 : return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
3642 : }
3643 :
3644 : GEN
3645 63 : znorder(GEN x, GEN o)
3646 : {
3647 63 : pari_sp av = avma;
3648 : GEN b, a;
3649 :
3650 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
3651 56 : b = gel(x,1); a = gel(x,2);
3652 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
3653 49 : if (!o)
3654 : {
3655 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
3656 35 : long i, l = lg(P);
3657 35 : o = gen_1;
3658 70 : for (i = 1; i < l; i++)
3659 : {
3660 35 : GEN p = gel(P,i);
3661 35 : long e = itos(gel(E,i));
3662 :
3663 35 : if (l == 2)
3664 35 : o = Zp_order(a, p, e, b);
3665 : else {
3666 0 : GEN pe = powiu(p,e);
3667 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
3668 : }
3669 : }
3670 35 : return gerepileuptoint(av, o);
3671 : }
3672 14 : return Fp_order(a, o, b);
3673 : }
3674 : GEN
3675 0 : order(GEN x) { return znorder(x, NULL); }
3676 :
3677 : /*********************************************************************/
3678 : /** **/
3679 : /** DISCRETE LOGARITHM in (Z/nZ)* **/
3680 : /** **/
3681 : /*********************************************************************/
3682 : static GEN
3683 61765 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
3684 : {
3685 61765 : pari_sp av = avma;
3686 : GEN h1, h2, F, G;
3687 61765 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return NULL;
3688 37074 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
3689 : {
3690 310 : GEN M = cgetg(3, t_MAT);
3691 310 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
3692 310 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
3693 310 : return gerepileupto(av, M);
3694 : }
3695 36764 : avma = av; return NULL;
3696 : }
3697 :
3698 : static GEN
3699 61765 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
3700 : {
3701 : GEN rel;
3702 : do
3703 : {
3704 61765 : (*e)++; *g = Fp_mul(*g, b, p);
3705 61765 : rel = Fp_log_halfgcd(bnd, C, *g, p);
3706 61765 : } while (!rel);
3707 310 : return rel;
3708 : }
3709 :
3710 : struct Fp_log_rel
3711 : {
3712 : GEN rel;
3713 : ulong prmax;
3714 : long nbrel, nbmax, nbgen;
3715 : };
3716 :
3717 : /* add u^e */
3718 : static void
3719 2583 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
3720 : {
3721 2583 : pari_sp av = avma;
3722 2583 : long off = r->prmax+1;
3723 2583 : GEN F = cgetg(3, t_MAT);
3724 2583 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
3725 2583 : gel(F,2) = vecsmall_append(gel(z,2), e);
3726 2583 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3727 2583 : }
3728 :
3729 : /* add u^-1 v^-1 */
3730 : static void
3731 99869 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
3732 : {
3733 99869 : pari_sp av = avma;
3734 99869 : long off = r->prmax+1;
3735 99869 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
3736 99869 : GEN F = cgetg(3, t_MAT);
3737 99869 : gel(F,1) = vecsmall_concat(gel(z,1), P);
3738 99869 : gel(F,2) = vecsmall_concat(gel(z,2), E);
3739 99869 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3740 99869 : }
3741 :
3742 : /*
3743 : Let p=C^2+c
3744 : Solve h = (C+x)*(C+a)-p = 0 [mod l]
3745 : h= -c+x*(C+a)+C*a = 0 [mod l]
3746 : x = (c-C*a)/(C+a) [mod l]
3747 : h = -c+C*(x+a)+a*x
3748 : */
3749 :
3750 : GEN
3751 38996 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
3752 : {
3753 38996 : pari_sp ltop = avma;
3754 38996 : long th = expi(mulis(C,a)), n = lg(pi)-1;
3755 : long i, j;
3756 38987 : GEN sieve = zero_zv(a+2)+1;
3757 39011 : GEN L = cgetg(1+a+2, t_VEC);
3758 39007 : pari_sp av = avma;
3759 39007 : long rel = 1;
3760 : GEN z, h;
3761 39007 : h = addis(C,a);
3762 38990 : if ((z = Z_issmooth_fact(h, prmax)))
3763 : {
3764 2415 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
3765 2414 : av = avma;
3766 : }
3767 17568477 : for(i=1; i<=n; i++)
3768 : {
3769 17529463 : ulong li = pi[i], s = sz[i], al = a % li;
3770 17529463 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
3771 18011890 : if (!iv) continue;
3772 17393168 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
3773 79505017 : for(j = u; j<=a; j+=li)
3774 62519416 : sieve[j] += s;
3775 : }
3776 39014 : th = th - expu(th)-1;
3777 27843694 : for(j=0; j<a; j++)
3778 27804656 : if (sieve[j]>=th)
3779 : {
3780 114221 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
3781 111601 : if ((z = Z_issmooth_fact(h, prmax)))
3782 : {
3783 104979 : gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
3784 105019 : av = avma;
3785 6657 : } else avma = av;
3786 : }
3787 : /* j = a */
3788 39038 : if (sieve[a]>=th)
3789 : {
3790 420 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
3791 420 : if ((z = Z_issmooth_fact(h, prmax)))
3792 : {
3793 343 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
3794 343 : av = avma;
3795 : }
3796 : }
3797 39038 : setlg(L, rel);
3798 39038 : return gerepilecopy(ltop, L);
3799 : }
3800 :
3801 : static long
3802 49 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
3803 : {
3804 : struct pari_mt pt;
3805 49 : long i, j, nb = 0;
3806 49 : GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
3807 : mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
3808 49 : long running, pending = 0;
3809 49 : GEN W = zerovec(r->nbgen);
3810 49 : mt_queue_start_lim(&pt, worker, r->nbgen);
3811 39193 : for (i = 0; (running = (i < r->nbgen)) || pending; i++)
3812 : {
3813 : GEN done;
3814 : long idx;
3815 39144 : mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
3816 39144 : done = mt_queue_get(&pt, &idx, &pending);
3817 39144 : if (!done || lg(done)==1) continue;
3818 34146 : gel(W, idx+1) = done;
3819 34146 : nb += lg(done)-1;
3820 34146 : if (DEBUGLEVEL && (i&127)==0)
3821 0 : err_printf("%ld%% ",100*nb/r->nbmax);
3822 : }
3823 49 : mt_queue_end(&pt);
3824 37870 : for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
3825 : {
3826 : long ll, m;
3827 37821 : GEN L = gel(W,j);
3828 37821 : if (isintzero(L)) continue;
3829 32942 : ll = lg(L);
3830 135394 : for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
3831 : {
3832 102452 : GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
3833 102452 : if (v[1] == 1)
3834 2583 : addifsmooth1(r, h, v[2], v[3]);
3835 : else
3836 99869 : addifsmooth2(r, h, v[2], v[3]);
3837 : }
3838 : }
3839 49 : return j;
3840 : }
3841 :
3842 : static GEN
3843 525 : ECP_psi(GEN x, GEN y)
3844 : {
3845 525 : long prec = realprec(x);
3846 525 : GEN lx = glog(x, prec), ly = glog(y, prec);
3847 525 : GEN u = gdiv(lx, ly);
3848 525 : return gpow(u, gneg(u),prec);
3849 : }
3850 :
3851 : struct computeG
3852 : {
3853 : GEN C;
3854 : long bnd, nbi;
3855 : };
3856 :
3857 : static GEN
3858 525 : _computeG(void *E, GEN gen)
3859 : {
3860 525 : struct computeG * d = (struct computeG *) E;
3861 525 : GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
3862 525 : return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
3863 : }
3864 :
3865 : static long
3866 49 : compute_nbgen(GEN C, long bnd, long nbi)
3867 : {
3868 : struct computeG d;
3869 49 : d.C = shifti(C, 1);
3870 49 : d.bnd = bnd;
3871 49 : d.nbi = nbi;
3872 49 : return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
3873 : }
3874 :
3875 : static GEN
3876 1367 : _psi(void*E, GEN y)
3877 : {
3878 1367 : GEN lx = (GEN) E;
3879 1367 : long prec = realprec(lx);
3880 1367 : GEN ly = glog(y, prec);
3881 1367 : GEN u = gdiv(lx, ly);
3882 1367 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
3883 : }
3884 :
3885 : static GEN
3886 49 : opt_param(GEN x, long prec)
3887 : {
3888 49 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
3889 : }
3890 :
3891 : static GEN
3892 49 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
3893 : {
3894 49 : pari_sp av = avma;
3895 49 : long lM = lg(M)-1, nbcol = lM;
3896 : for (;;)
3897 : {
3898 84 : GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
3899 84 : long i, f=0;
3900 84 : long l = lg(K), lm = lgefint(m);
3901 84 : GEN idx = diviiexact(subiu(p,1),m), g;
3902 : pari_timer ti;
3903 84 : if (DEBUGLEVEL) timer_start(&ti);
3904 140 : for(i=1; i<l; i++)
3905 140 : if (signe(gel(K,i)))
3906 84 : break;
3907 84 : g = Fp_pow(utoi(i), idx, p);
3908 84 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
3909 130438 : for(i=1; i<l; i++)
3910 : {
3911 130354 : GEN k = gel(K,i);
3912 130354 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
3913 130354 : if (signe(k)==0 || !equalii(Fp_pow(g, k, p),
3914 : Fp_pow(j, idx, p)))
3915 80542 : gel(K,i) = cgetineg(lm);
3916 : else
3917 49812 : f++;
3918 : }
3919 84 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld logs", f);
3920 133 : if(f > (nbg>>1)) return gerepileupto(av, K);
3921 9877 : for(i=1; i<=nbcol; i++)
3922 : {
3923 9842 : long a = 1+random_Fl(lM);
3924 9842 : swap(gel(M,a),gel(M,i));
3925 : }
3926 35 : if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
3927 35 : }
3928 : }
3929 :
3930 : static GEN
3931 98 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
3932 : {
3933 98 : pari_sp av=avma;
3934 98 : GEN aa = gen_1;
3935 98 : long AV = 0;
3936 : for(;;)
3937 : {
3938 310 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
3939 310 : GEN F = gel(A,1), E = gel(A,2);
3940 310 : GEN Ao = gen_0;
3941 310 : long i, l = lg(F);
3942 1518 : for(i=1; i<l; i++)
3943 : {
3944 1420 : GEN Ki = gel(K,F[i]);
3945 1420 : if (signe(Ki)<0) break;
3946 1208 : Ao = addii(Ao, mulis(Ki, E[i]));
3947 : }
3948 408 : if (i==l) return Fp_div(Ao, utoi(AV), m);
3949 212 : aa = gerepileuptoint(av, aa);
3950 212 : }
3951 : }
3952 :
3953 : static GEN
3954 49 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
3955 : {
3956 49 : pari_sp av = avma, av2;
3957 49 : long i, j, nbi, nbr = 0, nbrow, nbg;
3958 : GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
3959 : pari_timer ti;
3960 : struct Fp_log_rel r;
3961 49 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
3962 49 : ulong bnd = 4*bnds;
3963 49 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
3964 :
3965 49 : p_1 = subiu(p,1);
3966 49 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
3967 0 : m = diviiexact(p_1, Z_ppo(p_1, m));
3968 49 : pr = primes_upto_zv(bnd);
3969 49 : nbi = lg(pr)-1;
3970 49 : C = sqrtremi(p, &c);
3971 49 : av2 = avma;
3972 12236 : for (i = 1; i <= nbi; ++i)
3973 : {
3974 12187 : ulong lp = pr[i];
3975 37793 : while (lp <= bnd)
3976 : {
3977 13419 : nbr++;
3978 13419 : lp *= pr[i];
3979 : }
3980 : }
3981 49 : pi = cgetg(nbr+1,t_VECSMALL);
3982 49 : Ci = cgetg(nbr+1,t_VECSMALL);
3983 49 : ci = cgetg(nbr+1,t_VECSMALL);
3984 49 : sz = cgetg(nbr+1,t_VECSMALL);
3985 12236 : for (i = 1, j = 1; i <= nbi; ++i)
3986 : {
3987 12187 : ulong lp = pr[i], sp = expu(2*lp-1);
3988 37793 : while (lp <= bnd)
3989 : {
3990 13419 : pi[j] = lp;
3991 13419 : Ci[j] = umodiu(C, lp);
3992 13419 : ci[j] = umodiu(c, lp);
3993 13419 : sz[j] = sp;
3994 13419 : lp *= pr[i];
3995 13419 : j++;
3996 : }
3997 : }
3998 49 : r.nbrel = 0;
3999 49 : r.nbgen = compute_nbgen(C, bnd, nbi);
4000 49 : r.nbmax = 2*(nbi+r.nbgen);
4001 49 : r.rel = cgetg(r.nbmax+1,t_VEC);
4002 49 : r.prmax = pr[nbi];
4003 49 : if (DEBUGLEVEL)
4004 : {
4005 0 : err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
4006 0 : timer_start(&ti);
4007 : }
4008 49 : nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
4009 49 : nbrow = r.prmax + nbg;
4010 49 : if (DEBUGLEVEL)
4011 : {
4012 0 : err_printf("\n");
4013 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
4014 : }
4015 49 : setlg(r.rel,r.nbrel+1);
4016 49 : r.rel = gerepilecopy(av2, r.rel);
4017 49 : K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
4018 49 : if (DEBUGLEVEL) timer_start(&ti);
4019 49 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
4020 49 : if (DEBUGLEVEL) timer_printf(&ti," log element");
4021 49 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
4022 49 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
4023 49 : d = gcdii(Ao,Bo);
4024 49 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
4025 49 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
4026 49 : return gerepileuptoint(av, l);
4027 : }
4028 :
4029 : static int
4030 658132 : Fp_log_use_index(long e, long p)
4031 : {
4032 658132 : return (e >= 27 && 20*(p+6)<=e*e);
4033 : }
4034 :
4035 : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
4036 : static GEN
4037 989990 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
4038 : {
4039 989990 : pari_sp av = avma;
4040 989990 : GEN p = (GEN)E;
4041 : /* assume a reduced mod p, p not necessarily prime */
4042 989990 : if (equali1(a)) return gen_0;
4043 : /* p > 2 */
4044 608344 : if (equalii(subiu(p,1), a)) /* -1 */
4045 : {
4046 : pari_sp av2;
4047 : GEN t;
4048 188362 : ord = get_arith_Z(ord);
4049 188362 : if (mpodd(ord)) { avma = av; return cgetg(1, t_VEC); } /* no solution */
4050 188348 : t = shifti(ord,-1); /* only possible solution */
4051 188348 : av2 = avma;
4052 188348 : if (!equalii(Fp_pow(g, t, p), a)) { avma = av; return cgetg(1, t_VEC); }
4053 188152 : avma = av2; return gerepileuptoint(av, t);
4054 : }
4055 419982 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
4056 49 : return Fp_log_index(a, g, ord, p);
4057 419933 : avma = av; return NULL; /* not easy */
4058 : }
4059 :
4060 : GEN
4061 696384 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
4062 : {
4063 696384 : GEN v = get_arith_ZZM(ord);
4064 696356 : GEN F = gmael(v,2,1);
4065 696356 : long lF = lg(F)-1, lmax;
4066 696356 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
4067 580312 : lmax = expi(gel(F,lF));
4068 580312 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
4069 49 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
4070 580312 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
4071 : }
4072 :
4073 : static ulong
4074 0 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
4075 : {
4076 0 : ulong i, h=1;
4077 0 : for(i=0; i<ord; i++, h = Fl_mul(h, g, p))
4078 0 : if(a==h) return i;
4079 0 : return ~0UL;
4080 : }
4081 :
4082 : static ulong
4083 19088 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
4084 : {
4085 19088 : ulong i, h=1;
4086 48184 : for(i=0; i<ord; i++, h = Fl_mul_pre(h, g, p, pi))
4087 48184 : if(a==h) return i;
4088 0 : return ~0UL;
4089 : }
4090 :
4091 : static ulong
4092 0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
4093 : {
4094 0 : pari_sp av = avma;
4095 0 : GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
4096 0 : ulong z = typ(r)==t_INT ? itou(r): ~0L;
4097 0 : avma = av; return z;
4098 : }
4099 :
4100 : ulong
4101 19088 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
4102 : {
4103 19088 : if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
4104 0 : return Fl_log_Fp(a, g, ord, p);
4105 : }
4106 :
4107 : ulong
4108 0 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
4109 : {
4110 0 : if (ord <= 200)
4111 0 : return (p&HIGHMASK) ? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
4112 0 : : Fl_log_naive(a, g, ord, p);
4113 0 : return Fl_log_Fp(a, g, ord, p);
4114 : }
4115 :
4116 : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
4117 : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
4118 : static GEN
4119 112 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
4120 : {
4121 112 : long l = lg(P) - 1, e = E[l];
4122 112 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
4123 : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
4124 :
4125 112 : if (l == 1) {
4126 84 : hpe = h;
4127 84 : gpe = g;
4128 : } else {
4129 28 : hpe = modii(h, pe);
4130 28 : gpe = modii(g, pe);
4131 : }
4132 112 : if (e == 1) {
4133 28 : hp = hpe;
4134 28 : gp = gpe;
4135 : } else {
4136 84 : hp = remii(hpe, p);
4137 84 : gp = remii(gpe, p);
4138 : }
4139 112 : if (hp == gen_0 || gp == gen_0) return NULL;
4140 91 : if (absequaliu(p, 2))
4141 : {
4142 35 : GEN n = int2n(e);
4143 35 : ogpe = Zp_order(gpe, gen_2, e, n);
4144 35 : a = Fp_log(hpe, gpe, ogpe, n);
4145 35 : if (typ(a) != t_INT) return NULL;
4146 : }
4147 : else
4148 : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
4149 : is trivial */
4150 : /* [order(gp), factor(order(gp))] */
4151 56 : GEN v = Fp_factored_order(gp, subiu(p,1), p);
4152 56 : GEN ogp = gel(v,1);
4153 56 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
4154 56 : a = Fp_log(hp, gp, v, p);
4155 56 : if (typ(a) != t_INT) return NULL;
4156 56 : if (e == 1) ogpe = ogp;
4157 : else
4158 : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
4159 : /* use p-adic log: O(log p + e) mul*/
4160 : long vpogpe, vpohpe;
4161 :
4162 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
4163 28 : gpe = Fp_pow(gpe, ogp, pe);
4164 : /* g,h = 1 mod p; compute b s.t. h = g^b */
4165 :
4166 : /* v_p(order g mod pe) */
4167 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
4168 : /* v_p(order h mod pe) */
4169 28 : vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
4170 28 : if (vpohpe > vpogpe) return NULL;
4171 :
4172 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
4173 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
4174 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
4175 28 : a = addii(a, mulii(ogp, padic_to_Q(b)));
4176 : }
4177 : }
4178 : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
4179 77 : if (l == 1) return a;
4180 :
4181 28 : N = diviiexact(N, pe); /* make N coprime to p */
4182 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
4183 28 : g = Fp_pow(g, modii(ogpe, phi), N);
4184 28 : setlg(P, l); /* remove last element */
4185 28 : setlg(E, l);
4186 28 : b = znlog_rec(h, g, N, P, E, PHI);
4187 28 : if (!b) return NULL;
4188 28 : return addmulii(a, b, ogpe);
4189 : }
4190 :
4191 : static GEN
4192 84 : get_PHI(GEN P, GEN E)
4193 : {
4194 84 : long i, l = lg(P);
4195 84 : GEN PHI = cgetg(l, t_VEC);
4196 84 : gel(PHI,1) = gen_1;
4197 112 : for (i=1; i<l-1; i++)
4198 : {
4199 28 : GEN t, p = gel(P,i);
4200 28 : long e = E[i];
4201 28 : t = mulii(powiu(p, e-1), subiu(p,1));
4202 28 : if (i > 1) t = mulii(t, gel(PHI,i));
4203 28 : gel(PHI,i+1) = t;
4204 : }
4205 84 : return PHI;
4206 : }
4207 :
4208 : GEN
4209 224 : znlog(GEN h, GEN g, GEN o)
4210 : {
4211 224 : pari_sp av = avma;
4212 : GEN N, fa, P, E, x;
4213 224 : switch (typ(g))
4214 : {
4215 : case t_PADIC:
4216 : {
4217 28 : GEN p = gel(g,2);
4218 28 : long v = valp(g);
4219 28 : if (v < 0) pari_err_DIM("znlog");
4220 28 : if (v > 0) {
4221 0 : long k = gvaluation(h, p);
4222 0 : if (k % v) return cgetg(1,t_VEC);
4223 0 : k /= v;
4224 0 : if (!gequal(h, gpowgs(g,k))) { avma = av; return cgetg(1,t_VEC); }
4225 0 : avma = av; return stoi(k);
4226 : }
4227 28 : N = gel(g,3);
4228 28 : g = Rg_to_Fp(g, N);
4229 28 : break;
4230 : }
4231 : case t_INTMOD:
4232 189 : N = gel(g,1);
4233 189 : g = gel(g,2); break;
4234 7 : default: pari_err_TYPE("znlog", g);
4235 : return NULL; /* LCOV_EXCL_LINE */
4236 : }
4237 217 : if (equali1(N)) { avma = av; return gen_0; }
4238 217 : h = Rg_to_Fp(h, N);
4239 210 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
4240 84 : fa = Z_factor(N);
4241 84 : P = gel(fa,1);
4242 84 : E = vec_to_vecsmall(gel(fa,2));
4243 84 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
4244 84 : if (!x) { avma = av; return cgetg(1,t_VEC); }
4245 49 : return gerepileuptoint(av, x);
4246 : }
4247 :
4248 : GEN
4249 61301 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
4250 : {
4251 61301 : if (lgefint(p)==3)
4252 : {
4253 60909 : long nn = itos_or_0(n);
4254 60909 : if (nn)
4255 : {
4256 60909 : ulong pp = p[2];
4257 : ulong uz;
4258 60909 : ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
4259 60888 : if (r==ULONG_MAX) return NULL;
4260 60853 : if (zeta) *zeta = utoi(uz);
4261 60853 : return utoi(r);
4262 : }
4263 : }
4264 392 : a = modii(a,p);
4265 392 : if (!signe(a))
4266 : {
4267 0 : if (zeta) *zeta = gen_1;
4268 0 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
4269 0 : return gen_0;
4270 : }
4271 392 : if (absequaliu(n,2))
4272 : {
4273 224 : if (zeta) *zeta = subiu(p,1);
4274 224 : return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
4275 : }
4276 168 : return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
4277 : }
4278 :
4279 : /*********************************************************************/
4280 : /** **/
4281 : /** FUNDAMENTAL DISCRIMINANTS **/
4282 : /** **/
4283 : /*********************************************************************/
4284 : static int
4285 1407 : fa_isfundamental(GEN F)
4286 : {
4287 1407 : GEN P = gel(F,1), E = gel(F,2);
4288 1407 : long i, s, l = lg(P);
4289 :
4290 1407 : if (l == 1) return 1;
4291 1400 : s = signe(gel(P,1)); /* = signe(x) */
4292 1400 : if (!s) return 0;
4293 1393 : if (s < 0) { l--; P = vecslice(P,2,l); E = vecslice(E,2,l); }
4294 1393 : if (l == 1) return 0;
4295 1386 : if (!absequaliu(gel(P,1), 2))
4296 686 : i = 1; /* need x = 1 mod 4 */
4297 : else
4298 : {
4299 700 : i = 2;
4300 700 : switch(itou(gel(E,1)))
4301 : {
4302 182 : case 2: s = -s; break; /* need x/4 = 3 mod 4 */
4303 84 : case 3: s = 0; break; /* no condition mod 4 */
4304 434 : default: return 0;
4305 : }
4306 : }
4307 1974 : for(; i < l; i++)
4308 : {
4309 1190 : if (!equali1(gel(E,i))) return 0;
4310 1022 : if (s && Mod4(gel(P,i)) == 3) s = -s;
4311 : }
4312 784 : return s >= 0;
4313 : }
4314 : long
4315 17528 : isfundamental(GEN x)
4316 : {
4317 17528 : if (typ(x) != t_INT)
4318 : {
4319 1407 : pari_sp av = avma;
4320 1407 : int v = fa_isfundamental(check_arith_all(x,"isfundamental"));
4321 1407 : avma = av; return v;
4322 : }
4323 16121 : return Z_isfundamental(x);
4324 : }
4325 :
4326 : /* x fundamental ? */
4327 : long
4328 10583 : uposisfundamental(ulong x)
4329 : {
4330 10583 : ulong r = x & 15; /* x mod 16 */
4331 10583 : if (!r) return 0;
4332 10002 : switch(r & 3)
4333 : { /* x mod 4 */
4334 1947 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
4335 3109 : case 1: return uissquarefree(x);
4336 4946 : default: return 0;
4337 : }
4338 : }
4339 : /* -x fundamental ? */
4340 : long
4341 21148 : unegisfundamental(ulong x)
4342 : {
4343 21148 : ulong r = x & 15; /* x mod 16 */
4344 21148 : if (!r) return 0;
4345 20147 : switch(r & 3)
4346 : { /* x mod 4 */
4347 3669 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
4348 9642 : case 3: return uissquarefree(x);
4349 6836 : default: return 0;
4350 : }
4351 : }
4352 : long
4353 10353 : sisfundamental(long x)
4354 10353 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
4355 :
4356 : long
4357 16688 : Z_isfundamental(GEN x)
4358 : {
4359 : long r;
4360 16688 : switch(lgefint(x))
4361 : {
4362 7 : case 2: return 0;
4363 36776 : case 3: return signe(x) < 0? unegisfundamental(x[2])
4364 22105 : : uposisfundamental(x[2]);
4365 : }
4366 2010 : r = mod16(x);
4367 2010 : if (!r) return 0;
4368 1884 : if ((r & 3) == 0)
4369 : {
4370 : pari_sp av;
4371 376 : r >>= 2; /* |x|/4 mod 4 */
4372 376 : if (signe(x) < 0) r = 4-r;
4373 376 : if (r == 1) return 0;
4374 250 : av = avma;
4375 250 : r = Z_issquarefree( shifti(x,-2) );
4376 250 : avma = av; return r;
4377 : }
4378 1508 : r &= 3; /* |x| mod 4 */
4379 1508 : if (signe(x) < 0) r = 4-r;
4380 1508 : return (r==1) ? Z_issquarefree(x) : 0;
4381 : }
4382 :
4383 : static GEN
4384 2821 : fa_quaddisc(GEN f)
4385 : {
4386 2821 : GEN P = gel(f,1), E = gel(f,2), s = gen_1;
4387 2821 : long i, l = lg(P);
4388 9051 : for (i = 1; i < l; i++) /* possibly including -1 */
4389 6230 : if (mpodd(gel(E,i))) s = mulii(s, gel(P,i));
4390 2821 : if (Mod4(s) > 1) s = shifti(s,2);
4391 2821 : return s;
4392 : }
4393 :
4394 : GEN
4395 2821 : quaddisc(GEN x)
4396 : {
4397 2821 : const pari_sp av = avma;
4398 2821 : if (is_rational_t(typ(x))) x = factor(x);
4399 1407 : else x = check_arith_all(x,"quaddisc");
4400 2821 : return gerepileuptoint(av, fa_quaddisc(x));
4401 : }
4402 :
4403 : /*********************************************************************/
4404 : /** **/
4405 : /** FACTORIAL **/
4406 : /** **/
4407 : /*********************************************************************/
4408 : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
4409 : * first is slower ... ] */
4410 : GEN
4411 1023771 : mulu_interval(ulong a, ulong b)
4412 : {
4413 1023771 : pari_sp av = avma;
4414 : ulong k, l, N, n;
4415 : long lx;
4416 : GEN x;
4417 :
4418 1023771 : if (!a) return gen_0;
4419 1023771 : n = b - a + 1;
4420 1023771 : if (n < 61)
4421 : {
4422 1015636 : x = utoi(a);
4423 1015636 : for (k=a+1; k<=b; k++) x = mului(k,x);
4424 1015634 : return gerepileuptoint(av, x);
4425 : }
4426 8135 : lx = 1; x = cgetg(2 + n/2, t_VEC);
4427 8135 : N = b + a;
4428 982734 : for (k = a;; k++)
4429 : {
4430 982734 : l = N - k; if (l <= k) break;
4431 974599 : gel(x,lx++) = muluu(k,l);
4432 974599 : }
4433 8135 : if (l == k) gel(x,lx++) = utoipos(k);
4434 8135 : setlg(x, lx);
4435 8135 : return gerepileuptoint(av, ZV_prod(x));
4436 : }
4437 : GEN
4438 448 : muls_interval(long a, long b)
4439 : {
4440 448 : pari_sp av = avma;
4441 448 : long lx, k, l, N, n = b - a + 1;
4442 : GEN x;
4443 :
4444 448 : if (a <= 0 && b >= 0) return gen_0;
4445 287 : if (n < 61)
4446 : {
4447 287 : x = stoi(a);
4448 287 : for (k=a+1; k<=b; k++) x = mulsi(k,x);
4449 287 : return gerepileuptoint(av, x);
4450 : }
4451 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
4452 0 : N = b + a;
4453 0 : for (k = a;; k++)
4454 : {
4455 0 : l = N - k; if (l <= k) break;
4456 0 : gel(x,lx++) = mulss(k,l);
4457 0 : }
4458 0 : if (l == k) gel(x,lx++) = stoi(k);
4459 0 : setlg(x, lx);
4460 0 : return gerepileuptoint(av, ZV_prod(x));
4461 : }
4462 :
4463 : GEN
4464 2464121 : mpfact(long n)
4465 : {
4466 2464121 : if (n < 2)
4467 : {
4468 1561345 : if (n < 0) pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
4469 1561345 : return gen_1;
4470 : }
4471 902776 : return mulu_interval(2UL, (ulong)n);
4472 : }
4473 :
4474 : /*******************************************************************/
4475 : /** **/
4476 : /** LUCAS & FIBONACCI **/
4477 : /** **/
4478 : /*******************************************************************/
4479 : static void
4480 56 : lucas(ulong n, GEN *a, GEN *b)
4481 : {
4482 : GEN z, t, zt;
4483 112 : if (!n) { *a = gen_2; *b = gen_1; return; }
4484 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
4485 49 : switch(n & 3) {
4486 14 : case 0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
4487 14 : case 1: *a = subiu(zt,1); *b = addiu(sqri(t),2); break;
4488 7 : case 2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
4489 14 : case 3: *a = addiu(zt,1); *b = subiu(sqri(t),2);
4490 : }
4491 : }
4492 :
4493 : GEN
4494 7 : fibo(long n)
4495 : {
4496 7 : pari_sp av = avma;
4497 : GEN a, b;
4498 7 : if (!n) return gen_0;
4499 7 : lucas((ulong)(labs(n)-1), &a, &b);
4500 7 : a = diviuexact(addii(shifti(a,1),b), 5);
4501 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
4502 7 : return gerepileuptoint(av, a);
4503 : }
4504 :
4505 : /*******************************************************************/
4506 : /* */
4507 : /* CONTINUED FRACTIONS */
4508 : /* */
4509 : /*******************************************************************/
4510 : static GEN
4511 488070 : icopy_lg(GEN x, long l)
4512 : {
4513 488070 : long lx = lgefint(x);
4514 : GEN y;
4515 :
4516 488070 : if (lx >= l) return icopy(x);
4517 35 : y = cgeti(l); affii(x, y); return y;
4518 : }
4519 :
4520 : /* continued fraction of a/b. If y != NULL, stop when partial quotients
4521 : * differ from y */
4522 : static GEN
4523 488366 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
4524 : {
4525 : GEN z, c;
4526 488366 : ulong i, l, ly = lgefint(b);
4527 :
4528 : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
4529 488366 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
4530 488366 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
4531 488366 : if (l > LGBITS) l = LGBITS;
4532 :
4533 488366 : z = cgetg(l,t_VEC);
4534 488366 : l--;
4535 488366 : if (y) {
4536 296 : pari_sp av = avma;
4537 296 : if (l >= (ulong)lg(y)) l = lg(y)-1;
4538 19488 : for (i = 1; i <= l; i++)
4539 : {
4540 19303 : GEN q = gel(y,i);
4541 19303 : gel(z,i) = q;
4542 19303 : c = b; if (!gequal1(q)) c = mulii(q, b);
4543 19303 : c = subii(a, c);
4544 19303 : if (signe(c) < 0)
4545 : { /* partial quotient too large */
4546 110 : c = addii(c, b);
4547 110 : if (signe(c) >= 0) i++; /* by 1 */
4548 110 : break;
4549 : }
4550 19193 : if (cmpii(c, b) >= 0)
4551 : { /* partial quotient too small */
4552 1 : c = subii(c, b);
4553 1 : if (cmpii(c, b) < 0) {
4554 : /* by 1. If next quotient is 1 in y, add 1 */
4555 0 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
4556 0 : i++;
4557 : }
4558 1 : break;
4559 : }
4560 19192 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
4561 19192 : a = b; b = c;
4562 : }
4563 : } else {
4564 488070 : a = icopy_lg(a, ly);
4565 488070 : b = icopy(b);
4566 1737469 : for (i = 1; i <= l; i++)
4567 : {
4568 1737205 : gel(z,i) = truedvmdii(a,b,&c);
4569 1737205 : if (c == gen_0) { i++; break; }
4570 1249399 : affii(c, a); cgiv(c); c = a;
4571 1249399 : a = b; b = c;
4572 : }
4573 : }
4574 488366 : i--;
4575 488366 : if (i > 1 && gequal1(gel(z,i)))
4576 : {
4577 85 : cgiv(gel(z,i)); --i;
4578 85 : gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
4579 : }
4580 488366 : setlg(z,i+1); return z;
4581 : }
4582 :
4583 : static GEN
4584 0 : sersfcont(GEN a, GEN b, long k)
4585 : {
4586 0 : long i, l = typ(a) == t_POL? lg(a): 3;
4587 : GEN y, c;
4588 0 : if (lg(b) > l) l = lg(b);
4589 0 : if (k > 0 && l > k+1) l = k+1;
4590 0 : y = cgetg(l,t_VEC);
4591 0 : for (i=1; i<l; i++)
4592 : {
4593 0 : gel(y,i) = poldivrem(a,b,&c);
4594 0 : if (gequal0(c)) { i++; break; }
4595 0 : a = b; b = c;
4596 : }
4597 0 : setlg(y, i); return y;
4598 : }
4599 :
4600 : GEN
4601 488588 : gboundcf(GEN x, long k)
4602 : {
4603 : pari_sp av;
4604 488588 : long tx = typ(x), e;
4605 : GEN y, a, b, c;
4606 :
4607 488588 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
4608 488581 : if (is_scalar_t(tx))
4609 : {
4610 488581 : if (gequal0(x)) return mkvec(gen_0);
4611 488532 : switch(tx)
4612 : {
4613 455 : case t_INT: return mkveccopy(x);
4614 : case t_REAL:
4615 303 : av = avma;
4616 303 : c = mantissa_real(x,&e);
4617 303 : if (e < 0) pari_err_PREC("gboundcf");
4618 296 : y = int2n(e);
4619 296 : a = Qsfcont(c,y, NULL, k);
4620 296 : b = addsi(signe(x), c);
4621 296 : return gerepilecopy(av, Qsfcont(b,y, a, k));
4622 :
4623 : case t_FRAC:
4624 487774 : av = avma;
4625 487774 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
4626 : }
4627 0 : pari_err_TYPE("gboundcf",x);
4628 : }
4629 :
4630 0 : switch(tx)
4631 : {
4632 0 : case t_POL: return mkveccopy(x);
4633 : case t_SER:
4634 0 : av = avma;
4635 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
4636 : case t_RFRAC:
4637 0 : av = avma;
4638 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
4639 : }
4640 0 : pari_err_TYPE("gboundcf",x);
4641 : return NULL; /* LCOV_EXCL_LINE */
4642 : }
4643 :
4644 : static GEN
4645 14 : sfcont2(GEN b, GEN x, long k)
4646 : {
4647 14 : pari_sp av = avma;
4648 14 : long lb = lg(b), tx = typ(x), i;
4649 : GEN y,p1;
4650 :
4651 14 : if (k)
4652 : {
4653 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
4654 0 : lb = k+1;
4655 : }
4656 7 : y = cgetg(lb,t_VEC);
4657 7 : if (lb==1) return y;
4658 7 : if (is_scalar_t(tx))
4659 : {
4660 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
4661 : }
4662 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
4663 :
4664 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
4665 7 : for (i = 1;;)
4666 : {
4667 35 : if (tx == t_REAL)
4668 : {
4669 35 : long e = expo(x);
4670 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
4671 35 : gel(y,i) = floorr(x);
4672 35 : p1 = subri(x, gel(y,i));
4673 : }
4674 : else
4675 : {
4676 0 : gel(y,i) = gfloor(x);
4677 0 : p1 = gsub(x, gel(y,i));
4678 : }
4679 35 : if (++i >= lb) break;
4680 28 : if (gequal0(p1)) break;
4681 28 : x = gdiv(gel(b,i),p1);
4682 28 : }
4683 7 : setlg(y,i);
4684 7 : return gerepilecopy(av,y);
4685 : }
4686 :
4687 :
4688 : GEN
4689 91 : gcf(GEN x) { return gboundcf(x,0); }
4690 : GEN
4691 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
4692 : GEN
4693 49 : contfrac0(GEN x, GEN b, long nmax)
4694 : {
4695 : long tb;
4696 :
4697 49 : if (!b) return gboundcf(x,nmax);
4698 28 : tb = typ(b);
4699 28 : if (tb == t_INT) return gboundcf(x,itos(b));
4700 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
4701 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
4702 14 : return sfcont2(b,x,nmax);
4703 : }
4704 :
4705 : GEN
4706 231 : contfracpnqn(GEN x, long n)
4707 : {
4708 231 : pari_sp av = avma;
4709 231 : long i, lx = lg(x);
4710 : GEN M,A,B, p0,p1, q0,q1;
4711 :
4712 231 : if (lx == 1)
4713 : {
4714 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
4715 21 : if (n >= 0) return cgetg(1,t_MAT);
4716 7 : return matid(2);
4717 : }
4718 203 : switch(typ(x))
4719 : {
4720 161 : case t_VEC: case t_COL: A = x; B = NULL; break;
4721 : case t_MAT:
4722 42 : switch(lgcols(x))
4723 : {
4724 0 : case 2: A = row(x,1); B = NULL; break;
4725 35 : case 3: A = row(x,2); B = row(x,1); break;
4726 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
4727 : return NULL; /*LCOV_EXCL_LINE*/
4728 : }
4729 35 : break;
4730 0 : default: pari_err_TYPE("pnqn",x);
4731 : return NULL; /*LCOV_EXCL_LINE*/
4732 : }
4733 196 : p1 = gel(A,1);
4734 196 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
4735 196 : if (n >= 0)
4736 : {
4737 161 : lx = minss(lx, n+2);
4738 161 : if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
4739 : }
4740 35 : else if (lx == 2)
4741 7 : return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
4742 : /* lx >= 3 */
4743 98 : p0 = gen_1;
4744 98 : q0 = gen_0; /* p[-1], q[-1] */
4745 98 : M = cgetg(lx, t_MAT);
4746 98 : gel(M,1) = mkcol2(p1,q1);
4747 336 : for (i=2; i<lx; i++)
4748 : {
4749 238 : GEN a = gel(A,i), p2,q2;
4750 238 : if (B) {
4751 84 : GEN b = gel(B,i);
4752 84 : p0 = gmul(b,p0);
4753 84 : q0 = gmul(b,q0);
4754 : }
4755 238 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
4756 238 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
4757 238 : gel(M,i) = mkcol2(p1,q1);
4758 : }
4759 98 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
4760 98 : return gerepilecopy(av, M);
4761 : }
4762 : GEN
4763 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
4764 : /* x = [a0, ..., an] from gboundcf, n >= 0;
4765 : * return [[p0, ..., pn], [q0,...,qn]] */
4766 : GEN
4767 480550 : ZV_allpnqn(GEN x)
4768 : {
4769 480550 : long i, lx = lg(x);
4770 480550 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
4771 :
4772 480550 : gel(v,1) = P = cgetg(lx, t_VEC);
4773 480550 : gel(v,2) = Q = cgetg(lx, t_VEC);
4774 480550 : p0 = gen_1; q0 = gen_0;
4775 480550 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
4776 1684844 : for (i=2; i<lx; i++)
4777 : {
4778 1204294 : GEN a = gel(x,i);
4779 1204294 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
4780 1204294 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
4781 : }
4782 480550 : return v;
4783 : }
4784 :
4785 : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
4786 : static GEN
4787 42 : mod_to_frac(GEN x, GEN N, GEN B)
4788 : {
4789 : GEN a, b, A;
4790 42 : if (B) A = divii(shifti(N, -1), B);
4791 : else
4792 : {
4793 14 : A = sqrti(shifti(N, -1));
4794 14 : B = A;
4795 : }
4796 42 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
4797 28 : return equali1(b)? a: mkfrac(a,b);
4798 : }
4799 :
4800 : static GEN
4801 70 : mod_to_rfrac(GEN x, GEN N, long B)
4802 : {
4803 : GEN a, b;
4804 70 : long A, d = degpol(N);
4805 70 : if (B >= 0) A = d-1 - B;
4806 : else
4807 : {
4808 42 : B = d >> 1;
4809 42 : A = odd(d)? B : B-1;
4810 : }
4811 70 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
4812 70 : if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
4813 56 : return gdiv(a,b);
4814 : }
4815 :
4816 : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
4817 : * of the continued fraction of x with b <= k maximal */
4818 : static GEN
4819 0 : bestappr_frac(GEN x, GEN k)
4820 : {
4821 : pari_sp av;
4822 : GEN p0, p1, p, q0, q1, q, a, y;
4823 :
4824 0 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
4825 0 : av = avma; y = x;
4826 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
4827 0 : q1 = gen_0; q0 = gen_1;
4828 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
4829 : for(;;)
4830 : {
4831 0 : x = ginv(x); /* > 1 */
4832 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
4833 0 : if (cmpii(a,k) > 0)
4834 : { /* next partial quotient will overflow limits */
4835 : GEN n, d;
4836 0 : a = divii(subii(k, q1), q0);
4837 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4838 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4839 : /* compare |y-p0/q0|, |y-p1/q1| */
4840 0 : n = gel(y,1);
4841 0 : d = gel(y,2);
4842 0 : if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
4843 : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
4844 0 : { p1 = p0; q1 = q0; }
4845 0 : break;
4846 : }
4847 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4848 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4849 :
4850 0 : if (cmpii(q0,k) > 0) break;
4851 0 : x = gsub(x,a); /* 0 <= x < 1 */
4852 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
4853 :
4854 0 : }
4855 0 : return gerepileupto(av, gdiv(p1,q1));
4856 : }
4857 : /* bestappr(t_REAL != 0), to maximal accuracy */
4858 : static GEN
4859 287 : bestappr_real_max(GEN x)
4860 : {
4861 287 : pari_sp av = avma;
4862 : GEN p0, p1, p, q0, q1, q, a;
4863 : long e;
4864 287 : p1 = gen_1; a = p0 = floorr(x); q1 = gen_0; q0 = gen_1;
4865 287 : x = subri(x,a); /* 0 <= x < 1 */
4866 287 : e = bit_prec(x) - expo(x);
4867 : for(;;)
4868 : {
4869 : long d;
4870 1113 : if (!signe(x) || expi(q0) > e) { p1 = p0; q1 = q0; break; }
4871 1004 : x = invr(x); /* > 1 */
4872 1004 : d = nbits2prec(expo(x) + 1);
4873 1004 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4874 :
4875 826 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4876 826 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4877 826 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4878 826 : x = subri(x,a); /* 0 <= x < 1 */
4879 826 : }
4880 287 : return gerepileupto(av, gdiv(p1,q1));
4881 : }
4882 : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
4883 : * of the continued fraction of x with b <= k maximal */
4884 : static GEN
4885 295568 : bestappr_real(GEN x, GEN k)
4886 : {
4887 295568 : pari_sp av = avma;
4888 295568 : GEN kr, p0, p1, p, q0, q1, q, a, y = x;
4889 :
4890 295568 : p1 = gen_1; a = p0 = floorr(x);
4891 295568 : q1 = gen_0; q0 = gen_1;
4892 295568 : x = subri(x,a); /* 0 <= x < 1 */
4893 295568 : if (!signe(x)) { cgiv(x); return a; }
4894 288999 : kr = itor(k, realprec(x));
4895 : for(;;)
4896 : {
4897 : long d;
4898 757369 : x = invr(x); /* > 1 */
4899 757369 : if (cmprr(x,kr) > 0)
4900 : { /* next partial quotient will overflow limits */
4901 284421 : a = divii(subii(k, q1), q0);
4902 284421 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4903 284421 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4904 : /* compare |y-p0/q0|, |y-p1/q1| */
4905 284421 : if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
4906 : mulir(q0, subri(mulir(q1,y), p1))) < 0)
4907 5048 : { p1 = p0; q1 = q0; }
4908 284421 : break;
4909 : }
4910 472948 : d = nbits2prec(expo(x) + 1);
4911 472948 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4912 :
4913 472942 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4914 472942 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4915 472942 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4916 :
4917 472942 : if (cmpii(q0,k) > 0) break;
4918 468775 : x = subri(x,a); /* 0 <= x < 1 */
4919 468775 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
4920 468370 : }
4921 288999 : return gerepileupto(av, gdiv(p1,q1));
4922 : }
4923 :
4924 : /* k t_INT or NULL */
4925 : static GEN
4926 466580 : bestappr_Q(GEN x, GEN k)
4927 : {
4928 466580 : long lx, tx = typ(x), i;
4929 : GEN a, y;
4930 :
4931 466580 : switch(tx)
4932 : {
4933 77 : case t_INT: return icopy(x);
4934 0 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
4935 : case t_REAL:
4936 332463 : if (!signe(x)) return gen_0;
4937 295856 : if (bit_prec(x) < expo(x)) return NULL;
4938 295855 : return k? bestappr_real(x, k): bestappr_real_max(x);
4939 :
4940 : case t_INTMOD: {
4941 28 : pari_sp av = avma;
4942 28 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
4943 21 : return gerepilecopy(av, a);
4944 : }
4945 : case t_PADIC: {
4946 14 : pari_sp av = avma;
4947 14 : long v = valp(x);
4948 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
4949 7 : if (v) a = gmul(a, powis(gel(x,2), v));
4950 7 : return gerepilecopy(av, a);
4951 : }
4952 :
4953 : case t_COMPLEX: {
4954 126 : pari_sp av = avma;
4955 126 : y = cgetg(3, t_COMPLEX);
4956 126 : gel(y,2) = bestappr(gel(x,2), k);
4957 126 : gel(y,1) = bestappr(gel(x,1), k);
4958 126 : if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
4959 0 : return y;
4960 : }
4961 : case t_SER:
4962 0 : if (ser_isexactzero(x)) return gcopy(x);
4963 : /* fall through */
4964 : case t_POLMOD: case t_POL: case t_RFRAC:
4965 : case t_VEC: case t_COL: case t_MAT:
4966 133872 : y = cgetg_copy(x, &lx);
4967 133872 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
4968 584295 : for (; i<lx; i++)
4969 : {
4970 450425 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
4971 450423 : gel(y,i) = a;
4972 : }
4973 133870 : if (tx == t_POL) return normalizepol(y);
4974 133856 : if (tx == t_SER) return normalize(y);
4975 133856 : return y;
4976 : }
4977 0 : pari_err_TYPE("bestappr_Q",x);
4978 : return NULL; /* LCOV_EXCL_LINE */
4979 : }
4980 :
4981 : static GEN
4982 56 : bestappr_ser(GEN x, long B)
4983 : {
4984 56 : long dN, v = valp(x), lx = lg(x);
4985 : GEN t;
4986 56 : x = normalizepol(ser2pol_i(x, lx));
4987 56 : dN = lx-2;
4988 56 : if (v > 0)
4989 : {
4990 14 : x = RgX_shift_shallow(x, v);
4991 14 : dN += v;
4992 : }
4993 42 : else if (v < 0)
4994 : {
4995 7 : if (B >= 0) B = maxss(B+v, 0);
4996 : }
4997 56 : t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
4998 56 : if (!t) return NULL;
4999 42 : if (v < 0)
5000 : {
5001 : GEN a, b;
5002 : long vx;
5003 7 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
5004 : /* t_RFRAC */
5005 7 : vx = varn(x);
5006 7 : a = gel(t,1);
5007 7 : b = gel(t,2);
5008 7 : v -= RgX_valrem(b, &b);
5009 7 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
5010 7 : if (v < 0) b = RgX_shift(b, -v);
5011 0 : else if (v > 0) {
5012 0 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
5013 0 : a = RgX_shift(a, v);
5014 : }
5015 7 : t = mkrfraccopy(a, b);
5016 : }
5017 42 : return t;
5018 : }
5019 : static GEN bestappr_RgX(GEN x, long B);
5020 : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
5021 : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
5022 : static GEN
5023 77 : bestappr_RgX(GEN x, long B)
5024 : {
5025 77 : long i, lx, tx = typ(x);
5026 : GEN y, t;
5027 77 : switch(tx)
5028 : {
5029 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
5030 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
5031 0 : return gcopy(x);
5032 :
5033 : case t_RFRAC: {
5034 14 : pari_sp av = avma;
5035 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
5036 7 : x = rfrac_to_ser(x, 2*B+1);
5037 7 : t = bestappr_ser(x, B); if (!t) return NULL;
5038 0 : return gerepileupto(av, t);
5039 : }
5040 : case t_POLMOD: {
5041 14 : pari_sp av = avma;
5042 14 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
5043 14 : return gerepileupto(av, t);
5044 : }
5045 : case t_SER: {
5046 49 : pari_sp av = avma;
5047 49 : t = bestappr_ser(x, B); if (!t) return NULL;
5048 42 : return gerepileupto(av, t);
5049 : }
5050 :
5051 : case t_VEC: case t_COL: case t_MAT:
5052 0 : y = cgetg_copy(x, &lx);
5053 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
5054 0 : for (; i<lx; i++)
5055 : {
5056 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
5057 0 : gel(y,i) = t;
5058 : }
5059 0 : return y;
5060 : }
5061 0 : pari_err_TYPE("bestappr_RgX",x);
5062 : return NULL; /* LCOV_EXCL_LINE */
5063 : }
5064 :
5065 : /* allow k = NULL: maximal accuracy */
5066 : GEN
5067 16155 : bestappr(GEN x, GEN k)
5068 : {
5069 16155 : pari_sp av = avma;
5070 16155 : if (k) { /* replace by floor(k) */
5071 15889 : switch(typ(k))
5072 : {
5073 : case t_INT:
5074 1029 : break;
5075 : case t_REAL: case t_FRAC:
5076 14860 : k = floor_safe(k); /* left on stack for efficiency */
5077 14860 : if (!signe(k)) k = gen_1;
5078 14860 : break;
5079 : default:
5080 0 : pari_err_TYPE("bestappr [bound type]", k);
5081 0 : break;
5082 : }
5083 : }
5084 16155 : x = bestappr_Q(x, k);
5085 16155 : if (!x) { avma = av; return cgetg(1,t_VEC); }
5086 16140 : return x;
5087 : }
5088 : GEN
5089 77 : bestapprPade(GEN x, long B)
5090 : {
5091 77 : pari_sp av = avma;
5092 77 : GEN t = bestappr_RgX(x, B);
5093 77 : if (!t) { avma = av; return cgetg(1,t_VEC); }
5094 63 : return t;
5095 : }
5096 :
5097 : /***********************************************************************/
5098 : /** **/
5099 : /** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
5100 : /** **/
5101 : /***********************************************************************/
5102 :
5103 : static GEN
5104 14 : get_quad(GEN f, GEN pol, long r)
5105 : {
5106 14 : GEN p1 = gcoeff(f,1,2), q1 = gcoeff(f,2,2);
5107 14 : return mkquad(pol, r? subii(p1,q1): p1, q1);
5108 : }
5109 :
5110 : /* replace f by f * [a,1; 1,0] */
5111 : static void
5112 14 : update_f(GEN f, GEN a)
5113 : {
5114 : GEN p1;
5115 14 : p1 = gcoeff(f,1,1);
5116 14 : gcoeff(f,1,1) = addii(mulii(a,p1), gcoeff(f,1,2));
5117 14 : gcoeff(f,1,2) = p1;
5118 :
5119 14 : p1 = gcoeff(f,2,1);
5120 14 : gcoeff(f,2,1) = addii(mulii(a,p1), gcoeff(f,2,2));
5121 14 : gcoeff(f,2,2) = p1;
5122 14 : }
5123 :
5124 : GEN
5125 7 : quadunit(GEN x)
5126 : {
5127 7 : pari_sp av = avma, av2;
5128 : GEN pol, y, a, u, v, sqd, f;
5129 : long r;
5130 :
5131 7 : check_quaddisc_real(x, &r, "quadunit");
5132 7 : pol = quadpoly(x);
5133 7 : sqd = sqrti(x); av2 = avma;
5134 7 : a = shifti(addui(r,sqd),-1);
5135 7 : f = mkmat2(mkcol2(a, gen_1), mkcol2(gen_1, gen_0)); /* [a,0; 1,0] */
5136 7 : u = stoi(r); v = gen_2;
5137 : for(;;)
5138 : {
5139 : GEN u1, v1;
5140 14 : u1 = subii(mulii(a,v),u);
5141 14 : v1 = divii(subii(x,sqri(u1)),v);
5142 14 : if ( equalii(v,v1) ) {
5143 7 : y = get_quad(f,pol,r);
5144 7 : update_f(f,a);
5145 7 : y = gdiv(get_quad(f,pol,r), gconj(y));
5146 7 : break;
5147 : }
5148 7 : a = divii(addii(sqd,u1), v1);
5149 7 : if ( equalii(u,u1) ) {
5150 0 : y = get_quad(f,pol,r);
5151 0 : y = gdiv(y, gconj(y));
5152 0 : break;
5153 : }
5154 7 : update_f(f,a);
5155 7 : u = u1; v = v1;
5156 7 : if (gc_needed(av2,2))
5157 : {
5158 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit");
5159 0 : gerepileall(av2,4, &a,&f,&u,&v);
5160 : }
5161 7 : }
5162 7 : if (signe(gel(y,3)) < 0) y = gneg(y);
5163 7 : return gerepileupto(av, y);
5164 : }
5165 :
5166 : GEN
5167 7 : quadunit0(GEN x, long v)
5168 : {
5169 7 : GEN y = quadunit(x);
5170 7 : if (v==-1) v = fetch_user_var("w");
5171 7 : setvarn(gel(y,1), v);
5172 7 : return y;
5173 : }
5174 :
5175 : GEN
5176 21 : quadregulator(GEN x, long prec)
5177 : {
5178 21 : pari_sp av = avma, av2;
5179 : GEN R, rsqd, u, v, sqd;
5180 : long r, Rexpo;
5181 :
5182 21 : check_quaddisc_real(x, &r, "quadregulator");
5183 21 : sqd = sqrti(x);
5184 21 : rsqd = gsqrt(x,prec);
5185 21 : Rexpo = 0; R = real2n(1, prec); /* = 2 */
5186 21 : av2 = avma;
5187 21 : u = stoi(r); v = gen_2;
5188 : for(;;)
5189 : {
5190 70 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
5191 70 : GEN v1 = divii(subii(x,sqri(u1)),v);
5192 70 : if (equalii(v,v1))
5193 : {
5194 7 : R = sqrr(R); shiftr_inplace(R, -1);
5195 7 : R = mulrr(R, divri(addir(u1,rsqd),v));
5196 7 : break;
5197 : }
5198 63 : if (equalii(u,u1))
5199 : {
5200 14 : R = sqrr(R); shiftr_inplace(R, -1);
5201 14 : break;
5202 : }
5203 49 : R = mulrr(R, divri(addir(u1,rsqd),v));
5204 49 : Rexpo += expo(R); setexpo(R,0);
5205 49 : u = u1; v = v1;
5206 49 : if (Rexpo & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
5207 49 : if (gc_needed(av2,2))
5208 : {
5209 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
5210 0 : gerepileall(av2,3, &R,&u,&v);
5211 : }
5212 49 : }
5213 21 : R = logr_abs(divri(R,v));
5214 21 : if (Rexpo)
5215 : {
5216 21 : GEN t = mulsr(Rexpo, mplog2(prec));
5217 21 : shiftr_inplace(t, 1);
5218 21 : R = addrr(R,t);
5219 : }
5220 21 : return gerepileuptoleaf(av, R);
5221 : }
5222 :
5223 : /*************************************************************************/
5224 : /** **/
5225 : /** CLASS NUMBER **/
5226 : /** **/
5227 : /*************************************************************************/
5228 :
5229 : int
5230 12869702 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
5231 :
5232 18205555 : static GEN qfi_pow(void *E, GEN f, GEN n)
5233 18205555 : { return E? nupow(f,n,(GEN)E): powgi(f,n); }
5234 22882105 : static GEN qfi_comp(void *E, GEN f, GEN g)
5235 22882105 : { return E? nucomp(f,g,(GEN)E): qficomp(f,g); }
5236 : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
5237 : gidentical,qfb_equal1,NULL};
5238 :
5239 : GEN
5240 2926471 : qfi_order(GEN q, GEN o)
5241 2926471 : { return gen_order(q, o, NULL, &qfi_group); }
5242 :
5243 : GEN
5244 0 : qfi_log(GEN a, GEN g, GEN o)
5245 0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
5246 :
5247 : GEN
5248 623434 : qfi_Shanks(GEN a, GEN g, long n)
5249 : {
5250 623434 : pari_sp av = avma;
5251 : GEN T, X;
5252 : long rt_n, c;
5253 :
5254 623434 : a = redimag(a);
5255 623434 : g = redimag(g);
5256 :
5257 623434 : rt_n = sqrt((double)n);
5258 623434 : c = n / rt_n;
5259 623434 : c = (c * rt_n < n + 1) ? c + 1 : c;
5260 :
5261 623434 : T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
5262 623434 : X = gen_Shanks(T, a, c, NULL, &qfi_group);
5263 :
5264 623434 : if (!X) { avma = av; return X; }
5265 330659 : return gerepileuptoint(av, X);
5266 : }
5267 :
5268 : GEN
5269 140 : qfbclassno0(GEN x,long flag)
5270 : {
5271 140 : switch(flag)
5272 : {
5273 126 : case 0: return map_proto_G(classno,x);
5274 14 : case 1: return map_proto_G(classno2,x);
5275 0 : default: pari_err_FLAG("qfbclassno");
5276 : }
5277 : return NULL; /* LCOV_EXCL_LINE */
5278 : }
5279 :
5280 : /* f^h = 1, return order(f). Set *pfao to its factorization */
5281 : static GEN
5282 2756128 : find_order(void *E, GEN f, GEN h, GEN *pfao)
5283 : {
5284 2756128 : GEN v = gen_factored_order(f, h, E, &qfi_group);
5285 2756128 : *pfao = gel(v,2); return gel(v,1);
5286 : }
5287 :
5288 : static int
5289 6663 : ok_q(GEN q, GEN h, GEN d2, long r2)
5290 : {
5291 6663 : if (d2)
5292 : {
5293 7 : if (r2 <= 2 && !mpodd(q)) return 0;
5294 7 : return is_pm1(Z_ppo(q,d2));
5295 : }
5296 : else
5297 : {
5298 6656 : if (r2 <= 1 && !mpodd(q)) return 0;
5299 6656 : return is_pm1(Z_ppo(q,h));
5300 : }
5301 : }
5302 :
5303 : /* a,b given by their factorizations. Return factorization of lcm(a,b).
5304 : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
5305 : static GEN
5306 359404 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
5307 : {
5308 359404 : GEN P = ZV_union_shallow(gel(Fa,1), gel(Fb,1));
5309 359404 : GEN A = gen_1, B = gen_1;
5310 359404 : long i, l = lg(P);
5311 359404 : GEN E = cgetg(l, t_COL);
5312 1062920 : for (i=1; i<l; i++)
5313 : {
5314 703516 : GEN p = gel(P,i);
5315 703516 : long va = Z_pval(a,p);
5316 703516 : long vb = Z_pval(b,p);
5317 703516 : if (va < vb)
5318 : {
5319 361517 : B = mulii(B,powiu(p,vb));
5320 361517 : gel(E,i) = utoi(vb);
5321 : }
5322 : else
5323 : {
5324 341999 : A = mulii(A,powiu(p,va));
5325 341999 : gel(E,i) = utoi(va);
5326 : }
5327 : }
5328 359404 : *pA = A;
5329 359404 : *pB = B; return mkmat2(P,E);
5330 : }
5331 :
5332 : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
5333 : static void
5334 359404 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
5335 : {
5336 359404 : GEN A,B, g1 = *pg1, d1 = *pd1;
5337 359404 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
5338 359404 : *pg1 = gmul(powgi(g1, diviiexact(d1,A)), powgi(f, diviiexact(o,B)));
5339 359404 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
5340 359404 : }
5341 :
5342 : /* Write x = Df^2, where D = fundamental discriminant,
5343 : * P^E = factorisation of conductor f, with E[i] >= 0 */
5344 : static void
5345 2044969 : corediscfact(GEN x, long xmod4, GEN *ptD, GEN *ptP, GEN *ptE)
5346 : {
5347 2044969 : long s = signe(x), l, i;
5348 2044969 : GEN fa = absZ_factor(x);
5349 2044969 : GEN d, P = gel(fa,1), E = gtovecsmall(gel(fa,2));
5350 :
5351 2044969 : l = lg(P); d = gen_1;
5352 5323296 : for (i=1; i<l; i++)
5353 : {
5354 3278327 : if (E[i] & 1) d = mulii(d, gel(P,i));
5355 3278327 : E[i] >>= 1;
5356 : }
5357 2044969 : if (!xmod4 && mod4(d) != ((s < 0)? 3: 1)) { d = shifti(d,2); E[1]--; }
5358 2044969 : *ptD = (s < 0)? negi(d): d;
5359 2044969 : *ptP = P;
5360 2044969 : *ptE = E;
5361 2044969 : }
5362 :
5363 : static GEN
5364 2037166 : conductor_part(GEN x, long xmod4, GEN *ptD, GEN *ptreg)
5365 : {
5366 2037166 : long l, i, s = signe(x);
5367 : GEN E, H, D, P, reg;
5368 :
5369 2037166 : corediscfact(x, xmod4, &D, &P, &E);
5370 2037166 : H = gen_1; l = lg(P);
5371 : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
5372 5289649 : for (i=1; i<l; i++)
5373 : {
5374 3252483 : long e = E[i];
5375 3252483 : if (e)
5376 : {
5377 7 : GEN p = gel(P,i);
5378 7 : H = mulii(H, subis(p, kronecker(D,p)));
5379 7 : if (e >= 2) H = mulii(H, powiu(p,e-1));
5380 : }
5381 : }
5382 :
5383 : /* divide by [ O_K^* : O^* ] */
5384 2037166 : if (s < 0)
5385 : {
5386 2037152 : reg = NULL;
5387 2037152 : switch(itou_or_0(D))
5388 : {
5389 0 : case 4: H = shifti(H,-1); break;
5390 0 : case 3: H = divis(H,3); break;
5391 : }
5392 : } else {
5393 14 : reg = quadregulator(D,DEFAULTPREC);
5394 14 : if (!equalii(x,D))
5395 0 : H = divii(H, roundr(divrr(quadregulator(x,DEFAULTPREC), reg)));
5396 : }
5397 2037166 : if (ptreg) *ptreg = reg;
5398 2037166 : *ptD = D; return H;
5399 : }
5400 :
5401 : static long
5402 2037145 : two_rank(GEN x)
5403 : {
5404 2037145 : GEN p = gel(absZ_factor(x),1);
5405 2037145 : long l = lg(p)-1;
5406 : #if 0 /* positive disc not needed */
5407 : if (signe(x) > 0)
5408 : {
5409 : long i;
5410 : for (i=1; i<=l; i++)
5411 : if (mod4(gel(p,i)) == 3) { l--; break; }
5412 : }
5413 : #endif
5414 2037145 : return l-1;
5415 : }
5416 :
5417 : static GEN
5418 38703604 : sqr_primeform(GEN x, ulong p) { return redimag(qfisqr(primeform_u(x, p))); }
5419 : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
5420 : static GEN
5421 2037145 : get_forms(GEN D, GEN *pL)
5422 : {
5423 2037145 : const long MAXFORM = 20;
5424 2037145 : GEN L, sqrtD = gsqrt(absi(D),DEFAULTPREC), forms = vectrunc_init(MAXFORM+1);
5425 2037145 : long s, nforms = 0;
5426 : ulong p;
5427 : forprime_t S;
5428 2037145 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
5429 2037145 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
5430 2037145 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
5431 2037145 : if (s < 10) s = 200;
5432 1887131 : else if (s < 20) s = 1000;
5433 1477 : else if (s < 5000) s = 5000;
5434 2037145 : u_forprime_init(&S, 2, s);
5435 341758226 : while ( (p = u_forprime_next(&S)) )
5436 : {
5437 337683936 : long d, k = kroiu(D,p);
5438 : pari_sp av2;
5439 337683936 : if (!k) continue;
5440 335470522 : if (k > 0)
5441 : {
5442 168268695 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
5443 168268695 : d = p-1;
5444 : }
5445 : else
5446 167201827 : d = p+1;
5447 335470522 : av2 = avma; affrr(divru(mulur(p,L),d), L); avma = av2;
5448 : }
5449 2037145 : *pL = L; return forms;
5450 : }
5451 :
5452 : /* h ~ #G, return o = order of f, set fao = its factorization */
5453 : static GEN
5454 2037194 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
5455 : {
5456 2037194 : long s = minss(itos(sqrti(h)), 10000);
5457 2037194 : GEN T = gen_Shanks_init(f, s, E, &qfi_group);
5458 2037194 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
5459 2037194 : return find_order(E, f, addiu(v,1), pfao);
5460 : }
5461 :
5462 : /* if g = 1 in G/<f> ? */
5463 : static int
5464 518 : equal1(void *E, GEN T, ulong N, GEN g)
5465 518 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
5466 :
5467 : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
5468 : * FIXME: should be gen_order, but equal1 has the wrong prototype */
5469 : static GEN
5470 112 : relative_order(void *E, GEN a, GEN o, ulong N, GEN T)
5471 : {
5472 112 : pari_sp av = avma;
5473 : long i, l;
5474 : GEN m;
5475 :
5476 112 : m = get_arith_ZZM(o);
5477 112 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
5478 112 : o = gel(m,1);
5479 112 : m = gel(m,2); l = lgcols(m);
5480 322 : for (i = l-1; i; i--)
5481 : {
5482 210 : GEN t, y, p = gcoeff(m,i,1);
5483 210 : long j, e = itos(gcoeff(m,i,2));
5484 210 : if (l == 2) {
5485 35 : t = gen_1;
5486 35 : y = a;
5487 : } else {
5488 175 : t = diviiexact(o, powiu(p,e));
5489 175 : y = powgi(a, t);
5490 : }
5491 210 : if (equal1(E, T,N,y)) o = t;
5492 : else {
5493 126 : for (j = 1; j < e; j++)
5494 : {
5495 28 : y = powgi(y, p);
5496 28 : if (equal1(E, T,N,y)) break;
5497 : }
5498 119 : if (j < e) {
5499 21 : if (j > 1) p = powiu(p, j);
5500 21 : o = mulii(t, p);
5501 : }
5502 : }
5503 : }
5504 112 : return gerepilecopy(av, o);
5505 : }
5506 :
5507 : /* h(x) for x<0 using Baby Step/Giant Step.
5508 : * Assumes G is not too far from being cyclic.
5509 : *
5510 : * Compute G^2 instead of G so as to kill most of the non-cyclicity */
5511 : GEN
5512 2039521 : classno(GEN x)
5513 : {
5514 2039521 : pari_sp av = avma;
5515 : long r2, k, s, i, l;
5516 : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
5517 : void *E;
5518 :
5519 2039521 : if (signe(x) >= 0) return classno2(x);
5520 :
5521 2039514 : check_quaddisc(x, &s, &k, "classno");
5522 2039514 : if (abscmpiu(x,12) <= 0) return gen_1;
5523 :
5524 2037145 : Hf = conductor_part(x, k, &D, NULL);
5525 2037145 : if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
5526 2037145 : forms = get_forms(D, &L);
5527 2037145 : r2 = two_rank(D);
5528 2037145 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
5529 :
5530 2037145 : l = lg(forms);
5531 2037145 : order_bound = const_vec(l-1, NULL);
5532 2037145 : E = expi(D) > 60? (void*)sqrtnint(shifti(absi(D),-2),4): NULL;
5533 2037145 : g1 = gel(forms,1);
5534 2037145 : gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
5535 2037145 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
5536 2037145 : d2 = NULL; /* not computed yet */
5537 2037145 : if (is_pm1(q)) goto END;
5538 505347 : for (i=2; i < l; i++)
5539 : {
5540 498621 : GEN o, fao, a, F, fd, f = gel(forms,i);
5541 498621 : fd = powgi(f, d1); if (is_pm1(gel(fd,1))) continue;
5542 359404 : F = powgi(fd, q);
5543 359404 : a = gel(F,1);
5544 359404 : o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
5545 : /* f^(d1 q) = 1 */
5546 359404 : fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
5547 359404 : o = find_order(E, f, fao, &fao);
5548 359404 : gel(order_bound,i) = o;
5549 : /* o = order of f, fao = factor(o) */
5550 359404 : update_g1(&g1,&d1,&fad1, f,o,fao);
5551 359404 : q = diviiround(hin, d1);
5552 359404 : if (is_pm1(q)) goto END;
5553 : }
5554 : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
5555 6726 : if (expi(q) > 3)
5556 : { /* q large: compute d2, 2nd elt divisor */
5557 70 : ulong N, n = 2*itou(sqrti(d1));
5558 70 : GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
5559 70 : d2 = gen_1;
5560 70 : N = itou( gceil(gdivgs(d1,n)) ); /* order(g1) <= n*N */
5561 287 : for (i = 1; i < l; i++)
5562 : {
5563 280 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
5564 280 : if (!B) B = find_order(E, f, fad1, /*junk*/&d);
5565 280 : f = powgi(f,d2);
5566 280 : if (equal1(E, T,N,f)) continue;
5567 112 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
5568 : /* f^B = 1 */
5569 112 : d = relative_order(E, f, B, N,T);
5570 112 : d2= mulii(d,d2);
5571 112 : D = mulii(d1,d2);
5572 112 : q = diviiround(hin,D);
5573 112 : if (is_pm1(q)) { d1 = D; goto END; }
5574 : }
5575 : /* very probably, d2 is the 2nd elementary divisor */
5576 7 : d1 = D; /* product of first two elt divisors */
5577 : }
5578 : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
5579 : * 2-rank */
5580 6663 : if (!ok_q(q,d1,d2,r2))
5581 : {
5582 0 : GEN q0 = q;
5583 : long d;
5584 0 : if (cmpii(mulii(q,d1), hin) < 0)
5585 : { /* try q = q0+1,-1,+2,-2 */
5586 0 : d = 1;
5587 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
5588 : }
5589 : else
5590 : { /* q0-1,+1,-2,+2 */
5591 0 : d = -1;
5592 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
5593 : }
5594 : }
5595 6663 : d1 = mulii(d1,q);
5596 :
5597 : END:
5598 2037145 : return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
5599 : }
5600 :
5601 : GEN
5602 0 : quadclassno(GEN x)
5603 : {
5604 0 : pari_sp av = avma;
5605 : GEN Hf, D;
5606 : long s, r;
5607 0 : check_quaddisc(x, &s, &r, "quadclassno");
5608 0 : if (s < 0 && cmpiu(x,12) <= 0) return gen_1;
5609 0 : Hf = conductor_part(x, r, &D, NULL);
5610 0 : return gerepileuptoint(av, mulii(Hf, gel(quadclassunit0(D,0,NULL,0),1)));
5611 : }
5612 :
5613 : /* use Euler products */
5614 : GEN
5615 21 : classno2(GEN x)
5616 : {
5617 21 : pari_sp av = avma;
5618 21 : const long prec = DEFAULTPREC;
5619 : long n, i, r, s;
5620 : GEN p1, p2, S, p4, p5, p7, Hf, Pi, reg, logd, d, dr, D, half;
5621 :
5622 21 : check_quaddisc(x, &s, &r, "classno2");
5623 21 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
5624 :
5625 21 : Hf = conductor_part(x, r, &D, ®);
5626 21 : if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
5627 :
5628 21 : Pi = mppi(prec);
5629 21 : d = absi(D); dr = itor(d, prec);
5630 21 : logd = logr_abs(dr);
5631 21 : p1 = sqrtr(divrr(mulir(d,logd), gmul2n(Pi,1)));
5632 21 : if (s > 0)
5633 : {
5634 14 : GEN invlogd = invr(logd);
5635 14 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
5636 14 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
5637 : }
5638 21 : n = itos_or_0( mptrunc(p1) );
5639 21 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
5640 :
5641 21 : p4 = divri(Pi,d);
5642 21 : p7 = invr(sqrtr_abs(Pi));
5643 21 : half = real2n(-1, prec);
5644 21 : if (s > 0)
5645 : { /* i = 1, shortcut */
5646 14 : p1 = sqrtr_abs(dr);
5647 14 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5648 14 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
5649 546 : for (i=2; i<=n; i++)
5650 : {
5651 532 : long k = kroiu(D,i); if (!k) continue;
5652 434 : p2 = mulir(sqru(i), p4);
5653 434 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5654 434 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
5655 434 : S = (k>0)? addrr(S,p5): subrr(S,p5);
5656 : }
5657 14 : S = shiftr(divrr(S,reg),-1);
5658 : }
5659 : else
5660 : { /* i = 1, shortcut */
5661 7 : p1 = gdiv(sqrtr_abs(dr), Pi);
5662 7 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5663 7 : S = addrr(p5, divrr(p1, mpexp(p4)));
5664 952 : for (i=2; i<=n; i++)
5665 : {
5666 945 : long k = kroiu(D,i); if (!k) continue;
5667 945 : p2 = mulir(sqru(i), p4);
5668 945 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5669 945 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
5670 945 : S = (k>0)? addrr(S,p5): subrr(S,p5);
5671 : }
5672 : }
5673 21 : return gerepileuptoint(av, mulii(Hf, roundr(S)));
5674 : }
5675 :
5676 : /* 1 + q + ... + q^v, v > 0 */
5677 : static GEN
5678 120 : geomsumu(ulong q, long v)
5679 : {
5680 120 : GEN u = utoipos(1+q);
5681 120 : for (; v > 1; v--) u = addui(1, mului(q, u));
5682 120 : return u;
5683 : }
5684 : static GEN
5685 120 : geomsum(GEN q, long v)
5686 : {
5687 : GEN u;
5688 120 : if (lgefint(q) == 3) return geomsumu(q[2], v);
5689 0 : u = addiu(q,1);
5690 0 : for (; v > 1; v--) u = addui(1, mulii(q, u));
5691 0 : return u;
5692 : }
5693 :
5694 : static GEN
5695 7803 : hclassno6_large(GEN x)
5696 : {
5697 : long i, l, s, xmod4;
5698 : GEN Q, H, D, P, E;
5699 :
5700 7803 : x = negi(x);
5701 7803 : check_quaddisc(x, &s, &xmod4, "hclassno");
5702 7803 : corediscfact(x, xmod4, &D, &P, &E);
5703 :
5704 7803 : Q = quadclassunit0(D, 0, NULL, 0);
5705 7803 : H = gel(Q,1); l = lg(P);
5706 :
5707 : /* H \prod_{p^e||f} (1 + (p^e-1)/(p-1))[ p - (D/p) ] */
5708 33647 : for (i=1; i<l; i++)
5709 : {
5710 25844 : long e = E[i], s;
5711 : GEN p, t;
5712 25844 : if (!e) continue;
5713 4996 : p = gel(P,i); s = kronecker(D,p);
5714 4996 : if (e == 1) t = addiu(p, 1-s);
5715 1000 : else if (s == 1) t = powiu(p,e);
5716 120 : else t = addui(1, mulii(subis(p, s), geomsum(p,e-1)));
5717 4996 : H = mulii(H,t);
5718 : }
5719 7803 : switch( itou_or_0(D) )
5720 : {
5721 0 : case 3: H = shifti(H,1);break;
5722 0 : case 4: H = muliu(H,3); break;
5723 7803 : default:H = muliu(H,6); break;
5724 : }
5725 7803 : return H;
5726 : }
5727 :
5728 : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
5729 : GEN
5730 121520 : hclassno6(GEN x)
5731 : {
5732 121520 : ulong d = itou_or_0(x);
5733 121520 : if (!d || d > 500000) return hclassno6_large(x);
5734 113717 : return utoipos(hclassno6u(d));
5735 : }
5736 :
5737 : GEN
5738 45738 : hclassno(GEN x)
5739 : {
5740 : long a, s;
5741 45738 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
5742 45738 : s = signe(x);
5743 45738 : if (s < 0) return gen_0;
5744 45738 : if (!s) return gdivgs(gen_1, -12);
5745 45738 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
5746 45738 : return gdivgs(hclassno6(x), 6);
5747 : }
5748 : /******************************************************************/
5749 : /* */
5750 : /* RAMANUJAN's TAU FUNCTION */
5751 : /* */
5752 : /******************************************************************/
5753 : /* 4|N > 0, not fundamental at 2; 6 * Hurwitz class number in level 2,
5754 : * equal to 6*(H(N)+2H(N/4)), H=qfbhclassno */
5755 : static GEN
5756 36750 : Hspec(GEN N)
5757 : {
5758 36750 : long v2 = Z_lvalrem(N, 2, &N), v2f = v2 >> 1;
5759 : GEN t;
5760 36750 : if (odd(v2)) { v2f--; N = shifti(N,3); }
5761 32557 : else if (mod4(N)!=3) { v2f--; N = shifti(N,2); }
5762 : /* N fundamental at 2, v2f = v2(f) s.t. N = f^2 D, D fundamental */
5763 36750 : t = addui(3, muliu(subiu(int2n(v2f+1), 3), 2 - kroiu(N,2)));
5764 36750 : return mulii(t, hclassno6(N));
5765 : }
5766 :
5767 : /* Ramanujan tau function for p prime */
5768 : static GEN
5769 14903 : tauprime(GEN p)
5770 : {
5771 14903 : pari_sp av = avma, av2;
5772 : GEN s, p2, p2_7, p_9, T;
5773 : ulong lim, t, tin;
5774 :
5775 14903 : if (absequaliu(p, 2)) return utoineg(24);
5776 : /* p > 2 */
5777 11396 : p2 = sqri(p);
5778 11396 : p2_7 = mului(7, p2);
5779 11396 : p_9 = mului(9, p);
5780 11396 : av2 = avma;
5781 11396 : lim = itou(sqrtint(p));
5782 11396 : tin = mod4(p) == 3? 1: 0;
5783 11396 : s = gen_0;
5784 87178 : for (t = 1; t <= lim; ++t)
5785 : {
5786 75782 : GEN h, a, t2 = sqru(t), D = shifti(subii(p, t2), 2); /* 4(p-t^2) */
5787 : /* t mod 2 != tin <=> D not fundamental at 2 */
5788 75782 : h = ((t&1UL) == tin)? hclassno6(D): Hspec(D);
5789 75782 : a = mulii(powiu(t2,3), addii(p2_7, mulii(t2, subii(shifti(t2,2), p_9))));
5790 75782 : s = addii(s, mulii(a,h));
5791 75782 : if (!(t & 255)) s = gerepileuptoint(av2, s);
5792 : }
5793 : /* 28p^3 - 28p^2 - 90p - 35 */
5794 11396 : T = subii(shifti(mulii(p2_7, subiu(p,1)), 2), addiu(mului(90,p), 35));
5795 11396 : s = shifti(diviuexact(s, 3), 6);
5796 11396 : return gerepileuptoint(av, subii(mulii(mulii(p2,p),T), addui(1, s)));
5797 : }
5798 :
5799 : /* Ramanujan tau function, return 0 for <= 0 */
5800 : GEN
5801 7035 : ramanujantau(GEN n)
5802 : {
5803 7035 : pari_sp ltop = avma;
5804 : GEN T, F, P, E;
5805 : long j, lP;
5806 :
5807 7035 : if (!(F = check_arith_all(n,"ramanujantau")))
5808 : {
5809 7014 : if (signe(n) <= 0) return gen_0;
5810 7007 : F = Z_factor(n);
5811 : }
5812 : else
5813 : {
5814 21 : P = gel(F,1);
5815 21 : if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
5816 : }
5817 :
5818 7014 : P = gel(F,1);
5819 7014 : E = gel(F,2); lP = lg(P);
5820 7014 : T = gen_1;
5821 21917 : for (j = 1; j < lP; j++)
5822 : {
5823 14903 : GEN p = gel(P,j), tp = tauprime(p), t1 = tp, t0 = gen_1;
5824 14903 : long k, e = itou(gel(E,j));
5825 20160 : for (k = 1; k < e; k++)
5826 : {
5827 5257 : GEN t2 = subii(mulii(tp, t1), mulii(powiu(p, 11), t0));
5828 5257 : t0 = t1; t1 = t2;
5829 : }
5830 14903 : T = mulii(T, t1);
5831 : }
5832 7014 : return gerepileuptoint(ltop, T);
5833 : }
|
