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