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