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 41467 : 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 41466 : 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 34 : is_gener_expo(GEN p, GEN L0)
40 : {
41 34 : GEN L, q = shifti(p,-1);
42 : long i, l;
43 34 : if (L0) {
44 14 : l = lg(L0);
45 14 : L = cgetg(l, t_VEC);
46 : } else {
47 20 : L0 = L = odd_prime_divisors(q);
48 20 : l = lg(L);
49 : }
50 34 : for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
51 34 : return L;
52 : }
53 : static GEN
54 41549 : u_is_gener_expo(ulong p, GEN L0)
55 : {
56 41549 : const ulong q = p >> 1;
57 : long i, l;
58 : GEN L;
59 41549 : if (L0) {
60 2525 : l = lg(L0);
61 2525 : L = cgetg(l, t_VECSMALL);
62 : } else {
63 39024 : L0 = L = u_odd_prime_divisors(q);
64 39017 : l = lg(L);
65 : }
66 41548 : for (i=1; i<l; i++) L[i] = q / uel(L0,i);
67 41548 : return L;
68 : }
69 :
70 : int
71 110184 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
72 : {
73 : long i;
74 110184 : if (krouu(x, p) >= 0) return 0;
75 121688 : for (i=lg(L)-1; i; i--)
76 : {
77 77488 : ulong t = Fl_powu(x, uel(L,i), p);
78 77505 : if (t == p_1 || t == 1) return 0;
79 : }
80 44200 : return 1;
81 : }
82 : /* assume p prime */
83 : ulong
84 135358 : pgener_Fl_local(ulong p, GEN L0)
85 : {
86 135358 : const pari_sp av = avma;
87 135358 : const ulong p_1 = p-1;
88 : long x;
89 : GEN L;
90 135358 : if (p <= 19) switch(p)
91 : { /* quick trivial cases */
92 21 : case 2: return 1;
93 : case 7:
94 17170 : case 17: return 3;
95 76652 : default: return 2;
96 : }
97 41515 : L = u_is_gener_expo(p,L0);
98 41518 : for (x=2;;x++) { if (is_gener_Fl(x,p,p_1,L)) break; }
99 41520 : avma = av; return x;
100 : }
101 : ulong
102 112918 : 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 142 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
108 : {
109 142 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
110 142 : if (t >= 0) return 0;
111 373 : for (i = lg(L)-1; i; i--)
112 : {
113 258 : GEN t = Fp_pow(x, gel(L,i), p);
114 258 : if (equalii(t, p_1) || equali1(t)) return 0;
115 : }
116 115 : return 1;
117 : }
118 :
119 : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
120 : GEN
121 23975 : pgener_Fp_local(GEN p, GEN L0)
122 : {
123 23975 : pari_sp av0 = avma;
124 : GEN x, p_1, L;
125 23975 : if (lgefint(p) == 3)
126 : {
127 : ulong z;
128 23946 : if (p[2] == 2) return gen_1;
129 18017 : if (L0) L0 = ZV_to_nv(L0);
130 18017 : z = pgener_Fl_local(uel(p,2), L0);
131 18017 : avma = av0; return utoipos(z);
132 : }
133 29 : p_1 = subiu(p,1); L = is_gener_expo(p, L0);
134 29 : x = utoipos(2);
135 32 : for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
136 29 : avma = av0; return utoipos(uel(x,2));
137 : }
138 :
139 : GEN
140 23422 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
141 :
142 : ulong
143 70541 : pgener_Zl(ulong p)
144 : {
145 70541 : 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 70541 : if (p == 40487) return 10;
148 : #ifndef LONG_IS_64BIT
149 10073 : return pgener_Fl(p);
150 : #else
151 60468 : 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 70546 : pgener_Zp(GEN p)
168 : {
169 70546 : 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 7 : rootsof1u_Fp(ulong n, GEN p)
246 : {
247 7 : pari_sp av = avma;
248 7 : GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
249 7 : z = pgener_Fp_local(p, Flv_to_ZV(L));
250 7 : z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
251 7 : return gerepileuptoint(av, z);
252 : }
253 :
254 : ulong
255 2436 : rootsof1_Fl(ulong n, ulong p)
256 : {
257 2436 : pari_sp av = avma;
258 2436 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
259 2436 : ulong z = pgener_Fl_local(p, L);
260 2436 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
261 2436 : 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 251937 : ulogintall(ulong B, ulong y, ulong *ptq)
375 : {
376 : ulong r, r2;
377 : long e;
378 :
379 251937 : if (y == 2)
380 : {
381 4810 : long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
382 4810 : if (ptq) *ptq = 1UL << eB;
383 4810 : return eB;
384 : }
385 247127 : r = y, r2 = 1UL;
386 898696 : for (e=1;; e++)
387 : { /* here, r = y^e, r2 = y^(e-1) */
388 898696 : if (r >= B)
389 : {
390 247012 : if (r != B) { e--; r = r2; }
391 247012 : if (ptq) *ptq = r;
392 247012 : return e;
393 : }
394 651684 : r2 = r;
395 651684 : r = umuluu_or_0(y, r);
396 651684 : if (!r)
397 : {
398 115 : if (ptq) *ptq = r2;
399 115 : return e;
400 : }
401 651569 : }
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 259505 : logintall(GEN B, GEN y, GEN *ptq)
408 : {
409 : pari_sp av;
410 259505 : long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
411 : GEN q, pow2;
412 :
413 259505 : if (lgefint(B) == 3)
414 : {
415 : ulong q;
416 251937 : if (lgefint(y) > 3)
417 : {
418 0 : if (ptq) *ptq = gen_1;
419 0 : return 0;
420 : }
421 251937 : if (!ptq) return ulogintall(B[2], y[2], NULL);
422 46569 : e = ulogintall(B[2], y[2], &q);
423 46569 : *ptq = utoi(q); return e;
424 : }
425 7568 : if (equaliu(y,2))
426 : {
427 146 : if (ptq) *ptq = int2n(eB);
428 146 : return eB;
429 : }
430 7422 : av = avma;
431 7422 : ey = expi(y);
432 : /* eB/(ey+1) - 1 < e <= eB/ey */
433 7422 : emax = eB/ey;
434 7422 : if (emax <= 13) /* e small, be naive */
435 : {
436 1507 : GEN r = y, r2 = gen_1;
437 16022 : for (e=1;; e++)
438 : { /* here, r = y^e, r2 = y^(e-1) */
439 16022 : long fl = cmpii(r, B);
440 16022 : if (fl >= 0)
441 : {
442 1507 : if (fl) { e--; cgiv(r); r = r2; }
443 1507 : if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
444 1507 : return e;
445 : }
446 14515 : r2 = r; r = mulii(r,y);
447 14515 : }
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 5915 : pow2 = new_chunk((long)log2(eB)+2);
454 5915 : gel(pow2,0) = y;
455 5915 : for (i=0, q=y;; )
456 : {
457 34225 : GEN r = gel(pow2,i); /* r = y^2^i */
458 34225 : long fl = cmpii(r,B);
459 34225 : 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 34225 : if (fl > 0) { i--; break; }
466 32323 : q = r;
467 32323 : if (1L<<(i+1) > emax) break;
468 28310 : gel(pow2,++i) = sqri(q);
469 28310 : }
470 :
471 5915 : for (e = 1L<<i;;)
472 : { /* y^e = q < B < r = q * y^(2^i) */
473 32302 : pari_sp av2 = avma;
474 : long fl;
475 : GEN r;
476 32302 : if (--i < 0) break;
477 26394 : r = mulii(q, gel(pow2,i));
478 26394 : fl = cmpii(r, B);
479 26394 : if (fl > 0) avma = av2;
480 : else
481 : {
482 12704 : e += (1L<<i);
483 12704 : q = r;
484 12704 : if (!fl) break; /* B = r */
485 : }
486 26387 : }
487 5915 : if (ptq) *ptq = gerepileuptoint(av, q); else avma = av;
488 5915 : 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 39439 : sqrtint(GEN a)
508 : {
509 39439 : if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
510 39439 : switch (signe(a))
511 : {
512 39425 : 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 13525781 : carremod(ulong A)
526 : {
527 13525781 : 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 13525781 : 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 13525781 : 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 13525781 : const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
537 27051562 : return (carresmod64[A & 0x3fUL]
538 5364125 : && carresmod63[A % 63UL]
539 3192422 : && carresmod65[A % 65UL]
540 16171710 : && carresmod11[A % 11UL]);
541 : }
542 :
543 : /* emulate Z_issquareall on single-word integers */
544 : long
545 12036236 : uissquareall(ulong A, ulong *sqrtA)
546 : {
547 12036236 : if (!A) { *sqrtA = 0; return 1; }
548 12036236 : if (carremod(A))
549 : {
550 1796632 : ulong a = usqrt(A);
551 1796577 : if (a * a == A) { *sqrtA = a; return 1; }
552 : }
553 10330448 : return 0;
554 : }
555 : long
556 120651 : uissquare(ulong A)
557 : {
558 120651 : if (!A) return 1;
559 120651 : if (carremod(A))
560 : {
561 3469 : ulong a = usqrt(A);
562 3469 : if (a * a == A) return 1;
563 : }
564 117206 : return 0;
565 : }
566 :
567 : long
568 6180825 : Z_issquareall(GEN x, GEN *pt)
569 : {
570 : pari_sp av;
571 : GEN y, r;
572 :
573 6180825 : switch(signe(x))
574 : {
575 2101314 : case -1: return 0;
576 1064 : case 0: if (pt) *pt=gen_0; return 1;
577 : }
578 4078447 : if (lgefint(x) == 3)
579 : {
580 2709576 : ulong u = uel(x,2), a;
581 2709576 : if (!pt) return uissquare(u);
582 2588925 : if (!uissquareall(u, &a)) return 0;
583 1377119 : *pt = utoipos(a); return 1;
584 : }
585 1368871 : if (!carremod(umodiu(x, 64*63*65*11))) return 0;
586 607918 : av = avma; y = sqrtremi(x, &r);
587 607918 : if (r != gen_0) { avma = av; return 0; }
588 17119 : if (pt) { *pt = y; avma = (pari_sp)y; } else avma = av;
589 17119 : 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 2275 : polissquareall(GEN x, GEN *pt)
608 : {
609 : pari_sp av;
610 : long v;
611 : GEN y, a, b, p;
612 :
613 2275 : if (!signe(x))
614 : {
615 7 : if (pt) *pt = gcopy(x);
616 7 : return 1;
617 : }
618 2268 : if (odd(degpol(x))) return 0; /* odd degree */
619 1477 : av = avma;
620 1477 : v = RgX_valrem(x, &x);
621 1477 : if (v & 1) { avma = av; return 0; }
622 1470 : a = gel(x,2); /* test constant coeff */
623 1470 : if (!pt)
624 70 : { if (!issquare(a)) { avma = av; return 0; } }
625 : else
626 1400 : { if (!issquareall(a,&b)) { avma = av; return 0; } }
627 1470 : if (!degpol(x)) { /* constant polynomial */
628 77 : if (!pt) { avma = av; return 1; }
629 35 : y = scalarpol(b, varn(x)); goto END;
630 : }
631 1393 : p = characteristic(x);
632 1393 : 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 1358 : long m = 1;
656 1358 : x = RgX_Rg_div(x,a);
657 : /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
658 1358 : if (!signe(p)) x = RgX_deflate_max(x,&m);
659 1358 : y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
660 1939 : if (!RgX_equal(RgX_sqr(y), x)) { avma = av; return 0; }
661 784 : if (!pt) { avma = av; return 1; }
662 777 : if (!gequal1(a)) y = gmul(b, y);
663 777 : if (m != 1) y = RgX_inflate(y,m);
664 : }
665 : END:
666 826 : if (v) y = RgX_shift_shallow(y, v>>1);
667 826 : *pt = gerepilecopy(av, y); return 1;
668 : }
669 :
670 : /* b unit mod p */
671 : static int
672 350 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
673 : {
674 350 : if (d == 1)
675 : { /* mod p: faster */
676 266 : if (!Fp_ispower(b, K, p)) return 0;
677 266 : 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 329 : 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 490 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
691 : {
692 : GEN t, A;
693 490 : long v = Z_pvalrem(*pa, p, &A), d = e - v;
694 490 : if (d <= 0) t = gen_0;
695 : else
696 : {
697 : ulong r;
698 434 : v = udivui_rem(v, K, &r);
699 434 : if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
700 329 : if (L && v) t = mulii(t, powiu(p, v));
701 : }
702 385 : if (q) *pa = modii(*pa, q);
703 385 : if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
704 385 : return 1;
705 : }
706 : long
707 385 : 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 385 : 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 371 : av = avma;
725 :
726 371 : 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 371 : if (!mod2(K)
740 203 : && kronecker(a, shifti(q,-vali(q))) == -1) { avma = av; return 0; }
741 364 : L = pt? vectrunc_init(expi(q)+1): NULL;
742 364 : u_forprime_init(&S, 2, tridiv_bound(q));
743 364 : while ((pp = u_forprime_next(&S)))
744 : {
745 : int stop;
746 883428 : e = Z_lvalrem_stop(&q, pp, &stop);
747 883428 : if (!e) continue;
748 266 : if (!handle_pe(&a, q, L, K, utoipos(pp), e)) { avma = av; return 0; }
749 224 : if (stop)
750 : {
751 189 : if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
752 189 : 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 259 : if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
805 259 : return 1;
806 : }
807 :
808 : static long
809 126 : polmodispower(GEN x, GEN K, GEN *pt)
810 : {
811 126 : pari_sp av = avma;
812 126 : GEN p = NULL, T = NULL;
813 126 : if (Rg_is_FpXQ(x, &T,&p) && p)
814 : {
815 112 : x = liftall_shallow(x);
816 112 : if (!Fq_ispower(x, K, T, p)) { avma = av; return 0; }
817 105 : if (!pt) { avma = av; return 1; }
818 98 : x = Fq_sqrtn(x, K, T,p, NULL);
819 98 : if (typ(x) == t_INT)
820 7 : x = Fp_to_mod(x,p);
821 : else
822 91 : x = mkpolmod(FpX_to_mod(x,p), FpX_to_mod(T,p));
823 98 : *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 147406 : issquareall(GEN x, GEN *pt)
831 : {
832 147406 : long tx = typ(x);
833 : GEN F;
834 : pari_sp av;
835 :
836 147406 : if (!pt) return issquare(x);
837 4564 : switch(tx)
838 : {
839 1925 : case t_INT: return Z_issquareall(x, pt);
840 259 : case t_FRAC: av = avma;
841 259 : F = cgetg(3, t_FRAC);
842 259 : 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 2191 : 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 63 : if (!issquare(x)) return 0;
857 63 : *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 143122 : issquare(GEN x)
871 : {
872 : pari_sp av;
873 : GEN a, p;
874 : long i, v;
875 :
876 143122 : switch(typ(x))
877 : {
878 : case t_INT:
879 142716 : return Z_issquare(x);
880 :
881 : case t_REAL:
882 14 : 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 63 : return 1;
894 :
895 : case t_PADIC:
896 105 : a = gel(x,4); if (!signe(a)) return 1;
897 105 : if (valp(x)&1) return 0;
898 91 : p = gel(x,2);
899 91 : 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 497 : polispower(GEN x, GEN K, GEN *pt)
975 : {
976 : pari_sp av;
977 497 : long v, d, k = itos(K);
978 : GEN y, a, b;
979 :
980 497 : if (!signe(x))
981 : {
982 7 : if (pt) *pt = gcopy(x);
983 7 : return 1;
984 : }
985 490 : if (degpol(x) % k) return 0; /* degree not multiple of k */
986 483 : av = avma;
987 483 : y = NULL; /*-Wall*/
988 483 : v = RgX_valrem(x, &x);
989 483 : if (v % k) return 0;
990 476 : v /= k;
991 476 : a = gel(x,2); b = NULL;
992 476 : if (!ispower(a, K, &b)) { avma = av; return 0; }
993 462 : d = degpol(x);
994 462 : if (d)
995 : {
996 329 : GEN p = characteristic(x);
997 329 : a = leading_coeff(x);
998 357 : if (!ispower(a, K, &b)) { avma = av; return 0; }
999 329 : x = RgX_normalize(x);
1000 329 : if (signe(p))
1001 : {
1002 168 : GEN T0, T = NULL;
1003 168 : if (!BPSW_isprime(p))
1004 0 : pari_err_IMPL("ispower in non-prime characteristic");
1005 168 : if (RgX_is_FpXQX(x,&T,&p))
1006 : { /* over Fq */
1007 168 : T0 = T;
1008 168 : if (T && typ(T) == t_FFELT) T = FF_mod(T);
1009 168 : x = RgX_to_FqX(x,T,p);
1010 168 : if (!FqX_ispower(x, k, T,p, pt)) { avma = av; return 0; }
1011 140 : if (pt)
1012 : {
1013 140 : y = *pt;
1014 140 : if (!T) y = FpX_to_mod(y, p);
1015 105 : else if (typ(T0) == t_FFELT)
1016 70 : y = FqX_to_FFX(y, T0);
1017 : else
1018 : {
1019 35 : T = FpX_to_mod(T, p);
1020 35 : y = gmul(y, gmodulsg(1,T));
1021 : }
1022 : }
1023 140 : goto END;
1024 : }
1025 0 : if (cmpii(p,K) <= 0)
1026 0 : pari_err_IMPL("ispower(general t_POL) in small characteristic");
1027 : }
1028 161 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
1029 161 : if (!RgX_equal(powgi(y, K), x)) { avma = av; return 0; }
1030 : }
1031 : else
1032 133 : y = pol_1(varn(x));
1033 : END:
1034 434 : if (pt)
1035 : {
1036 434 : if (!gequal1(a))
1037 : {
1038 14 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
1039 14 : y = gmul(b,y);
1040 : }
1041 434 : if (v) y = RgX_shift_shallow(y, v);
1042 434 : *pt = gerepilecopy(av, y);
1043 : }
1044 0 : else avma = av;
1045 434 : return 1;
1046 : }
1047 :
1048 : long
1049 1400 : Z_ispowerall(GEN x, ulong k, GEN *pt)
1050 : {
1051 1400 : long s = signe(x);
1052 : ulong mask;
1053 1400 : if (!s) { if (pt) *pt = gen_0; return 1; }
1054 1400 : if (s > 0) {
1055 1246 : if (k == 2) return Z_issquareall(x, pt);
1056 875 : if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
1057 161 : if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
1058 147 : if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
1059 140 : return is_kth_power(x, k, pt);
1060 : }
1061 154 : if (!odd(k)) return 0;
1062 140 : if (Z_ispowerall(absi(x), k, pt))
1063 : {
1064 140 : if (pt) *pt = negi(*pt);
1065 140 : return 1;
1066 : };
1067 0 : return 0;
1068 : }
1069 :
1070 : /* is x a K-th power mod p ? Assume p prime. */
1071 : int
1072 287 : Fp_ispower(GEN x, GEN K, GEN p)
1073 : {
1074 287 : pari_sp av = avma;
1075 : GEN p_1;
1076 : long r;
1077 287 : x = modii(x, p);
1078 287 : if (!signe(x) || equali1(x)) { avma = av; return 1; }
1079 : /* implies p > 2 */
1080 119 : p_1 = subiu(p,1);
1081 119 : K = gcdii(K, p_1);
1082 119 : if (absequaliu(K, 2)) { r = kronecker(x,p); avma = av; return (r > 0); }
1083 42 : x = Fp_pow(x, diviiexact(p_1,K), p);
1084 42 : avma = av; return equali1(x);
1085 : }
1086 :
1087 : /* x unit defined modulo 2^e, e > 0, p prime */
1088 : static int
1089 2373 : U2_issquare(GEN x, long e)
1090 : {
1091 2373 : long r = signe(x)>=0?mod8(x):8-mod8(x);
1092 2373 : if (e==1) return 1;
1093 2373 : if (e==2) return (r&3L) == 1;
1094 2009 : return r == 1;
1095 : }
1096 : /* x unit defined modulo p^e, e > 0, p prime */
1097 : static int
1098 4690 : Up_issquare(GEN x, GEN p, long e)
1099 4690 : { return (absequaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
1100 :
1101 : long
1102 2548 : Zn_issquare(GEN d, GEN fn)
1103 : {
1104 : long j, np;
1105 2548 : if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
1106 2548 : if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
1107 : /* integer factorization */
1108 2548 : np = nbrows(fn);
1109 5320 : for (j = 1; j <= np; ++j)
1110 : {
1111 4970 : GEN r, p = gcoeff(fn, j, 1);
1112 4970 : long e = itos(gcoeff(fn, j, 2));
1113 4970 : long v = Z_pvalrem(d,p,&r);
1114 4970 : if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
1115 : }
1116 350 : return 1;
1117 : }
1118 :
1119 : static long
1120 1113 : Qp_ispower(GEN x, GEN K, GEN *pt)
1121 : {
1122 1113 : pari_sp av = avma;
1123 1113 : GEN z = Qp_sqrtn(x, K, NULL);
1124 1113 : if (!z) { avma = av; return 0; }
1125 819 : if (pt) *pt = z;
1126 819 : return 1;
1127 : }
1128 :
1129 : long
1130 7003003 : ispower(GEN x, GEN K, GEN *pt)
1131 : {
1132 : GEN z;
1133 :
1134 7003003 : if (!K) return gisanypower(x, pt);
1135 2821 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
1136 2821 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
1137 2821 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
1138 2786 : switch(typ(x)) {
1139 : case t_INT:
1140 616 : return Z_ispowerall(x, itou(K), pt);
1141 : case t_FRAC:
1142 : {
1143 21 : GEN a = gel(x,1), b = gel(x,2);
1144 21 : ulong k = itou(K);
1145 21 : if (pt) {
1146 14 : z = cgetg(3, t_FRAC);
1147 14 : if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
1148 14 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1149 : }
1150 0 : avma = (pari_sp)(z + 3); return 0;
1151 : }
1152 7 : return Z_ispower(a, k) && Z_ispower(b, k);
1153 : }
1154 : case t_INTMOD:
1155 245 : return Zn_ispower(gel(x,2), gel(x,1), K, pt);
1156 : case t_FFELT:
1157 182 : return FF_ispower(x, K, pt);
1158 :
1159 : case t_PADIC:
1160 1113 : return Qp_ispower(x, K, pt);
1161 : case t_POLMOD:
1162 91 : return polmodispower(x, K, pt);
1163 : case t_POL:
1164 490 : return polispower(x, K, pt);
1165 : case t_RFRAC: {
1166 7 : GEN a = gel(x,1), b = gel(x,2);
1167 7 : if (pt) {
1168 7 : z = cgetg(3, t_RFRAC);
1169 7 : if (ispower(a, K, &a) && polispower(b, K, &b)) {
1170 7 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1171 : }
1172 0 : avma = (pari_sp)(z + 3); return 0;
1173 : }
1174 0 : return (ispower(a, K, NULL) && polispower(b, K, NULL));
1175 : }
1176 : case t_REAL:
1177 7 : if (signe(x) < 0 && !mpodd(K)) return 0;
1178 : case t_COMPLEX:
1179 14 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1180 14 : return 1;
1181 :
1182 : case t_SER:
1183 7 : if (signe(x) && (!dvdsi(valp(x), K) || !ispower(gel(x,2), K, NULL)))
1184 0 : return 0;
1185 7 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1186 7 : return 1;
1187 : }
1188 0 : pari_err_TYPE("ispower",x);
1189 : return 0; /* LCOV_EXCL_LINE */
1190 : }
1191 :
1192 : long
1193 7000182 : gisanypower(GEN x, GEN *pty)
1194 : {
1195 7000182 : long tx = typ(x);
1196 : ulong k, h;
1197 7000182 : if (tx == t_INT) return Z_isanypower(x, pty);
1198 14 : if (tx == t_FRAC)
1199 : {
1200 14 : pari_sp av = avma;
1201 14 : GEN fa, P, E, a = gel(x,1), b = gel(x,2);
1202 : long i, j, p, e;
1203 14 : int sw = (abscmpii(a, b) > 0);
1204 :
1205 14 : if (sw) swap(a, b);
1206 14 : k = Z_isanypower(a, pty? &a: NULL);
1207 14 : if (!k)
1208 : { /* a = -1,1 or not a pure power */
1209 7 : if (!is_pm1(a)) { avma = av; return 0; }
1210 7 : if (signe(a) < 0) b = negi(b);
1211 7 : k = Z_isanypower(b, pty? &b: NULL);
1212 7 : if (!k || !pty) { avma = av; return k; }
1213 7 : *pty = gerepileupto(av, ginv(b));
1214 7 : return k;
1215 : }
1216 7 : fa = factoru(k);
1217 7 : P = gel(fa,1);
1218 7 : E = gel(fa,2); h = k;
1219 14 : for (i = lg(P) - 1; i > 0; i--)
1220 : {
1221 7 : p = P[i];
1222 7 : e = E[i];
1223 21 : for (j = 0; j < e; j++)
1224 14 : if (!is_kth_power(b, p, &b)) break;
1225 7 : if (j < e) k /= upowuu(p, e - j);
1226 : }
1227 7 : if (k == 1) { avma = av; return 0; }
1228 7 : if (!pty) { avma = av; return k; }
1229 0 : if (k != h) a = powiu(a, h/k);
1230 0 : *pty = gerepilecopy(av, mkfrac(a, b));
1231 0 : return k;
1232 : }
1233 0 : pari_err_TYPE("gisanypower", x);
1234 : return 0; /* LCOV_EXCL_LINE */
1235 : }
1236 :
1237 : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
1238 : * No need to optimize for 2,3,5,7 powers (done before) */
1239 : static long
1240 505729 : split_exponent(ulong e, GEN *x)
1241 : {
1242 : GEN fa, P, E;
1243 505729 : long i, j, l, k = 1;
1244 505729 : if (e == 1) return 1;
1245 14 : fa = factoru(e);
1246 14 : P = gel(fa,1);
1247 14 : E = gel(fa,2); l = lg(P);
1248 28 : for (i = 1; i < l; i++)
1249 : {
1250 14 : ulong p = P[i];
1251 28 : for (j = 0; j < E[i]; j++)
1252 : {
1253 : GEN y;
1254 14 : if (!is_kth_power(*x, p, &y)) break;
1255 14 : k *= p; *x = y;
1256 : }
1257 : }
1258 14 : return k;
1259 : }
1260 :
1261 : static long
1262 864570 : Z_isanypower_nosmalldiv(GEN *px)
1263 : { /* any prime divisor of x is > 102 */
1264 864570 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
1265 864570 : const double LOG103 = 4.6347; /* lower bound for log(103) */
1266 : forprime_t T;
1267 864570 : ulong mask = 7, e2;
1268 : long k, ex;
1269 864570 : GEN y, x = *px;
1270 :
1271 864570 : k = 1;
1272 864570 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
1273 864570 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
1274 864570 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
1275 864570 : if (u_forprime_init(&T, 11, e2))
1276 : {
1277 16905 : GEN logx = NULL;
1278 16905 : const ulong Q = 30011; /* prime */
1279 : ulong p, xmodQ;
1280 16905 : double dlogx = 0;
1281 : /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
1282 : * for large p the modular checks are no longer competitively fast */
1283 33852 : while ( (ex = is_pth_power(x, &y, &T, 30)) )
1284 : {
1285 42 : k *= ex; x = y;
1286 42 : e2 = (ulong)((expi(x) + 1) / LOG2_103);
1287 42 : u_forprime_restrict(&T, e2);
1288 : }
1289 16905 : if (DEBUGLEVEL>4) err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x));
1290 16905 : xmodQ = umodiu(x, Q);
1291 : /* test Q | x, just in case */
1292 16905 : if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
1293 : /* x^(1/p) < 2^31 */
1294 16891 : p = T.p;
1295 16891 : if (p <= e2)
1296 : {
1297 16877 : logx = logr_abs( itor(x, DEFAULTPREC) );
1298 16877 : dlogx = rtodbl(logx);
1299 16877 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1300 : }
1301 152621 : while (p && p <= e2)
1302 : { /* is x a p-th power ? By computing y = round(x^(1/p)).
1303 : * Check whether y^p = x, first mod Q, then exactly. */
1304 118839 : pari_sp av = avma;
1305 : long e;
1306 118839 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
1307 118839 : ulong ymodQ = umodiu(y,Q);
1308 118839 : if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
1309 21 : || !equalii(powiu(y, p), x)) avma = av;
1310 : else
1311 : {
1312 21 : k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
1313 21 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1314 21 : u_forprime_restrict(&T, e2);
1315 21 : continue; /* if success, retry same p */
1316 : }
1317 118818 : p = u_forprime_next(&T);
1318 : }
1319 : }
1320 864556 : *px = x; return k;
1321 : }
1322 :
1323 : static ulong tinyprimes[] = {
1324 : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
1325 : 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
1326 : 157, 163, 167, 173, 179, 181, 191, 193, 197, 199
1327 : };
1328 :
1329 : /* disregard the sign of x, caller will take care of x < 0 */
1330 : static long
1331 7000819 : Z_isanypower_aux(GEN x, GEN *pty)
1332 : {
1333 : long ex, v, i, l, k;
1334 : GEN y, P, E;
1335 7000819 : ulong mask, e = 0;
1336 :
1337 7000819 : if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
1338 :
1339 7000805 : if (signe(x) < 0) x = negi(x);
1340 7000805 : k = l = 1;
1341 7000805 : P = cgetg(26 + 1, t_VECSMALL);
1342 7000805 : E = cgetg(26 + 1, t_VECSMALL);
1343 : /* trial division */
1344 122938942 : for(i = 0; i < 26; i++)
1345 : {
1346 60135033 : ulong p = tinyprimes[i];
1347 : int stop;
1348 60135033 : v = Z_lvalrem_stop(&x, p, &stop);
1349 60135033 : if (v)
1350 : {
1351 7922439 : P[l] = p;
1352 7922439 : E[l] = v; l++;
1353 13588806 : e = ugcd(e, v); if (e == 1) goto END;
1354 : }
1355 54716718 : if (stop) {
1356 248052 : if (is_pm1(x)) k = e;
1357 248052 : goto END;
1358 : }
1359 : }
1360 :
1361 1334438 : if (e)
1362 : { /* Bingo. Result divides e */
1363 : long v3, v5, v7;
1364 505715 : ulong e2 = e;
1365 505715 : v = u_lvalrem(e2, 2, &e2);
1366 505715 : if (v)
1367 : {
1368 375277 : for (i = 0; i < v; i++)
1369 : {
1370 374185 : if (!Z_issquareall(x, &y)) break;
1371 1302 : k <<= 1; x = y;
1372 : }
1373 : }
1374 505715 : mask = 0;
1375 505715 : v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
1376 505715 : v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
1377 505715 : v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
1378 1011507 : while ( (ex = is_357_power(x, &y, &mask)) ) {
1379 77 : x = y;
1380 77 : switch(ex)
1381 : {
1382 28 : case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
1383 28 : case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
1384 21 : case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
1385 : }
1386 : }
1387 505715 : k *= split_exponent(e2, &x);
1388 : }
1389 : else
1390 828723 : k = Z_isanypower_nosmalldiv(&x);
1391 : END:
1392 7000805 : if (pty && k != 1)
1393 : {
1394 8001 : if (e)
1395 : { /* add missing small factors */
1396 6881 : y = powuu(P[1], E[1] / k);
1397 6881 : for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
1398 6881 : x = equali1(x)? y: mulii(x,y);
1399 : }
1400 8001 : *pty = x;
1401 : }
1402 7000805 : return k == 1? 0: k;
1403 : }
1404 :
1405 : long
1406 7000819 : Z_isanypower(GEN x, GEN *pty)
1407 : {
1408 7000819 : pari_sp av = avma;
1409 7000819 : long k = Z_isanypower_aux(x, pty);
1410 7000819 : if (!k) { avma = av; return 0; }
1411 8064 : if (signe(x) < 0)
1412 : {
1413 42 : long v = vals(k);
1414 42 : if (v)
1415 : {
1416 : GEN y;
1417 28 : k >>= v;
1418 28 : if (k == 1) { avma = av; return 0; }
1419 21 : if (!pty) { avma = av; return k; }
1420 14 : y = *pty;
1421 14 : y = powiu(y, 1<<v);
1422 14 : togglesign(y);
1423 14 : *pty = gerepileuptoint(av, y);
1424 14 : return k;
1425 : }
1426 14 : if (pty) togglesign_safe(pty);
1427 : }
1428 8036 : if (pty) *pty = gerepilecopy(av, *pty); else avma = av;
1429 8036 : return k;
1430 : }
1431 :
1432 : /* Faster than */
1433 : /* !cmpii(n, int2n(vali(n))) */
1434 : /* !cmpis(shifti(n, -vali(n)), 1) */
1435 : /* expi(n) == vali(n) */
1436 : /* hamming(n) == 1 */
1437 : /* even for single-word values, and much faster for multiword values. */
1438 : /* If all you have is a word, you can just use n & !(n & (n-1)). */
1439 : long
1440 32933 : Z_ispow2(GEN n)
1441 : {
1442 : GEN xp;
1443 : long i, lx;
1444 : ulong u;
1445 32933 : if (signe(n) != 1) return 0;
1446 32926 : xp = int_LSW(n);
1447 32926 : lx = lgefint(n);
1448 32926 : u = *xp;
1449 32958 : for (i = 3; i < lx; ++i)
1450 : {
1451 30349 : if (u) return 0;
1452 32 : xp = int_nextW(xp);
1453 32 : u = *xp;
1454 : }
1455 2609 : return !(u & (u-1));
1456 : }
1457 :
1458 : static long
1459 842006 : isprimepower_i(GEN n, GEN *pt, long flag)
1460 : {
1461 842006 : pari_sp av = avma;
1462 : long i, v;
1463 :
1464 842006 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1465 842006 : if (signe(n) <= 0) return 0;
1466 :
1467 842006 : if (lgefint(n) == 3)
1468 : {
1469 : ulong p;
1470 541114 : v = uisprimepower(n[2], &p);
1471 541114 : if (v)
1472 : {
1473 54902 : if (pt) *pt = utoipos(p);
1474 54902 : return v;
1475 : }
1476 486212 : return 0;
1477 : }
1478 1662590 : for (i = 0; i < 26; i++)
1479 : {
1480 1626743 : ulong p = tinyprimes[i];
1481 1626743 : v = Z_lvalrem(n, p, &n);
1482 1626744 : if (v)
1483 : {
1484 265046 : avma = av;
1485 265046 : if (!is_pm1(n)) return 0;
1486 650 : if (pt) *pt = utoipos(p);
1487 647 : return v;
1488 : }
1489 : }
1490 : /* p | n => p >= 103 */
1491 35847 : v = Z_isanypower_nosmalldiv(&n); /* expensive */
1492 35847 : if (!(flag? isprime(n): BPSW_psp(n))) { avma = av; return 0; }
1493 5535 : if (pt) *pt = gerepilecopy(av, n); else avma = av;
1494 5535 : return v;
1495 : }
1496 : long
1497 840098 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
1498 : long
1499 1908 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1500 :
1501 : long
1502 541751 : uisprimepower(ulong n, ulong *pp)
1503 : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1504 : * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1505 541751 : const ulong CUTOFF = 200UL;
1506 541751 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1507 541751 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1508 : #ifdef LONG_IS_64BIT
1509 : /* primes preceeding the appropriate root of ULONG_MAX. */
1510 481512 : const ulong ROOT9 = 137;
1511 481512 : const ulong ROOT8 = 251;
1512 481512 : const ulong ROOT7 = 563;
1513 481512 : const ulong ROOT5 = 7129;
1514 481512 : const ulong ROOT4 = 65521;
1515 : #else
1516 60239 : const ulong ROOT9 = 11;
1517 60239 : const ulong ROOT8 = 13;
1518 60239 : const ulong ROOT7 = 23;
1519 60239 : const ulong ROOT5 = 83;
1520 60239 : const ulong ROOT4 = 251;
1521 : #endif
1522 : ulong mask;
1523 : long v, i;
1524 : int e;
1525 541751 : if (n < 2) return 0;
1526 541737 : if (!odd(n)) {
1527 270718 : if (n & (n-1)) return 0;
1528 914 : *pp = 2; return vals(n);
1529 : }
1530 271019 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1531 3654186 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1532 : {
1533 3595105 : ulong p = tinyprimes[i];
1534 3595105 : if (n % p == 0)
1535 : {
1536 211532 : v = u_lvalrem(n, p, &n);
1537 211532 : if (n == 1) { *pp = p; return v; }
1538 209438 : return 0;
1539 : }
1540 : }
1541 : /* p | n => p >= CUTOFF */
1542 :
1543 59081 : if (n < CUTOFF3)
1544 : {
1545 46354 : if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1546 0 : if (uissquareall(n, &n)) { *pp = n; return 2; }
1547 0 : return 0;
1548 : }
1549 :
1550 : /* Check for squares, fourth powers, and eighth powers as appropriate. */
1551 12727 : v = 1;
1552 12727 : if (uissquareall(n, &n)) {
1553 0 : v <<= 1;
1554 0 : if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1555 0 : v <<= 1;
1556 0 : if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1557 : }
1558 : }
1559 :
1560 12727 : if (CUTOFF > ROOT5) mask = 1;
1561 : else
1562 : {
1563 12726 : const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1564 12726 : if (n < CUTOFF5) mask = 1; else mask = 3;
1565 12726 : if (CUTOFF <= ROOT7)
1566 : {
1567 12726 : const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1568 12726 : if (n >= CUTOFF7) mask = 7;
1569 : }
1570 : }
1571 :
1572 12727 : if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1573 12727 : if ((e = uis_357_power(n, &n, &mask))) v *= e;
1574 :
1575 12727 : if (uisprime_101(n)) { *pp = n; return v; }
1576 6984 : return 0;
1577 : }
1578 :
1579 : /*********************************************************************/
1580 : /** **/
1581 : /** KRONECKER SYMBOL **/
1582 : /** **/
1583 : /*********************************************************************/
1584 : /* t = 3,5 mod 8 ? (= 2 not a square mod t) */
1585 : static int
1586 650459712 : ome(long t)
1587 : {
1588 650459712 : switch(t & 7)
1589 : {
1590 : case 3:
1591 373775996 : case 5: return 1;
1592 276683716 : default: return 0;
1593 : }
1594 : }
1595 : /* t a t_INT, is t = 3,5 mod 8 ? */
1596 : static int
1597 5056202 : gome(GEN t)
1598 5056202 : { return signe(t)? ome( mod2BIL(t) ): 0; }
1599 :
1600 : /* assume y odd, return kronecker(x,y) * s */
1601 : static long
1602 509073688 : krouu_s(ulong x, ulong y, long s)
1603 : {
1604 509073688 : ulong x1 = x, y1 = y, z;
1605 2661791045 : while (x1)
1606 : {
1607 1643721552 : long r = vals(x1);
1608 1643827999 : if (r)
1609 : {
1610 892808629 : if (odd(r) && ome(y1)) s = -s;
1611 892624299 : x1 >>= r;
1612 : }
1613 1643643669 : if (x1 & y1 & 2) s = -s;
1614 1643643669 : z = y1 % x1; y1 = x1; x1 = z;
1615 : }
1616 508995805 : return (y1 == 1)? s: 0;
1617 : }
1618 :
1619 : long
1620 5052006 : kronecker(GEN x, GEN y)
1621 : {
1622 5052006 : pari_sp av = avma;
1623 5052006 : long s = 1, r;
1624 : ulong xu, yu;
1625 :
1626 5052006 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
1627 5052006 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
1628 5052006 : switch (signe(y))
1629 : {
1630 0 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
1631 0 : case 0: return is_pm1(x);
1632 : }
1633 5052006 : r = vali(y);
1634 5052006 : if (r)
1635 : {
1636 12111 : if (!mpodd(x)) { avma = av; return 0; }
1637 10396 : if (odd(r) && gome(x)) s = -s;
1638 10396 : y = shifti(y,-r);
1639 : }
1640 5050291 : x = modii(x,y);
1641 10946723 : while (lgefint(x) > 3) /* x < y */
1642 : {
1643 : GEN z;
1644 846141 : r = vali(x);
1645 846141 : if (r)
1646 : {
1647 461687 : if (odd(r) && gome(y)) s = -s;
1648 461687 : x = shifti(x,-r);
1649 : }
1650 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1651 846141 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
1652 846141 : z = remii(y,x); y = x; x = z;
1653 846141 : if (gc_needed(av,2))
1654 : {
1655 0 : if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
1656 0 : gerepileall(av, 2, &x, &y);
1657 : }
1658 : }
1659 5050291 : xu = itou(x);
1660 5050291 : if (!xu) return is_pm1(y)? s: 0;
1661 5029273 : r = vals(xu);
1662 5029273 : if (r)
1663 : {
1664 3496280 : if (odd(r) && gome(y)) s = -s;
1665 3496280 : xu >>= r;
1666 : }
1667 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1668 5029273 : if (xu & mod2BIL(y) & 2) s = -s;
1669 5029273 : yu = umodiu(y, xu);
1670 5029273 : avma = av; return krouu_s(yu, xu, s);
1671 : }
1672 :
1673 : long
1674 29246 : krois(GEN x, long y)
1675 : {
1676 : ulong yu;
1677 29246 : long s = 1;
1678 :
1679 29246 : if (y <= 0)
1680 : {
1681 14 : if (y == 0) return is_pm1(x);
1682 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
1683 : }
1684 : else
1685 29232 : yu = (ulong)y;
1686 29232 : if (!odd(yu))
1687 : {
1688 : long r;
1689 13461 : if (!mpodd(x)) return 0;
1690 10073 : r = vals(yu); yu >>= r;
1691 10073 : if (odd(r) && gome(x)) s = -s;
1692 : }
1693 25844 : return krouu_s(umodiu(x, yu), yu, s);
1694 : }
1695 : /* assume y != 0 */
1696 : long
1697 325884784 : kroiu(GEN x, ulong y)
1698 : {
1699 : long r;
1700 325884784 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
1701 2026493 : if (!mpodd(x)) return 0;
1702 2004443 : r = vals(y); y >>= r;
1703 2004443 : return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
1704 : }
1705 :
1706 : /* assume y > 0, odd, return s * kronecker(x,y) */
1707 : static long
1708 632471 : krouodd(ulong x, GEN y, long s)
1709 : {
1710 : long r;
1711 632471 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
1712 553882 : if (!x) return 0; /* y != 1 */
1713 553882 : r = vals(x);
1714 553882 : if (r)
1715 : {
1716 218652 : if (odd(r) && gome(y)) s = -s;
1717 218652 : x >>= r;
1718 : }
1719 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1720 553882 : if (x & mod2BIL(y) & 2) s = -s;
1721 553882 : return krouu_s(umodiu(y,x), x, s);
1722 : }
1723 :
1724 : long
1725 632347 : krosi(long x, GEN y)
1726 : {
1727 632347 : const pari_sp av = avma;
1728 632347 : long s = 1, r;
1729 632347 : switch (signe(y))
1730 : {
1731 0 : case -1: y = negi(y); if (x < 0) s = -1; break;
1732 0 : case 0: return (x==1 || x==-1);
1733 : }
1734 632347 : r = vali(y);
1735 632347 : if (r)
1736 : {
1737 8421 : if (!odd(x)) { avma = av; return 0; }
1738 8421 : if (odd(r) && ome(x)) s = -s;
1739 8421 : y = shifti(y,-r);
1740 : }
1741 632347 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
1742 632347 : s = krouodd((ulong)x, y, s);
1743 632347 : avma = av; return s;
1744 : }
1745 :
1746 : long
1747 124 : kroui(ulong x, GEN y)
1748 : {
1749 124 : const pari_sp av = avma;
1750 124 : long s = 1, r;
1751 124 : switch (signe(y))
1752 : {
1753 0 : case -1: y = negi(y); break;
1754 0 : case 0: return x==1UL;
1755 : }
1756 124 : r = vali(y);
1757 124 : if (r)
1758 : {
1759 0 : if (!odd(x)) { avma = av; return 0; }
1760 0 : if (odd(r) && ome(x)) s = -s;
1761 0 : y = shifti(y,-r);
1762 : }
1763 124 : s = krouodd(x, y, s);
1764 124 : avma = av; return s;
1765 : }
1766 :
1767 : long
1768 80030012 : kross(long x, long y)
1769 : {
1770 : ulong yu;
1771 80030012 : long s = 1;
1772 :
1773 80030012 : if (y <= 0)
1774 : {
1775 385 : if (y == 0) return (labs(x)==1);
1776 385 : yu = (ulong)-y; if (x < 0) s = -1;
1777 : }
1778 : else
1779 80029627 : yu = (ulong)y;
1780 80030012 : if (!odd(yu))
1781 : {
1782 : long r;
1783 24080029 : if (!odd(x)) return 0;
1784 15645204 : r = vals(yu); yu >>= r;
1785 15645204 : if (odd(r) && ome(x)) s = -s;
1786 : }
1787 71595187 : x %= (long)yu; if (x < 0) x += yu;
1788 71595187 : return krouu_s((ulong)x, yu, s);
1789 : }
1790 :
1791 : long
1792 105819567 : krouu(ulong x, ulong y)
1793 : {
1794 : long r;
1795 105819567 : if (odd(y)) return krouu_s(x, y, 1);
1796 2040 : if (!odd(x)) return 0;
1797 2040 : r = vals(y); y >>= r;
1798 2040 : return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
1799 : }
1800 :
1801 : /*********************************************************************/
1802 : /** **/
1803 : /** HILBERT SYMBOL **/
1804 : /** **/
1805 : /*********************************************************************/
1806 : /* x,y are t_INT or t_REAL */
1807 : static long
1808 9982 : mphilbertoo(GEN x, GEN y)
1809 : {
1810 9982 : long sx = signe(x), sy = signe(y);
1811 9982 : if (!sx || !sy) return 0;
1812 9982 : return (sx < 0 && sy < 0)? -1: 1;
1813 : }
1814 :
1815 : long
1816 53144 : hilbertii(GEN x, GEN y, GEN p)
1817 : {
1818 : pari_sp av;
1819 : long oddvx, oddvy, z;
1820 :
1821 53144 : if (!p) return mphilbertoo(x,y);
1822 43183 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
1823 43183 : if (!signe(x) || !signe(y)) return 0;
1824 43162 : av = avma;
1825 43162 : oddvx = odd(Z_pvalrem(x,p,&x));
1826 43162 : oddvy = odd(Z_pvalrem(y,p,&y));
1827 : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
1828 43162 : if (absequaliu(p, 2))
1829 : {
1830 10689 : z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
1831 10689 : if (oddvx && gome(y)) z = -z;
1832 10689 : if (oddvy && gome(x)) z = -z;
1833 : }
1834 : else
1835 : {
1836 32473 : z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
1837 32473 : if (oddvx && kronecker(y,p) < 0) z = -z;
1838 32473 : if (oddvy && kronecker(x,p) < 0) z = -z;
1839 : }
1840 43162 : avma = av; return z;
1841 : }
1842 :
1843 : static void
1844 196 : err_prec(void) { pari_err_PREC("hilbert"); }
1845 : static void
1846 161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
1847 : static void
1848 56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
1849 :
1850 : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
1851 : * Return lift(x) provided it's p-adic accuracy is large enough to decide
1852 : * hilbert()'s value [ problem at p = 2 ] */
1853 : static GEN
1854 420 : lift_intmod(GEN x, GEN *pp)
1855 : {
1856 420 : GEN p = *pp, N = gel(x,1);
1857 420 : x = gel(x,2);
1858 420 : if (!p)
1859 : {
1860 266 : *pp = p = N;
1861 266 : switch(itos_or_0(p))
1862 : {
1863 : case 2:
1864 126 : case 4: err_prec();
1865 : }
1866 140 : return x;
1867 : }
1868 154 : if (!signe(p)) err_oo(N);
1869 112 : if (absequaliu(p,2))
1870 42 : { if (vali(N) <= 2) err_prec(); }
1871 : else
1872 70 : { if (!dvdii(N,p)) err_p(N,p); }
1873 28 : if (!signe(x)) err_prec();
1874 21 : return x;
1875 : }
1876 : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
1877 : * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
1878 : * to decide hilbert()'s value [ problem at p = 2 ]*/
1879 : static GEN
1880 210 : lift_padic(GEN x, GEN *pp)
1881 : {
1882 210 : GEN p = *pp, q = gel(x,2), y = gel(x,4);
1883 210 : if (!p) *pp = p = q;
1884 147 : else if (!equalii(p,q)) err_p(p, q);
1885 105 : if (absequaliu(p,2) && precp(x) <= 2) err_prec();
1886 70 : if (!signe(y)) err_prec();
1887 70 : return odd(valp(x))? mulii(p,y): y;
1888 : }
1889 :
1890 : long
1891 658 : hilbert(GEN x, GEN y, GEN p)
1892 : {
1893 658 : pari_sp av = avma;
1894 658 : long tx = typ(x), ty = typ(y), z;
1895 :
1896 658 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
1897 658 : if (tx == t_REAL)
1898 : {
1899 77 : if (p && signe(p)) err_oo(p);
1900 63 : switch (ty)
1901 : {
1902 : case t_INT:
1903 7 : case t_REAL: return mphilbertoo(x,y);
1904 0 : case t_FRAC: return mphilbertoo(x,gel(y,1));
1905 56 : default: pari_err_TYPE2("hilbert",x,y);
1906 : }
1907 : }
1908 581 : if (ty == t_REAL)
1909 : {
1910 14 : if (p && signe(p)) err_oo(p);
1911 14 : switch (tx)
1912 : {
1913 : case t_INT:
1914 14 : case t_REAL: return mphilbertoo(x,y);
1915 0 : case t_FRAC: return mphilbertoo(gel(x,1),y);
1916 0 : default: pari_err_TYPE2("hilbert",x,y);
1917 : }
1918 : }
1919 567 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
1920 364 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
1921 :
1922 308 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
1923 245 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
1924 :
1925 168 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
1926 168 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
1927 :
1928 168 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
1929 168 : if (p && !signe(p)) p = NULL;
1930 168 : z = hilbertii(x,y,p); avma = av; return z;
1931 : }
1932 :
1933 : /*******************************************************************/
1934 : /* */
1935 : /* SQUARE ROOT MODULO p */
1936 : /* */
1937 : /*******************************************************************/
1938 :
1939 : static ulong
1940 26125149 : Fl_2gener_pre_all(long e, ulong p, ulong pi)
1941 : {
1942 : ulong y, m;
1943 : long k, i;
1944 26125149 : ulong q = (p-1) >> e; /* q = (p-1)/2^oo is odd */
1945 45706887 : for (k=2; ; k++)
1946 : { /* loop terminates for k < p (even if p composite) */
1947 45706887 : i = krouu(k, p);
1948 45706931 : if (i >= 0)
1949 : {
1950 19581739 : if (i) continue;
1951 7 : pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1952 : }
1953 26125192 : y = m = Fl_powu_pre(k, q, p, pi);
1954 70544895 : for (i=1; i<e; i++)
1955 44419558 : if ((m = Fl_sqr_pre(m, p, pi)) == 1) break;
1956 26125337 : if (i == e) break; /* success */
1957 19581738 : }
1958 26125331 : return y;
1959 : }
1960 :
1961 : ulong
1962 145045 : Fl_2gener_pre(ulong p, ulong pi)
1963 145045 : { return Fl_2gener_pre_all(vals(p-1), p, pi); }
1964 :
1965 : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
1966 : ulong
1967 62756392 : Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
1968 : {
1969 : long i, e, k;
1970 : ulong p1, q, v, w;
1971 :
1972 62756392 : if (!a) return 0;
1973 61349097 : p1 = p - 1; e = vals(p1);
1974 61354445 : if (e == 0) /* p = 2 */
1975 : {
1976 418831 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1977 418824 : return ((a & 1) == 0)? 0: 1;
1978 : }
1979 60935614 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
1980 60935614 : if (e == 1) y = p1;
1981 25983533 : else if (y==0) y = Fl_2gener_pre_all(e, p, pi);
1982 60935607 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
1983 60910861 : if (!p1) return 0;
1984 60910861 : v = Fl_mul_pre(a, p1, p, pi);
1985 60917216 : w = Fl_mul_pre(v, p1, p, pi);
1986 143366034 : while (w != 1)
1987 : { /* a*w = v^2, y primitive 2^e-th root of 1
1988 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1989 21642346 : p1 = Fl_sqr_pre(w,p,pi);
1990 21642346 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
1991 21642346 : if (k == e) return ~0UL;
1992 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1993 21534608 : p1 = y;
1994 21534608 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
1995 21534608 : y = Fl_sqr_pre(p1, p, pi); e = k;
1996 21534608 : w = Fl_mul_pre(y, w, p, pi);
1997 21534608 : v = Fl_mul_pre(v, p1, p, pi);
1998 : }
1999 60807673 : p1 = p - v; if (v > p1) v = p1;
2000 60807673 : return v;
2001 : }
2002 :
2003 : ulong
2004 58268667 : Fl_sqrt(ulong a, ulong p)
2005 : {
2006 58268667 : ulong pi = get_Fl_red(p);
2007 58266122 : return Fl_sqrt_pre_i(a, 0, p, pi);
2008 : }
2009 :
2010 : ulong
2011 4461997 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
2012 : {
2013 4461997 : return Fl_sqrt_pre_i(a, 0, p, pi);
2014 : }
2015 :
2016 : static ulong
2017 69774 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
2018 : {
2019 : ulong x, y, m;
2020 69774 : ulong le1 = upowuu(l, e-1);
2021 111569 : for (x = 2; ; x++)
2022 : {
2023 111569 : y = Fl_powu_pre(x, r, p, pi);
2024 111569 : if (y==1) continue;
2025 86371 : m = Fl_powu_pre(y, le1, p, pi);
2026 86371 : if (m != 1) break;
2027 41795 : }
2028 69774 : *pt_m = m;
2029 69774 : return y;
2030 : }
2031 :
2032 : /* solve x^l = a , l prime in G of order q.
2033 : *
2034 : * q = (l^e)*r, e >= 1, (r,l) = 1
2035 : * y generates the l-Sylow of G
2036 : * m = y^(l^(e-1)) != 1 */
2037 : static ulong
2038 132501 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
2039 : {
2040 : ulong p1, v, w, z, dl, zm;
2041 : ulong r, e, u2;
2042 132501 : if (a==0) return a;
2043 132495 : e = u_lvalrem(p-1, l, &r);
2044 132495 : u2 = Fl_inv(l%r, r);
2045 132495 : v = Fl_powu_pre(a, u2, p,pi);
2046 132495 : w = Fl_powu_pre(v, l, p,pi);
2047 132494 : w = Fl_mul_pre(w, Fl_inv(a, p),p,pi);
2048 132496 : if (w==1) return v;
2049 69774 : if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
2050 160319 : while (w!=1)
2051 : {
2052 74605 : ulong k = 0;
2053 74605 : p1 = w;
2054 : do
2055 : {
2056 111001 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
2057 111001 : k++;
2058 111001 : } while (p1!=1);
2059 74605 : if (k==e) return ~0UL;
2060 20771 : dl = 0; zm = 1;
2061 73187 : while (z!=zm)
2062 : {
2063 31645 : zm = Fl_mul_pre(zm, m, p, pi); dl++;
2064 : }
2065 20771 : dl = Fl_neg(dl, l);
2066 20771 : p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
2067 20771 : m = Fl_powu_pre(m, dl, p, pi);
2068 20771 : e = k;
2069 20771 : v = Fl_mul_pre(p1,v,p,pi);
2070 20771 : y = Fl_powu_pre(p1,l,p,pi);
2071 20771 : w = Fl_mul_pre(y,w,p,pi);
2072 : }
2073 15940 : return v;
2074 : }
2075 :
2076 : ulong
2077 132502 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
2078 : {
2079 132502 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2080 : }
2081 :
2082 : ulong
2083 0 : Fl_sqrtl(ulong a, ulong l, ulong p)
2084 : {
2085 0 : ulong pi = get_Fl_red(p);
2086 0 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2087 : }
2088 :
2089 : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
2090 : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
2091 : * and in average is about 2 or 2.5 times worse. But try both algorithms for
2092 : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
2093 : *
2094 : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
2095 : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
2096 : * If (a|p)=1, then sqrt(a) is in F_p.
2097 : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
2098 :
2099 : /* compute y^2, y = y[1] + y[2] X */
2100 : static GEN
2101 449 : sqrt_Cipolla_sqr(void *data, GEN y)
2102 : {
2103 449 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
2104 449 : GEN u2 = sqri(u), v2 = sqri(v);
2105 449 : v = subii(sqri(addii(v,u)), addii(u2,v2));
2106 449 : u = addii(u2, mulii(v2,n));
2107 : /* NOT mkvec2: must be gerepileupto-able */
2108 449 : retmkvec2(modii(u,p), modii(v,p));
2109 : }
2110 : /* compute (t+X) y^2 */
2111 : static GEN
2112 23 : sqrt_Cipolla_msqr(void *data, GEN y)
2113 : {
2114 23 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
2115 23 : ulong t = gt[2];
2116 23 : GEN d = addii(u, mului(t,v)), d2= sqri(d);
2117 23 : GEN b = remii(mulii(a,v), p);
2118 23 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
2119 23 : v = subii(d2, mulii(b,v));
2120 : /* NOT mkvec2: must be gerepileupto-able */
2121 23 : retmkvec2(modii(u,p), modii(v,p));
2122 : }
2123 : /* assume a reduced mod p [ otherwise correct but inefficient ] */
2124 : static GEN
2125 8 : sqrt_Cipolla(GEN a, GEN p)
2126 : {
2127 : pari_sp av1;
2128 : GEN u, v, n, y, pov2;
2129 : ulong t;
2130 :
2131 8 : if (kronecker(a, p) < 0) return NULL;
2132 8 : pov2 = shifti(p,-1);
2133 8 : if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
2134 :
2135 8 : av1 = avma;
2136 41 : for(t=1; ; t++)
2137 : {
2138 41 : n = subsi((long)(t*t), a);
2139 41 : if (kronecker(n, p) < 0) break;
2140 33 : avma = av1;
2141 33 : }
2142 :
2143 : /* compute (t+X)^((p-1)/2) =: u+vX */
2144 8 : u = utoipos(t);
2145 8 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
2146 : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
2147 : /* Now u+vX = (t+X)^((p-1)/2); thus
2148 : * (u+vX)(t+X) = sqrt(a) + 0 X
2149 : * Whence,
2150 : * sqrt(a) = (u+vt)t - v*a
2151 : * 0 = (u+vt)
2152 : * Thus a square root is v*a */
2153 :
2154 8 : v = Fp_mul(gel(y, 2), a, p);
2155 8 : if (cmpii(v,pov2) > 0) v = subii(p,v);
2156 8 : return v;
2157 : }
2158 :
2159 : /* Return NULL if p is found to be composite */
2160 : static GEN
2161 6627 : Fp_2gener_all(long e, GEN p)
2162 : {
2163 : GEN y, m;
2164 : long k;
2165 6627 : GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
2166 6627 : if (e==0 && !equaliu(p,2)) return NULL;
2167 11561 : for (k=2; ; k++)
2168 : {
2169 11561 : long i = krosi(k, p);
2170 11561 : if (i >= 0)
2171 : {
2172 4934 : if (i) continue;
2173 0 : return NULL;
2174 : }
2175 6627 : y = m = Fp_pow(utoi(k), q, p);
2176 16632 : for (i=1; i<e; i++)
2177 10005 : if (equali1(m = Fp_sqr(m, p))) break;
2178 6627 : if (i == e) break; /* success */
2179 4934 : }
2180 6627 : return y;
2181 : }
2182 :
2183 : /* Return NULL if p is found to be composite */
2184 : GEN
2185 686 : Fp_2gener(GEN p)
2186 686 : { return Fp_2gener_all(vali(subis(p,1)),p); }
2187 :
2188 : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
2189 : GEN
2190 2349422 : Fp_sqrt_i(GEN a, GEN y, GEN p)
2191 : {
2192 2349422 : pari_sp av = avma, av1;
2193 : long i, k, e;
2194 : GEN p1, q, v, w;
2195 :
2196 2349422 : if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
2197 2349422 : if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
2198 2349422 : if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
2199 2349422 : if (lgefint(p) == 3)
2200 : {
2201 2321698 : ulong pp = uel(p,2), u = Fl_sqrt(umodiu(a, pp), pp);
2202 2321684 : if (u == ~0UL) return NULL;
2203 2321635 : return utoi(u);
2204 : }
2205 :
2206 27724 : p1 = subiu(p,1); e = vali(p1);
2207 27724 : a = modii(a, p);
2208 :
2209 : /* On average, the algorithm of Cipolla is better than the algorithm of
2210 : * Tonelli and Shanks if and only if e(e-1)>8*log2(n)+20
2211 : * see LNCS 2286 pp 430 [GTL] */
2212 27724 : if (e*(e-1) > 20 + 8 * expi(p))
2213 : {
2214 8 : v = sqrt_Cipolla(a,p);
2215 8 : if (!v) { avma = av; return NULL; }
2216 8 : return gerepileuptoint(av,v);
2217 : }
2218 :
2219 27716 : if (e == 0) /* p = 2 */
2220 : {
2221 0 : avma = av;
2222 0 : if (!absequaliu(p,2)) pari_err_PRIME("Fp_sqrt [modulus]",p);
2223 0 : if (!signe(a) || !mod2(a)) return gen_0;
2224 0 : return gen_1;
2225 : }
2226 27716 : q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
2227 27716 : if (e == 1) y = p1;
2228 12409 : else if (!y) y = Fp_2gener_all(e, p);
2229 27716 : if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
2230 27716 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ [(q-1)/2] */
2231 27716 : if (!signe(p1)) { avma=av; return gen_0; }
2232 27709 : v = Fp_mul(a, p1, p);
2233 27709 : w = Fp_mul(v, p1, p);
2234 66682 : while (!equali1(w))
2235 : { /* a*w = v^2, y primitive 2^e-th root of 1
2236 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
2237 11274 : p1 = Fp_sqr(w,p);
2238 11274 : for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
2239 11274 : if (k == e) { avma=av; return NULL; } /* p composite or (a/p) != 1 */
2240 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
2241 11264 : p1 = y;
2242 11264 : for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
2243 11264 : y = Fp_sqr(p1, p); e = k;
2244 11264 : w = Fp_mul(y, w, p);
2245 11264 : v = Fp_mul(v, p1, p);
2246 11264 : if (gc_needed(av,1))
2247 : {
2248 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
2249 0 : gerepileall(av,3, &y,&w,&v);
2250 : }
2251 : }
2252 27699 : av1 = avma;
2253 27699 : p1 = subii(p,v); if (cmpii(v,p1) > 0) v = p1; else avma = av1;
2254 27699 : return gerepileuptoint(av, v);
2255 : }
2256 :
2257 : GEN
2258 2333413 : Fp_sqrt(GEN a, GEN p)
2259 : {
2260 2333413 : return Fp_sqrt_i(a, NULL, p);
2261 : }
2262 :
2263 : /*********************************************************************/
2264 : /** **/
2265 : /** GCD & BEZOUT **/
2266 : /** **/
2267 : /*********************************************************************/
2268 :
2269 : GEN
2270 14185488 : lcmii(GEN x, GEN y)
2271 : {
2272 : pari_sp av;
2273 : GEN a, b;
2274 14185488 : if (!signe(x) || !signe(y)) return gen_0;
2275 14185488 : av = avma;
2276 14185488 : a = gcdii(x,y); if (!equali1(a)) y = diviiexact(y,a);
2277 14185488 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
2278 : }
2279 :
2280 : /*********************************************************************/
2281 : /** **/
2282 : /** CHINESE REMAINDERS **/
2283 : /** **/
2284 : /*********************************************************************/
2285 :
2286 : /* Chinese Remainder Theorem. x and y must have the same type (integermod,
2287 : * polymod, or polynomial/vector/matrix recursively constructed with these
2288 : * as coefficients). Creates (with the same type) a z in the same residue
2289 : * class as x and the same residue class as y, if it is possible.
2290 : *
2291 : * We also allow (during recursion) two identical objects even if they are
2292 : * not integermod or polymod. For example:
2293 : *
2294 : * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
2295 : * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
2296 : * ? chinese(x, y)
2297 : * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
2298 :
2299 : static GEN
2300 708665 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
2301 : {
2302 708665 : GEN z = gassoc_proto(f,x,NULL);
2303 708658 : if (z == gen_1) retmkintmod(gen_0,gen_1);
2304 708630 : return z;
2305 : }
2306 :
2307 : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
2308 : * call chinese: makes Mod(0,1) a better "neutral" element */
2309 : static GEN
2310 21 : chinese_intpol(GEN x,GEN y)
2311 : {
2312 21 : pari_sp av = avma;
2313 21 : GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
2314 21 : return gerepileupto(av, chinese(z, y));
2315 : }
2316 :
2317 : GEN
2318 203 : chinese1(GEN x) { return gen_chinese(x,chinese); }
2319 :
2320 : GEN
2321 16520 : chinese(GEN x, GEN y)
2322 : {
2323 : pari_sp av;
2324 16520 : long tx = typ(x), ty;
2325 : GEN z,p1,p2,d,u,v;
2326 :
2327 16520 : if (!y) return chinese1(x);
2328 16471 : if (gidentical(x,y)) return gcopy(x);
2329 16464 : ty = typ(y);
2330 16464 : if (tx == ty) switch(tx)
2331 : {
2332 : case t_POLMOD:
2333 : {
2334 28 : GEN A = gel(x,1), B = gel(y,1);
2335 28 : GEN a = gel(x,2), b = gel(y,2);
2336 28 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
2337 28 : if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
2338 28 : av = avma;
2339 28 : d = RgX_extgcd(A,B,&u,&v);
2340 28 : p2 = gsub(b, a);
2341 28 : if (!gequal0(gmod(p2, d))) break;
2342 28 : p1 = gdiv(A,d);
2343 28 : p2 = gadd(a, gmul(gmul(u,p1), p2));
2344 :
2345 28 : z = cgetg(3, t_POLMOD);
2346 28 : gel(z,1) = gmul(p1,B);
2347 28 : gel(z,2) = gmod(p2,gel(z,1));
2348 28 : return gerepileupto(av, z);
2349 : }
2350 : case t_INTMOD:
2351 : {
2352 16401 : GEN A = gel(x,1), B = gel(y,1);
2353 16401 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
2354 16401 : z = cgetg(3,t_INTMOD);
2355 16401 : Z_chinese_pre(A, B, &C, &U, &d);
2356 16401 : c = Z_chinese_post(a, b, C, U, d);
2357 16401 : if (!c) pari_err_OP("chinese", x,y);
2358 16401 : gel(z,1) = icopy_avma(C, (pari_sp)z);
2359 16401 : gel(z,2) = icopy_avma(c, (pari_sp)gel(z,1));
2360 16401 : avma = (pari_sp)gel(z,2); return z;
2361 : }
2362 : case t_POL:
2363 : {
2364 7 : long i, lx = lg(x), ly = lg(y);
2365 7 : if (varn(x) != varn(y)) break;
2366 7 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
2367 7 : z = cgetg(lx, t_POL); z[1] = x[1];
2368 7 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2369 7 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
2370 7 : return z;
2371 : }
2372 :
2373 : case t_VEC: case t_COL: case t_MAT:
2374 : {
2375 : long i, lx;
2376 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
2377 7 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2378 7 : return z;
2379 : }
2380 : }
2381 21 : if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
2382 7 : if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
2383 0 : pari_err_OP("chinese",x,y);
2384 : return NULL; /* LCOV_EXCL_LINE */
2385 : }
2386 :
2387 : /* init chinese(Mod(.,A), Mod(.,B)) */
2388 : void
2389 218589 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
2390 : {
2391 218589 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
2392 218589 : GEN t = diviiexact(A,d);
2393 218589 : *pU = mulii(u, t);
2394 218589 : *pC = mulii(t, B);
2395 218589 : if (pd) *pd = d;
2396 218589 : }
2397 : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
2398 : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
2399 : * If d not NULL, check wether a = b mod d. */
2400 : GEN
2401 798319 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
2402 : {
2403 : GEN b_a;
2404 798319 : if (!signe(a))
2405 : {
2406 430564 : if (d && remii(b, d) != gen_0) return NULL;
2407 430564 : return Fp_mul(b, U, C);
2408 : }
2409 367755 : b_a = subii(b,a);
2410 367755 : if (d && remii(b_a, d) != gen_0) return NULL;
2411 367755 : return modii(addii(a, mulii(U, b_a)), C);
2412 : }
2413 : static ulong
2414 79289 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
2415 : {
2416 79289 : if (!a) return Fl_mul(b, U, C);
2417 79289 : return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
2418 : }
2419 :
2420 : GEN
2421 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
2422 : {
2423 2142 : pari_sp av = avma;
2424 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
2425 2142 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
2426 : }
2427 : GEN
2428 199990 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
2429 : {
2430 199990 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
2431 199990 : return Z_chinese_post(a,b, *pC, U, NULL);
2432 : }
2433 :
2434 : /* return lift(chinese(a mod A, b mod B))
2435 : * assume(A,B)=1, a,b,A,B integers and C = A*B */
2436 : GEN
2437 1050 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
2438 : {
2439 1050 : pari_sp av = avma;
2440 1050 : GEN U = mulii(Fp_inv(A,B), A);
2441 1050 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2442 : }
2443 : ulong
2444 79289 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
2445 79289 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
2446 :
2447 : /* chinese1 for coprime moduli in Z */
2448 : static GEN
2449 578414 : chinese1_coprime_Z_aux(GEN x, GEN y)
2450 : {
2451 578414 : GEN z = cgetg(3, t_INTMOD);
2452 578414 : GEN A = gel(x,1), a = gel(x, 2);
2453 578414 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
2454 578414 : pari_sp av = avma;
2455 578414 : GEN U = mulii(Fp_inv(A,B), A);
2456 578414 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2457 578414 : gel(z,1) = C; return z;
2458 : }
2459 : GEN
2460 708462 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
2461 :
2462 : /*********************************************************************/
2463 : /** **/
2464 : /** MODULAR EXPONENTIATION **/
2465 : /** **/
2466 : /*********************************************************************/
2467 :
2468 : /* xa, ya = t_VECSMALL */
2469 : GEN
2470 75395 : ZV_producttree(GEN xa)
2471 : {
2472 75395 : long n = lg(xa)-1;
2473 75395 : long m = n==1 ? 1: expu(n-1)+1;
2474 75396 : GEN T = cgetg(m+1, t_VEC), t;
2475 : long i, j, k;
2476 75395 : t = cgetg(((n+1)>>1)+1, t_VEC);
2477 75394 : if (typ(xa)==t_VECSMALL)
2478 : {
2479 7250090 : for (j=1, k=1; k<n; j++, k+=2)
2480 7200146 : gel(t, j) = muluu(xa[k], xa[k+1]);
2481 49944 : if (k==n) gel(t, j) = utoi(xa[k]);
2482 : } else {
2483 123388 : for (j=1, k=1; k<n; j++, k+=2)
2484 97937 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
2485 25451 : if (k==n) gel(t, j) = icopy(gel(xa,k));
2486 : }
2487 75395 : gel(T,1) = t;
2488 176807 : for (i=2; i<=m; i++)
2489 : {
2490 101412 : GEN u = gel(T, i-1);
2491 101412 : long n = lg(u)-1;
2492 101412 : t = cgetg(((n+1)>>1)+1, t_VEC);
2493 7354433 : for (j=1, k=1; k<n; j++, k+=2)
2494 7253021 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
2495 101412 : if (k==n) gel(t, j) = gel(u, k);
2496 101412 : gel(T, i) = t;
2497 : }
2498 75395 : return T;
2499 : }
2500 :
2501 : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
2502 : GEN
2503 842576 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
2504 : {
2505 : long i,j,k;
2506 842576 : long m = lg(T)-1, n = lg(P)-1;
2507 : GEN t;
2508 842576 : GEN Tp = cgetg(m+1, t_VEC);
2509 842574 : gel(Tp, m) = mkvec(A);
2510 1560153 : for (i=m-1; i>=1; i--)
2511 : {
2512 717581 : GEN u = gel(T, i);
2513 717581 : GEN v = gel(Tp, i+1);
2514 717581 : long n = lg(u)-1;
2515 717581 : t = cgetg(n+1, t_VEC);
2516 1781403 : for (j=1, k=1; k<n; j++, k+=2)
2517 : {
2518 1063805 : gel(t, k) = modii(gel(v, j), gel(u, k));
2519 1063804 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
2520 : }
2521 717598 : if (k==n) gel(t, k) = gel(v, j);
2522 717598 : gel(Tp, i) = t;
2523 : }
2524 : {
2525 842572 : GEN u = gel(T, i+1);
2526 842572 : GEN v = gel(Tp, i+1);
2527 842572 : long l = lg(u)-1;
2528 842572 : if (typ(P)==t_VECSMALL)
2529 : {
2530 767214 : GEN R = cgetg(n+1, t_VECSMALL);
2531 2401780 : for (j=1, k=1; j<=l; j++, k+=2)
2532 : {
2533 1634553 : uel(R,k) = umodiu(gel(v, j), P[k]);
2534 1634567 : if (k < n)
2535 1343292 : uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
2536 : }
2537 767227 : return R;
2538 : }
2539 : else
2540 : {
2541 75358 : GEN R = cgetg(n+1, t_VEC);
2542 347175 : for (j=1, k=1; j<=l; j++, k+=2)
2543 : {
2544 271812 : gel(R,k) = modii(gel(v, j), gel(P,k));
2545 271815 : if (k < n)
2546 241516 : gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
2547 : }
2548 75363 : return R;
2549 : }
2550 : }
2551 : }
2552 :
2553 : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
2554 : GEN
2555 2354116 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
2556 : {
2557 2354116 : long m = lg(T)-1, n = lg(A)-1;
2558 : long i,j,k;
2559 2354116 : GEN Tp = cgetg(m+1, t_VEC);
2560 2345933 : GEN M = gel(T, 1);
2561 2345933 : GEN t = cgetg(lg(M), t_VEC);
2562 2390965 : if (typ(P)==t_VECSMALL)
2563 : {
2564 18823222 : for (j=1, k=1; k<n; j++, k+=2)
2565 : {
2566 16509184 : pari_sp av = avma;
2567 16509184 : GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
2568 16354760 : GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
2569 16386810 : gel(t, j) = gerepileuptoint(av, tj);
2570 : }
2571 2314038 : if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
2572 : } else
2573 : {
2574 111110 : for (j=1, k=1; k<n; j++, k+=2)
2575 : {
2576 87010 : pari_sp av = avma;
2577 87010 : GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
2578 87010 : GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
2579 87010 : gel(t, j) = gerepileuptoint(av, tj);
2580 : }
2581 24100 : if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
2582 : }
2583 2338389 : gel(Tp, 1) = t;
2584 8803337 : for (i=2; i<=m; i++)
2585 : {
2586 6469547 : GEN u = gel(T, i-1), M = gel(T, i);
2587 6469547 : GEN t = cgetg(lg(M), t_VEC);
2588 6562579 : GEN v = gel(Tp, i-1);
2589 6562579 : long n = lg(v)-1;
2590 21465609 : for (j=1, k=1; k<n; j++, k+=2)
2591 : {
2592 15000661 : pari_sp av = avma;
2593 60002644 : gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
2594 45001983 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
2595 : }
2596 6464948 : if (k==n) gel(t, j) = gel(v, k);
2597 6464948 : gel(Tp, i) = t;
2598 : }
2599 2333790 : return gmael(Tp,m,1);
2600 : }
2601 :
2602 : static GEN
2603 46922 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2604 : {
2605 46922 : long i, l = lg(gel(vA,1)), n = lg(P);
2606 46922 : GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
2607 2275984 : for (i=1; i < l; i++)
2608 : {
2609 2228674 : pari_sp av = avma;
2610 2228674 : GEN c, A = cgetg(n, t_VECSMALL);
2611 : long j;
2612 2247475 : for (j=1; j < n; j++) A[j] = mael(vA,j,i);
2613 2247475 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
2614 2229482 : gel(V,i) = gerepileuptoint(av, c);
2615 : }
2616 47310 : return V;
2617 : }
2618 :
2619 : static GEN
2620 4200 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
2621 : {
2622 4200 : long i, l = lg(gel(vA,1)), n = lg(P);
2623 4200 : GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_POL);
2624 4200 : V[1] = evalsigne(1) | evalvarn(0);
2625 58163 : for (i=2; i < l; i++)
2626 : {
2627 53963 : pari_sp av = avma;
2628 53963 : GEN c, A = cgetg(n, t_VECSMALL);
2629 : long j;
2630 53963 : for (j=1; j < n; j++) A[j] = mael(vA,j,i);
2631 53963 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
2632 53963 : gel(V,i) = gerepileuptoint(av, c);
2633 : }
2634 4200 : return V;
2635 : }
2636 :
2637 : GEN
2638 43393 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
2639 : {
2640 43393 : return ncV_polint_center_tree(vA, P, T, R, m2);
2641 : }
2642 :
2643 : static GEN
2644 2433 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
2645 : {
2646 2433 : long i, j, l = lg(gel(mA,1)), n = lg(P);
2647 2433 : long pending = 0, workid, cnt = 0;
2648 : struct pari_mt pt;
2649 : GEN worker, done, va, M;
2650 2433 : GEN A = cgetg(n, t_VEC);
2651 2433 : worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
2652 2433 : va = mkvec(gen_0);
2653 2433 : M = cgetg(l, t_MAT);
2654 2433 : if (DEBUGLEVEL>2) err_printf("Start parallel Chinese remainder: ");
2655 2433 : mt_queue_start_lim(&pt, worker, l-1);
2656 50078 : for (i=1; i<l || pending; i++)
2657 : {
2658 47645 : for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
2659 47645 : gel(va, 1) = A;
2660 47645 : mt_queue_submit(&pt, i, i<l? va: NULL);
2661 47645 : done = mt_queue_get(&pt, &workid, &pending);
2662 47645 : if (done)
2663 : {
2664 43835 : gel(M,workid) = done;
2665 43835 : if (DEBUGLEVEL>2) err_printf("%ld%% ",(++cnt)*100/(l-1));
2666 : }
2667 : }
2668 2433 : if (DEBUGLEVEL>2) err_printf("\n");
2669 2433 : mt_queue_end(&pt);
2670 2433 : return M;
2671 : }
2672 :
2673 : /* return [A mod P[i], i=1..#P] */
2674 : GEN
2675 0 : Z_ZV_mod(GEN A, GEN P)
2676 : {
2677 0 : pari_sp av = avma;
2678 0 : return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
2679 : }
2680 : /* P a t_VECSMALL */
2681 : GEN
2682 0 : Z_nv_mod(GEN A, GEN P)
2683 : {
2684 0 : pari_sp av = avma;
2685 0 : return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
2686 : }
2687 : /* B a ZX, T = ZV_producttree(P) */
2688 : GEN
2689 78038 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
2690 : {
2691 78038 : long i, j, l = lg(B), n = lg(A)-1;
2692 78038 : GEN V = cgetg(n+1, t_VEC);
2693 335557 : for (j=1; j <= n; j++)
2694 : {
2695 257517 : gel(V, j) = cgetg(l, t_VECSMALL);
2696 257515 : mael(V, j, 1) = B[1]&VARNBITS;
2697 : }
2698 845273 : for (i=2; i < l; i++)
2699 : {
2700 767235 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
2701 3745079 : for (j=1; j <= n; j++)
2702 2977846 : mael(V, j, i) = v[j];
2703 : }
2704 335557 : for (j=1; j <= n; j++)
2705 257520 : (void) Flx_renormalize(gel(V, j), l);
2706 78037 : return V;
2707 : }
2708 :
2709 : static GEN
2710 29582 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
2711 :
2712 : GEN
2713 4200 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
2714 : {
2715 4200 : pari_sp av = avma;
2716 4200 : long i, j, l = lg(P), n = lg(xa)-1, vP = varn(P);
2717 4200 : GEN V = cgetg(n+1, t_VEC);
2718 13834 : for (j=1; j <= n; j++)
2719 : {
2720 9634 : gel(V, j) = cgetg(l, t_POL);
2721 9634 : mael(V, j, 1) = vP;
2722 : }
2723 33782 : for (i=2; i < l; i++)
2724 : {
2725 29582 : GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
2726 102139 : for (j=1; j <= n; j++)
2727 72557 : gmael(V, j, i) = gel(v,j);
2728 : }
2729 13834 : for (j=1; j <= n; j++)
2730 9634 : (void) FlxX_renormalize(gel(V, j), l);
2731 4200 : return gerepilecopy(av, V);
2732 : }
2733 :
2734 : static GEN
2735 74011 : ZV_sqr(GEN z)
2736 : {
2737 74011 : long i,l = lg(z);
2738 74011 : GEN x = cgetg(l, t_VEC);
2739 74011 : if (typ(z)==t_VECSMALL)
2740 49911 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
2741 : else
2742 24100 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
2743 74009 : return x;
2744 : }
2745 :
2746 : static GEN
2747 723732 : ZT_sqr(GEN z)
2748 : {
2749 723732 : if (typ(z) == t_INT)
2750 479166 : return sqri(z);
2751 : else
2752 : {
2753 244566 : long i,l = lg(z);
2754 244566 : GEN x = cgetg(l, t_VEC);
2755 244568 : for (i=1; i<l; i++) gel(x,i) = ZT_sqr(gel(z,i));
2756 244570 : return x;
2757 : }
2758 : }
2759 :
2760 : static GEN
2761 74010 : ZV_invdivexact(GEN y, GEN x)
2762 : {
2763 74010 : long i, l = lg(y);
2764 74010 : GEN z = cgetg(l,t_VEC);
2765 74012 : if (typ(x)==t_VECSMALL)
2766 357791 : for (i=1; i<l; i++)
2767 : {
2768 307881 : pari_sp av = avma;
2769 307881 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
2770 307881 : avma = av;
2771 307881 : gel(z,i) = utoipos(a);
2772 : }
2773 : else
2774 207693 : for (i=1; i<l; i++)
2775 183593 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
2776 74010 : return z;
2777 : }
2778 :
2779 : /* P t_VECSMALL or t_VEC of t_INT */
2780 : GEN
2781 74009 : ZV_chinesetree(GEN P, GEN T)
2782 : {
2783 74009 : GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
2784 74009 : GEN mod = gmael(T,lg(T)-1,1);
2785 74009 : return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
2786 : }
2787 :
2788 : static GEN
2789 30083 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
2790 : {
2791 30083 : if (!pt_mod)
2792 2433 : return gerepileupto(av, a);
2793 : else
2794 : {
2795 27650 : GEN mod = gmael(T, lg(T)-1, 1);
2796 27650 : gerepileall(av, 2, &a, &mod);
2797 27650 : *pt_mod = mod;
2798 27650 : return a;
2799 : }
2800 : }
2801 :
2802 : GEN
2803 24100 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
2804 : {
2805 24100 : pari_sp av = avma;
2806 24100 : GEN T = ZV_producttree(P);
2807 24100 : GEN R = ZV_chinesetree(P, T);
2808 24100 : GEN a = ZV_chinese_tree(A, P, T, R);
2809 24100 : return gc_chinese(av, T, a, pt_mod);
2810 : }
2811 :
2812 : GEN
2813 4200 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2814 : {
2815 4200 : pari_sp av = avma;
2816 4200 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2817 4200 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
2818 4200 : return gerepileupto(av, a);
2819 : }
2820 :
2821 : GEN
2822 0 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2823 : {
2824 0 : pari_sp av = avma;
2825 0 : GEN T = ZV_producttree(P);
2826 0 : GEN R = ZV_chinesetree(P, T);
2827 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2828 0 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
2829 0 : return gc_chinese(av, T, a, pt_mod);
2830 : }
2831 :
2832 : GEN
2833 3550 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2834 : {
2835 3550 : pari_sp av = avma;
2836 3550 : GEN T = ZV_producttree(P);
2837 3550 : GEN R = ZV_chinesetree(P, T);
2838 3550 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2839 3550 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
2840 3550 : return gc_chinese(av, T, a, pt_mod);
2841 : }
2842 :
2843 : GEN
2844 2433 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2845 : {
2846 2433 : pari_sp av = avma;
2847 2433 : GEN T = ZV_producttree(P);
2848 2433 : GEN R = ZV_chinesetree(P, T);
2849 2433 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2850 2433 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
2851 2433 : return gc_chinese(av, T, a, pt_mod);
2852 : }
2853 :
2854 : /**********************************************************************
2855 : ** **
2856 : ** Powering over (Z/NZ)^*, small N **
2857 : ** **
2858 : **********************************************************************/
2859 :
2860 : /* 2^n mod p; assume n > 1 */
2861 : static ulong
2862 15695522 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
2863 : {
2864 15695522 : ulong y = 2;
2865 15695522 : int j = 1+bfffo(n);
2866 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
2867 15695522 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
2868 184104207 : for (; j; n<<=1,j--)
2869 : {
2870 168408723 : y = Fl_sqr_pre(y,p,pi);
2871 168408966 : if (n & HIGHBIT) y = Fl_double(y, p);
2872 : }
2873 15695484 : return y;
2874 : }
2875 :
2876 : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
2877 : static ulong
2878 4853151 : Fl_2powu(ulong n, ulong p)
2879 : {
2880 4853151 : ulong y = 2;
2881 4853151 : int j = 1+bfffo(n);
2882 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
2883 4853151 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
2884 59799448 : for (; j; n<<=1,j--)
2885 : {
2886 54961841 : y = (y*y) % p;
2887 54961841 : if (n & HIGHBIT) y = Fl_double(y, p);
2888 : }
2889 4837607 : return y;
2890 : }
2891 :
2892 : ulong
2893 90459524 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
2894 : {
2895 : ulong y, z, n;
2896 90459524 : if (n0 <= 1)
2897 : { /* frequent special cases */
2898 12604741 : if (n0 == 1) return x;
2899 4749786 : if (n0 == 0) return 1;
2900 : }
2901 77854783 : if (x <= 2)
2902 : {
2903 16861235 : if (x == 2) return Fl_2powu_pre(n0, p, pi);
2904 1165824 : return x; /* 0 or 1 */
2905 : }
2906 60993548 : y = 1; z = x; n = n0;
2907 : for(;;)
2908 : {
2909 552332877 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
2910 552921269 : n>>=1; if (!n) return y;
2911 491941090 : z = Fl_sqr_pre(z,p,pi);
2912 491339329 : }
2913 : }
2914 :
2915 : ulong
2916 24005530 : Fl_powu(ulong x, ulong n0, ulong p)
2917 : {
2918 : ulong y, z, n;
2919 24005530 : if (n0 <= 2)
2920 : { /* frequent special cases */
2921 12948853 : if (n0 == 2) return Fl_sqr(x,p);
2922 4463593 : if (n0 == 1) return x;
2923 2617 : if (n0 == 0) return 1;
2924 : }
2925 11056677 : if (x <= 1) return x; /* 0 or 1 */
2926 11021370 : if (p & HIGHMASK)
2927 2109460 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
2928 8911910 : if (x == 2) return Fl_2powu(n0, p);
2929 4079894 : y = 1; z = x; n = n0;
2930 : for(;;)
2931 : {
2932 29856947 : if (n&1) y = (y*z) % p;
2933 29856947 : n>>=1; if (!n) return y;
2934 25777053 : z = (z*z) % p;
2935 25777053 : }
2936 : }
2937 :
2938 : /* Reduce data dependency to maximize internal parallelism */
2939 : GEN
2940 8682403 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
2941 : {
2942 : long i, k;
2943 8682403 : GEN powers = cgetg(n + 2, t_VECSMALL);
2944 8682763 : powers[1] = 1; if (n == 0) return powers;
2945 8682763 : powers[2] = x;
2946 45571542 : for (i = 3, k=2; i <= n; i+=2, k++)
2947 : {
2948 36893548 : powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
2949 36893066 : powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
2950 : }
2951 8677994 : if (i==n+1)
2952 7633773 : powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
2953 8678271 : return powers;
2954 : }
2955 :
2956 : GEN
2957 1520 : Fl_powers(ulong x, long n, ulong p)
2958 : {
2959 1520 : return Fl_powers_pre(x, n, p, get_Fl_red(p));
2960 : }
2961 :
2962 : /**********************************************************************
2963 : ** **
2964 : ** Powering over (Z/NZ)^*, large N **
2965 : ** **
2966 : **********************************************************************/
2967 :
2968 : static GEN
2969 4525948 : Fp_dblsqr(GEN x, GEN N)
2970 : {
2971 4525948 : GEN z = shifti(Fp_sqr(x, N), 1);
2972 4525948 : return cmpii(z, N) >= 0? subii(z, N): z;
2973 : }
2974 :
2975 : typedef struct muldata {
2976 : GEN (*sqr)(void * E, GEN x);
2977 : GEN (*mul)(void * E, GEN x, GEN y);
2978 : GEN (*mul2)(void * E, GEN x);
2979 : } muldata;
2980 :
2981 : /* modified Barrett reduction with one fold */
2982 : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
2983 :
2984 : static GEN
2985 9957 : Fp_invmBarrett(GEN p, long s)
2986 : {
2987 9957 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
2988 9957 : return mkvec2(Q,R);
2989 : }
2990 :
2991 : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
2992 : * a = r (mod N) */
2993 : static GEN
2994 4744824 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
2995 : {
2996 4744824 : pari_sp av = avma;
2997 4744824 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
2998 4744824 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
2999 4744824 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
3000 4744824 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
3001 4744824 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
3002 4744824 : GEN r = subii(A, mulii(q, N));
3003 4744824 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
3004 4744824 : if (signe(sr)<0) return gerepileuptoint(av, r);
3005 2842388 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
3006 2842388 : if (signe(sr)<0) return gerepileuptoint(av, r);
3007 123474 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
3008 123474 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
3009 : }
3010 :
3011 : /* Montgomery reduction */
3012 :
3013 : INLINE ulong
3014 325438 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
3015 :
3016 : struct montred
3017 : {
3018 : GEN N;
3019 : ulong inv;
3020 : };
3021 :
3022 : /* Montgomery reduction */
3023 : static GEN
3024 32761189 : _sqr_montred(void * E, GEN x)
3025 : {
3026 32761189 : struct montred * D = (struct montred *) E;
3027 32761189 : return red_montgomery(sqri(x), D->N, D->inv);
3028 : }
3029 :
3030 : /* Montgomery reduction */
3031 : static GEN
3032 3133385 : _mul_montred(void * E, GEN x, GEN y)
3033 : {
3034 3133385 : struct montred * D = (struct montred *) E;
3035 3133385 : return red_montgomery(mulii(x, y), D->N, D->inv);
3036 : }
3037 :
3038 : static GEN
3039 5586926 : _mul2_montred(void * E, GEN x)
3040 : {
3041 5586926 : struct montred * D = (struct montred *) E;
3042 5586926 : GEN z = shifti(_sqr_montred(E, x), 1);
3043 5586923 : long l = lgefint(D->N);
3044 5586923 : while (lgefint(z) > l) z = subii(z, D->N);
3045 5586924 : return z;
3046 : }
3047 :
3048 : static GEN
3049 12964288 : _sqr_remii(void* N, GEN x)
3050 12964288 : { return remii(sqri(x), (GEN) N); }
3051 :
3052 : static GEN
3053 1347469 : _mul_remii(void* N, GEN x, GEN y)
3054 1347469 : { return remii(mulii(x, y), (GEN) N); }
3055 :
3056 : static GEN
3057 3299485 : _mul2_remii(void* N, GEN x)
3058 3299485 : { return Fp_dblsqr(x, (GEN) N); }
3059 :
3060 : struct redbarrett
3061 : {
3062 : GEN iM, N;
3063 : long s;
3064 : };
3065 :
3066 : static GEN
3067 4258621 : _sqr_remiibar(void *E, GEN x)
3068 : {
3069 4258621 : struct redbarrett * D = (struct redbarrett *) E;
3070 4258621 : return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
3071 : }
3072 :
3073 : static GEN
3074 486203 : _mul_remiibar(void *E, GEN x, GEN y)
3075 : {
3076 486203 : struct redbarrett * D = (struct redbarrett *) E;
3077 486203 : return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
3078 : }
3079 :
3080 : static GEN
3081 1226463 : _mul2_remiibar(void *E, GEN x)
3082 : {
3083 1226463 : struct redbarrett * D = (struct redbarrett *) E;
3084 1226463 : return Fp_dblsqr(x, D->N);
3085 : }
3086 :
3087 : static long
3088 423257 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
3089 : {
3090 423257 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
3091 : {
3092 9957 : struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
3093 9957 : D->sqr = &_sqr_remiibar;
3094 9957 : D->mul = &_mul_remiibar;
3095 9957 : D->mul2 = &_mul2_remiibar;
3096 9957 : E->N = N;
3097 9957 : E->s = 1+(expi(N)>>1);
3098 9957 : E->iM = Fp_invmBarrett(N, E->s);
3099 9957 : *pt_E = (void*) E;
3100 9957 : return 0;
3101 : }
3102 413300 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
3103 : {
3104 325440 : struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
3105 325444 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
3106 325440 : D->sqr = &_sqr_montred;
3107 325440 : D->mul = &_mul_montred;
3108 325440 : D->mul2 = &_mul2_montred;
3109 325440 : E->N = N;
3110 325440 : E->inv = init_montdata(N);
3111 325434 : *pt_E = (void*) E;
3112 325434 : return 1;
3113 : }
3114 : else
3115 : {
3116 87860 : D->sqr = &_sqr_remii;
3117 87860 : D->mul = &_mul_remii;
3118 87860 : D->mul2 = &_mul2_remii;
3119 87860 : *pt_E = (void*) N;
3120 87860 : return 0;
3121 : }
3122 : }
3123 :
3124 : GEN
3125 983152 : Fp_powu(GEN A, ulong k, GEN N)
3126 : {
3127 983152 : long lN = lgefint(N), sA;
3128 : int base_is_2, use_montgomery;
3129 : muldata D;
3130 : void *E;
3131 : pari_sp av;
3132 :
3133 983152 : if (lN == 3) {
3134 76530 : ulong n = uel(N,2);
3135 76530 : return utoi( Fl_powu(umodiu(A, n), k, n) );
3136 : }
3137 906622 : if (k <= 2)
3138 : { /* frequent special cases */
3139 621620 : if (k == 2) return Fp_sqr(A,N);
3140 170354 : if (k == 1) return A;
3141 0 : if (k == 0) return gen_1;
3142 : }
3143 285002 : sA = signe(A)==-1 && odd(k);
3144 285002 : base_is_2 = 0;
3145 285002 : if (lgefint(A) == 3) switch(A[2])
3146 : {
3147 495 : case 1: return sA ? gen_m1 : gen_1;
3148 37408 : case 2: base_is_2 = 1; break;
3149 : }
3150 :
3151 : /* TODO: Move this out of here and use for general modular computations */
3152 284507 : av = avma;
3153 284507 : use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
3154 284507 : if (base_is_2)
3155 37408 : A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
3156 : else
3157 247099 : A = gen_powu_i(A, k, E, D.sqr, D.mul);
3158 284508 : if (use_montgomery)
3159 : {
3160 253364 : A = red_montgomery(A, N, ((struct montred *) E)->inv);
3161 253363 : if (cmpii(A, N) >= 0) A = subii(A,N);
3162 253363 : if (sA) A = subii(N, A);
3163 : }
3164 284507 : return gerepileuptoint(av, A);
3165 : }
3166 :
3167 : GEN
3168 21070 : Fp_pows(GEN A, long k, GEN N)
3169 : {
3170 21070 : if (lgefint(N) == 3) {
3171 7043 : ulong n = N[2];
3172 7043 : ulong a = umodiu(A, n);
3173 7043 : if (k < 0) {
3174 14 : a = Fl_inv(a, n);
3175 14 : k = -k;
3176 : }
3177 7043 : return utoi( Fl_powu(a, (ulong)k, n) );
3178 : }
3179 14027 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
3180 14027 : return Fp_powu(A, (ulong)k, N);
3181 : }
3182 :
3183 : /* A^K mod N */
3184 : GEN
3185 10935350 : Fp_pow(GEN A, GEN K, GEN N)
3186 : {
3187 10935350 : pari_sp av = avma;
3188 10935350 : long t,s, lN = lgefint(N), sA;
3189 : int base_is_2, use_montgomery;
3190 : GEN y;
3191 : muldata D;
3192 : void *E;
3193 :
3194 10935350 : s = signe(K);
3195 10935350 : if (!s)
3196 : {
3197 138426 : t = signe(remii(A,N)); avma = av;
3198 138426 : return t? gen_1: gen_0;
3199 : }
3200 10796924 : if (lN == 3)
3201 : {
3202 10500336 : ulong k, n = N[2], a = umodiu(A, n);
3203 10500336 : if (s < 0) a = Fl_inv(a, n);
3204 10500322 : if (a <= 1) return utoi(a); /* 0 or 1 */
3205 9730342 : if (lgefint(K) > 3)
3206 : { /* silly case : huge exponent, small modulus */
3207 21 : pari_warn(warner, "Mod(a,b)^n with n >> b : wasteful");
3208 21 : if (s > 0)
3209 : {
3210 14 : ulong d = ugcd(a, n);
3211 14 : if (d != 1)
3212 : { /* write n = n1 n2, with n2 maximal such that (n1,a) = 1 */
3213 7 : ulong n1 = u_ppo(n, d), n2 = n/n1;
3214 :
3215 7 : k = umodiu(K, eulerphiu(n1));
3216 : /* CRT: = a^K (mod n1), = 0 (mod n2)*/
3217 7 : return utoi( Fl_mul(Fl_powu(a, k, n1), n2 * Fl_inv(n2,n1), n) );
3218 : }
3219 : /* gcd(a,n) = 1 */
3220 7 : k = umodiu(K, eulerphiu(n));
3221 : }
3222 : else
3223 7 : k = umodiu(negi(K), eulerphiu(n));
3224 : }
3225 : else
3226 9730321 : k = uel(K,2);
3227 9730335 : return utoi(Fl_powu(a, k, n));
3228 : }
3229 :
3230 296588 : if (s < 0) y = Fp_inv(A,N);
3231 : else
3232 : {
3233 296209 : y = modii(A,N);
3234 296204 : if (!signe(y)) { avma = av; return gen_0; }
3235 : }
3236 296576 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
3237 :
3238 138811 : base_is_2 = 0;
3239 138811 : sA = signe(y)==-1 && mod2(K);
3240 138811 : if (lgefint(y) == 3) switch(y[2])
3241 : {
3242 61 : case 1: return sA ? gen_m1 : gen_1;
3243 90246 : case 2: base_is_2 = 1; break;
3244 : }
3245 :
3246 : /* TODO: Move this out of here and use for general modular computations */
3247 138750 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
3248 138744 : if (base_is_2)
3249 90246 : y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
3250 : else
3251 48498 : y = gen_pow_i(y, K, E, D.sqr, D.mul);
3252 138751 : if (use_montgomery)
3253 : {
3254 72078 : y = red_montgomery(y, N, ((struct montred *) E)->inv);
3255 72077 : if (cmpii(y,N) >= 0) y = subii(y,N);
3256 72079 : if (sA) y = subii(N, y);
3257 : }
3258 138752 : return gerepileuptoint(av,y);
3259 : }
3260 :
3261 : static GEN
3262 1301976 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
3263 :
3264 : static GEN
3265 948 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
3266 :
3267 : static GEN
3268 104 : _Fp_one(void *E) { (void) E; return gen_1; }
3269 :
3270 : GEN
3271 1624 : Fp_powers(GEN x, long n, GEN p)
3272 : {
3273 1624 : if (lgefint(p) == 3)
3274 1520 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
3275 104 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
3276 : }
3277 :
3278 : static GEN
3279 8196880 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
3280 :
3281 : static GEN
3282 946 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
3283 :
3284 : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
3285 :
3286 : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
3287 : equalii,equali1,Fp_easylog};
3288 :
3289 : static GEN
3290 971871 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
3291 :
3292 : static GEN
3293 997162 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
3294 :
3295 : static GEN
3296 107441 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
3297 :
3298 : static GEN
3299 1236143 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
3300 :
3301 : static GEN
3302 18501 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
3303 :
3304 : static int
3305 521974 : _Fp_equal0(GEN x) { return signe(x)==0; }
3306 :
3307 : static GEN
3308 82036 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
3309 :
3310 : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
3311 : _Fp_inv,_Fp_equal0,_Fp_s};
3312 :
3313 4578 : const struct bb_field *get_Fp_field(void **E, GEN p)
3314 : {
3315 4578 : *E = (void*)p; return &Fp_field;
3316 : }
3317 :
3318 : /*********************************************************************/
3319 : /** **/
3320 : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
3321 : /** **/
3322 : /*********************************************************************/
3323 : ulong
3324 5341 : Fl_order(ulong a, ulong o, ulong p)
3325 : {
3326 5341 : pari_sp av = avma;
3327 : GEN m, P, E;
3328 : long i;
3329 5341 : if (!o) o = p-1;
3330 5341 : m = factoru(o);
3331 5341 : P = gel(m,1);
3332 5341 : E = gel(m,2);
3333 13412 : for (i = lg(P)-1; i; i--)
3334 : {
3335 8071 : ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
3336 8071 : if (y == 1) o = t;
3337 7504 : else for (j = 1; j < e; j++)
3338 : {
3339 2303 : y = Fl_powu(y, l, p);
3340 2303 : if (y == 1) { o = t * upowuu(l, j); break; }
3341 : }
3342 : }
3343 5341 : avma = av; return o;
3344 : }
3345 :
3346 : /*Find the exact order of a assuming a^o==1*/
3347 : GEN
3348 10739 : Fp_order(GEN a, GEN o, GEN p) {
3349 10739 : if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
3350 : {
3351 21 : ulong pp = p[2], oo = (o && lgefint(o)==3)? o[2]: pp-1;
3352 21 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
3353 : }
3354 10718 : return gen_order(a, o, (void*)p, &Fp_star);
3355 : }
3356 : GEN
3357 56 : Fp_factored_order(GEN a, GEN o, GEN p)
3358 56 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
3359 :
3360 : /* return order of a mod p^e, e > 0, pe = p^e */
3361 : static GEN
3362 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
3363 : {
3364 : GEN ap, op;
3365 70 : if (absequaliu(p, 2))
3366 : {
3367 56 : if (e == 1) return gen_1;
3368 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
3369 49 : if (mod4(a) == 1)
3370 14 : op = gen_1;
3371 : else {
3372 35 : op = gen_2;
3373 35 : a = Fp_sqr(a, pe);
3374 : }
3375 : } else {
3376 14 : ap = (e == 1)? a: remii(a,p);
3377 14 : op = Fp_order(ap, subiu(p,1), p);
3378 14 : if (e == 1) return op;
3379 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
3380 : }
3381 49 : if (equali1(a)) return op;
3382 7 : return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
3383 : }
3384 :
3385 : GEN
3386 63 : znorder(GEN x, GEN o)
3387 : {
3388 63 : pari_sp av = avma;
3389 : GEN b, a;
3390 :
3391 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
3392 56 : b = gel(x,1); a = gel(x,2);
3393 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
3394 49 : if (!o)
3395 : {
3396 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
3397 35 : long i, l = lg(P);
3398 35 : o = gen_1;
3399 70 : for (i = 1; i < l; i++)
3400 : {
3401 35 : GEN p = gel(P,i);
3402 35 : long e = itos(gel(E,i));
3403 :
3404 35 : if (l == 2)
3405 35 : o = Zp_order(a, p, e, b);
3406 : else {
3407 0 : GEN pe = powiu(p,e);
3408 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
3409 : }
3410 : }
3411 35 : return gerepileuptoint(av, o);
3412 : }
3413 14 : return Fp_order(a, o, b);
3414 : }
3415 : GEN
3416 0 : order(GEN x) { return znorder(x, NULL); }
3417 :
3418 : /*********************************************************************/
3419 : /** **/
3420 : /** DISCRETE LOGARITHM in (Z/nZ)* **/
3421 : /** **/
3422 : /*********************************************************************/
3423 : static GEN
3424 67189 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
3425 : {
3426 67189 : pari_sp av = avma;
3427 : GEN h1, h2, F, G;
3428 67189 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return NULL;
3429 40421 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
3430 : {
3431 256 : GEN M = cgetg(3, t_MAT);
3432 256 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
3433 256 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
3434 256 : return gerepileupto(av, M);
3435 : }
3436 40165 : avma = av; return NULL;
3437 : }
3438 :
3439 : static GEN
3440 67189 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
3441 : {
3442 : GEN rel;
3443 : do
3444 : {
3445 67189 : (*e)++; *g = Fp_mul(*g, b, p);
3446 67189 : rel = Fp_log_halfgcd(bnd, C, *g, p);
3447 67189 : } while (!rel);
3448 256 : return rel;
3449 : }
3450 :
3451 : struct Fp_log_rel
3452 : {
3453 : GEN rel;
3454 : ulong prmax;
3455 : long nbrel, nbmax, nbgen;
3456 : };
3457 :
3458 : /* add u^e */
3459 : static void
3460 2583 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
3461 : {
3462 2583 : pari_sp av = avma;
3463 2583 : long off = r->prmax+1;
3464 2583 : GEN F = cgetg(3, t_MAT);
3465 2583 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
3466 2583 : gel(F,2) = vecsmall_append(gel(z,2), e);
3467 2583 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3468 2583 : }
3469 :
3470 : /* add u^-1 v^-1 */
3471 : static void
3472 99869 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
3473 : {
3474 99869 : pari_sp av = avma;
3475 99869 : long off = r->prmax+1;
3476 99869 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
3477 99869 : GEN F = cgetg(3, t_MAT);
3478 99869 : gel(F,1) = vecsmall_concat(gel(z,1), P);
3479 99869 : gel(F,2) = vecsmall_concat(gel(z,2), E);
3480 99869 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3481 99869 : }
3482 :
3483 : /*
3484 : Let p=C^2+c
3485 : Solve h = (C+x)*(C+a)-p = 0 [mod l]
3486 : h= -c+x*(C+a)+C*a = 0 [mod l]
3487 : x = (c-C*a)/(C+a) [mod l]
3488 : h = -c+C*(x+a)+a*x
3489 : */
3490 :
3491 : GEN
3492 39020 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
3493 : {
3494 39020 : pari_sp ltop = avma;
3495 39020 : long th = expi(mulis(C,a)), n = lg(pi)-1;
3496 : long i, j;
3497 39032 : GEN sieve = zero_zv(a+2)+1;
3498 39037 : GEN L = cgetg(1+a+2, t_VEC);
3499 39032 : pari_sp av = avma;
3500 39032 : long rel = 1;
3501 : GEN z, h;
3502 39032 : h = addis(C,a);
3503 39008 : if ((z = Z_issmooth_fact(h, prmax)))
3504 : {
3505 2415 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
3506 2414 : av = avma;
3507 : }
3508 17502744 : for(i=1; i<=n; i++)
3509 : {
3510 17463728 : ulong li = pi[i], s = sz[i], al = a % li;
3511 17463728 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
3512 18053486 : if (!iv) continue;
3513 17445955 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
3514 79432853 : for(j = u; j<=a; j+=li)
3515 62480370 : sieve[j] += s;
3516 : }
3517 39016 : th = th - expu(th)-1;
3518 27910680 : for(j=0; j<a; j++)
3519 27871642 : if (sieve[j]>=th)
3520 : {
3521 114617 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
3522 111607 : if ((z = Z_issmooth_fact(h, prmax)))
3523 : {
3524 104953 : gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
3525 104973 : av = avma;
3526 6660 : } else avma = av;
3527 : }
3528 : /* j = a */
3529 39038 : if (sieve[a]>=th)
3530 : {
3531 420 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
3532 420 : if ((z = Z_issmooth_fact(h, prmax)))
3533 : {
3534 343 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
3535 343 : av = avma;
3536 : }
3537 : }
3538 39038 : setlg(L, rel);
3539 39039 : return gerepilecopy(ltop, L);
3540 : }
3541 :
3542 : static long
3543 49 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
3544 : {
3545 : struct pari_mt pt;
3546 49 : long i, j, nb = 0;
3547 49 : GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
3548 : mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
3549 49 : long running, pending = 0;
3550 49 : GEN W = zerovec(r->nbgen);
3551 49 : mt_queue_start_lim(&pt, worker, r->nbgen);
3552 39193 : for (i = 0; (running = (i < r->nbgen)) || pending; i++)
3553 : {
3554 : GEN done;
3555 : long idx;
3556 39144 : mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
3557 39144 : done = mt_queue_get(&pt, &idx, &pending);
3558 39144 : if (!done || lg(done)==1) continue;
3559 34146 : gel(W, idx+1) = done;
3560 34146 : nb += lg(done)-1;
3561 34146 : if (DEBUGLEVEL && (i&127)==0)
3562 0 : err_printf("%ld%% ",100*nb/r->nbmax);
3563 : }
3564 49 : mt_queue_end(&pt);
3565 37870 : for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
3566 : {
3567 : long ll, m;
3568 37821 : GEN L = gel(W,j);
3569 37821 : if (isintzero(L)) continue;
3570 32942 : ll = lg(L);
3571 135394 : for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
3572 : {
3573 102452 : GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
3574 102452 : if (v[1] == 1)
3575 2583 : addifsmooth1(r, h, v[2], v[3]);
3576 : else
3577 99869 : addifsmooth2(r, h, v[2], v[3]);
3578 : }
3579 : }
3580 49 : return j;
3581 : }
3582 :
3583 : static GEN
3584 525 : ECP_psi(GEN x, GEN y)
3585 : {
3586 525 : long prec = realprec(x);
3587 525 : GEN lx = glog(x, prec), ly = glog(y, prec);
3588 525 : GEN u = gdiv(lx, ly);
3589 525 : return gpow(u, gneg(u),prec);
3590 : }
3591 :
3592 : struct computeG
3593 : {
3594 : GEN C;
3595 : long bnd, nbi;
3596 : };
3597 :
3598 : static GEN
3599 525 : _computeG(void *E, GEN gen)
3600 : {
3601 525 : struct computeG * d = (struct computeG *) E;
3602 525 : GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
3603 525 : return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
3604 : }
3605 :
3606 : static long
3607 49 : compute_nbgen(GEN C, long bnd, long nbi)
3608 : {
3609 : struct computeG d;
3610 49 : d.C = shifti(C, 1);
3611 49 : d.bnd = bnd;
3612 49 : d.nbi = nbi;
3613 49 : return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
3614 : }
3615 :
3616 : static GEN
3617 1367 : _psi(void*E, GEN y)
3618 : {
3619 1367 : GEN lx = (GEN) E;
3620 1367 : long prec = realprec(lx);
3621 1367 : GEN ly = glog(y, prec);
3622 1367 : GEN u = gdiv(lx, ly);
3623 1367 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
3624 : }
3625 :
3626 : static GEN
3627 49 : opt_param(GEN x, long prec)
3628 : {
3629 49 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
3630 : }
3631 :
3632 : static GEN
3633 49 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
3634 : {
3635 49 : pari_sp av = avma;
3636 49 : long lM = lg(M)-1, nbcol = lM;
3637 : for (;;)
3638 : {
3639 84 : GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
3640 84 : long i, f=0;
3641 84 : long l = lg(K), lm = lgefint(m);
3642 84 : GEN idx = diviiexact(subiu(p,1),m), g;
3643 : pari_timer ti;
3644 84 : if (DEBUGLEVEL) timer_start(&ti);
3645 140 : for(i=1; i<l; i++)
3646 140 : if (signe(gel(K,i)))
3647 84 : break;
3648 84 : g = Fp_pow(utoi(i), idx, p);
3649 84 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
3650 130438 : for(i=1; i<l; i++)
3651 : {
3652 130354 : GEN k = gel(K,i);
3653 130354 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
3654 130354 : if (signe(k)==0 || !equalii(Fp_pow(g, k, p),
3655 : Fp_pow(j, idx, p)))
3656 80542 : gel(K,i) = cgetineg(lm);
3657 : else
3658 49812 : f++;
3659 : }
3660 84 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld logs", f);
3661 133 : if(f > (nbg>>1)) return gerepileupto(av, K);
3662 9877 : for(i=1; i<=nbcol; i++)
3663 : {
3664 9842 : long a = 1+random_Fl(lM);
3665 9842 : swap(gel(M,a),gel(M,i));
3666 : }
3667 35 : if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
3668 35 : }
3669 : }
3670 :
3671 : static GEN
3672 98 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
3673 : {
3674 98 : pari_sp av=avma;
3675 98 : GEN aa = gen_1;
3676 98 : long AV = 0;
3677 : for(;;)
3678 : {
3679 256 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
3680 256 : GEN F = gel(A,1), E = gel(A,2);
3681 256 : GEN Ao = gen_0;
3682 256 : long i, l = lg(F);
3683 1304 : for(i=1; i<l; i++)
3684 : {
3685 1206 : GEN Ki = gel(K,F[i]);
3686 1206 : if (signe(Ki)<0) break;
3687 1048 : Ao = addii(Ao, mulis(Ki, E[i]));
3688 : }
3689 354 : if (i==l) return Fp_div(Ao, utoi(AV), m);
3690 158 : aa = gerepileuptoint(av, aa);
3691 158 : }
3692 : }
3693 :
3694 : static GEN
3695 49 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
3696 : {
3697 49 : pari_sp av = avma, av2;
3698 49 : long i, j, nbi, nbr = 0, nbrow, nbg;
3699 : GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
3700 : pari_timer ti;
3701 : struct Fp_log_rel r;
3702 49 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
3703 49 : ulong bnd = 4*bnds;
3704 49 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
3705 :
3706 49 : p_1 = subiu(p,1);
3707 49 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
3708 0 : m = diviiexact(p_1, Z_ppo(p_1, m));
3709 49 : pr = primes_upto_zv(bnd);
3710 49 : nbi = lg(pr)-1;
3711 49 : C = sqrtremi(p, &c);
3712 49 : av2 = avma;
3713 12236 : for (i = 1; i <= nbi; ++i)
3714 : {
3715 12187 : ulong lp = pr[i];
3716 37793 : while (lp <= bnd)
3717 : {
3718 13419 : nbr++;
3719 13419 : lp *= pr[i];
3720 : }
3721 : }
3722 49 : pi = cgetg(nbr+1,t_VECSMALL);
3723 49 : Ci = cgetg(nbr+1,t_VECSMALL);
3724 49 : ci = cgetg(nbr+1,t_VECSMALL);
3725 49 : sz = cgetg(nbr+1,t_VECSMALL);
3726 12236 : for (i = 1, j = 1; i <= nbi; ++i)
3727 : {
3728 12187 : ulong lp = pr[i], sp = expu(2*lp-1);
3729 37793 : while (lp <= bnd)
3730 : {
3731 13419 : pi[j] = lp;
3732 13419 : Ci[j] = umodiu(C, lp);
3733 13419 : ci[j] = umodiu(c, lp);
3734 13419 : sz[j] = sp;
3735 13419 : lp *= pr[i];
3736 13419 : j++;
3737 : }
3738 : }
3739 49 : r.nbrel = 0;
3740 49 : r.nbgen = compute_nbgen(C, bnd, nbi);
3741 49 : r.nbmax = 2*(nbi+r.nbgen);
3742 49 : r.rel = cgetg(r.nbmax+1,t_VEC);
3743 49 : r.prmax = pr[nbi];
3744 49 : if (DEBUGLEVEL)
3745 : {
3746 0 : err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
3747 0 : timer_start(&ti);
3748 : }
3749 49 : nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
3750 49 : nbrow = r.prmax + nbg;
3751 49 : if (DEBUGLEVEL)
3752 : {
3753 0 : err_printf("\n");
3754 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
3755 : }
3756 49 : setlg(r.rel,r.nbrel+1);
3757 49 : r.rel = gerepilecopy(av2, r.rel);
3758 49 : K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
3759 49 : if (DEBUGLEVEL) timer_start(&ti);
3760 49 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
3761 49 : if (DEBUGLEVEL) timer_printf(&ti," log element");
3762 49 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
3763 49 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
3764 49 : d = gcdii(Ao,Bo);
3765 49 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
3766 49 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
3767 49 : return gerepileuptoint(av, l);
3768 : }
3769 :
3770 : static int
3771 1807805 : Fp_log_use_index(long e, long p)
3772 : {
3773 1807805 : return (e >= 27 && 20*(p+6)<=e*e);
3774 : }
3775 :
3776 : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
3777 : static GEN
3778 2885731 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
3779 : {
3780 2885731 : pari_sp av = avma;
3781 2885731 : GEN p = (GEN)E;
3782 : /* assume a reduced mod p, p not necessarily prime */
3783 2885731 : if (equali1(a)) return gen_0;
3784 : /* p > 2 */
3785 1798288 : if (equalii(subiu(p,1), a)) /* -1 */
3786 : {
3787 : pari_sp av2;
3788 : GEN t;
3789 571828 : ord = get_arith_Z(ord);
3790 571828 : if (mpodd(ord)) { avma = av; return cgetg(1, t_VEC); } /* no solution */
3791 571814 : t = shifti(ord,-1); /* only possible solution */
3792 571814 : av2 = avma;
3793 571814 : if (!equalii(Fp_pow(g, t, p), a)) { avma = av; return cgetg(1, t_VEC); }
3794 571632 : avma = av2; return gerepileuptoint(av, t);
3795 : }
3796 1226460 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
3797 49 : return Fp_log_index(a, g, ord, p);
3798 1226411 : avma = av; return NULL; /* not easy */
3799 : }
3800 :
3801 : GEN
3802 1707967 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
3803 : {
3804 1707967 : GEN v = get_arith_ZZM(ord);
3805 1707939 : GEN F = gmael(v,2,1);
3806 1707939 : long lF = lg(F)-1, lmax;
3807 1707939 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
3808 1591748 : lmax = expi(gel(F,lF));
3809 1591748 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
3810 49 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
3811 1591748 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
3812 : }
3813 :
3814 : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
3815 : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
3816 : static GEN
3817 112 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
3818 : {
3819 112 : long l = lg(P) - 1, e = E[l];
3820 112 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
3821 : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
3822 :
3823 112 : if (l == 1) {
3824 84 : hpe = h;
3825 84 : gpe = g;
3826 : } else {
3827 28 : hpe = modii(h, pe);
3828 28 : gpe = modii(g, pe);
3829 : }
3830 112 : if (e == 1) {
3831 28 : hp = hpe;
3832 28 : gp = gpe;
3833 : } else {
3834 84 : hp = remii(hpe, p);
3835 84 : gp = remii(gpe, p);
3836 : }
3837 112 : if (hp == gen_0 || gp == gen_0) return NULL;
3838 91 : if (absequaliu(p, 2))
3839 : {
3840 35 : GEN n = int2n(e);
3841 35 : ogpe = Zp_order(gpe, gen_2, e, n);
3842 35 : a = Fp_log(hpe, gpe, ogpe, n);
3843 35 : if (typ(a) != t_INT) return NULL;
3844 : }
3845 : else
3846 : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
3847 : is trivial */
3848 : /* [order(gp), factor(order(gp))] */
3849 56 : GEN v = Fp_factored_order(gp, subiu(p,1), p);
3850 56 : GEN ogp = gel(v,1);
3851 56 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
3852 56 : a = Fp_log(hp, gp, v, p);
3853 56 : if (typ(a) != t_INT) return NULL;
3854 56 : if (e == 1) ogpe = ogp;
3855 : else
3856 : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
3857 : /* use p-adic log: O(log p + e) mul*/
3858 : long vpogpe, vpohpe;
3859 :
3860 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
3861 28 : gpe = Fp_pow(gpe, ogp, pe);
3862 : /* g,h = 1 mod p; compute b s.t. h = g^b */
3863 :
3864 : /* v_p(order g mod pe) */
3865 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
3866 : /* v_p(order h mod pe) */
3867 28 : vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
3868 28 : if (vpohpe > vpogpe) return NULL;
3869 :
3870 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
3871 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
3872 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
3873 28 : a = addii(a, mulii(ogp, padic_to_Q(b)));
3874 : }
3875 : }
3876 : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
3877 77 : if (l == 1) return a;
3878 :
3879 28 : N = diviiexact(N, pe); /* make N coprime to p */
3880 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
3881 28 : g = Fp_pow(g, modii(ogpe, phi), N);
3882 28 : setlg(P, l); /* remove last element */
3883 28 : setlg(E, l);
3884 28 : b = znlog_rec(h, g, N, P, E, PHI);
3885 28 : if (!b) return NULL;
3886 28 : return addmulii(a, b, ogpe);
3887 : }
3888 :
3889 : static GEN
3890 84 : get_PHI(GEN P, GEN E)
3891 : {
3892 84 : long i, l = lg(P);
3893 84 : GEN PHI = cgetg(l, t_VEC);
3894 84 : gel(PHI,1) = gen_1;
3895 112 : for (i=1; i<l-1; i++)
3896 : {
3897 28 : GEN t, p = gel(P,i);
3898 28 : long e = E[i];
3899 28 : t = mulii(powiu(p, e-1), subiu(p,1));
3900 28 : if (i > 1) t = mulii(t, gel(PHI,i));
3901 28 : gel(PHI,i+1) = t;
3902 : }
3903 84 : return PHI;
3904 : }
3905 :
3906 : GEN
3907 224 : znlog(GEN h, GEN g, GEN o)
3908 : {
3909 224 : pari_sp av = avma;
3910 : GEN N, fa, P, E, x;
3911 224 : switch (typ(g))
3912 : {
3913 : case t_PADIC:
3914 : {
3915 28 : GEN p = gel(g,2);
3916 28 : long v = valp(g);
3917 28 : if (v < 0) pari_err_DIM("znlog");
3918 28 : if (v > 0) {
3919 0 : long k = gvaluation(h, p);
3920 0 : if (k % v) return cgetg(1,t_VEC);
3921 0 : k /= v;
3922 0 : if (!gequal(h, gpowgs(g,k))) { avma = av; return cgetg(1,t_VEC); }
3923 0 : avma = av; return stoi(k);
3924 : }
3925 28 : N = gel(g,3);
3926 28 : g = Rg_to_Fp(g, N);
3927 28 : break;
3928 : }
3929 : case t_INTMOD:
3930 189 : N = gel(g,1);
3931 189 : g = gel(g,2); break;
3932 7 : default: pari_err_TYPE("znlog", g);
3933 : return NULL; /* LCOV_EXCL_LINE */
3934 : }
3935 217 : if (equali1(N)) { avma = av; return gen_0; }
3936 217 : h = Rg_to_Fp(h, N);
3937 210 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
3938 84 : fa = Z_factor(N);
3939 84 : P = gel(fa,1);
3940 84 : E = vec_to_vecsmall(gel(fa,2));
3941 84 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
3942 84 : if (!x) { avma = av; return cgetg(1,t_VEC); }
3943 49 : return gerepileuptoint(av, x);
3944 : }
3945 :
3946 : GEN
3947 63491 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
3948 : {
3949 63491 : a = modii(a,p);
3950 63491 : if (!signe(a))
3951 : {
3952 48118 : if (zeta) *zeta = gen_1;
3953 48118 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
3954 48111 : return gen_0;
3955 : }
3956 15373 : if (absequaliu(n,2))
3957 : {
3958 2282 : if (zeta) *zeta = subiu(p,1);
3959 2282 : return Fp_sqrt(a,p);
3960 : }
3961 13091 : return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
3962 : }
3963 :
3964 : /*********************************************************************/
3965 : /** **/
3966 : /** FUNDAMENTAL DISCRIMINANTS **/
3967 : /** **/
3968 : /*********************************************************************/
3969 : static int
3970 1407 : fa_isfundamental(GEN F)
3971 : {
3972 1407 : GEN P = gel(F,1), E = gel(F,2);
3973 1407 : long i, s, l = lg(P);
3974 :
3975 1407 : if (l == 1) return 1;
3976 1400 : s = signe(gel(P,1)); /* = signe(x) */
3977 1400 : if (!s) return 0;
3978 1393 : if (s < 0) { l--; P = vecslice(P,2,l); E = vecslice(E,2,l); }
3979 1393 : if (l == 1) return 0;
3980 1386 : if (!absequaliu(gel(P,1), 2))
3981 686 : i = 1; /* need x = 1 mod 4 */
3982 : else
3983 : {
3984 700 : i = 2;
3985 700 : switch(itou(gel(E,1)))
3986 : {
3987 182 : case 2: s = -s; break; /* need x/4 = 3 mod 4 */
3988 84 : case 3: s = 0; break; /* no condition mod 4 */
3989 434 : default: return 0;
3990 : }
3991 : }
3992 1974 : for(; i < l; i++)
3993 : {
3994 1190 : if (!equali1(gel(E,i))) return 0;
3995 1022 : if (s && Mod4(gel(P,i)) == 3) s = -s;
3996 : }
3997 784 : return s >= 0;
3998 : }
3999 : long
4000 17185 : isfundamental(GEN x)
4001 : {
4002 17185 : if (typ(x) != t_INT)
4003 : {
4004 1407 : pari_sp av = avma;
4005 1407 : int v = fa_isfundamental(check_arith_all(x,"isfundamental"));
4006 1407 : avma = av; return v;
4007 : }
4008 15778 : return Z_isfundamental(x);
4009 : }
4010 :
4011 : /* x fundamental ? */
4012 : long
4013 10240 : uposisfundamental(ulong x)
4014 : {
4015 10240 : ulong r = x & 15; /* x mod 16 */
4016 10240 : if (!r) return 0;
4017 9666 : switch(r & 3)
4018 : { /* x mod 4 */
4019 1870 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
4020 3018 : case 1: return uissquarefree(x);
4021 4778 : default: return 0;
4022 : }
4023 : }
4024 : /* -x fundamental ? */
4025 : long
4026 21043 : unegisfundamental(ulong x)
4027 : {
4028 21043 : ulong r = x & 15; /* x mod 16 */
4029 21043 : if (!r) return 0;
4030 20042 : switch(r & 3)
4031 : { /* x mod 4 */
4032 3578 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
4033 9460 : case 3: return uissquarefree(x);
4034 7004 : default: return 0;
4035 : }
4036 : }
4037 : long
4038 10451 : sisfundamental(long x)
4039 10451 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
4040 :
4041 : long
4042 16338 : Z_isfundamental(GEN x)
4043 : {
4044 : long r;
4045 16338 : switch(lgefint(x))
4046 : {
4047 7 : case 2: return 0;
4048 35908 : case 3: return signe(x) < 0? unegisfundamental(x[2])
4049 21587 : : uposisfundamental(x[2]);
4050 : }
4051 2010 : r = mod16(x);
4052 2010 : if (!r) return 0;
4053 1884 : if ((r & 3) == 0)
4054 : {
4055 : pari_sp av;
4056 376 : r >>= 2; /* |x|/4 mod 4 */
4057 376 : if (signe(x) < 0) r = 4-r;
4058 376 : if (r == 1) return 0;
4059 250 : av = avma;
4060 250 : r = Z_issquarefree( shifti(x,-2) );
4061 250 : avma = av; return r;
4062 : }
4063 1508 : r &= 3; /* |x| mod 4 */
4064 1508 : if (signe(x) < 0) r = 4-r;
4065 1508 : return (r==1) ? Z_issquarefree(x) : 0;
4066 : }
4067 :
4068 : GEN
4069 7 : quaddisc(GEN x)
4070 : {
4071 7 : const pari_sp av = avma;
4072 7 : long i,r,tx=typ(x);
4073 : GEN P,E,f,s;
4074 :
4075 7 : if (!is_rational_t(tx)) pari_err_TYPE("quaddisc",x);
4076 7 : f = factor(x);
4077 7 : P = gel(f,1);
4078 7 : E = gel(f,2); s = gen_1;
4079 35 : for (i=1; i<lg(P); i++)
4080 28 : if (odd(mael(E,i,2))) s = mulii(s,gel(P,i));
4081 7 : r = mod4(s); if (gsigne(x) < 0) r = 4-r;
4082 7 : if (r>1) s = shifti(s,2);
4083 7 : return gerepileuptoint(av, s);
4084 : }
4085 :
4086 : /*********************************************************************/
4087 : /** **/
4088 : /** FACTORIAL **/
4089 : /** **/
4090 : /*********************************************************************/
4091 : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
4092 : * first is slower ... ] */
4093 : GEN
4094 983921 : mulu_interval(ulong a, ulong b)
4095 : {
4096 983921 : pari_sp av = avma;
4097 : ulong k, l, N, n;
4098 : long lx;
4099 : GEN x;
4100 :
4101 983921 : if (!a) return gen_0;
4102 983921 : n = b - a + 1;
4103 983921 : if (n < 61)
4104 : {
4105 975153 : x = utoi(a);
4106 975158 : for (k=a+1; k<=b; k++) x = mului(k,x);
4107 975158 : return gerepileuptoint(av, x);
4108 : }
4109 8768 : lx = 1; x = cgetg(2 + n/2, t_VEC);
4110 8768 : N = b + a;
4111 1007781 : for (k = a;; k++)
4112 : {
4113 1007781 : l = N - k; if (l <= k) break;
4114 999013 : gel(x,lx++) = muluu(k,l);
4115 999013 : }
4116 8768 : if (l == k) gel(x,lx++) = utoipos(k);
4117 8768 : setlg(x, lx);
4118 8768 : return gerepileuptoint(av, ZV_prod(x));
4119 : }
4120 : GEN
4121 448 : muls_interval(long a, long b)
4122 : {
4123 448 : pari_sp av = avma;
4124 448 : long lx, k, l, N, n = b - a + 1;
4125 : GEN x;
4126 :
4127 448 : if (a <= 0 && b >= 0) return gen_0;
4128 287 : if (n < 61)
4129 : {
4130 287 : x = stoi(a);
4131 287 : for (k=a+1; k<=b; k++) x = mulsi(k,x);
4132 287 : return gerepileuptoint(av, x);
4133 : }
4134 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
4135 0 : N = b + a;
4136 0 : for (k = a;; k++)
4137 : {
4138 0 : l = N - k; if (l <= k) break;
4139 0 : gel(x,lx++) = mulss(k,l);
4140 0 : }
4141 0 : if (l == k) gel(x,lx++) = stoi(k);
4142 0 : setlg(x, lx);
4143 0 : return gerepileuptoint(av, ZV_prod(x));
4144 : }
4145 :
4146 : GEN
4147 2333859 : mpfact(long n)
4148 : {
4149 2333859 : if (n < 2)
4150 : {
4151 1470338 : if (n < 0) pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
4152 1470338 : return gen_1;
4153 : }
4154 863521 : return mulu_interval(2UL, (ulong)n);
4155 : }
4156 :
4157 : /*******************************************************************/
4158 : /** **/
4159 : /** LUCAS & FIBONACCI **/
4160 : /** **/
4161 : /*******************************************************************/
4162 : static void
4163 56 : lucas(ulong n, GEN *a, GEN *b)
4164 : {
4165 : GEN z, t, zt;
4166 112 : if (!n) { *a = gen_2; *b = gen_1; return; }
4167 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
4168 49 : switch(n & 3) {
4169 14 : case 0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
4170 14 : case 1: *a = subiu(zt,1); *b = addiu(sqri(t),2); break;
4171 7 : case 2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
4172 14 : case 3: *a = addiu(zt,1); *b = subiu(sqri(t),2);
4173 : }
4174 : }
4175 :
4176 : GEN
4177 7 : fibo(long n)
4178 : {
4179 7 : pari_sp av = avma;
4180 : GEN a, b;
4181 7 : if (!n) return gen_0;
4182 7 : lucas((ulong)(labs(n)-1), &a, &b);
4183 7 : a = diviuexact(addii(shifti(a,1),b), 5);
4184 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
4185 7 : return gerepileuptoint(av, a);
4186 : }
4187 :
4188 : /*******************************************************************/
4189 : /* */
4190 : /* CONTINUED FRACTIONS */
4191 : /* */
4192 : /*******************************************************************/
4193 : static GEN
4194 359403 : icopy_lg(GEN x, long l)
4195 : {
4196 359403 : long lx = lgefint(x);
4197 : GEN y;
4198 :
4199 359403 : if (lx >= l) return icopy(x);
4200 35 : y = cgeti(l); affii(x, y); return y;
4201 : }
4202 :
4203 : /* continued fraction of a/b. If y != NULL, stop when partial quotients
4204 : * differ from y */
4205 : static GEN
4206 359692 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
4207 : {
4208 : GEN z, c;
4209 359692 : ulong i, l, ly = lgefint(b);
4210 :
4211 : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
4212 359692 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
4213 359692 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
4214 359692 : if (l > LGBITS) l = LGBITS;
4215 :
4216 359692 : z = cgetg(l,t_VEC);
4217 359692 : l--;
4218 359692 : if (y) {
4219 289 : pari_sp av = avma;
4220 289 : if (l >= (ulong)lg(y)) l = lg(y)-1;
4221 19467 : for (i = 1; i <= l; i++)
4222 : {
4223 19289 : GEN q = gel(y,i);
4224 19289 : gel(z,i) = q;
4225 19289 : c = b; if (!gequal1(q)) c = mulii(q, b);
4226 19289 : c = subii(a, c);
4227 19289 : if (signe(c) < 0)
4228 : { /* partial quotient too large */
4229 110 : c = addii(c, b);
4230 110 : if (signe(c) >= 0) i++; /* by 1 */
4231 110 : break;
4232 : }
4233 19179 : if (cmpii(c, b) >= 0)
4234 : { /* partial quotient too small */
4235 1 : c = subii(c, b);
4236 1 : if (cmpii(c, b) < 0) {
4237 : /* by 1. If next quotient is 1 in y, add 1 */
4238 0 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
4239 0 : i++;
4240 : }
4241 1 : break;
4242 : }
4243 19178 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
4244 19178 : a = b; b = c;
4245 : }
4246 : } else {
4247 359403 : a = icopy_lg(a, ly);
4248 359403 : b = icopy(b);
4249 1255197 : for (i = 1; i <= l; i++)
4250 : {
4251 1254933 : gel(z,i) = truedvmdii(a,b,&c);
4252 1254933 : if (c == gen_0) { i++; break; }
4253 895794 : affii(c, a); cgiv(c); c = a;
4254 895794 : a = b; b = c;
4255 : }
4256 : }
4257 359692 : i--;
4258 359692 : if (i > 1 && gequal1(gel(z,i)))
4259 : {
4260 85 : cgiv(gel(z,i)); --i;
4261 85 : gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
4262 : }
4263 359692 : setlg(z,i+1); return z;
4264 : }
4265 :
4266 : static GEN
4267 0 : sersfcont(GEN a, GEN b, long k)
4268 : {
4269 0 : long i, l = typ(a) == t_POL? lg(a): 3;
4270 : GEN y, c;
4271 0 : if (lg(b) > l) l = lg(b);
4272 0 : if (k > 0 && l > k+1) l = k+1;
4273 0 : y = cgetg(l,t_VEC);
4274 0 : for (i=1; i<l; i++)
4275 : {
4276 0 : gel(y,i) = poldivrem(a,b,&c);
4277 0 : if (gequal0(c)) { i++; break; }
4278 0 : a = b; b = c;
4279 : }
4280 0 : setlg(y, i); return y;
4281 : }
4282 :
4283 : GEN
4284 359718 : gboundcf(GEN x, long k)
4285 : {
4286 : pari_sp av;
4287 359718 : long tx = typ(x), e;
4288 : GEN y, a, b, c;
4289 :
4290 359718 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
4291 359711 : if (is_scalar_t(tx))
4292 : {
4293 359711 : if (gequal0(x)) return mkvec(gen_0);
4294 359410 : switch(tx)
4295 : {
4296 0 : case t_INT: return mkveccopy(x);
4297 : case t_REAL:
4298 296 : av = avma;
4299 296 : c = mantissa_real(x,&e);
4300 296 : if (e < 0) pari_err_PREC("gboundcf");
4301 289 : y = int2n(e);
4302 289 : a = Qsfcont(c,y, NULL, k);
4303 289 : b = addsi(signe(x), c);
4304 289 : return gerepilecopy(av, Qsfcont(b,y, a, k));
4305 :
4306 : case t_FRAC:
4307 359114 : av = avma;
4308 359114 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
4309 : }
4310 0 : pari_err_TYPE("gboundcf",x);
4311 : }
4312 :
4313 0 : switch(tx)
4314 : {
4315 0 : case t_POL: return mkveccopy(x);
4316 : case t_SER:
4317 0 : av = avma;
4318 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
4319 : case t_RFRAC:
4320 0 : av = avma;
4321 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
4322 : }
4323 0 : pari_err_TYPE("gboundcf",x);
4324 : return NULL; /* LCOV_EXCL_LINE */
4325 : }
4326 :
4327 : static GEN
4328 14 : sfcont2(GEN b, GEN x, long k)
4329 : {
4330 14 : pari_sp av = avma;
4331 14 : long lb = lg(b), tx = typ(x), i;
4332 : GEN y,p1;
4333 :
4334 14 : if (k)
4335 : {
4336 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
4337 0 : lb = k+1;
4338 : }
4339 7 : y = cgetg(lb,t_VEC);
4340 7 : if (lb==1) return y;
4341 7 : if (is_scalar_t(tx))
4342 : {
4343 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
4344 : }
4345 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
4346 :
4347 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
4348 7 : for (i = 1;;)
4349 : {
4350 35 : if (tx == t_REAL)
4351 : {
4352 35 : long e = expo(x);
4353 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
4354 35 : gel(y,i) = floorr(x);
4355 35 : p1 = subri(x, gel(y,i));
4356 : }
4357 : else
4358 : {
4359 0 : gel(y,i) = gfloor(x);
4360 0 : p1 = gsub(x, gel(y,i));
4361 : }
4362 35 : if (++i >= lb) break;
4363 28 : if (gequal0(p1)) break;
4364 28 : x = gdiv(gel(b,i),p1);
4365 28 : }
4366 7 : setlg(y,i);
4367 7 : return gerepilecopy(av,y);
4368 : }
4369 :
4370 :
4371 : GEN
4372 0 : gcf(GEN x) { return gboundcf(x,0); }
4373 : GEN
4374 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
4375 : GEN
4376 49 : contfrac0(GEN x, GEN b, long nmax)
4377 : {
4378 : long tb;
4379 :
4380 49 : if (!b) return gboundcf(x,nmax);
4381 28 : tb = typ(b);
4382 28 : if (tb == t_INT) return gboundcf(x,itos(b));
4383 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
4384 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
4385 14 : return sfcont2(b,x,nmax);
4386 : }
4387 :
4388 : GEN
4389 126 : contfracpnqn(GEN x, long n)
4390 : {
4391 126 : pari_sp av = avma;
4392 126 : long i, lx = lg(x);
4393 : GEN M,A,B, p0,p1, q0,q1;
4394 :
4395 126 : if (lx == 1)
4396 : {
4397 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
4398 21 : if (n >= 0) return cgetg(1,t_MAT);
4399 7 : return matid(2);
4400 : }
4401 98 : switch(typ(x))
4402 : {
4403 56 : case t_VEC: case t_COL: A = x; B = NULL; break;
4404 : case t_MAT:
4405 42 : switch(lgcols(x))
4406 : {
4407 0 : case 2: A = row(x,1); B = NULL; break;
4408 35 : case 3: A = row(x,2); B = row(x,1); break;
4409 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
4410 : return NULL; /*LCOV_EXCL_LINE*/
4411 : }
4412 35 : break;
4413 0 : default: pari_err_TYPE("pnqn",x);
4414 : return NULL; /*LCOV_EXCL_LINE*/
4415 : }
4416 91 : p1 = gel(A,1);
4417 91 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
4418 91 : if (n >= 0)
4419 : {
4420 56 : lx = minss(lx, n+2);
4421 56 : if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
4422 : }
4423 35 : else if (lx == 2)
4424 7 : return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
4425 : /* lx >= 3 */
4426 56 : p0 = gen_1;
4427 56 : q0 = gen_0; /* p[-1], q[-1] */
4428 56 : M = cgetg(lx, t_MAT);
4429 56 : gel(M,1) = mkcol2(p1,q1);
4430 217 : for (i=2; i<lx; i++)
4431 : {
4432 161 : GEN a = gel(A,i), p2,q2;
4433 161 : if (B) {
4434 84 : GEN b = gel(B,i);
4435 84 : p0 = gmul(b,p0);
4436 84 : q0 = gmul(b,q0);
4437 : }
4438 161 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
4439 161 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
4440 161 : gel(M,i) = mkcol2(p1,q1);
4441 : }
4442 56 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
4443 56 : return gerepilecopy(av, M);
4444 : }
4445 : GEN
4446 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
4447 : /* x = [a0, ..., an] from gboundcf, n >= 0;
4448 : * return [[p0, ..., pn], [q0,...,qn]] */
4449 : GEN
4450 352863 : ZV_allpnqn(GEN x)
4451 : {
4452 352863 : long i, lx = lg(x);
4453 352863 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
4454 :
4455 352863 : gel(v,1) = P = cgetg(lx, t_VEC);
4456 352863 : gel(v,2) = Q = cgetg(lx, t_VEC);
4457 352863 : p0 = gen_1; q0 = gen_0;
4458 352863 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
4459 1205029 : for (i=2; i<lx; i++)
4460 : {
4461 852166 : GEN a = gel(x,i);
4462 852166 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
4463 852166 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
4464 : }
4465 352863 : return v;
4466 : }
4467 :
4468 : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
4469 : static GEN
4470 35 : mod_to_frac(GEN x, GEN N, GEN B)
4471 : {
4472 : GEN a, b, A;
4473 35 : if (B)
4474 : {
4475 21 : A = divii(shifti(N, -1), B);
4476 : /* denominator bound useless, don't use it */
4477 21 : if (cmpii(A, B) < 0) B = NULL;
4478 : }
4479 35 : if (!B)
4480 : {
4481 14 : A = sqrti(shifti(N, -1));
4482 14 : B = A;
4483 : }
4484 35 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
4485 21 : return equali1(b)? a: mkfrac(a,b);
4486 : }
4487 :
4488 : static GEN
4489 63 : mod_to_rfrac(GEN x, GEN N, long B)
4490 : {
4491 : GEN a, b;
4492 63 : long A, d = degpol(N);
4493 63 : if (B >= 0)
4494 : {
4495 21 : A = d-1 - B;
4496 : /* denominator bound useless, don't use it */
4497 21 : if (A < B) B = -1;
4498 : }
4499 63 : if (B < 0)
4500 : {
4501 56 : B = d >> 1;
4502 56 : A = odd(d)? B : B-1;
4503 : }
4504 63 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
4505 63 : if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
4506 56 : return gdiv(a,b);
4507 : }
4508 :
4509 : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
4510 : * of the continued fraction of x with b <= k maximal */
4511 : static GEN
4512 0 : bestappr_frac(GEN x, GEN k)
4513 : {
4514 : pari_sp av;
4515 : GEN p0, p1, p, q0, q1, q, a, y;
4516 :
4517 0 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
4518 0 : av = avma; y = x;
4519 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
4520 0 : q1 = gen_0; q0 = gen_1;
4521 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
4522 : for(;;)
4523 : {
4524 0 : x = ginv(x); /* > 1 */
4525 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
4526 0 : if (cmpii(a,k) > 0)
4527 : { /* next partial quotient will overflow limits */
4528 : GEN n, d;
4529 0 : a = divii(subii(k, q1), q0);
4530 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4531 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4532 : /* compare |y-p0/q0|, |y-p1/q1| */
4533 0 : n = gel(y,1);
4534 0 : d = gel(y,2);
4535 0 : if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
4536 : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
4537 0 : { p1 = p0; q1 = q0; }
4538 0 : break;
4539 : }
4540 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4541 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4542 :
4543 0 : if (cmpii(q0,k) > 0) break;
4544 0 : x = gsub(x,a); /* 0 <= x < 1 */
4545 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
4546 :
4547 0 : }
4548 0 : return gerepileupto(av, gdiv(p1,q1));
4549 : }
4550 : /* bestappr(t_REAL != 0), to maximal accuracy */
4551 : static GEN
4552 287 : bestappr_real_max(GEN x)
4553 : {
4554 287 : pari_sp av = avma;
4555 : GEN p0, p1, p, q0, q1, q, a;
4556 : long e;
4557 287 : p1 = gen_1; a = p0 = floorr(x); q1 = gen_0; q0 = gen_1;
4558 287 : x = subri(x,a); /* 0 <= x < 1 */
4559 287 : e = bit_prec(x) - expo(x);
4560 : for(;;)
4561 : {
4562 : long d;
4563 1092 : if (!signe(x) || expi(q0) > e) { p1 = p0; q1 = q0; break; }
4564 974 : x = invr(x); /* > 1 */
4565 974 : d = nbits2prec(expo(x) + 1);
4566 974 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4567 :
4568 805 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4569 805 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4570 805 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4571 805 : x = subri(x,a); /* 0 <= x < 1 */
4572 805 : }
4573 287 : return gerepileupto(av, gdiv(p1,q1));
4574 : }
4575 : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
4576 : * of the continued fraction of x with b <= k maximal */
4577 : static GEN
4578 259704 : bestappr_real(GEN x, GEN k)
4579 : {
4580 259704 : pari_sp av = avma;
4581 : GEN kr, p0, p1, p, q0, q1, q, a, y;
4582 :
4583 259704 : y = x;
4584 259704 : p1 = gen_1; a = p0 = floorr(x);
4585 259703 : q1 = gen_0; q0 = gen_1;
4586 259703 : x = subri(x,a); /* 0 <= x < 1 */
4587 259703 : if (!signe(x)) { cgiv(x); return a; }
4588 253614 : kr = itor(k, realprec(x));
4589 : for(;;)
4590 : {
4591 : long d;
4592 684028 : x = invr(x); /* > 1 */
4593 684028 : if (cmprr(x,kr) > 0)
4594 : { /* next partial quotient will overflow limits */
4595 249124 : a = divii(subii(k, q1), q0);
4596 249124 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4597 249124 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4598 : /* compare |y-p0/q0|, |y-p1/q1| */
4599 249124 : if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
4600 : mulir(q0, subri(mulir(q1,y), p1))) < 0)
4601 4978 : { p1 = p0; q1 = q0; }
4602 249124 : break;
4603 : }
4604 434904 : d = nbits2prec(expo(x) + 1);
4605 434904 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4606 :
4607 434898 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4608 434898 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4609 434898 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4610 :
4611 434898 : if (cmpii(q0,k) > 0) break;
4612 430731 : x = subri(x,a); /* 0 <= x < 1 */
4613 430731 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
4614 430414 : }
4615 253614 : return gerepileupto(av, gdiv(p1,q1));
4616 : }
4617 :
4618 : /* k t_INT or NULL */
4619 : static GEN
4620 413636 : bestappr_Q(GEN x, GEN k)
4621 : {
4622 413636 : long lx, tx = typ(x), i;
4623 : GEN a, y;
4624 :
4625 413636 : switch(tx)
4626 : {
4627 77 : case t_INT: return icopy(x);
4628 0 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
4629 : case t_REAL:
4630 292389 : if (!signe(x)) return gen_0;
4631 259991 : return k? bestappr_real(x, k): bestappr_real_max(x);
4632 :
4633 : case t_INTMOD: {
4634 21 : pari_sp av = avma;
4635 21 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
4636 14 : return gerepilecopy(av, a);
4637 : }
4638 : case t_PADIC: {
4639 14 : pari_sp av = avma;
4640 14 : long v = valp(x);
4641 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
4642 7 : if (v) a = gmul(a, powis(gel(x,2), v));
4643 7 : return gerepilecopy(av, a);
4644 : }
4645 :
4646 : case t_COMPLEX: {
4647 7 : pari_sp av = avma;
4648 7 : y = cgetg(3, t_COMPLEX);
4649 7 : gel(y,2) = bestappr(gel(x,2), k);
4650 7 : gel(y,1) = bestappr(gel(x,1), k);
4651 7 : if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
4652 0 : return y;
4653 : }
4654 : case t_SER:
4655 0 : if (ser_isexactzero(x)) return gcopy(x);
4656 : /* fall through */
4657 : case t_POLMOD: case t_POL: case t_RFRAC:
4658 : case t_VEC: case t_COL: case t_MAT:
4659 121128 : y = cgetg_copy(x, &lx);
4660 121128 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
4661 519717 : for (; i<lx; i++)
4662 : {
4663 398591 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
4664 398589 : gel(y,i) = a;
4665 : }
4666 121126 : if (tx == t_POL) return normalizepol(y);
4667 121112 : if (tx == t_SER) return normalize(y);
4668 121112 : return y;
4669 : }
4670 0 : pari_err_TYPE("bestappr_Q",x);
4671 : return NULL; /* LCOV_EXCL_LINE */
4672 : }
4673 :
4674 : static GEN
4675 56 : bestappr_ser(GEN x, long B)
4676 : {
4677 56 : long dN, v = valp(x), lx = lg(x);
4678 : GEN t;
4679 56 : x = normalizepol(ser2pol_i(x, lx));
4680 56 : dN = lx-2;
4681 56 : if (v > 0)
4682 : {
4683 14 : x = RgX_shift_shallow(x, v);
4684 14 : dN += v;
4685 : }
4686 42 : else if (v < 0)
4687 : {
4688 7 : if (B >= 0) B = maxss(B+v, 0);
4689 : }
4690 56 : t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
4691 56 : if (!t) return NULL;
4692 49 : if (v < 0)
4693 : {
4694 : GEN a, b;
4695 : long vx;
4696 7 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
4697 : /* t_RFRAC */
4698 7 : vx = varn(x);
4699 7 : a = gel(t,1);
4700 7 : b = gel(t,2);
4701 7 : v -= RgX_valrem(b, &b);
4702 7 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
4703 7 : if (v < 0) b = RgX_shift(b, -v);
4704 0 : else if (v > 0) {
4705 0 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
4706 0 : a = RgX_shift(a, v);
4707 : }
4708 7 : t = mkrfraccopy(a, b);
4709 : }
4710 49 : return t;
4711 : }
4712 : static GEN bestappr_RgX(GEN x, long B);
4713 : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
4714 : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
4715 : static GEN
4716 70 : bestappr_RgX(GEN x, long B)
4717 : {
4718 70 : long i, lx, tx = typ(x);
4719 : GEN y, t;
4720 70 : switch(tx)
4721 : {
4722 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
4723 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
4724 0 : return gcopy(x);
4725 :
4726 : case t_RFRAC: {
4727 14 : pari_sp av = avma;
4728 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
4729 7 : x = rfrac_to_ser(x, 2*B+1);
4730 7 : t = bestappr_ser(x, B); if (!t) return NULL;
4731 7 : return gerepileupto(av, t);
4732 : }
4733 : case t_POLMOD: {
4734 7 : pari_sp av = avma;
4735 7 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
4736 7 : return gerepileupto(av, t);
4737 : }
4738 : case t_SER: {
4739 49 : pari_sp av = avma;
4740 49 : t = bestappr_ser(x, B); if (!t) return NULL;
4741 42 : return gerepileupto(av, t);
4742 : }
4743 :
4744 : case t_VEC: case t_COL: case t_MAT:
4745 0 : y = cgetg_copy(x, &lx);
4746 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
4747 0 : for (; i<lx; i++)
4748 : {
4749 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
4750 0 : gel(y,i) = t;
4751 : }
4752 0 : return y;
4753 : }
4754 0 : pari_err_TYPE("bestappr_RgX",x);
4755 : return NULL; /* LCOV_EXCL_LINE */
4756 : }
4757 :
4758 : /* allow k = NULL: maximal accuracy */
4759 : GEN
4760 15045 : bestappr(GEN x, GEN k)
4761 : {
4762 15045 : pari_sp av = avma;
4763 15045 : if (k) { /* replace by floor(k) */
4764 14779 : switch(typ(k))
4765 : {
4766 : case t_INT:
4767 168 : break;
4768 : case t_REAL: case t_FRAC:
4769 14611 : k = floor_safe(k); /* left on stack for efficiency */
4770 14611 : if (!signe(k)) k = gen_1;
4771 14611 : break;
4772 : default:
4773 0 : pari_err_TYPE("bestappr [bound type]", k);
4774 0 : break;
4775 : }
4776 : }
4777 15045 : x = bestappr_Q(x, k);
4778 15044 : if (!x) { avma = av; return cgetg(1,t_VEC); }
4779 15030 : return x;
4780 : }
4781 : GEN
4782 70 : bestapprPade(GEN x, long B)
4783 : {
4784 70 : pari_sp av = avma;
4785 70 : GEN t = bestappr_RgX(x, B);
4786 70 : if (!t) { avma = av; return cgetg(1,t_VEC); }
4787 63 : return t;
4788 : }
4789 :
4790 : /***********************************************************************/
4791 : /** **/
4792 : /** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
4793 : /** **/
4794 : /***********************************************************************/
4795 :
4796 : static GEN
4797 14 : get_quad(GEN f, GEN pol, long r)
4798 : {
4799 14 : GEN p1 = gcoeff(f,1,2), q1 = gcoeff(f,2,2);
4800 14 : return mkquad(pol, r? subii(p1,q1): p1, q1);
4801 : }
4802 :
4803 : /* replace f by f * [a,1; 1,0] */
4804 : static void
4805 14 : update_f(GEN f, GEN a)
4806 : {
4807 : GEN p1;
4808 14 : p1 = gcoeff(f,1,1);
4809 14 : gcoeff(f,1,1) = addii(mulii(a,p1), gcoeff(f,1,2));
4810 14 : gcoeff(f,1,2) = p1;
4811 :
4812 14 : p1 = gcoeff(f,2,1);
4813 14 : gcoeff(f,2,1) = addii(mulii(a,p1), gcoeff(f,2,2));
4814 14 : gcoeff(f,2,2) = p1;
4815 14 : }
4816 :
4817 : GEN
4818 7 : quadunit(GEN x)
4819 : {
4820 7 : pari_sp av = avma, av2;
4821 : GEN pol, y, a, u, v, sqd, f;
4822 : long r;
4823 :
4824 7 : check_quaddisc_real(x, &r, "quadunit");
4825 7 : pol = quadpoly(x);
4826 7 : sqd = sqrti(x); av2 = avma;
4827 7 : a = shifti(addui(r,sqd),-1);
4828 7 : f = mkmat2(mkcol2(a, gen_1), mkcol2(gen_1, gen_0)); /* [a,0; 1,0] */
4829 7 : u = stoi(r); v = gen_2;
4830 : for(;;)
4831 : {
4832 : GEN u1, v1;
4833 14 : u1 = subii(mulii(a,v),u);
4834 14 : v1 = divii(subii(x,sqri(u1)),v);
4835 14 : if ( equalii(v,v1) ) {
4836 7 : y = get_quad(f,pol,r);
4837 7 : update_f(f,a);
4838 7 : y = gdiv(get_quad(f,pol,r), gconj(y));
4839 7 : break;
4840 : }
4841 7 : a = divii(addii(sqd,u1), v1);
4842 7 : if ( equalii(u,u1) ) {
4843 0 : y = get_quad(f,pol,r);
4844 0 : y = gdiv(y, gconj(y));
4845 0 : break;
4846 : }
4847 7 : update_f(f,a);
4848 7 : u = u1; v = v1;
4849 7 : if (gc_needed(av2,2))
4850 : {
4851 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit");
4852 0 : gerepileall(av2,4, &a,&f,&u,&v);
4853 : }
4854 7 : }
4855 7 : if (signe(gel(y,3)) < 0) y = gneg(y);
4856 7 : return gerepileupto(av, y);
4857 : }
4858 :
4859 : GEN
4860 7 : quadunit0(GEN x, long v)
4861 : {
4862 7 : GEN y = quadunit(x);
4863 7 : if (v==-1) v = fetch_user_var("w");
4864 7 : setvarn(gel(y,1), v);
4865 7 : return y;
4866 : }
4867 :
4868 : GEN
4869 21 : quadregulator(GEN x, long prec)
4870 : {
4871 21 : pari_sp av = avma, av2;
4872 : GEN R, rsqd, u, v, sqd;
4873 : long r, Rexpo;
4874 :
4875 21 : check_quaddisc_real(x, &r, "quadregulator");
4876 21 : sqd = sqrti(x);
4877 21 : rsqd = gsqrt(x,prec);
4878 21 : Rexpo = 0; R = real2n(1, prec); /* = 2 */
4879 21 : av2 = avma;
4880 21 : u = stoi(r); v = gen_2;
4881 : for(;;)
4882 : {
4883 70 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
4884 70 : GEN v1 = divii(subii(x,sqri(u1)),v);
4885 70 : if (equalii(v,v1))
4886 : {
4887 7 : R = sqrr(R); shiftr_inplace(R, -1);
4888 7 : R = mulrr(R, divri(addir(u1,rsqd),v));
4889 7 : break;
4890 : }
4891 63 : if (equalii(u,u1))
4892 : {
4893 14 : R = sqrr(R); shiftr_inplace(R, -1);
4894 14 : break;
4895 : }
4896 49 : R = mulrr(R, divri(addir(u1,rsqd),v));
4897 49 : Rexpo += expo(R); setexpo(R,0);
4898 49 : u = u1; v = v1;
4899 49 : if (Rexpo & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
4900 49 : if (gc_needed(av2,2))
4901 : {
4902 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
4903 0 : gerepileall(av2,3, &R,&u,&v);
4904 : }
4905 49 : }
4906 21 : R = logr_abs(divri(R,v));
4907 21 : if (Rexpo)
4908 : {
4909 21 : GEN t = mulsr(Rexpo, mplog2(prec));
4910 21 : shiftr_inplace(t, 1);
4911 21 : R = addrr(R,t);
4912 : }
4913 21 : return gerepileuptoleaf(av, R);
4914 : }
4915 :
4916 : /*************************************************************************/
4917 : /** **/
4918 : /** CLASS NUMBER **/
4919 : /** **/
4920 : /*************************************************************************/
4921 :
4922 : int
4923 12264017 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
4924 :
4925 17358976 : static GEN qfi_pow(void *E, GEN f, GEN n)
4926 17358976 : { return E? nupow(f,n,(GEN)E): powgi(f,n); }
4927 21814307 : static GEN qfi_comp(void *E, GEN f, GEN g)
4928 21814307 : { return E? nucomp(f,g,(GEN)E): qficomp(f,g); }
4929 : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
4930 : gidentical,qfb_equal1,NULL};
4931 :
4932 : GEN
4933 2789255 : qfi_order(GEN q, GEN o)
4934 2789255 : { return gen_order(q, o, NULL, &qfi_group); }
4935 :
4936 : GEN
4937 0 : qfi_log(GEN a, GEN g, GEN o)
4938 0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
4939 :
4940 : GEN
4941 596537 : qfi_Shanks(GEN a, GEN g, long n)
4942 : {
4943 596537 : pari_sp av = avma;
4944 : GEN T, X;
4945 : long rt_n, c;
4946 :
4947 596537 : a = redimag(a);
4948 596537 : g = redimag(g);
4949 :
4950 596537 : rt_n = sqrt((double)n);
4951 596537 : c = n / rt_n;
4952 596537 : c = (c * rt_n < n + 1) ? c + 1 : c;
4953 :
4954 596537 : T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
4955 596537 : X = gen_Shanks(T, a, c, NULL, &qfi_group);
4956 :
4957 596537 : if (!X) { avma = av; return X; }
4958 316928 : return gerepileuptoint(av, X);
4959 : }
4960 :
4961 : GEN
4962 140 : qfbclassno0(GEN x,long flag)
4963 : {
4964 140 : switch(flag)
4965 : {
4966 126 : case 0: return map_proto_G(classno,x);
4967 14 : case 1: return map_proto_G(classno2,x);
4968 0 : default: pari_err_FLAG("qfbclassno");
4969 : }
4970 : return NULL; /* LCOV_EXCL_LINE */
4971 : }
4972 :
4973 : /* f^h = 1, return order(f). Set *pfao to its factorization */
4974 : static GEN
4975 2635497 : find_order(void *E, GEN f, GEN h, GEN *pfao)
4976 : {
4977 2635497 : GEN v = gen_factored_order(f, h, E, &qfi_group);
4978 2635497 : *pfao = gel(v,2); return gel(v,1);
4979 : }
4980 :
4981 : static int
4982 6327 : ok_q(GEN q, GEN h, GEN d2, long r2)
4983 : {
4984 6327 : if (d2)
4985 : {
4986 7 : if (r2 <= 2 && !mpodd(q)) return 0;
4987 7 : return is_pm1(Z_ppo(q,d2));
4988 : }
4989 : else
4990 : {
4991 6320 : if (r2 <= 1 && !mpodd(q)) return 0;
4992 6320 : return is_pm1(Z_ppo(q,h));
4993 : }
4994 : }
4995 :
4996 : /* a,b given by their factorizations. Return factorization of lcm(a,b).
4997 : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
4998 : static GEN
4999 344004 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
5000 : {
5001 344004 : GEN P = ZV_union_shallow(gel(Fa,1), gel(Fb,1));
5002 344004 : GEN A = gen_1, B = gen_1;
5003 344004 : long i, l = lg(P);
5004 344004 : GEN E = cgetg(l, t_COL);
5005 1015957 : for (i=1; i<l; i++)
5006 : {
5007 671953 : GEN p = gel(P,i);
5008 671953 : long va = Z_pval(a,p);
5009 671953 : long vb = Z_pval(b,p);
5010 671953 : if (va < vb)
5011 : {
5012 345970 : B = mulii(B,powiu(p,vb));
5013 345970 : gel(E,i) = utoi(vb);
5014 : }
5015 : else
5016 : {
5017 325983 : A = mulii(A,powiu(p,va));
5018 325983 : gel(E,i) = utoi(va);
5019 : }
5020 : }
5021 344004 : *pA = A;
5022 344004 : *pB = B; return mkmat2(P,E);
5023 : }
5024 :
5025 : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
5026 : static void
5027 344004 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
5028 : {
5029 344004 : GEN A,B, g1 = *pg1, d1 = *pd1;
5030 344004 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
5031 344004 : *pg1 = gmul(powgi(g1, diviiexact(d1,A)), powgi(f, diviiexact(o,B)));
5032 344004 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
5033 344004 : }
5034 :
5035 : /* Write x = Df^2, where D = fundamental discriminant,
5036 : * P^E = factorisation of conductor f, with E[i] >= 0 */
5037 : static void
5038 1955208 : corediscfact(GEN x, long xmod4, GEN *ptD, GEN *ptP, GEN *ptE)
5039 : {
5040 1955208 : long s = signe(x), l, i;
5041 1955208 : GEN fa = absZ_factor(x);
5042 1955208 : GEN d, P = gel(fa,1), E = gtovecsmall(gel(fa,2));
5043 :
5044 1955208 : l = lg(P); d = gen_1;
5045 5088425 : for (i=1; i<l; i++)
5046 : {
5047 3133217 : if (E[i] & 1) d = mulii(d, gel(P,i));
5048 3133217 : E[i] >>= 1;
5049 : }
5050 1955208 : if (!xmod4 && mod4(d) != ((s < 0)? 3: 1)) { d = shifti(d,2); E[1]--; }
5051 1955208 : *ptD = (s < 0)? negi(d): d;
5052 1955208 : *ptP = P;
5053 1955208 : *ptE = E;
5054 1955208 : }
5055 :
5056 : static GEN
5057 1947335 : conductor_part(GEN x, long xmod4, GEN *ptD, GEN *ptreg)
5058 : {
5059 1947335 : long l, i, s = signe(x);
5060 : GEN E, H, D, P, reg;
5061 :
5062 1947335 : corediscfact(x, xmod4, &D, &P, &E);
5063 1947335 : H = gen_1; l = lg(P);
5064 : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
5065 5054526 : for (i=1; i<l; i++)
5066 : {
5067 3107191 : long e = E[i];
5068 3107191 : if (e)
5069 : {
5070 7 : GEN p = gel(P,i);
5071 7 : H = mulii(H, subis(p, kronecker(D,p)));
5072 7 : if (e >= 2) H = mulii(H, powiu(p,e-1));
5073 : }
5074 : }
5075 :
5076 : /* divide by [ O_K^* : O^* ] */
5077 1947335 : if (s < 0)
5078 : {
5079 1947321 : reg = NULL;
5080 1947321 : switch(itou_or_0(D))
5081 : {
5082 0 : case 4: H = shifti(H,-1); break;
5083 0 : case 3: H = divis(H,3); break;
5084 : }
5085 : } else {
5086 14 : reg = quadregulator(D,DEFAULTPREC);
5087 14 : if (!equalii(x,D))
5088 0 : H = divii(H, roundr(divrr(quadregulator(x,DEFAULTPREC), reg)));
5089 : }
5090 1947335 : if (ptreg) *ptreg = reg;
5091 1947335 : *ptD = D; return H;
5092 : }
5093 :
5094 : static long
5095 1947314 : two_rank(GEN x)
5096 : {
5097 1947314 : GEN p = gel(absZ_factor(x),1);
5098 1947314 : long l = lg(p)-1;
5099 : #if 0 /* positive disc not needed */
5100 : if (signe(x) > 0)
5101 : {
5102 : long i;
5103 : for (i=1; i<=l; i++)
5104 : if (mod4(gel(p,i)) == 3) { l--; break; }
5105 : }
5106 : #endif
5107 1947314 : return l-1;
5108 : }
5109 :
5110 : static GEN
5111 36996962 : sqr_primeform(GEN x, ulong p) { return redimag(qfisqr(primeform_u(x, p))); }
5112 : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
5113 : static GEN
5114 1947314 : get_forms(GEN D, GEN *pL)
5115 : {
5116 1947314 : const long MAXFORM = 20;
5117 1947314 : GEN L, sqrtD = gsqrt(absi(D),DEFAULTPREC), forms = vectrunc_init(MAXFORM+1);
5118 1947314 : long s, nforms = 0;
5119 : ulong p;
5120 : forprime_t S;
5121 1947314 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
5122 1947314 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
5123 1947314 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
5124 1947314 : if (s < 10) s = 200;
5125 1804118 : else if (s < 20) s = 1000;
5126 1477 : else if (s < 5000) s = 5000;
5127 1947314 : u_forprime_init(&S, 2, s);
5128 327318752 : while ( (p = u_forprime_next(&S)) )
5129 : {
5130 323424124 : long d, k = kroiu(D,p);
5131 : pari_sp av2;
5132 323424124 : if (!k) continue;
5133 321311195 : if (k > 0)
5134 : {
5135 161155757 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
5136 161155757 : d = p-1;
5137 : }
5138 : else
5139 160155438 : d = p+1;
5140 321311195 : av2 = avma; affrr(divru(mulur(p,L),d), L); avma = av2;
5141 : }
5142 1947314 : *pL = L; return forms;
5143 : }
5144 :
5145 : /* h ~ #G, return o = order of f, set fao = its factorization */
5146 : static GEN
5147 1947363 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
5148 : {
5149 1947363 : long s = minss(itos(sqrti(h)), 10000);
5150 1947363 : GEN T = gen_Shanks_init(f, s, E, &qfi_group);
5151 1947363 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
5152 1947363 : return find_order(E, f, addiu(v,1), pfao);
5153 : }
5154 :
5155 : /* if g = 1 in G/<f> ? */
5156 : static int
5157 518 : equal1(void *E, GEN T, ulong N, GEN g)
5158 518 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
5159 :
5160 : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
5161 : * FIXME: should be gen_order, but equal1 has the wrong prototype */
5162 : static GEN
5163 112 : relative_order(void *E, GEN a, GEN o, ulong N, GEN T)
5164 : {
5165 112 : pari_sp av = avma;
5166 : long i, l;
5167 : GEN m;
5168 :
5169 112 : m = get_arith_ZZM(o);
5170 112 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
5171 112 : o = gel(m,1);
5172 112 : m = gel(m,2); l = lgcols(m);
5173 322 : for (i = l-1; i; i--)
5174 : {
5175 210 : GEN t, y, p = gcoeff(m,i,1);
5176 210 : long j, e = itos(gcoeff(m,i,2));
5177 210 : if (l == 2) {
5178 35 : t = gen_1;
5179 35 : y = a;
5180 : } else {
5181 175 : t = diviiexact(o, powiu(p,e));
5182 175 : y = powgi(a, t);
5183 : }
5184 210 : if (equal1(E, T,N,y)) o = t;
5185 : else {
5186 126 : for (j = 1; j < e; j++)
5187 : {
5188 28 : y = powgi(y, p);
5189 28 : if (equal1(E, T,N,y)) break;
5190 : }
5191 119 : if (j < e) {
5192 21 : if (j > 1) p = powiu(p, j);
5193 21 : o = mulii(t, p);
5194 : }
5195 : }
5196 : }
5197 112 : return gerepilecopy(av, o);
5198 : }
5199 :
5200 : /* h(x) for x<0 using Baby Step/Giant Step.
5201 : * Assumes G is not too far from being cyclic.
5202 : *
5203 : * Compute G^2 instead of G so as to kill most of the non-cyclicity */
5204 : GEN
5205 1949557 : classno(GEN x)
5206 : {
5207 1949557 : pari_sp av = avma;
5208 : long r2, k, s, i, l;
5209 : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
5210 : void *E;
5211 :
5212 1949557 : if (signe(x) >= 0) return classno2(x);
5213 :
5214 1949550 : check_quaddisc(x, &s, &k, "classno");
5215 1949550 : if (abscmpiu(x,12) <= 0) return gen_1;
5216 :
5217 1947314 : Hf = conductor_part(x, k, &D, NULL);
5218 1947314 : if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
5219 1947314 : forms = get_forms(D, &L);
5220 1947314 : r2 = two_rank(D);
5221 1947314 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
5222 :
5223 1947314 : l = lg(forms);
5224 1947314 : order_bound = const_vec(l-1, NULL);
5225 1947314 : E = expi(D) > 60? (void*)sqrtnint(shifti(absi(D),-2),4): NULL;
5226 1947314 : g1 = gel(forms,1);
5227 1947314 : gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
5228 1947314 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
5229 1947314 : d2 = NULL; /* not computed yet */
5230 1947314 : if (is_pm1(q)) goto END;
5231 482765 : for (i=2; i < l; i++)
5232 : {
5233 476375 : GEN o, fao, a, F, fd, f = gel(forms,i);
5234 476375 : fd = powgi(f, d1); if (is_pm1(gel(fd,1))) continue;
5235 344004 : F = powgi(fd, q);
5236 344004 : a = gel(F,1);
5237 344004 : o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
5238 : /* f^(d1 q) = 1 */
5239 344004 : fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
5240 344004 : o = find_order(E, f, fao, &fao);
5241 344004 : gel(order_bound,i) = o;
5242 : /* o = order of f, fao = factor(o) */
5243 344004 : update_g1(&g1,&d1,&fad1, f,o,fao);
5244 344004 : q = diviiround(hin, d1);
5245 344004 : if (is_pm1(q)) goto END;
5246 : }
5247 : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
5248 6390 : if (expi(q) > 3)
5249 : { /* q large: compute d2, 2nd elt divisor */
5250 70 : ulong N, n = 2*itou(sqrti(d1));
5251 70 : GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
5252 70 : d2 = gen_1;
5253 70 : N = itou( gceil(gdivgs(d1,n)) ); /* order(g1) <= n*N */
5254 287 : for (i = 1; i < l; i++)
5255 : {
5256 280 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
5257 280 : if (!B) B = find_order(E, f, fad1, /*junk*/&d);
5258 280 : f = powgi(f,d2);
5259 280 : if (equal1(E, T,N,f)) continue;
5260 112 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
5261 : /* f^B = 1 */
5262 112 : d = relative_order(E, f, B, N,T);
5263 112 : d2= mulii(d,d2);
5264 112 : D = mulii(d1,d2);
5265 112 : q = diviiround(hin,D);
5266 112 : if (is_pm1(q)) { d1 = D; goto END; }
5267 : }
5268 : /* very probably, d2 is the 2nd elementary divisor */
5269 7 : d1 = D; /* product of first two elt divisors */
5270 : }
5271 : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
5272 : * 2-rank */
5273 6327 : if (!ok_q(q,d1,d2,r2))
5274 : {
5275 0 : GEN q0 = q;
5276 : long d;
5277 0 : if (cmpii(mulii(q,d1), hin) < 0)
5278 : { /* try q = q0+1,-1,+2,-2 */
5279 0 : d = 1;
5280 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
5281 : }
5282 : else
5283 : { /* q0-1,+1,-2,+2 */
5284 0 : d = -1;
5285 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
5286 : }
5287 : }
5288 6327 : d1 = mulii(d1,q);
5289 :
5290 : END:
5291 1947314 : return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
5292 : }
5293 :
5294 : GEN
5295 0 : quadclassno(GEN x)
5296 : {
5297 0 : pari_sp av = avma;
5298 : GEN Hf, D;
5299 : long s, r;
5300 0 : check_quaddisc(x, &s, &r, "quadclassno");
5301 0 : if (s < 0 && cmpiu(x,12) <= 0) return gen_1;
5302 0 : Hf = conductor_part(x, r, &D, NULL);
5303 0 : return gerepileuptoint(av, mulii(Hf, gel(quadclassunit0(D,0,NULL,0),1)));
5304 : }
5305 :
5306 : /* use Euler products */
5307 : GEN
5308 21 : classno2(GEN x)
5309 : {
5310 21 : pari_sp av = avma;
5311 21 : const long prec = DEFAULTPREC;
5312 : long n, i, r, s;
5313 : GEN p1, p2, S, p4, p5, p7, Hf, Pi, reg, logd, d, dr, D, half;
5314 :
5315 21 : check_quaddisc(x, &s, &r, "classno2");
5316 21 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
5317 :
5318 21 : Hf = conductor_part(x, r, &D, ®);
5319 21 : if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
5320 :
5321 21 : Pi = mppi(prec);
5322 21 : d = absi(D); dr = itor(d, prec);
5323 21 : logd = logr_abs(dr);
5324 21 : p1 = sqrtr(divrr(mulir(d,logd), gmul2n(Pi,1)));
5325 21 : if (s > 0)
5326 : {
5327 14 : GEN invlogd = invr(logd);
5328 14 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
5329 14 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
5330 : }
5331 21 : n = itos_or_0( mptrunc(p1) );
5332 21 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
5333 :
5334 21 : p4 = divri(Pi,d);
5335 21 : p7 = invr(sqrtr_abs(Pi));
5336 21 : half = real2n(-1, prec);
5337 21 : if (s > 0)
5338 : { /* i = 1, shortcut */
5339 14 : p1 = sqrtr_abs(dr);
5340 14 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5341 14 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
5342 546 : for (i=2; i<=n; i++)
5343 : {
5344 532 : long k = kroiu(D,i); if (!k) continue;
5345 434 : p2 = mulir(sqru(i), p4);
5346 434 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5347 434 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
5348 434 : S = (k>0)? addrr(S,p5): subrr(S,p5);
5349 : }
5350 14 : S = shiftr(divrr(S,reg),-1);
5351 : }
5352 : else
5353 : { /* i = 1, shortcut */
5354 7 : p1 = gdiv(sqrtr_abs(dr), Pi);
5355 7 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5356 7 : S = addrr(p5, divrr(p1, mpexp(p4)));
5357 952 : for (i=2; i<=n; i++)
5358 : {
5359 945 : long k = kroiu(D,i); if (!k) continue;
5360 945 : p2 = mulir(sqru(i), p4);
5361 945 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5362 945 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
5363 945 : S = (k>0)? addrr(S,p5): subrr(S,p5);
5364 : }
5365 : }
5366 21 : return gerepileuptoint(av, mulii(Hf, roundr(S)));
5367 : }
5368 :
5369 : /* 1 + q + ... + q^v, v > 0 */
5370 : static GEN
5371 120 : geomsumu(ulong q, long v)
5372 : {
5373 120 : GEN u = utoipos(1+q);
5374 120 : for (; v > 1; v--) u = addui(1, mului(q, u));
5375 120 : return u;
5376 : }
5377 : static GEN
5378 120 : geomsum(GEN q, long v)
5379 : {
5380 : GEN u;
5381 120 : if (lgefint(q) == 3) return geomsumu(q[2], v);
5382 0 : u = addiu(q,1);
5383 0 : for (; v > 1; v--) u = addui(1, mulii(q, u));
5384 0 : return u;
5385 : }
5386 :
5387 : static GEN
5388 7873 : hclassno6_large(GEN x)
5389 : {
5390 : long i, l, s, xmod4;
5391 : GEN Q, H, D, P, E;
5392 :
5393 7873 : x = negi(x);
5394 7873 : check_quaddisc(x, &s, &xmod4, "hclassno");
5395 7873 : corediscfact(x, xmod4, &D, &P, &E);
5396 :
5397 7873 : Q = quadclassunit0(D, 0, NULL, 0);
5398 7873 : H = gel(Q,1); l = lg(P);
5399 :
5400 : /* H \prod_{p^e||f} (1 + (p^e-1)/(p-1))[ p - (D/p) ] */
5401 33899 : for (i=1; i<l; i++)
5402 : {
5403 26026 : long e = E[i], s;
5404 : GEN p, t;
5405 26026 : if (!e) continue;
5406 5108 : p = gel(P,i); s = kronecker(D,p);
5407 5108 : if (e == 1) t = addiu(p, 1-s);
5408 1028 : else if (s == 1) t = powiu(p,e);
5409 120 : else t = addui(1, mulii(subis(p, s), geomsum(p,e-1)));
5410 5108 : H = mulii(H,t);
5411 : }
5412 7873 : switch( itou_or_0(D) )
5413 : {
5414 0 : case 3: H = shifti(H,1);break;
5415 0 : case 4: H = muliu(H,3); break;
5416 7873 : default:H = muliu(H,6); break;
5417 : }
5418 7873 : return H;
5419 : }
5420 :
5421 : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
5422 : GEN
5423 120505 : hclassno6(GEN x)
5424 : {
5425 120505 : ulong d = itou_or_0(x);
5426 120505 : if (!d || d > 500000) return hclassno6_large(x);
5427 112632 : return utoipos(hclassno6u(d));
5428 : }
5429 :
5430 : GEN
5431 44576 : hclassno(GEN x)
5432 : {
5433 : long a, s;
5434 44576 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
5435 44576 : s = signe(x);
5436 44576 : if (s < 0) return gen_0;
5437 44576 : if (!s) return gdivgs(gen_1, -12);
5438 44576 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
5439 44576 : return gdivgs(hclassno6(x), 6);
5440 : }
5441 : /******************************************************************/
5442 : /* */
5443 : /* RAMANUJAN's TAU FUNCTION */
5444 : /* */
5445 : /******************************************************************/
5446 : /* 4|N > 0, not fundamental at 2; 6 * Hurwitz class number in level 2,
5447 : * equal to 6*(H(N)+2H(N/4)), H=qfbhclassno */
5448 : static GEN
5449 36806 : Hspec(GEN N)
5450 : {
5451 36806 : long v2 = Z_lvalrem(N, 2, &N), v2f = v2 >> 1;
5452 : GEN t;
5453 36806 : if (odd(v2)) { v2f--; N = shifti(N,3); }
5454 32613 : else if (mod4(N)!=3) { v2f--; N = shifti(N,2); }
5455 : /* N fundamental at 2, v2f = v2(f) s.t. N = f^2 D, D fundamental */
5456 36806 : t = addui(3, muliu(subiu(int2n(v2f+1), 3), 2 - kroiu(N,2)));
5457 36806 : return mulii(t, hclassno6(N));
5458 : }
5459 :
5460 : /* Ramanujan tau function for p prime */
5461 : static GEN
5462 15043 : tauprime(GEN p)
5463 : {
5464 15043 : pari_sp av = avma, av2;
5465 : GEN s, p2, p2_7, p_9, T;
5466 : ulong lim, t, tin;
5467 :
5468 15043 : if (absequaliu(p, 2)) return utoineg(24);
5469 : /* p > 2 */
5470 11480 : p2 = sqri(p);
5471 11480 : p2_7 = mului(7, p2);
5472 11480 : p_9 = mului(9, p);
5473 11480 : av2 = avma;
5474 11480 : lim = itou(sqrtint(p));
5475 11480 : tin = mod4(p) == 3? 1: 0;
5476 11480 : s = gen_0;
5477 87409 : for (t = 1; t <= lim; ++t)
5478 : {
5479 75929 : GEN h, a, t2 = sqru(t), D = shifti(subii(p, t2), 2); /* 4(p-t^2) */
5480 : /* t mod 2 != tin <=> D not fundamental at 2 */
5481 75929 : h = ((t&1UL) == tin)? hclassno6(D): Hspec(D);
5482 75929 : a = mulii(powiu(t2,3), addii(p2_7, mulii(t2, subii(shifti(t2,2), p_9))));
5483 75929 : s = addii(s, mulii(a,h));
5484 75929 : if (!(t & 255)) s = gerepileuptoint(av2, s);
5485 : }
5486 : /* 28p^3 - 28p^2 - 90p - 35 */
5487 11480 : T = subii(shifti(mulii(p2_7, subiu(p,1)), 2), addiu(mului(90,p), 35));
5488 11480 : s = shifti(diviuexact(s, 3), 6);
5489 11480 : return gerepileuptoint(av, subii(mulii(mulii(p2,p),T), addui(1, s)));
5490 : }
5491 :
5492 : /* Ramanujan tau function, return 0 for <= 0 */
5493 : GEN
5494 7168 : ramanujantau(GEN n)
5495 : {
5496 7168 : pari_sp ltop = avma;
5497 : GEN T, F, P, E;
5498 : long j, lP;
5499 :
5500 7168 : if (!(F = check_arith_all(n,"ramanujantau")))
5501 : {
5502 7147 : if (signe(n) <= 0) return gen_0;
5503 7126 : F = Z_factor(n);
5504 : }
5505 : else
5506 : {
5507 21 : P = gel(F,1);
5508 21 : if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
5509 : }
5510 :
5511 7133 : P = gel(F,1);
5512 7133 : E = gel(F,2); lP = lg(P);
5513 7133 : T = gen_1;
5514 22176 : for (j = 1; j < lP; j++)
5515 : {
5516 15043 : GEN p = gel(P,j), tp = tauprime(p), t1 = tp, t0 = gen_1;
5517 15043 : long k, e = itou(gel(E,j));
5518 20335 : for (k = 1; k < e; k++)
5519 : {
5520 5292 : GEN t2 = subii(mulii(tp, t1), mulii(powiu(p, 11), t0));
5521 5292 : t0 = t1; t1 = t2;
5522 : }
5523 15043 : T = mulii(T, t1);
5524 : }
5525 7133 : return gerepileuptoint(ltop, T);
5526 : }
|
