Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; 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 : /** ARITHMETIC FUNCTIONS **/
17 : /** (first part) **/
18 : /*********************************************************************/
19 : #include "pari.h"
20 : #include "paripriv.h"
21 :
22 : #define DEBUGLEVEL DEBUGLEVEL_arith
23 :
24 : /******************************************************************/
25 : /* GENERATOR of (Z/mZ)* */
26 : /******************************************************************/
27 : static GEN
28 1812 : remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
29 : static ulong
30 464062 : u_remove2(ulong q) { return q >> vals(q); }
31 : GEN
32 1812 : odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
33 : static GEN
34 464062 : u_odd_prime_divisors(ulong q) { return gel(factoru(u_remove2(q)), 1); }
35 : /* p odd prime, q=(p-1)/2; L0 list of (some) divisors of q = (p-1)/2 or NULL
36 : * (all prime divisors of q); return the q/l, l in L0 */
37 : static GEN
38 4909 : is_gener_expo(GEN p, GEN L0)
39 : {
40 4909 : GEN L, q = shifti(p,-1);
41 : long i, l;
42 4909 : if (L0) {
43 3134 : l = lg(L0);
44 3134 : L = cgetg(l, t_VEC);
45 : } else {
46 1775 : L0 = L = odd_prime_divisors(q);
47 1775 : l = lg(L);
48 : }
49 14199 : for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
50 4909 : return L;
51 : }
52 : static GEN
53 545055 : u_is_gener_expo(ulong p, GEN L0)
54 : {
55 545055 : const ulong q = p >> 1;
56 : long i;
57 : GEN L;
58 545055 : if (!L0) L0 = u_odd_prime_divisors(q);
59 545055 : L = cgetg_copy(L0,&i);
60 1177593 : while (--i) L[i] = q / uel(L0,i);
61 545055 : return L;
62 : }
63 :
64 : int
65 1662762 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
66 : {
67 : long i;
68 1662762 : if (krouu(x, p) >= 0) return 0;
69 1402770 : for (i=lg(L)-1; i; i--)
70 : {
71 852934 : ulong t = Fl_powu(x, uel(L,i), p);
72 852934 : if (t == p_1 || t == 1) return 0;
73 : }
74 549836 : return 1;
75 : }
76 : /* assume p prime */
77 : ulong
78 1081451 : pgener_Fl_local(ulong p, GEN L0)
79 : {
80 1081451 : const pari_sp av = avma;
81 1081451 : const ulong p_1 = p-1;
82 : long x;
83 : GEN L;
84 1081451 : if (p <= 19) switch(p)
85 : { /* quick trivial cases */
86 63 : case 2: return 1;
87 116213 : case 7:
88 116213 : case 17: return 3;
89 420150 : default: return 2;
90 : }
91 545025 : L = u_is_gener_expo(p,L0);
92 545025 : for (x = 2;; x++)
93 1655012 : if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
94 : }
95 : ulong
96 575299 : pgener_Fl(ulong p) { return pgener_Fl_local(p, NULL); }
97 :
98 : /* L[i] = set of (p-1)/2l, l ODD prime divisor of p-1 (l=2 can be included,
99 : * but wasteful) */
100 : int
101 13712 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
102 : {
103 13712 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
104 13712 : if (t >= 0) return 0;
105 21765 : for (i = lg(L)-1; i; i--)
106 : {
107 14241 : GEN t = Fp_pow(x, gel(L,i), p);
108 14241 : if (equalii(t, p_1) || equali1(t)) return 0;
109 : }
110 7524 : return 1;
111 : }
112 :
113 : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
114 : GEN
115 412365 : pgener_Fp_local(GEN p, GEN L0)
116 : {
117 412365 : pari_sp av0 = avma;
118 : GEN x, p_1, L;
119 412365 : if (lgefint(p) == 3)
120 : {
121 : ulong z;
122 407461 : if (p[2] == 2) return gen_1;
123 288398 : if (L0) L0 = ZV_to_nv(L0);
124 288398 : z = pgener_Fl_local(uel(p,2), L0);
125 288398 : return gc_utoipos(av0, z);
126 : }
127 4904 : p_1 = subiu(p,1); L = is_gener_expo(p, L0);
128 4904 : x = utoipos(2);
129 9927 : for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
130 4904 : return gc_utoipos(av0, uel(x,2));
131 : }
132 :
133 : GEN
134 44282 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
135 :
136 : ulong
137 205857 : pgener_Zl(ulong p)
138 : {
139 205857 : if (p == 2) pari_err_DOMAIN("pgener_Zl","p","=",gen_2,gen_2);
140 : /* only p < 2^32 such that znprimroot(p) != znprimroot(p^2) */
141 205857 : if (p == 40487) return 10;
142 : #ifndef LONG_IS_64BIT
143 29829 : return pgener_Fl(p);
144 : #else
145 176028 : if (p < (1UL<<32)) return pgener_Fl(p);
146 : else
147 : {
148 30 : const pari_sp av = avma;
149 30 : const ulong p_1 = p-1;
150 : long x ;
151 30 : GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
152 30 : for (x=2;;x++)
153 102 : if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2)))
154 30 : return gc_ulong(av, x);
155 : }
156 : #endif
157 : }
158 :
159 : /* p prime. Return a primitive root modulo p^e, e > 1 */
160 : GEN
161 171101 : pgener_Zp(GEN p)
162 : {
163 171101 : if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
164 : else
165 : {
166 5 : const pari_sp av = avma;
167 5 : GEN p_1 = subiu(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
168 5 : GEN x = utoipos(2);
169 12 : for (;; x[2]++)
170 17 : if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
171 5 : return gc_utoipos(av, uel(x,2));
172 : }
173 : }
174 :
175 : static GEN
176 259 : gener_Zp(GEN q, GEN F)
177 : {
178 259 : GEN p = NULL;
179 259 : long e = 0;
180 259 : if (F)
181 : {
182 14 : GEN P = gel(F,1), E = gel(F,2);
183 14 : long i, l = lg(P);
184 42 : for (i = 1; i < l; i++)
185 : {
186 28 : p = gel(P,i);
187 28 : if (absequaliu(p, 2)) continue;
188 14 : if (i < l-1) pari_err_DOMAIN("znprimroot", "n","=",F,F);
189 14 : e = itos(gel(E,i));
190 : }
191 14 : if (!p) pari_err_DOMAIN("znprimroot", "n","=",F,F);
192 : }
193 : else
194 245 : e = Z_isanypower(q, &p);
195 259 : if (!BPSW_psp(e? p: q)) pari_err_DOMAIN("znprimroot", "n","=", q,q);
196 245 : return e > 1? pgener_Zp(p): pgener_Fp(q);
197 : }
198 :
199 : GEN
200 329 : znprimroot(GEN N)
201 : {
202 329 : pari_sp av = avma;
203 : GEN x, n, F;
204 :
205 329 : if ((F = check_arith_non0(N,"znprimroot")))
206 : {
207 14 : F = clean_Z_factor(F);
208 14 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
209 : }
210 322 : N = absi_shallow(N);
211 322 : if (abscmpiu(N, 4) <= 0) { set_avma(av); return mkintmodu(N[2]-1,N[2]); }
212 273 : switch(mod4(N))
213 : {
214 14 : case 0: /* N = 0 mod 4 */
215 14 : pari_err_DOMAIN("znprimroot", "n","=",N,N);
216 0 : x = NULL; break;
217 28 : case 2: /* N = 2 mod 4 */
218 28 : n = shifti(N,-1); /* becomes odd */
219 28 : x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
220 21 : break;
221 231 : default: /* N odd */
222 231 : x = gener_Zp(N,F);
223 224 : break;
224 : }
225 245 : return gc_GEN(av, mkintmod(x, N));
226 : }
227 :
228 : /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
229 : GEN
230 0 : rootsof1_Fp(GEN n, GEN p)
231 : {
232 0 : pari_sp av = avma;
233 0 : GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
234 0 : GEN z = pgener_Fp_local(p, L);
235 0 : z = Fp_pow(z, diviiexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
236 0 : return gc_INT(av, z);
237 : }
238 :
239 : GEN
240 3033 : rootsof1u_Fp(ulong n, GEN p)
241 : {
242 3033 : pari_sp av = avma;
243 3033 : GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
244 3033 : z = pgener_Fp_local(p, Flv_to_ZV(L));
245 3033 : z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
246 3033 : return gc_INT(av, z);
247 : }
248 :
249 : ulong
250 215577 : rootsof1_Fl(ulong n, ulong p)
251 : {
252 215577 : pari_sp av = avma;
253 215577 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
254 215577 : ulong z = pgener_Fl_local(p, L);
255 215577 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
256 215577 : return gc_ulong(av,z);
257 : }
258 :
259 : /*********************************************************************/
260 : /** INVERSE TOTIENT FUNCTION **/
261 : /*********************************************************************/
262 : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
263 : * primes. Return factor(N) */
264 : GEN
265 350651 : Z_factor_listP(GEN N, GEN L)
266 : {
267 350651 : long i, k, l = lg(L);
268 350651 : GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
269 1346688 : for (i = k = 1; i < l; i++)
270 : {
271 996037 : GEN p = gel(L,i);
272 996037 : long v = Z_pvalrem(N, p, &N);
273 996037 : if (v)
274 : {
275 792176 : gel(P,k) = p;
276 792176 : gel(E,k) = utoipos(v);
277 792176 : k++;
278 : }
279 : }
280 350651 : setlg(P, k);
281 350651 : setlg(E, k); return mkmat2(P,E);
282 : }
283 :
284 : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
285 : * L is a list of primes containing all prime divisors of n. */
286 : static long
287 621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
288 : {
289 621565 : pari_sp av = avma, av2;
290 : GEN k, D;
291 : long i, v;
292 621565 : if (m && mod2(n))
293 : {
294 270914 : if (!equali1(n)) return 0;
295 69986 : if (px) *px = gen_1;
296 69986 : return 1;
297 : }
298 350651 : D = divisors(Z_factor_listP(shifti(n, -1), L));
299 : /* loop through primes p > m, d = p-1 | n */
300 350651 : av2 = avma;
301 350651 : if (!m)
302 : { /* special case p = 2, d = 1 */
303 69986 : k = n;
304 69986 : for (v = 1;; v++) {
305 69986 : if (istotient_i(k, gen_2, L, px)) {
306 69986 : if (px) *px = shifti(*px, v);
307 69986 : return 1;
308 : }
309 0 : if (mod2(k)) break;
310 0 : k = shifti(k,-1);
311 : }
312 0 : set_avma(av2);
313 : }
314 1099462 : for (i = 1; i < lg(D); ++i)
315 : {
316 1001588 : GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
317 1001588 : if (m && cmpii(d, m) < 0) continue;
318 677782 : p = addiu(d, 1);
319 677782 : if (!isprime(p)) continue;
320 442064 : k = diviiexact(n, d);
321 481593 : for (v = 1;; v++) {
322 : GEN r;
323 481593 : if (istotient_i(k, p, L, px)) {
324 182791 : if (px) *px = mulii(*px, powiu(p, v));
325 182791 : return 1;
326 : }
327 298802 : k = dvmdii(k, p, &r);
328 298802 : if (r != gen_0) break;
329 : }
330 259273 : set_avma(av2);
331 : }
332 97874 : return gc_long(av,0);
333 : }
334 :
335 : /* find x such that phi(x) = n */
336 : long
337 70000 : istotient(GEN n, GEN *px)
338 : {
339 70000 : pari_sp av = avma;
340 70000 : if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
341 70000 : if (signe(n) < 1) return 0;
342 70000 : if (mod2(n))
343 : {
344 14 : if (!equali1(n)) return 0;
345 14 : if (px) *px = gen_1;
346 14 : return 1;
347 : }
348 69986 : if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
349 : {
350 69986 : if (!px) set_avma(av);
351 : else
352 69986 : *px = gc_INT(av, *px);
353 69986 : return 1;
354 : }
355 0 : return gc_long(av,0);
356 : }
357 :
358 : /*********************************************************************/
359 : /** KRONECKER SYMBOL **/
360 : /*********************************************************************/
361 : /* t = 3,5 mod 8 ? (= 2 not a square mod t) */
362 : static int
363 361519238 : ome(long t)
364 : {
365 361519238 : switch(t & 7)
366 : {
367 204239588 : case 3:
368 204239588 : case 5: return 1;
369 157279650 : default: return 0;
370 : }
371 : }
372 : /* t a t_INT, is t = 3,5 mod 8 ? */
373 : static int
374 5975452 : gome(GEN t)
375 5975452 : { return signe(t)? ome( mod2BIL(t) ): 0; }
376 :
377 : /* assume y odd, return kronecker(x,y) * s */
378 : static long
379 250635026 : krouu_s(ulong x, ulong y, long s)
380 : {
381 250635026 : ulong x1 = x, y1 = y, z;
382 1163425922 : while (x1)
383 : {
384 912790896 : long r = vals(x1);
385 912790896 : if (r)
386 : {
387 485677141 : if (odd(r) && ome(y1)) s = -s;
388 485677141 : x1 >>= r;
389 : }
390 912790896 : if (x1 & y1 & 2) s = -s;
391 912790896 : z = y1 % x1; y1 = x1; x1 = z;
392 : }
393 250635026 : return (y1 == 1)? s: 0;
394 : }
395 :
396 : long
397 12966820 : kronecker(GEN x, GEN y)
398 : {
399 12966820 : pari_sp av = avma;
400 12966820 : long s = 1, r;
401 : ulong xu;
402 :
403 12966820 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
404 12966820 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
405 12966820 : switch (signe(y))
406 : {
407 63 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
408 133 : case 0: return is_pm1(x);
409 : }
410 12966687 : r = vali(y);
411 12966687 : if (r)
412 : {
413 1384660 : if (!mpodd(x)) return gc_long(av,0);
414 336361 : if (odd(r) && gome(x)) s = -s;
415 336361 : y = shifti(y,-r);
416 : }
417 11918388 : x = modii(x,y);
418 14288224 : while (lgefint(x) > 3) /* x < y */
419 : {
420 : GEN z;
421 2369836 : r = vali(x);
422 2369836 : if (r)
423 : {
424 1292950 : if (odd(r) && gome(y)) s = -s;
425 1292950 : x = shifti(x,-r);
426 : }
427 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
428 2369836 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
429 2369836 : z = remii(y,x); y = x; x = z;
430 2369836 : if (gc_needed(av,2))
431 : {
432 0 : if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
433 0 : (void)gc_all(av, 2, &x, &y);
434 : }
435 : }
436 11918388 : xu = itou(x);
437 11918388 : if (!xu) return is_pm1(y)? s: 0;
438 11805820 : r = vals(xu);
439 11805820 : if (r)
440 : {
441 6242134 : if (odd(r) && gome(y)) s = -s;
442 6242134 : xu >>= r;
443 : }
444 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
445 11805820 : if (xu & mod2BIL(y) & 2) s = -s;
446 11805820 : return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
447 : }
448 :
449 : long
450 40089 : krois(GEN x, long y)
451 : {
452 : ulong yu;
453 40089 : long s = 1;
454 :
455 40089 : if (y <= 0)
456 : {
457 35 : if (y == 0) return is_pm1(x);
458 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
459 : }
460 : else
461 40054 : yu = (ulong)y;
462 40054 : if (!odd(yu))
463 : {
464 : long r;
465 18578 : if (!mpodd(x)) return 0;
466 12600 : r = vals(yu); yu >>= r;
467 12600 : if (odd(r) && gome(x)) s = -s;
468 : }
469 34076 : return krouu_s(umodiu(x, yu), yu, s);
470 : }
471 : /* assume y != 0 */
472 : long
473 37769593 : kroiu(GEN x, ulong y)
474 : {
475 : long r;
476 37769593 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
477 390219 : if (!mpodd(x)) return 0;
478 266947 : r = vals(y); y >>= r;
479 266947 : return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
480 : }
481 :
482 : /* assume y > 0, odd, return s * kronecker(x,y) */
483 : static long
484 195202 : krouodd(ulong x, GEN y, long s)
485 : {
486 : long r;
487 195202 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
488 28326 : if (!x) return 0; /* y != 1 */
489 28326 : r = vals(x);
490 28326 : if (r)
491 : {
492 14672 : if (odd(r) && gome(y)) s = -s;
493 14672 : x >>= r;
494 : }
495 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
496 28326 : if (x & mod2BIL(y) & 2) s = -s;
497 28326 : return krouu_s(umodiu(y,x), x, s);
498 : }
499 :
500 : long
501 159773 : krosi(long x, GEN y)
502 : {
503 159773 : const pari_sp av = avma;
504 159773 : long s = 1, r;
505 159773 : switch (signe(y))
506 : {
507 0 : case -1: y = negi(y); if (x < 0) s = -1; break;
508 0 : case 0: return (x==1 || x==-1);
509 : }
510 159773 : r = vali(y);
511 159773 : if (r)
512 : {
513 16884 : if (!odd(x)) return gc_long(av,0);
514 16884 : if (odd(r) && ome(x)) s = -s;
515 16884 : y = shifti(y,-r);
516 : }
517 159773 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
518 159773 : return gc_long(av, krouodd((ulong)x, y, s));
519 : }
520 :
521 : long
522 35429 : kroui(ulong x, GEN y)
523 : {
524 35429 : const pari_sp av = avma;
525 35429 : long s = 1, r;
526 35429 : switch (signe(y))
527 : {
528 0 : case -1: y = negi(y); break;
529 0 : case 0: return x==1UL;
530 : }
531 35429 : r = vali(y);
532 35429 : if (r)
533 : {
534 0 : if (!odd(x)) return gc_long(av,0);
535 0 : if (odd(r) && ome(x)) s = -s;
536 0 : y = shifti(y,-r);
537 : }
538 35429 : return gc_long(av, krouodd(x, y, s));
539 : }
540 :
541 : long
542 100889841 : kross(long x, long y)
543 : {
544 : ulong yu;
545 100889841 : long s = 1;
546 :
547 100889841 : if (y <= 0)
548 : {
549 67494 : if (y == 0) return (labs(x)==1);
550 67466 : yu = (ulong)-y; if (x < 0) s = -1;
551 : }
552 : else
553 100822347 : yu = (ulong)y;
554 100889813 : if (!odd(yu))
555 : {
556 : long r;
557 23976382 : if (!odd(x)) return 0;
558 17040405 : r = vals(yu); yu >>= r;
559 17040405 : if (odd(r) && ome(x)) s = -s;
560 : }
561 93953836 : x %= (long)yu; if (x < 0) x += yu;
562 93953836 : return krouu_s((ulong)x, yu, s);
563 : }
564 :
565 : long
566 106999771 : krouu(ulong x, ulong y)
567 : {
568 : long r;
569 106999771 : if (odd(y)) return krouu_s(x, y, 1);
570 17282 : if (!odd(x)) return 0;
571 17282 : r = vals(y); y >>= r;
572 17282 : return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
573 : }
574 :
575 : /*********************************************************************/
576 : /** HILBERT SYMBOL **/
577 : /*********************************************************************/
578 : /* x,y are t_INT or t_REAL */
579 : static long
580 7329 : mphilbertoo(GEN x, GEN y)
581 : {
582 7329 : long sx = signe(x), sy = signe(y);
583 7329 : if (!sx || !sy) return 0;
584 7329 : return (sx < 0 && sy < 0)? -1: 1;
585 : }
586 :
587 : long
588 140826 : hilbertii(GEN x, GEN y, GEN p)
589 : {
590 : pari_sp av;
591 : long oddvx, oddvy, z;
592 :
593 140826 : if (!p) return mphilbertoo(x,y);
594 133518 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
595 133518 : if (!signe(x) || !signe(y)) return 0;
596 133497 : av = avma;
597 133497 : oddvx = odd(Z_pvalrem(x,p,&x));
598 133497 : oddvy = odd(Z_pvalrem(y,p,&y));
599 : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
600 133497 : if (absequaliu(p, 2))
601 : {
602 12355 : z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
603 12355 : if (oddvx && gome(y)) z = -z;
604 12355 : if (oddvy && gome(x)) z = -z;
605 : }
606 : else
607 : {
608 121142 : z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
609 121142 : if (oddvx && kronecker(y,p) < 0) z = -z;
610 121142 : if (oddvy && kronecker(x,p) < 0) z = -z;
611 : }
612 133497 : return gc_long(av, z);
613 : }
614 :
615 : static void
616 196 : err_prec(void) { pari_err_PREC("hilbert"); }
617 : static void
618 161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
619 : static void
620 56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
621 :
622 : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
623 : * Return lift(x) provided it's p-adic accuracy is large enough to decide
624 : * hilbert()'s value [ problem at p = 2 ] */
625 : static GEN
626 420 : lift_intmod(GEN x, GEN *pp)
627 : {
628 420 : GEN p = *pp, N = gel(x,1);
629 420 : x = gel(x,2);
630 420 : if (!p)
631 : {
632 266 : *pp = p = N;
633 266 : switch(itos_or_0(p))
634 : {
635 126 : case 2:
636 126 : case 4: err_prec();
637 : }
638 140 : return x;
639 : }
640 154 : if (!signe(p)) err_oo(N);
641 112 : if (absequaliu(p,2))
642 42 : { if (vali(N) <= 2) err_prec(); }
643 : else
644 70 : { if (!dvdii(N,p)) err_p(N,p); }
645 28 : if (!signe(x)) err_prec();
646 21 : return x;
647 : }
648 : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
649 : * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
650 : * to decide hilbert()'s value [ problem at p = 2 ]*/
651 : static GEN
652 210 : lift_padic(GEN x, GEN *pp)
653 : {
654 210 : GEN p = *pp, q = padic_p(x), u = padic_u(x);
655 210 : if (!p) *pp = p = q;
656 147 : else if (!equalii(p,q)) err_p(p, q);
657 105 : if (absequaliu(p,2) && precp(x) <= 2) err_prec();
658 70 : if (!signe(u)) err_prec();
659 70 : return odd(valp(x))? mulii(p,u): u;
660 : }
661 :
662 : long
663 62314 : hilbert(GEN x, GEN y, GEN p)
664 : {
665 62314 : pari_sp av = avma;
666 62314 : long tx = typ(x), ty = typ(y);
667 :
668 62314 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
669 62314 : if (tx == t_REAL)
670 : {
671 77 : if (p && signe(p)) err_oo(p);
672 63 : switch (ty)
673 : {
674 7 : case t_INT:
675 7 : case t_REAL: return mphilbertoo(x,y);
676 0 : case t_FRAC: return mphilbertoo(x,gel(y,1));
677 56 : default: pari_err_TYPE2("hilbert",x,y);
678 : }
679 : }
680 62237 : if (ty == t_REAL)
681 : {
682 14 : if (p && signe(p)) err_oo(p);
683 14 : switch (tx)
684 : {
685 14 : case t_INT:
686 14 : case t_REAL: return mphilbertoo(x,y);
687 0 : case t_FRAC: return mphilbertoo(gel(x,1),y);
688 0 : default: pari_err_TYPE2("hilbert",x,y);
689 : }
690 : }
691 62223 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
692 62020 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
693 :
694 61964 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
695 61901 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
696 :
697 61824 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
698 61824 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
699 :
700 61824 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
701 61824 : if (p && !signe(p)) p = NULL;
702 61824 : return gc_long(av, hilbertii(x,y,p));
703 : }
704 :
705 : /*******************************************************************/
706 : /* SQUARE ROOT MODULO p */
707 : /*******************************************************************/
708 : static void
709 3355369 : checkp(ulong q, ulong p)
710 3355369 : { if (!q) pari_err_PRIME("Fl_nonsquare",utoipos(p)); }
711 : /* p = 1 (mod 4) prime, return the first quadratic nonresidue, a prime */
712 : static ulong
713 14510041 : nonsquare1_Fl(ulong p)
714 : {
715 : forprime_t S;
716 : ulong q;
717 14510041 : if ((p & 7UL) != 1) return 2UL;
718 5482388 : q = p % 3; if (q == 2) return 3UL;
719 2069871 : checkp(q, p);
720 2069864 : q = p % 5; if (q == 2 || q == 3) return 5UL;
721 751793 : checkp(q, p);
722 751793 : q = p % 7; if (q != 4 && q >= 3) return 7UL;
723 305113 : checkp(q, p);
724 : /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
725 305113 : u_forprime_init(&S, 11, 1967);
726 533705 : while ((q = u_forprime_next(&S)))
727 : {
728 533705 : if (krouu(q, p) < 0) return q;
729 228592 : checkp(q, p);
730 : }
731 0 : checkp(0, p);
732 : return 0; /*LCOV_EXCL_LINE*/
733 : }
734 : /* p > 2 a prime */
735 : ulong
736 7935 : nonsquare_Fl(ulong p)
737 7935 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
738 :
739 : /* allow pi = 0 */
740 : ulong
741 192246 : Fl_2gener_pre(ulong p, ulong pi)
742 : {
743 192246 : ulong p1 = p-1;
744 192246 : long e = vals(p1);
745 192246 : if (e == 1) return p1;
746 70209 : return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
747 : }
748 :
749 : ulong
750 65964 : Fl_2gener_pre_i(ulong ns, ulong p, ulong pi)
751 : {
752 65964 : ulong p1 = p-1;
753 65964 : long e = vals(p1);
754 65964 : if (e == 1) return p1;
755 25255 : return Fl_powu_pre(ns, p1 >> e, p, pi);
756 : }
757 :
758 : static ulong
759 16076830 : Fl_sqrt_i(ulong a, ulong y, ulong p)
760 : {
761 : long i, e, k;
762 : ulong p1, q, v, w;
763 :
764 16076830 : if (!a) return 0;
765 14524366 : p1 = p - 1; e = vals(p1);
766 14524366 : if (e == 0) /* p = 2 */
767 : {
768 651342 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
769 651335 : return ((a & 1) == 0)? 0: 1;
770 : }
771 13873024 : if (e == 1)
772 : {
773 6141097 : v = Fl_powu(a, (p+1) >> 2, p);
774 6141097 : if (Fl_sqr(v, p) != a) return ~0UL;
775 6136207 : p1 = p - v; if (v > p1) v = p1;
776 6136207 : return v;
777 : }
778 7731927 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
779 7731927 : p1 = Fl_powu(a, q >> 1, p); /* a ^ [(q-1)/2] */
780 7731927 : if (!p1) return 0;
781 7731927 : v = Fl_mul(a, p1, p);
782 7731927 : w = Fl_mul(v, p1, p);
783 7731927 : if (!y) y = Fl_powu(nonsquare1_Fl(p), q, p);
784 13007657 : while (w != 1)
785 : { /* a*w = v^2, y primitive 2^e-th root of 1
786 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
787 5277814 : p1 = Fl_sqr(w, p);
788 8859546 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr(p1, p);
789 5277814 : if (k == e) return ~0UL;
790 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
791 5275737 : p1 = y;
792 7296473 : for (i=1; i < e-k; i++) p1 = Fl_sqr(p1, p);
793 5275737 : y = Fl_sqr(p1, p); e = k;
794 5275737 : w = Fl_mul(y, w, p);
795 5275737 : v = Fl_mul(v, p1, p);
796 : }
797 7729843 : p1 = p - v; if (v > p1) v = p1;
798 7729843 : return v;
799 : }
800 :
801 : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. Allow pi = 0 */
802 : ulong
803 41545939 : Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
804 : {
805 : long i, e, k;
806 : ulong p1, q, v, w;
807 :
808 41545939 : if (!pi) return Fl_sqrt_i(a, y, p);
809 25469109 : if (!a) return 0;
810 25340409 : p1 = p - 1; e = vals(p1);
811 25340409 : if (e == 0) /* p = 2 */
812 : {
813 0 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
814 0 : return ((a & 1) == 0)? 0: 1;
815 : }
816 25340409 : if (e == 1)
817 : {
818 18585204 : v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
819 18585204 : if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
820 18585170 : p1 = p - v; if (v > p1) v = p1;
821 18585170 : return v;
822 : }
823 6755205 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
824 6755205 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
825 6755205 : if (!p1) return 0;
826 6755205 : v = Fl_mul_pre(a, p1, p, pi);
827 6755205 : w = Fl_mul_pre(v, p1, p, pi);
828 6755205 : if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
829 12769909 : while (w != 1)
830 : { /* a*w = v^2, y primitive 2^e-th root of 1
831 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
832 6014796 : p1 = Fl_sqr_pre(w,p,pi);
833 11158258 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
834 6014796 : if (k == e) return ~0UL;
835 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
836 6014704 : p1 = y;
837 7905459 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
838 6014704 : y = Fl_sqr_pre(p1, p, pi); e = k;
839 6014704 : w = Fl_mul_pre(y, w, p, pi);
840 6014704 : v = Fl_mul_pre(v, p1, p, pi);
841 : }
842 6755113 : p1 = p - v; if (v > p1) v = p1;
843 6755113 : return v;
844 : }
845 :
846 : ulong
847 16039063 : Fl_sqrt(ulong a, ulong p)
848 16039063 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0; return Fl_sqrt_pre_i(a, 0, p, pi); }
849 :
850 : ulong
851 25318771 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
852 25318771 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
853 :
854 : /* allow pi = 0 */
855 : static ulong
856 220460 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
857 : {
858 : ulong m, m1;
859 : long i;
860 : for (;;)
861 : {
862 220460 : m = m1 = Fl_powu_pre(random_Fl(p-1)+1UL, r, p, pi);
863 220460 : if (m==1) continue;
864 245552 : for (i=1; i<e; i++)
865 : {
866 97482 : m = Fl_powu_pre(m, l, p, pi);
867 97482 : if (m == 1) break;
868 : }
869 166631 : if (i==e) break;
870 : }
871 148070 : *pt_m = m; return m1;
872 : }
873 :
874 : /* Solve x^l = a in G = Fp^* of order p-1 = (l^e)*r; l prime, (r,l) = 1, e >= 1
875 : * y generates the l-Sylow of G, m = y^(l^(e-1)) != 1. Allow y = 0, in which
876 : * case m is ignored and (y,m) are computed from scratch */
877 : static ulong
878 234965 : Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
879 : {
880 : ulong u2, v, w, z, dl;
881 234965 : u2 = Fl_inv(l%r, r);
882 234965 : v = Fl_powu_pre(a, u2, p, pi);
883 234965 : w = Fl_powu_pre(v, l, p, pi); if (w == a) return v;
884 145218 : w = pi? Fl_mul_pre(w, Fl_inv(a, p), p, pi): Fl_div(w, a, p);
885 145204 : if (!y) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
886 175277 : while (w != 1)
887 : {
888 152656 : ulong k = 0, p1 = w;
889 : do
890 : {
891 199904 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
892 199904 : if (++k == e) return ULONG_MAX;
893 77321 : } while (p1 != 1);
894 : /* z = w^(l^(k-1)) has order l */
895 30073 : dl = Fl_log_pre(z, m, l, p, pi);
896 30073 : dl = Fl_neg(dl, l);
897 30073 : p1 = Fl_powu_pre(y, dl*upowuu(l,e-k-1), p, pi);
898 30073 : m = Fl_powu_pre(m, dl, p, pi);
899 30073 : e = k;
900 30073 : v = pi? Fl_mul_pre(p1,v,p,pi): Fl_mul(p1,v,p);
901 30073 : y = Fl_powu_pre(p1,l,p,pi);
902 30073 : w = pi? Fl_mul_pre(y,w,p,pi): Fl_mul(y,w,p);
903 : }
904 22621 : return v;
905 : }
906 :
907 : /* allow pi = 0 */
908 : ulong
909 232010 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
910 : {
911 : ulong r, e;
912 232010 : if (!a) return 0;
913 231998 : e = u_lvalrem(p-1, l, &r);
914 231998 : return Fl_sqrtl_raw(a, l, e, r, p, pi, 0, 0);
915 : }
916 : ulong
917 0 : Fl_sqrtl(ulong a, ulong l, ulong p)
918 0 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
919 0 : return Fl_sqrtl_pre(a, l, p, pi); }
920 :
921 : /* allow pi = 0 */
922 : ulong
923 233471 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
924 : {
925 233471 : ulong m, q = p-1, z;
926 233471 : ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
927 233471 : if (a==0)
928 : {
929 116389 : if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
930 116382 : if (zetan) *zetan = 1UL;
931 116382 : return 0;
932 : }
933 : /* a != 0 */
934 117082 : if (n==1)
935 : {
936 420 : if (zetan) *zetan = 1;
937 420 : return n < 0? Fl_inv(a,p): a;
938 : }
939 116662 : if (n==2)
940 : {
941 42839 : if (zetan) *zetan = p-1;
942 42839 : return Fl_sqrt_pre_i(a, 0, p, pi);
943 : }
944 73823 : if (a == 1 && !zetan) return a;
945 44346 : m = ugcd(nn, q);
946 44346 : z = 1;
947 44346 : if (m!=1)
948 : {
949 2878 : GEN F = factoru(m);
950 : long i, j, e;
951 : ulong r, zeta, y, l;
952 6153 : for (i = nbrows(F); i; i--)
953 : {
954 3324 : l = ucoeff(F,i,1);
955 3324 : j = ucoeff(F,i,2);
956 3324 : e = u_lvalrem(q,l, &r);
957 3324 : y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
958 3324 : if (zetan)
959 : {
960 1585 : ulong Y = Fl_powu_pre(y, upowuu(l,e-j), p, pi);
961 1585 : z = pi? Fl_mul_pre(z, Y, p, pi): Fl_mul(z, Y, p);
962 : }
963 3324 : if (a!=1)
964 : do
965 : {
966 2967 : a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
967 2953 : if (a==ULONG_MAX) return ULONG_MAX;
968 2918 : } while (--j);
969 : }
970 : }
971 44297 : if (m != nn)
972 : {
973 41489 : ulong qm = q/m, nm = (nn/m) % qm;
974 41489 : a = Fl_powu_pre(a, Fl_inv(nm, qm), p, pi);
975 : }
976 44297 : if (n < 0) a = Fl_inv(a, p);
977 44297 : if (zetan) *zetan = z;
978 44297 : return a;
979 : }
980 :
981 : ulong
982 233471 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
983 : {
984 233471 : ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
985 233471 : return Fl_sqrtn_pre(a, n, p, pi, zetan);
986 : }
987 :
988 : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
989 : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
990 : * and in average is about 2 or 2.5 times worse. But try both algorithms for
991 : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
992 : *
993 : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
994 : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
995 : * If (a|p)=1, then sqrt(a) is in F_p.
996 : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
997 :
998 : /* compute y^2, y = y[1] + y[2] X */
999 : static GEN
1000 0 : sqrt_Cipolla_sqr(void *data, GEN y)
1001 : {
1002 0 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
1003 0 : GEN u2 = sqri(u), v2 = sqri(v);
1004 0 : v = subii(sqri(addii(v,u)), addii(u2,v2));
1005 0 : u = addii(u2, mulii(v2,n));
1006 0 : retmkvec2(modii(u,p), modii(v,p));
1007 : }
1008 : /* compute (t+X) y^2 */
1009 : static GEN
1010 0 : sqrt_Cipolla_msqr(void *data, GEN y)
1011 : {
1012 0 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2);
1013 0 : ulong t = gel(data,4)[2];
1014 0 : GEN d = addii(u, mului(t,v)), d2 = sqri(d);
1015 0 : GEN b = remii(mulii(a,v), p);
1016 0 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
1017 0 : v = subii(d2, mulii(b,v));
1018 0 : retmkvec2(modii(u,p), modii(v,p));
1019 : }
1020 : /* assume a reduced mod p [ otherwise correct but inefficient ] */
1021 : static GEN
1022 0 : sqrt_Cipolla(GEN a, GEN p)
1023 : {
1024 : pari_sp av;
1025 : GEN u, n, y, pov2;
1026 : ulong t;
1027 :
1028 0 : if (kronecker(a, p) < 0) return NULL;
1029 0 : pov2 = shifti(p,-1); /* center to avoid multiplying by huge base*/
1030 0 : if (cmpii(a,pov2) > 0) a = subii(a,p);
1031 0 : av = avma;
1032 0 : for (t=1; ; t++, set_avma(av))
1033 : {
1034 0 : n = subsi((long)(t*t), a);
1035 0 : if (kronecker(n, p) < 0) break;
1036 : }
1037 :
1038 : /* compute (t+X)^((p-1)/2) =: u+vX */
1039 0 : u = utoipos(t);
1040 0 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
1041 : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
1042 : /* Now u+vX = (t+X)^((p-1)/2); thus
1043 : * (u+vX)(t+X) = sqrt(a) + 0 X
1044 : * Whence,
1045 : * sqrt(a) = (u+vt)t - v*a
1046 : * 0 = (u+vt)
1047 : * Thus a square root is v*a */
1048 0 : return Fp_mul(gel(y,2), a, p);
1049 : }
1050 :
1051 : /* Return NULL if p is found to be composite.
1052 : * p odd, q = (p-1)/2^oo is odd */
1053 : static GEN
1054 5909 : Fp_2gener_all(GEN q, GEN p)
1055 : {
1056 : long k;
1057 5909 : for (k = 2;; k++)
1058 12261 : {
1059 18170 : long i = kroui(k, p);
1060 18170 : if (i < 0) return Fp_pow(utoipos(k), q, p);
1061 12261 : if (i == 0) return NULL;
1062 : }
1063 : }
1064 :
1065 : /* Return NULL if p is found to be composite */
1066 : GEN
1067 3222 : Fp_2gener(GEN p)
1068 : {
1069 3222 : GEN q = subiu(p, 1);
1070 3222 : long e = Z_lvalrem(q, 2, &q);
1071 3222 : if (e == 0 && !equaliu(p,2)) return NULL;
1072 3222 : return Fp_2gener_all(q, p);
1073 : }
1074 :
1075 : GEN
1076 19392 : Fp_2gener_i(GEN ns, GEN p)
1077 : {
1078 19392 : GEN q = subiu(p,1);
1079 19392 : long e = vali(q);
1080 19392 : if (e == 1) return q;
1081 18129 : return Fp_pow(ns, shifti(q,-e), p);
1082 : }
1083 :
1084 : static GEN
1085 1458 : nonsquare_Fp(GEN p)
1086 : {
1087 : forprime_t T;
1088 : ulong a;
1089 1458 : if (mod4(p)==3) return gen_m1;
1090 1458 : if (mod8(p)==5) return gen_2;
1091 712 : u_forprime_init(&T, 3, ULONG_MAX);
1092 1397 : while((a = u_forprime_next(&T)))
1093 1397 : if (kroui(a,p) < 0) return utoi(a);
1094 0 : pari_err_PRIME("Fp_sqrt [modulus]",p);
1095 : return NULL; /* LCOV_EXCL_LINE */
1096 : }
1097 :
1098 : static GEN
1099 820 : Fp_rootsof1(ulong l, GEN p)
1100 : {
1101 820 : GEN z, pl = diviuexact(subis(p,1),l);
1102 : ulong a;
1103 : forprime_t T;
1104 820 : u_forprime_init(&T, 3, ULONG_MAX);
1105 1066 : while((a = u_forprime_next(&T)))
1106 : {
1107 1066 : z = Fp_pow(utoi(a), pl, p);
1108 1066 : if (!equali1(z)) return z;
1109 : }
1110 0 : pari_err_PRIME("Fp_sqrt [modulus]",p);
1111 : return NULL; /* LCOV_EXCL_LINE */
1112 : }
1113 :
1114 : static GEN
1115 351 : Fp_gausssum(long D, GEN p)
1116 : {
1117 351 : long i, l = labs(D);
1118 351 : GEN z = Fp_rootsof1(l, p);
1119 351 : GEN s = z, x = z;
1120 3436 : for(i = 2; i < l; i++)
1121 : {
1122 3085 : long k = kross(i,l);
1123 3085 : x = mulii(x, z);
1124 3085 : if (k==1) s = addii(s, x);
1125 1718 : else if (k==-1) s = subii(s, x);
1126 : }
1127 351 : return s;
1128 : }
1129 :
1130 : static GEN
1131 18957 : Fp_sqrts(long a, GEN p)
1132 : {
1133 18957 : long v = vals(a)>>1;
1134 18957 : GEN r = gen_0;
1135 18957 : a >>= v << 1;
1136 18957 : switch(a)
1137 : {
1138 1 : case 1:
1139 1 : r = gen_1;
1140 1 : break;
1141 1110 : case -1:
1142 1110 : if (mod4(p)==1)
1143 1110 : r = Fp_pow(nonsquare_Fp(p), shifti(p,-2),p);
1144 : else
1145 0 : r = NULL;
1146 1110 : break;
1147 151 : case 2:
1148 151 : if (mod8(p)==1)
1149 : {
1150 151 : GEN z = Fp_pow(nonsquare_Fp(p), shifti(p,-3),p);
1151 151 : r = Fp_mul(z,Fp_sub(gen_1,Fp_sqr(z,p),p),p);
1152 0 : } else if (mod8(p)==7)
1153 0 : r = Fp_pow(gen_2, shifti(addiu(p,1),-2),p);
1154 : else
1155 0 : return NULL;
1156 151 : break;
1157 197 : case -2:
1158 197 : if (mod8(p)==1)
1159 : {
1160 197 : GEN z = Fp_pow(nonsquare_Fp(p), shifti(p,-3),p);
1161 197 : r = Fp_mul(z,Fp_add(gen_1,Fp_sqr(z,p),p),p);
1162 0 : } else if (mod8(p)==3)
1163 0 : r = Fp_pow(gen_m2, shifti(addiu(p,1),-2),p);
1164 : else
1165 0 : return NULL;
1166 197 : break;
1167 469 : case -3:
1168 469 : if (umodiu(p,3)==1)
1169 : {
1170 469 : GEN z = Fp_rootsof1(3, p);
1171 469 : r = Fp_sub(z,Fp_sqr(z,p),p);
1172 : }
1173 : else
1174 0 : return NULL;
1175 469 : break;
1176 2201 : case 5: case 13: case 17: case 21: case 29: case 33:
1177 : case -7: case -11: case -15: case -19: case -23:
1178 2201 : if (umodiu(p,labs(a))==1)
1179 351 : r = Fp_gausssum(a,p);
1180 : else
1181 1850 : return gen_0;
1182 351 : break;
1183 14828 : default:
1184 14828 : return gen_0;
1185 : }
1186 2279 : return remii(shifti(r, v), p);
1187 : }
1188 :
1189 : static GEN
1190 75211 : Fp_sqrt_ii(GEN a, GEN y, GEN p)
1191 : {
1192 75211 : pari_sp av = avma;
1193 75211 : GEN q, v, w, p1 = subiu(p,1);
1194 75211 : long i, k, e = vali(p1), as;
1195 :
1196 : /* direct formulas more efficient */
1197 75211 : if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
1198 75211 : if (e == 1)
1199 : {
1200 17646 : q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
1201 17646 : v = Fp_pow(a, q, p);
1202 : /* must check equality in case (a/p) = -1 or p not prime */
1203 17646 : av = avma; e = equalii(Fp_sqr(v,p), a); set_avma(av);
1204 17646 : return e? v: NULL;
1205 : }
1206 57565 : as = itos_or_0(a);
1207 57565 : if (!as) as = itos_or_0(subii(a,p));
1208 57565 : if (as)
1209 : {
1210 18957 : GEN res = Fp_sqrts(as, p);
1211 18957 : if (!res) return gc_NULL(av);
1212 18957 : if (signe(res)) return gc_upto(av, res);
1213 : }
1214 55286 : if (e == 2)
1215 : { /* Atkin's formula */
1216 17924 : GEN I, a2 = shifti(a,1);
1217 17924 : if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
1218 17924 : q = shifti(p1, -3); /* (p-5)/8 */
1219 17924 : v = Fp_pow(a2, q, p);
1220 17924 : I = Fp_mul(a2, Fp_sqr(v,p), p); /* I^2 = -1 */
1221 17924 : v = Fp_mul(a, Fp_mul(v, subiu(I,1), p), p);
1222 : /* must check equality in case (a/p) = -1 or p not prime */
1223 17924 : av = avma; e = equalii(Fp_sqr(v,p), a); set_avma(av);
1224 17924 : return e? v: NULL;
1225 : }
1226 : /* On average, Cipolla is better than Tonelli/Shanks if and only if
1227 : * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
1228 37362 : if (e*(e-1) > 20 + 8 * expi(p)) return sqrt_Cipolla(a,p);
1229 : /* Tonelli-Shanks */
1230 37362 : av = avma; q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
1231 37362 : if (!y)
1232 : {
1233 2687 : y = Fp_2gener_all(q, p);
1234 2687 : if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
1235 : }
1236 37362 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
1237 37362 : v = Fp_mul(a, p1, p);
1238 37362 : w = Fp_mul(v, p1, p);
1239 88490 : while (!equali1(w))
1240 : { /* a*w = v^2, y primitive 2^e-th root of 1
1241 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1242 51171 : p1 = Fp_sqr(w,p);
1243 105888 : for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
1244 51171 : if (k == e) return NULL; /* p composite or (a/p) != 1 */
1245 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1246 51128 : p1 = y;
1247 73292 : for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
1248 51128 : y = Fp_sqr(p1, p); e = k;
1249 51128 : w = Fp_mul(y, w, p);
1250 51128 : v = Fp_mul(v, p1, p);
1251 51128 : if (gc_needed(av,1))
1252 : {
1253 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
1254 0 : (void)gc_all(av,3, &y,&w,&v);
1255 : }
1256 : }
1257 37319 : return v;
1258 : }
1259 :
1260 : /* Assume p is prime and return NULL if (a,p) = -1; y = NULL or generator
1261 : * of Fp^* 2-Sylow */
1262 : GEN
1263 6140894 : Fp_sqrt_i(GEN a, GEN y, GEN p)
1264 : {
1265 6140894 : pari_sp av = avma, av2;
1266 : GEN q;
1267 :
1268 6140894 : if (lgefint(p) == 3)
1269 : {
1270 6065571 : ulong pp = uel(p,2), u = umodiu(a, pp);
1271 6065571 : if (!u) return gen_0;
1272 4853719 : u = Fl_sqrt(u, pp);
1273 4853705 : return (u == ~0UL)? NULL: utoipos(u);
1274 : }
1275 75323 : a = modii(a, p); if (!signe(a)) return gen_0;
1276 75211 : a = Fp_sqrt_ii(a, y, p); if (!a) return gc_NULL(av);
1277 : /* smallest square root */
1278 74793 : av2 = avma; q = subii(p, a);
1279 74793 : if (cmpii(a, q) > 0) a = q; else set_avma(av2);
1280 74793 : return gc_INT(av, a);
1281 : }
1282 : GEN
1283 6087200 : Fp_sqrt(GEN a, GEN p) { return Fp_sqrt_i(a, NULL, p); }
1284 :
1285 : /*********************************************************************/
1286 : /** GCD & BEZOUT **/
1287 : /*********************************************************************/
1288 :
1289 : GEN
1290 55319625 : lcmii(GEN x, GEN y)
1291 : {
1292 : pari_sp av;
1293 : GEN a, b;
1294 55319625 : if (!signe(x) || !signe(y)) return gen_0;
1295 55319625 : av = avma; a = gcdii(x,y);
1296 55319625 : if (absequalii(a,y)) { set_avma(av); return absi(x); }
1297 11926415 : if (!equali1(a)) y = diviiexact(y,a);
1298 11926415 : b = mulii(x,y); setabssign(b); return gc_INT(av, b);
1299 : }
1300 :
1301 : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
1302 : * set *pd = gcd(x,N) */
1303 : GEN
1304 6110397 : Fp_invgen(GEN x, GEN N, GEN *pd)
1305 : {
1306 : GEN d, d0, e, v;
1307 6110397 : if (lgefint(N) == 3)
1308 : {
1309 5326687 : ulong dd, NN = N[2], xx = umodiu(x,NN);
1310 5326687 : if (!xx) { *pd = N; return gen_0; }
1311 5326687 : xx = Fl_invgen(xx, NN, &dd);
1312 5326687 : *pd = utoi(dd); return utoi(xx);
1313 : }
1314 783710 : *pd = d = bezout(x, N, &v, NULL);
1315 783710 : if (equali1(d)) return v;
1316 : /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
1317 686283 : e = diviiexact(N,d);
1318 686283 : d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
1319 686283 : if (equali1(d0)) return v;
1320 544565 : if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
1321 544565 : return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
1322 : }
1323 :
1324 : /*********************************************************************/
1325 : /** CHINESE REMAINDERS **/
1326 : /*********************************************************************/
1327 :
1328 : /* Chinese Remainder Theorem. x and y must have the same type (integermod,
1329 : * polymod, or polynomial/vector/matrix recursively constructed with these
1330 : * as coefficients). Creates (with the same type) a z in the same residue
1331 : * class as x and the same residue class as y, if it is possible.
1332 : *
1333 : * We also allow (during recursion) two identical objects even if they are
1334 : * not integermod or polymod. For example:
1335 : *
1336 : * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
1337 : * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
1338 : * ? chinese(x, y)
1339 : * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
1340 :
1341 : static GEN
1342 2462764 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
1343 : {
1344 2462764 : GEN z = gassoc_proto(f,x,NULL);
1345 2462757 : if (z == gen_1) retmkintmod(gen_0,gen_1);
1346 2462708 : return z;
1347 : }
1348 :
1349 : GEN
1350 2415 : chinese1(GEN x) { return gen_chinese(x,chinese); }
1351 :
1352 : static GEN
1353 21 : padic2mod(GEN x)
1354 : {
1355 21 : pari_sp av = avma;
1356 21 : GEN pd = padic_pd(x), p = padic_p(x), u = padic_u(x);
1357 21 : long v = valp(x);
1358 21 : if (v < 0) pari_err_INV("chinese", mkintmod(gen_0, p));
1359 21 : if (v)
1360 : {
1361 0 : GEN pv = powiu(p, v);
1362 0 : pd = mulii(pd, pv);
1363 0 : u = mulii(u, pv);
1364 : }
1365 21 : return gc_GEN(av, mkintmod(u, pd));
1366 :
1367 : }
1368 : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod): makes Mod(0,1)
1369 : * a better "neutral" element */
1370 : static GEN
1371 21 : intmod2polmod(GEN x,GEN y)
1372 21 : { retmkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1)))); }
1373 :
1374 : GEN
1375 5495 : chinese(GEN x, GEN y)
1376 : {
1377 5495 : pari_sp av = avma;
1378 : long tx, ty;
1379 : GEN z;
1380 :
1381 5495 : if (!y) return chinese1(x);
1382 5446 : if (gidentical(x,y)) return gcopy(x);
1383 : /* allows GC optimization for this most frequent case */
1384 5439 : z = cgetg(3,t_INTMOD);
1385 5439 : tx = typ(x); if (tx == t_PADIC) { x = padic2mod(x); tx = t_INTMOD; }
1386 5439 : ty = typ(y); if (ty == t_PADIC) { y = padic2mod(y); ty = t_INTMOD; }
1387 5439 : if (tx == t_POLMOD && ty == t_INTMOD)
1388 14 : { y = intmod2polmod(y, x); ty = t_POLMOD; }
1389 5439 : if (ty == t_POLMOD && tx == t_INTMOD)
1390 7 : { x = intmod2polmod(x, y); tx = t_POLMOD; }
1391 5439 : if (tx == ty) switch(tx)
1392 : {
1393 3892 : case t_POLMOD:
1394 : {
1395 3892 : GEN A = gel(x,1), B = gel(y,1);
1396 3892 : GEN a = gel(x,2), b = gel(y,2), t, d, e, u, v;
1397 3892 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
1398 3892 : if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
1399 3892 : d = RgX_extgcd(A,B,&u,&v);
1400 3892 : e = gsub(b, a);
1401 3892 : if (!gequal0(gmod(e, d))) pari_err_OP("chinese",x,y);
1402 3892 : t = gdiv(A, d);
1403 3892 : e = gadd(a, gmul(gmul(u,t), e));
1404 :
1405 3892 : z = cgetg(3, t_POLMOD);
1406 3892 : gel(z,1) = RgX_mul(t, B);
1407 3892 : gel(z,2) = gmod(e, gel(z,1));
1408 3892 : return gc_upto(av, z);
1409 : }
1410 1519 : case t_INTMOD:
1411 : {
1412 1519 : GEN A = gel(x,1), B = gel(y,1);
1413 1519 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
1414 1519 : Z_chinese_pre(A, B, &C, &U, &d);
1415 1519 : c = Z_chinese_post(a, b, C, U, d);
1416 1519 : if (!c) pari_err_OP("chinese", x,y);
1417 1519 : set_avma((pari_sp)z); /* GC optimization */
1418 1519 : gel(z,1) = icopy(C);
1419 1519 : gel(z,2) = icopy(c); return z;
1420 : }
1421 14 : case t_POL:
1422 : {
1423 14 : long i, lx = lg(x), ly = lg(y);
1424 14 : if (varn(x) != varn(y)) pari_err_OP("chinese",x,y);
1425 14 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
1426 14 : set_avma(av);
1427 14 : z = cgetg(lx, t_POL); z[1] = x[1];
1428 42 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
1429 14 : if (i < lx)
1430 : {
1431 14 : GEN _0 = Rg_get_0(y);
1432 28 : for ( ; i<lx; i++) gel(z,i) = chinese(gel(x,i),_0);
1433 : }
1434 14 : return z;
1435 : }
1436 14 : case t_VEC: case t_COL: case t_MAT:
1437 : {
1438 : long i, lx;
1439 14 : set_avma(av);
1440 14 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) pari_err_OP("chinese",x,y);
1441 42 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
1442 14 : return z;
1443 : }
1444 : }
1445 0 : pari_err_OP("chinese",x,y);
1446 : return NULL; /* LCOV_EXCL_LINE */
1447 : }
1448 :
1449 : /* init chinese(Mod(.,A), Mod(.,B)) */
1450 : void
1451 451996 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
1452 : {
1453 451996 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
1454 451996 : GEN t = diviiexact(A,d);
1455 451996 : *pU = mulii(u, t);
1456 451996 : *pC = mulii(t, B); if (pd) *pd = d;
1457 451996 : }
1458 : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
1459 : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
1460 : * If d not NULL, check whether a = b mod d. */
1461 : GEN
1462 3251814 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
1463 : {
1464 : GEN e;
1465 3251814 : if (!signe(a))
1466 : {
1467 1010985 : if (d && !dvdii(b, d)) return NULL;
1468 1010985 : return Fp_mul(b, U, C);
1469 : }
1470 2240829 : e = subii(b,a);
1471 2240829 : if (d && !dvdii(e, d)) return NULL;
1472 2240829 : return modii(addii(a, mulii(U, e)), C);
1473 : }
1474 : static ulong
1475 1645000 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
1476 : {
1477 1645000 : if (!a) return Fl_mul(b, U, C);
1478 1645000 : return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
1479 : }
1480 :
1481 : GEN
1482 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
1483 : {
1484 2142 : pari_sp av = avma;
1485 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
1486 2142 : return gc_INT(av, Z_chinese_post(a,b, C, U, NULL));
1487 : }
1488 : GEN
1489 448279 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
1490 : {
1491 448279 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
1492 448279 : return Z_chinese_post(a,b, *pC, U, NULL);
1493 : }
1494 :
1495 : /* return lift(chinese(a mod A, b mod B))
1496 : * assume(A,B)=1, a,b,A,B integers and C = A*B */
1497 : GEN
1498 546014 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
1499 : {
1500 546014 : pari_sp av = avma;
1501 546014 : GEN U = mulii(Fp_inv(A,B), A);
1502 546014 : return gc_INT(av, Z_chinese_post(a,b,C,U, NULL));
1503 : }
1504 : ulong
1505 1645000 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
1506 1645000 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
1507 :
1508 : /* chinese1 for coprime moduli in Z */
1509 : static GEN
1510 2253538 : chinese1_coprime_Z_aux(GEN x, GEN y)
1511 : {
1512 2253538 : GEN z = cgetg(3, t_INTMOD);
1513 2253538 : GEN A = gel(x,1), a = gel(x, 2);
1514 2253538 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
1515 2253538 : pari_sp av = avma;
1516 2253538 : GEN U = mulii(Fp_inv(A,B), A);
1517 2253538 : gel(z,2) = gc_INT(av, Z_chinese_post(a,b,C,U, NULL));
1518 2253538 : gel(z,1) = C; return z;
1519 : }
1520 : GEN
1521 2460349 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
1522 :
1523 : /*********************************************************************/
1524 : /** MODULAR EXPONENTIATION **/
1525 : /*********************************************************************/
1526 : /* xa ZV or nv */
1527 : GEN
1528 2631145 : ZV_producttree(GEN xa)
1529 : {
1530 2631145 : long n = lg(xa)-1;
1531 2631145 : long m = n==1 ? 1: expu(n-1)+1;
1532 2631145 : GEN T = cgetg(m+1, t_VEC), t;
1533 : long i, j, k;
1534 2631145 : t = cgetg(((n+1)>>1)+1, t_VEC);
1535 2631145 : if (typ(xa)==t_VECSMALL)
1536 : {
1537 3598772 : for (j=1, k=1; k<n; j++, k+=2)
1538 2324708 : gel(t, j) = muluu(xa[k], xa[k+1]);
1539 1274064 : if (k==n) gel(t, j) = utoi(xa[k]);
1540 : } else {
1541 2856730 : for (j=1, k=1; k<n; j++, k+=2)
1542 1499649 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
1543 1357081 : if (k==n) gel(t, j) = icopy(gel(xa,k));
1544 : }
1545 2631145 : gel(T,1) = t;
1546 4243507 : for (i=2; i<=m; i++)
1547 : {
1548 1612362 : GEN u = gel(T, i-1);
1549 1612362 : long n = lg(u)-1;
1550 1612362 : t = cgetg(((n+1)>>1)+1, t_VEC);
1551 3622288 : for (j=1, k=1; k<n; j++, k+=2)
1552 2009926 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
1553 1612362 : if (k==n) gel(t, j) = gel(u, k);
1554 1612362 : gel(T, i) = t;
1555 : }
1556 2631145 : return T;
1557 : }
1558 :
1559 : /* not GC-clean */
1560 : GEN
1561 0 : ZMV_producttree(GEN xa)
1562 : {
1563 0 : long i, n = lg(xa)-1;
1564 : GEN T, worker;
1565 0 : long m = n==1 ? 1: expu(n-1)+1;
1566 : pari_timer ti;
1567 0 : if (DEBUGLEVEL>4) timer_start(&ti);
1568 0 : worker = snm_closure(is_entry("_ZM_mulrev"),NULL);
1569 0 : m = expu(n-1)+1; T = cgetg(m+1, t_VEC);
1570 0 : if (DEBUGLEVEL>5) err_printf("start ZMV Product tree:\nlevel 1: ");
1571 0 : gel(T,1) = gen_parpairwiseop_percent(worker,xa,DEBUGLEVEL>5);
1572 0 : if (DEBUGLEVEL>5) err_printf("\n");
1573 0 : if (m > 1)
1574 : {
1575 0 : for (i = 2; i < m-1; i++)
1576 : {
1577 0 : if (DEBUGLEVEL>5) err_printf("level %ld:",i);
1578 0 : gel(T, i) = gen_parpairwiseop_percent(worker,gel(T,i-1),DEBUGLEVEL>5);
1579 0 : if (DEBUGLEVEL>5) err_printf("\n");
1580 : }
1581 0 : if (m > 2)
1582 : {
1583 0 : if (DEBUGLEVEL>5) err_printf("level %ld:",m-1);
1584 0 : gel(T, m-1) = odd(lg(gel(T,m-2)))
1585 0 : ? mkvec2(RgM_ZM_mul(gmael(T,m-2,2), gmael(T,m-2,1)), RgM_ZM_mul(gmael(T,m-2,4), gmael(T,m-2,3)))
1586 0 : : mkvec2(RgM_ZM_mul(gmael(T,m-2,2), gmael(T,m-2,1)), gmael(T,m-2,3));
1587 : }
1588 0 : if (DEBUGLEVEL>5) err_printf("\nlevel %ld:",m);
1589 0 : gel(T, m) = mkvec(RgM_ZM_mul(gmael(T,m-1,2), gmael(T,m-1,1)));
1590 0 : if (DEBUGLEVEL>5) err_printf("\n");
1591 : }
1592 0 : if (DEBUGLEVEL>4) timer_printf(&ti,"ZMV_producttree");
1593 0 : return T;
1594 : }
1595 :
1596 : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
1597 : GEN
1598 58715462 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
1599 : {
1600 : long i,j,k;
1601 58715462 : long m = lg(T)-1, n = lg(P)-1;
1602 : GEN t;
1603 58715462 : GEN Tp = cgetg(m+1, t_VEC);
1604 58715462 : gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
1605 123156054 : for (i=m-1; i>=1; i--)
1606 : {
1607 64440592 : GEN u = gel(T, i);
1608 64440592 : GEN v = gel(Tp, i+1);
1609 64440592 : long n = lg(u)-1;
1610 64440592 : t = cgetg(n+1, t_VEC);
1611 155386351 : for (j=1, k=1; k<n; j++, k+=2)
1612 : {
1613 90945759 : gel(t, k) = modii(gel(v, j), gel(u, k));
1614 90945759 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
1615 : }
1616 64440592 : if (k==n) gel(t, k) = gel(v, j);
1617 64440592 : gel(Tp, i) = t;
1618 : }
1619 : {
1620 58715462 : GEN u = gel(T, i+1);
1621 58715462 : GEN v = gel(Tp, i+1);
1622 58715462 : long l = lg(u)-1;
1623 58715462 : if (typ(P)==t_VECSMALL)
1624 : {
1625 56084672 : GEN R = cgetg(n+1, t_VECSMALL);
1626 201108133 : for (j=1, k=1; j<=l; j++, k+=2)
1627 : {
1628 145023461 : uel(R,k) = umodiu(gel(v, j), P[k]);
1629 145023461 : if (k < n)
1630 114762501 : uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
1631 : }
1632 56084672 : return R;
1633 : }
1634 : else
1635 : {
1636 2630790 : GEN R = cgetg(n+1, t_VEC);
1637 7268550 : for (j=1, k=1; j<=l; j++, k+=2)
1638 : {
1639 4637760 : gel(R,k) = modii(gel(v, j), gel(P,k));
1640 4637760 : if (k < n)
1641 3821133 : gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
1642 : }
1643 2630790 : return R;
1644 : }
1645 : }
1646 : }
1647 :
1648 : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
1649 : GEN
1650 42246768 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
1651 : {
1652 42246768 : long m = lg(T)-1, n = lg(A)-1;
1653 : long i,j,k;
1654 42246768 : GEN Tp = cgetg(m+1, t_VEC);
1655 42246768 : GEN M = gel(T, 1);
1656 42246768 : GEN t = cgetg(lg(M), t_VEC);
1657 42246768 : if (typ(P)==t_VECSMALL)
1658 : {
1659 91883689 : for (j=1, k=1; k<n; j++, k+=2)
1660 : {
1661 67413232 : pari_sp av = avma;
1662 67413232 : GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
1663 67413232 : GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
1664 67413232 : gel(t, j) = gc_INT(av, tj);
1665 : }
1666 24470457 : if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
1667 : } else
1668 : {
1669 37892957 : for (j=1, k=1; k<n; j++, k+=2)
1670 : {
1671 20116646 : pari_sp av = avma;
1672 20116646 : GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
1673 20116646 : GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
1674 20116646 : gel(t, j) = gc_INT(av, tj);
1675 : }
1676 17776311 : if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
1677 : }
1678 42246768 : gel(Tp, 1) = t;
1679 79112547 : for (i=2; i<=m; i++)
1680 : {
1681 36865779 : GEN u = gel(T, i-1), M = gel(T, i);
1682 36865779 : GEN t = cgetg(lg(M), t_VEC);
1683 36865779 : GEN v = gel(Tp, i-1);
1684 36865779 : long n = lg(v)-1;
1685 98401069 : for (j=1, k=1; k<n; j++, k+=2)
1686 : {
1687 61535290 : pari_sp av = avma;
1688 61535290 : gel(t, j) = gc_INT(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
1689 61535290 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
1690 : }
1691 36865779 : if (k==n) gel(t, j) = gel(v, k);
1692 36865779 : gel(Tp, i) = t;
1693 : }
1694 42246768 : return gmael(Tp,m,1);
1695 : }
1696 :
1697 : static GEN
1698 1538431 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1699 : {
1700 1538431 : long i, l = lg(gel(vA,1)), n = lg(P);
1701 1538431 : GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
1702 35561613 : for (i=1; i < l; i++)
1703 : {
1704 34023182 : pari_sp av = avma;
1705 34023182 : GEN c, A = cgetg(n, typ(P));
1706 : long j;
1707 198984916 : for (j=1; j < n; j++) A[j] = mael(vA,j,i);
1708 34023182 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
1709 34023182 : gel(V,i) = gc_INT(av, c);
1710 : }
1711 1538431 : return V;
1712 : }
1713 :
1714 : static GEN
1715 875333 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1716 : {
1717 875333 : long i, j, l, n = lg(P);
1718 875333 : GEN mod = gmael(T, lg(T)-1, 1), V, w;
1719 875333 : w = cgetg(n, t_VECSMALL);
1720 3263087 : for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
1721 875333 : l = vecsmall_max(w);
1722 875333 : V = cgetg(l, t_POL);
1723 875333 : V[1] = mael(vA,1,1);
1724 6257721 : for (i=2; i < l; i++)
1725 : {
1726 5382388 : pari_sp av = avma;
1727 5382388 : GEN c, A = cgetg(n, typ(P));
1728 5382388 : if (typ(P)==t_VECSMALL)
1729 15841942 : for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
1730 : else
1731 6315112 : for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
1732 5382388 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
1733 5382388 : gel(V,i) = gc_INT(av, c);
1734 : }
1735 875333 : return ZX_renormalize(V, l);
1736 : }
1737 :
1738 : static GEN
1739 6585 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1740 : {
1741 6585 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1742 6585 : GEN A = cgetg(n, t_VEC);
1743 6585 : GEN V = cgetg(l, t_COL);
1744 149190 : for (i=1; i < l; i++)
1745 : {
1746 732416 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1747 142605 : gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
1748 : }
1749 6585 : return V;
1750 : }
1751 :
1752 : static GEN
1753 377545 : polint_chinese(GEN worker, GEN mA, GEN P)
1754 : {
1755 377545 : long cnt, pending, n, i, j, l = lg(gel(mA,1));
1756 : struct pari_mt pt;
1757 : GEN done, va, M, A;
1758 : pari_timer ti;
1759 :
1760 377545 : if (l == 1) return cgetg(1, t_MAT);
1761 369937 : cnt = pending = 0;
1762 369937 : n = lg(P);
1763 369937 : A = cgetg(n, t_VEC);
1764 369937 : va = mkvec(A);
1765 369937 : M = cgetg(l, t_MAT);
1766 369937 : if (DEBUGLEVEL>4) timer_start(&ti);
1767 369937 : if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
1768 369937 : mt_queue_start_lim(&pt, worker, l-1);
1769 1399258 : for (i=1; i<l || pending; i++)
1770 : {
1771 : long workid;
1772 4027808 : for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
1773 1029321 : mt_queue_submit(&pt, i, i<l? va: NULL);
1774 1029321 : done = mt_queue_get(&pt, &workid, &pending);
1775 1029321 : if (done)
1776 : {
1777 988366 : gel(M,workid) = done;
1778 988366 : if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
1779 : }
1780 : }
1781 369937 : if (DEBUGLEVEL>5) err_printf("\n");
1782 369937 : if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
1783 369937 : mt_queue_end(&pt);
1784 369937 : return M;
1785 : }
1786 :
1787 : GEN
1788 1225 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
1789 : {
1790 1225 : return nxCV_polint_center_tree(vA, P, T, R, m2);
1791 : }
1792 :
1793 : static GEN
1794 491 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1795 : {
1796 491 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1797 491 : GEN A = cgetg(n, t_VEC);
1798 491 : GEN V = cgetg(l, t_MAT);
1799 5851 : for (i=1; i < l; i++)
1800 : {
1801 26720 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1802 5360 : gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
1803 : }
1804 491 : return V;
1805 : }
1806 :
1807 : static GEN
1808 104 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
1809 : {
1810 104 : GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
1811 104 : return polint_chinese(worker, mA, P);
1812 : }
1813 :
1814 : static GEN
1815 120512 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1816 : {
1817 120512 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1818 120512 : GEN A = cgetg(n, t_VEC);
1819 120512 : GEN V = cgetg(l, t_MAT);
1820 655988 : for (i=1; i < l; i++)
1821 : {
1822 3054717 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1823 535476 : gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
1824 : }
1825 120512 : return V;
1826 : }
1827 :
1828 : GEN
1829 987141 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
1830 : {
1831 987141 : return ncV_polint_center_tree(vA, P, T, R, m2);
1832 : }
1833 :
1834 : static GEN
1835 377441 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
1836 : {
1837 377441 : GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
1838 377441 : return polint_chinese(worker, mA, P);
1839 : }
1840 :
1841 : /* return [A mod P[i], i=1..#P] */
1842 : GEN
1843 0 : Z_ZV_mod(GEN A, GEN P)
1844 : {
1845 0 : pari_sp av = avma;
1846 0 : return gc_GEN(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
1847 : }
1848 : /* P a t_VECSMALL */
1849 : GEN
1850 0 : Z_nv_mod(GEN A, GEN P)
1851 : {
1852 0 : pari_sp av = avma;
1853 0 : return gc_leaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
1854 : }
1855 : /* B a ZX, T = ZV_producttree(P) */
1856 : GEN
1857 3078666 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
1858 : {
1859 : pari_sp av;
1860 3078666 : long i, j, l = lg(B), n = lg(A)-1;
1861 3078666 : GEN V = cgetg(n+1, t_VEC);
1862 13886535 : for (j=1; j <= n; j++)
1863 : {
1864 10807869 : gel(V, j) = cgetg(l, t_VECSMALL);
1865 10807869 : mael(V, j, 1) = B[1]&VARNBITS;
1866 : }
1867 3078666 : av = avma;
1868 17172082 : for (i=2; i < l; i++)
1869 : {
1870 14093416 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
1871 92311435 : for (j=1; j <= n; j++)
1872 78218019 : mael(V, j, i) = v[j];
1873 14093416 : set_avma(av);
1874 : }
1875 13886535 : for (j=1; j <= n; j++)
1876 10807869 : (void) Flx_renormalize(gel(V, j), l);
1877 3078666 : return V;
1878 : }
1879 :
1880 : static GEN
1881 1890880 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
1882 :
1883 : GEN
1884 254703 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
1885 : {
1886 254703 : pari_sp av = avma;
1887 254703 : long i, j, l = lg(P), n = lg(xa)-1;
1888 254703 : GEN V = cgetg(n+1, t_VEC);
1889 914845 : for (j=1; j <= n; j++)
1890 : {
1891 660142 : gel(V, j) = cgetg(l, t_POL);
1892 660142 : mael(V, j, 1) = P[1]&VARNBITS;
1893 : }
1894 2018756 : for (i=2; i < l; i++)
1895 : {
1896 1764053 : GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
1897 6978394 : for (j=1; j <= n; j++)
1898 5214341 : gmael(V, j, i) = gel(v,j);
1899 : }
1900 914845 : for (j=1; j <= n; j++)
1901 660142 : (void) FlxX_renormalize(gel(V, j), l);
1902 254703 : return gc_GEN(av, V);
1903 : }
1904 :
1905 : GEN
1906 5634 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
1907 : {
1908 5634 : pari_sp av = avma;
1909 5634 : long i, j, l = lg(C), n = lg(xa)-1;
1910 5634 : GEN V = cgetg(n+1, t_VEC);
1911 28373 : for (j = 1; j <= n; j++)
1912 22739 : gel(V, j) = cgetg(l, t_COL);
1913 132461 : for (i = 1; i < l; i++)
1914 : {
1915 126827 : GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
1916 701415 : for (j = 1; j <= n; j++)
1917 574588 : gmael(V, j, i) = gel(v,j);
1918 : }
1919 5634 : return gc_GEN(av, V);
1920 : }
1921 :
1922 : GEN
1923 491 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
1924 : {
1925 491 : pari_sp av = avma;
1926 491 : long i, j, l = lg(M), n = lg(xa)-1;
1927 491 : GEN V = cgetg(n+1, t_VEC);
1928 2484 : for (j=1; j <= n; j++)
1929 1993 : gel(V, j) = cgetg(l, t_MAT);
1930 5851 : for (i=1; i < l; i++)
1931 : {
1932 5360 : GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
1933 26720 : for (j=1; j <= n; j++)
1934 21360 : gmael(V, j, i) = gel(v,j);
1935 : }
1936 491 : return gc_GEN(av, V);
1937 : }
1938 :
1939 : GEN
1940 1246310 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
1941 : {
1942 : pari_sp av;
1943 1246310 : long i, j, l = lg(B), n = lg(A)-1;
1944 1246310 : GEN V = cgetg(n+1, t_VEC);
1945 6478756 : for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
1946 1246310 : av = avma;
1947 43175092 : for (i=1; i < l; i++)
1948 : {
1949 41928782 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
1950 223156076 : for (j=1; j <= n; j++) mael(V, j, i) = v[j];
1951 41928782 : set_avma(av);
1952 : }
1953 1246310 : return V;
1954 : }
1955 :
1956 : static GEN
1957 220534 : ZM_nv_mod_tree_t(GEN M, GEN xa, GEN T, long t)
1958 : {
1959 220534 : pari_sp av = avma;
1960 220534 : long i, j, l = lg(M), n = lg(xa)-1;
1961 220534 : GEN V = cgetg(n+1, t_VEC);
1962 1261614 : for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t);
1963 1466615 : for (i=1; i < l; i++)
1964 : {
1965 1246081 : GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
1966 6477255 : for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
1967 : }
1968 220534 : return gc_GEN(av, V);
1969 : }
1970 :
1971 : GEN
1972 215077 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
1973 215077 : { return ZM_nv_mod_tree_t(M, xa, T, t_MAT); }
1974 :
1975 : GEN
1976 5457 : ZVV_nv_mod_tree(GEN M, GEN xa, GEN T)
1977 5457 : { return ZM_nv_mod_tree_t(M, xa, T, t_VEC); }
1978 :
1979 : static GEN
1980 2627290 : ZV_sqr(GEN z)
1981 : {
1982 2627290 : long i,l = lg(z);
1983 2627290 : GEN x = cgetg(l, t_VEC);
1984 2627290 : if (typ(z)==t_VECSMALL)
1985 6381790 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
1986 : else
1987 4681307 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
1988 2627290 : return x;
1989 : }
1990 :
1991 : static GEN
1992 13754836 : ZT_sqr(GEN x)
1993 : {
1994 13754836 : if (typ(x) == t_INT) return sqri(x);
1995 17988619 : pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
1996 : }
1997 :
1998 : static GEN
1999 2627290 : ZV_invdivexact(GEN y, GEN x)
2000 : {
2001 2627290 : long i, l = lg(y);
2002 2627290 : GEN z = cgetg(l,t_VEC);
2003 2627290 : if (typ(x)==t_VECSMALL)
2004 6381790 : for (i=1; i<l; i++)
2005 : {
2006 5108081 : pari_sp av = avma;
2007 5108081 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
2008 5108081 : set_avma(av); gel(z,i) = utoi(a);
2009 : }
2010 : else
2011 4681307 : for (i=1; i<l; i++)
2012 3327726 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
2013 2627290 : return z;
2014 : }
2015 :
2016 : /* P t_VECSMALL or t_VEC of t_INT */
2017 : GEN
2018 2627290 : ZV_chinesetree(GEN P, GEN T)
2019 : {
2020 2627290 : GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
2021 2627290 : GEN mod = gmael(T,lg(T)-1,1);
2022 2627290 : return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
2023 : }
2024 :
2025 : static GEN
2026 972154 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
2027 : {
2028 972154 : if (!pt_mod)
2029 12585 : return gc_upto(av, a);
2030 : else
2031 : {
2032 959569 : GEN mod = gmael(T, lg(T)-1, 1);
2033 959569 : (void)gc_all(av, 2, &a, &mod);
2034 959569 : *pt_mod = mod;
2035 959569 : return a;
2036 : }
2037 : }
2038 :
2039 : GEN
2040 150488 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2041 : {
2042 150488 : pari_sp av = avma;
2043 150488 : GEN T = ZV_producttree(P);
2044 150488 : GEN R = ZV_chinesetree(P, T);
2045 150488 : GEN a = ZV_chinese_tree(A, P, T, R);
2046 150488 : GEN mod = gmael(T, lg(T)-1, 1);
2047 150488 : GEN ca = Fp_center(a, mod, shifti(mod,-1));
2048 150488 : return gc_chinese(av, T, ca, pt_mod);
2049 : }
2050 :
2051 : GEN
2052 5141 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
2053 : {
2054 5141 : pari_sp av = avma;
2055 5141 : GEN T = ZV_producttree(P);
2056 5141 : GEN R = ZV_chinesetree(P, T);
2057 5141 : GEN a = ZV_chinese_tree(A, P, T, R);
2058 5141 : return gc_chinese(av, T, a, pt_mod);
2059 : }
2060 :
2061 : GEN
2062 304105 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2063 : {
2064 304105 : pari_sp av = avma;
2065 304105 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2066 304105 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
2067 304105 : return gc_upto(av, a);
2068 : }
2069 :
2070 : GEN
2071 428623 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2072 : {
2073 428623 : pari_sp av = avma;
2074 428623 : GEN T = ZV_producttree(P);
2075 428623 : GEN R = ZV_chinesetree(P, T);
2076 428623 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2077 428623 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
2078 428623 : return gc_chinese(av, T, a, pt_mod);
2079 : }
2080 :
2081 : GEN
2082 10357 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2083 : {
2084 10357 : pari_sp av = avma;
2085 10357 : GEN T = ZV_producttree(P);
2086 10357 : GEN R = ZV_chinesetree(P, T);
2087 10357 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2088 10357 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
2089 10357 : return gc_chinese(av, T, a, pt_mod);
2090 : }
2091 :
2092 : GEN
2093 5457 : ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2094 : {
2095 5457 : pari_sp av = avma;
2096 5457 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2097 5457 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
2098 5457 : return gc_upto(av, a);
2099 : }
2100 :
2101 : GEN
2102 0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2103 : {
2104 0 : pari_sp av = avma;
2105 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2106 0 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
2107 0 : return gc_upto(av, a);
2108 : }
2109 :
2110 : GEN
2111 120512 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
2112 : {
2113 120512 : pari_sp av = avma;
2114 120512 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2115 120512 : GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
2116 120512 : return gc_upto(av, a);
2117 : }
2118 :
2119 : GEN
2120 377441 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2121 : {
2122 377441 : pari_sp av = avma;
2123 377441 : GEN T = ZV_producttree(P);
2124 377441 : GEN R = ZV_chinesetree(P, T);
2125 377441 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2126 377441 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
2127 377441 : return gc_chinese(av, T, a, pt_mod);
2128 : }
2129 :
2130 : GEN
2131 0 : nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2132 : {
2133 0 : pari_sp av = avma;
2134 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2135 0 : GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
2136 0 : return gc_upto(av, a);
2137 : }
2138 :
2139 : GEN
2140 0 : nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2141 : {
2142 0 : pari_sp av = avma;
2143 0 : GEN T = ZV_producttree(P);
2144 0 : GEN R = ZV_chinesetree(P, T);
2145 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2146 0 : GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
2147 0 : return gc_chinese(av, T, a, pt_mod);
2148 : }
2149 :
2150 : GEN
2151 491 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
2152 : {
2153 491 : pari_sp av = avma;
2154 491 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2155 491 : GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
2156 491 : return gc_upto(av, a);
2157 : }
2158 :
2159 : GEN
2160 104 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2161 : {
2162 104 : pari_sp av = avma;
2163 104 : GEN T = ZV_producttree(P);
2164 104 : GEN R = ZV_chinesetree(P, T);
2165 104 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2166 104 : GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
2167 104 : return gc_chinese(av, T, a, pt_mod);
2168 : }
2169 :
2170 : /**********************************************************************
2171 : ** Powering over (Z/NZ)^*, small N **
2172 : **********************************************************************/
2173 :
2174 : /* 2^n mod p; assume n > 1 */
2175 : static ulong
2176 13231965 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
2177 : {
2178 13231965 : ulong y = 2;
2179 13231965 : int j = 1+bfffo(n);
2180 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
2181 13231965 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
2182 592969817 : for (; j; n<<=1,j--)
2183 : {
2184 579737852 : y = Fl_sqr_pre(y,p,pi);
2185 579737852 : if (n & HIGHBIT) y = Fl_double(y, p);
2186 : }
2187 13231965 : return y;
2188 : }
2189 :
2190 : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
2191 : static ulong
2192 5552334 : Fl_2powu(ulong n, ulong p)
2193 : {
2194 5552334 : ulong y = 2;
2195 5552334 : int j = 1+bfffo(n);
2196 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
2197 5552334 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
2198 39084477 : for (; j; n<<=1,j--)
2199 : {
2200 33532143 : y = (y*y) % p;
2201 33532143 : if (n & HIGHBIT) y = Fl_double(y, p);
2202 : }
2203 5552334 : return y;
2204 : }
2205 :
2206 : /* allow pi = 0 */
2207 : ulong
2208 169557291 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
2209 : {
2210 : ulong y, z, n;
2211 169557291 : if (!pi) return Fl_powu(x, n0, p);
2212 167111854 : if (n0 <= 1)
2213 : { /* frequent special cases */
2214 14362746 : if (n0 == 1) return x;
2215 117882 : if (n0 == 0) return 1;
2216 : }
2217 152749108 : if (x <= 2)
2218 : {
2219 13503678 : if (x == 2) return Fl_2powu_pre(n0, p, pi);
2220 271713 : return x; /* 0 or 1 */
2221 : }
2222 139245430 : y = 1; z = x; n = n0;
2223 : for(;;)
2224 : {
2225 781931074 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
2226 781931074 : n>>=1; if (!n) return y;
2227 642685644 : z = Fl_sqr_pre(z,p,pi);
2228 : }
2229 : }
2230 :
2231 : ulong
2232 151990533 : Fl_powu(ulong x, ulong n0, ulong p)
2233 : {
2234 : ulong y, z, n;
2235 151990533 : if (n0 <= 2)
2236 : { /* frequent special cases */
2237 68327927 : if (n0 == 2) return Fl_sqr(x,p);
2238 34062566 : if (n0 == 1) return x;
2239 2225104 : if (n0 == 0) return 1;
2240 : }
2241 83662606 : if (x <= 1) return x; /* 0 or 1 */
2242 83101272 : if (p & HIGHMASK)
2243 7998507 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
2244 75102765 : if (x == 2) return Fl_2powu(n0, p);
2245 69550431 : y = 1; z = x; n = n0;
2246 : for(;;)
2247 : {
2248 331521076 : if (n&1) y = (y*z) % p;
2249 331521076 : n>>=1; if (!n) return y;
2250 261970645 : z = (z*z) % p;
2251 : }
2252 : }
2253 :
2254 : /* Reduce data dependency to maximize internal parallelism; allow pi = 0 */
2255 : GEN
2256 13223567 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
2257 : {
2258 : long i, k;
2259 13223567 : GEN z = cgetg(n + 2, t_VECSMALL);
2260 13223567 : z[1] = 1; if (n == 0) return z;
2261 13223567 : z[2] = x;
2262 13223567 : if (pi)
2263 : {
2264 90239933 : for (i = 3, k=2; i <= n; i+=2, k++)
2265 : {
2266 77229474 : z[i] = Fl_sqr_pre(z[k], p, pi);
2267 77229474 : z[i+1] = Fl_mul_pre(z[k], z[k+1], p, pi);
2268 : }
2269 13010459 : if (i==n+1) z[i] = Fl_sqr_pre(z[k], p, pi);
2270 : }
2271 213108 : else if (p & HIGHMASK)
2272 : {
2273 0 : for (i = 3, k=2; i <= n; i+=2, k++)
2274 : {
2275 0 : z[i] = Fl_sqr(z[k], p);
2276 0 : z[i+1] = Fl_mul(z[k], z[k+1], p);
2277 : }
2278 0 : if (i==n+1) z[i] = Fl_sqr(z[k], p);
2279 : }
2280 : else
2281 400596287 : for (i = 2; i <= n; i++) z[i+1] = (z[i] * x) % p;
2282 13223567 : return z;
2283 : }
2284 :
2285 : GEN
2286 296057 : Fl_powers(ulong x, long n, ulong p)
2287 : {
2288 296057 : return Fl_powers_pre(x, n, p, (p & HIGHMASK)? get_Fl_red(p): 0);
2289 : }
2290 :
2291 : /**********************************************************************
2292 : ** Powering over (Z/NZ)^*, large N **
2293 : **********************************************************************/
2294 : typedef struct muldata {
2295 : GEN (*sqr)(void * E, GEN x);
2296 : GEN (*mul)(void * E, GEN x, GEN y);
2297 : GEN (*mul2)(void * E, GEN x);
2298 : } muldata;
2299 :
2300 : /* modified Barrett reduction with one fold */
2301 : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
2302 :
2303 : static GEN
2304 14033 : Fp_invmBarrett(GEN p, long s)
2305 : {
2306 14033 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
2307 14033 : return mkvec2(Q,R);
2308 : }
2309 :
2310 : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
2311 : * a = r (mod N) */
2312 : static GEN
2313 8295405 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
2314 : {
2315 8295405 : pari_sp av = avma;
2316 8295405 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
2317 8295405 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
2318 8295405 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
2319 8295405 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
2320 8295405 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
2321 8295405 : GEN r = subii(A, mulii(q, N));
2322 8295405 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
2323 8295405 : if (signe(sr)<0) return gc_INT(av, r);
2324 4485746 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
2325 4485746 : if (signe(sr)<0) return gc_INT(av, r);
2326 153748 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
2327 153748 : return gc_INT(av, signe(sr)>=0 ? sr:r);
2328 : }
2329 :
2330 : /* Montgomery reduction */
2331 :
2332 : INLINE ulong
2333 694611 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
2334 :
2335 : struct montred
2336 : {
2337 : GEN N;
2338 : ulong inv;
2339 : };
2340 :
2341 : /* Montgomery reduction */
2342 : static GEN
2343 63544565 : _sqr_montred(void * E, GEN x)
2344 : {
2345 63544565 : struct montred * D = (struct montred *) E;
2346 63544565 : return red_montgomery(sqri(x), D->N, D->inv);
2347 : }
2348 :
2349 : /* Montgomery reduction */
2350 : static GEN
2351 6705660 : _mul_montred(void * E, GEN x, GEN y)
2352 : {
2353 6705660 : struct montred * D = (struct montred *) E;
2354 6705660 : return red_montgomery(mulii(x, y), D->N, D->inv);
2355 : }
2356 :
2357 : static GEN
2358 9376494 : _mul2_montred(void * E, GEN x)
2359 : {
2360 9376494 : struct montred * D = (struct montred *) E;
2361 9376494 : GEN z = shifti(_sqr_montred(E, x), 1);
2362 9376494 : long l = lgefint(D->N);
2363 9892130 : while (lgefint(z) > l) z = subii(z, D->N);
2364 9376494 : return z;
2365 : }
2366 :
2367 : static GEN
2368 24540931 : _sqr_remii(void* N, GEN x)
2369 24540931 : { return remii(sqri(x), (GEN) N); }
2370 :
2371 : static GEN
2372 1505175 : _mul_remii(void* N, GEN x, GEN y)
2373 1505175 : { return remii(mulii(x, y), (GEN) N); }
2374 :
2375 : static GEN
2376 3700801 : _mul2_remii(void* N, GEN x)
2377 3700801 : { return Fp_double(_sqr_remii(N, x), (GEN)N); }
2378 :
2379 : struct redbarrett
2380 : {
2381 : GEN iM, N;
2382 : long s;
2383 : };
2384 :
2385 : static GEN
2386 7578139 : _sqr_remiibar(void *E, GEN x)
2387 : {
2388 7578139 : struct redbarrett * D = (struct redbarrett *) E;
2389 7578139 : return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
2390 : }
2391 :
2392 : static GEN
2393 717266 : _mul_remiibar(void *E, GEN x, GEN y)
2394 : {
2395 717266 : struct redbarrett * D = (struct redbarrett *) E;
2396 717266 : return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
2397 : }
2398 :
2399 : static GEN
2400 1797278 : _mul2_remiibar(void *E, GEN x)
2401 : {
2402 1797278 : struct redbarrett * D = (struct redbarrett *) E;
2403 1797278 : return Fp_double(_sqr_remiibar(E, x), D->N);
2404 : }
2405 :
2406 : static long
2407 903800 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
2408 : {
2409 903800 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
2410 : {
2411 14033 : struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
2412 14033 : D->sqr = &_sqr_remiibar;
2413 14033 : D->mul = &_mul_remiibar;
2414 14033 : D->mul2 = &_mul2_remiibar;
2415 14033 : E->N = N;
2416 14033 : E->s = 1+(expi(N)>>1);
2417 14033 : E->iM = Fp_invmBarrett(N, E->s);
2418 14033 : *pt_E = (void*) E;
2419 14033 : return 0;
2420 : }
2421 889767 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
2422 : {
2423 694611 : struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
2424 694611 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
2425 694611 : D->sqr = &_sqr_montred;
2426 694611 : D->mul = &_mul_montred;
2427 694611 : D->mul2 = &_mul2_montred;
2428 694611 : E->N = N;
2429 694611 : E->inv = init_montdata(N);
2430 694611 : *pt_E = (void*) E;
2431 694611 : return 1;
2432 : }
2433 : else
2434 : {
2435 195156 : D->sqr = &_sqr_remii;
2436 195156 : D->mul = &_mul_remii;
2437 195156 : D->mul2 = &_mul2_remii;
2438 195156 : *pt_E = (void*) N;
2439 195156 : return 0;
2440 : }
2441 : }
2442 :
2443 : GEN
2444 1901394 : Fp_powu(GEN A, ulong k, GEN N)
2445 : {
2446 1901394 : long lN = lgefint(N);
2447 : int base_is_2, use_montgomery;
2448 : muldata D;
2449 : void *E;
2450 : pari_sp av;
2451 :
2452 1901394 : if (lN == 3) {
2453 312300 : ulong n = uel(N,2);
2454 312300 : return utoi( Fl_powu(umodiu(A, n), k, n) );
2455 : }
2456 1589094 : if (k <= 2)
2457 : { /* frequent special cases */
2458 957464 : if (k == 2) return Fp_sqr(A,N);
2459 375403 : if (k == 1) return A;
2460 0 : if (k == 0) return gen_1;
2461 : }
2462 631630 : av = avma; A = modii(A,N);
2463 631630 : base_is_2 = 0;
2464 631630 : if (lgefint(A) == 3) switch(A[2])
2465 : {
2466 908 : case 1: set_avma(av); return gen_1;
2467 34122 : case 2: base_is_2 = 1; break;
2468 : }
2469 :
2470 : /* TODO: Move this out of here and use for general modular computations */
2471 630722 : use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
2472 630722 : if (base_is_2)
2473 34122 : A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
2474 : else
2475 596600 : A = gen_powu_i(A, k, E, D.sqr, D.mul);
2476 630722 : if (use_montgomery)
2477 : {
2478 531534 : A = red_montgomery(A, N, ((struct montred *) E)->inv);
2479 531534 : if (cmpii(A, N) >= 0) A = subii(A,N);
2480 : }
2481 630722 : return gc_INT(av, A);
2482 : }
2483 :
2484 : GEN
2485 1346239 : Fp_pows(GEN A, long k, GEN N)
2486 : {
2487 1346239 : if (lgefint(N) == 3) {
2488 1322361 : ulong n = N[2];
2489 1322361 : ulong a = umodiu(A, n);
2490 1322361 : if (k < 0) {
2491 58634 : a = Fl_inv(a, n);
2492 58634 : k = -k;
2493 : }
2494 1322361 : return utoi( Fl_powu(a, (ulong)k, n) );
2495 : }
2496 23878 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
2497 23878 : return Fp_powu(A, (ulong)k, N);
2498 : }
2499 :
2500 : /* A^K mod N */
2501 : GEN
2502 38935240 : Fp_pow(GEN A, GEN K, GEN N)
2503 : {
2504 : pari_sp av;
2505 38935240 : long s, lN = lgefint(N), sA, sy;
2506 : int base_is_2, use_montgomery;
2507 : GEN y;
2508 : muldata D;
2509 : void *E;
2510 :
2511 38935240 : s = signe(K);
2512 38935240 : if (!s) return dvdii(A,N)? gen_0: gen_1;
2513 37879344 : if (lN == 3 && lgefint(K) == 3)
2514 : {
2515 37158823 : ulong n = N[2], a = umodiu(A, n);
2516 37158823 : if (s < 0) a = Fl_inv(a, n);
2517 37158823 : if (a <= 1) return utoi(a); /* 0 or 1 */
2518 33454794 : return utoi(Fl_powu(a, uel(K,2), n));
2519 : }
2520 :
2521 720521 : av = avma;
2522 720521 : if (s < 0) y = Fp_inv(A,N);
2523 : else
2524 : {
2525 718567 : y = modii(A,N);
2526 718567 : if (!signe(y)) { set_avma(av); return gen_0; }
2527 : }
2528 720521 : if (lgefint(K) == 3) return gc_INT(av, Fp_powu(y, K[2], N));
2529 :
2530 273298 : base_is_2 = 0;
2531 273298 : sy = abscmpii(y, shifti(N,-1)) > 0;
2532 273298 : if (sy) y = subii(N,y);
2533 273298 : sA = sy && mod2(K);
2534 273298 : if (lgefint(y) == 3) switch(y[2])
2535 : {
2536 220 : case 1: set_avma(av); return sA ? subis(N,1): gen_1;
2537 152679 : case 2: base_is_2 = 1; break;
2538 : }
2539 :
2540 : /* TODO: Move this out of here and use for general modular computations */
2541 273078 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
2542 273078 : if (base_is_2)
2543 152679 : y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
2544 : else
2545 120399 : y = gen_pow_i(y, K, E, D.sqr, D.mul);
2546 273078 : if (use_montgomery)
2547 : {
2548 163077 : y = red_montgomery(y, N, ((struct montred *) E)->inv);
2549 163077 : if (cmpii(y,N) >= 0) y = subii(y,N);
2550 : }
2551 273078 : if (sA) y = subii(N, y);
2552 273078 : return gc_INT(av,y);
2553 : }
2554 :
2555 : static GEN
2556 14231484 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
2557 : static GEN
2558 8134253 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
2559 : static GEN
2560 47162 : _Fp_one(void *E) { (void) E; return gen_1; }
2561 :
2562 : GEN
2563 105 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
2564 105 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
2565 :
2566 : GEN
2567 43694 : Fp_pow_table(GEN R, GEN n, GEN p)
2568 43694 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
2569 :
2570 : GEN
2571 5931 : Fp_powers(GEN x, long n, GEN p)
2572 : {
2573 5931 : if (lgefint(p) == 3)
2574 2463 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
2575 3468 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
2576 : }
2577 :
2578 : GEN
2579 504 : FpV_prod(GEN V, GEN p) { return gen_product(V, (void *)p, &_Fp_mul); }
2580 :
2581 : static GEN
2582 28601216 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
2583 : static GEN
2584 160 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
2585 :
2586 : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
2587 : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
2588 : equalii,equali1,Fp_easylog};
2589 :
2590 : static GEN
2591 889913 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
2592 : static GEN
2593 1175565 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
2594 : static GEN
2595 1086846 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
2596 : static GEN
2597 575346 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
2598 : static GEN
2599 34307 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
2600 : static int
2601 260724 : _Fp_equal0(GEN x) { return signe(x)==0; }
2602 : static GEN
2603 19075 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
2604 :
2605 : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
2606 : _Fp_inv,_Fp_equal0,_Fp_s};
2607 :
2608 6963 : const struct bb_field *get_Fp_field(void **E, GEN p)
2609 6963 : { *E = (void*)p; return &Fp_field; }
2610 :
2611 : /*********************************************************************/
2612 : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
2613 : /*********************************************************************/
2614 : ulong
2615 546648 : Fl_order(ulong a, ulong o, ulong p)
2616 : {
2617 546648 : pari_sp av = avma;
2618 : GEN m, P, E;
2619 : long i;
2620 546648 : if (a==1) return 1;
2621 447518 : if (!o) o = p-1;
2622 447518 : m = factoru(o);
2623 447518 : P = gel(m,1);
2624 447518 : E = gel(m,2);
2625 1270965 : for (i = lg(P)-1; i; i--)
2626 : {
2627 823447 : ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
2628 823447 : if (y == 1) o = t;
2629 782513 : else for (j = 1; j < e; j++)
2630 : {
2631 386450 : y = Fl_powu(y, l, p);
2632 386450 : if (y == 1) { o = t * upowuu(l, j); break; }
2633 : }
2634 : }
2635 447518 : return gc_ulong(av, o);
2636 : }
2637 :
2638 : /*Find the exact order of a assuming a^o==1*/
2639 : GEN
2640 136429 : Fp_order(GEN a, GEN o, GEN p) {
2641 136429 : if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
2642 : {
2643 63559 : ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
2644 63559 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
2645 : }
2646 72870 : return gen_order(a, o, (void*)p, &Fp_star);
2647 : }
2648 : GEN
2649 70 : Fp_factored_order(GEN a, GEN o, GEN p)
2650 70 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
2651 :
2652 : /* return order of a mod p^e, e > 0, pe = p^e */
2653 : static GEN
2654 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
2655 : {
2656 : GEN ap, op;
2657 70 : if (absequaliu(p, 2))
2658 : {
2659 56 : if (e == 1) return gen_1;
2660 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
2661 49 : if (mod4(a) == 1) op = gen_1; else { op = gen_2; a = Fp_sqr(a, pe); }
2662 : } else {
2663 14 : ap = (e == 1)? a: remii(a,p);
2664 14 : op = Fp_order(ap, subiu(p,1), p);
2665 14 : if (e == 1) return op;
2666 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
2667 : }
2668 49 : if (equali1(a)) return op;
2669 7 : return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
2670 : }
2671 :
2672 : GEN
2673 63 : znorder(GEN x, GEN o)
2674 : {
2675 63 : pari_sp av = avma;
2676 : GEN b, a;
2677 :
2678 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
2679 56 : b = gel(x,1); a = gel(x,2);
2680 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
2681 49 : if (!o)
2682 : {
2683 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
2684 35 : long i, l = lg(P);
2685 35 : o = gen_1;
2686 70 : for (i = 1; i < l; i++)
2687 : {
2688 35 : GEN p = gel(P,i);
2689 35 : long e = itos(gel(E,i));
2690 :
2691 35 : if (l == 2)
2692 35 : o = Zp_order(a, p, e, b);
2693 : else {
2694 0 : GEN pe = powiu(p,e);
2695 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
2696 : }
2697 : }
2698 35 : return gc_INT(av, o);
2699 : }
2700 14 : return Fp_order(a, o, b);
2701 : }
2702 :
2703 : /*********************************************************************/
2704 : /** DISCRETE LOGARITHM in (Z/nZ)* **/
2705 : /*********************************************************************/
2706 : static GEN
2707 56566 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
2708 : {
2709 56566 : pari_sp av = avma;
2710 : GEN h1, h2, F, G;
2711 56566 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
2712 33987 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
2713 : {
2714 126 : GEN M = cgetg(3, t_MAT);
2715 126 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
2716 126 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
2717 126 : return gc_upto(av, M);
2718 : }
2719 33861 : return gc_NULL(av);
2720 : }
2721 :
2722 : static GEN
2723 56566 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
2724 : {
2725 : GEN rel;
2726 56566 : do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
2727 56566 : while (!rel);
2728 126 : return rel;
2729 : }
2730 :
2731 : struct Fp_log_rel
2732 : {
2733 : GEN rel;
2734 : ulong prmax;
2735 : long nbrel, nbmax, nbgen;
2736 : };
2737 :
2738 : static long
2739 59731 : tr(long i) { return odd(i) ? (i+1)>>1: -(i>>1); }
2740 :
2741 : static long
2742 169813 : rt(long i) { return i>0 ? 2*i-1: -2*i; }
2743 :
2744 : /* add u^e */
2745 : static void
2746 2163 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
2747 : {
2748 2163 : pari_sp av = avma;
2749 2163 : long off = r->prmax+1;
2750 2163 : GEN F = cgetg(3, t_MAT);
2751 2163 : gel(F,1) = vecsmall_append(gel(z,1), off+rt(u));
2752 2163 : gel(F,2) = vecsmall_append(gel(z,2), e);
2753 2163 : gel(r->rel,++r->nbrel) = gc_upto(av, F);
2754 2163 : }
2755 :
2756 : /* add u^-1 v^-1 */
2757 : static void
2758 83825 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
2759 : {
2760 83825 : pari_sp av = avma;
2761 83825 : long off = r->prmax+1;
2762 83825 : GEN P = mkvecsmall2(off+rt(u),off+rt(v)), E = mkvecsmall2(-1,-1);
2763 83825 : GEN F = cgetg(3, t_MAT);
2764 83825 : gel(F,1) = vecsmall_concat(gel(z,1), P);
2765 83825 : gel(F,2) = vecsmall_concat(gel(z,2), E);
2766 83825 : gel(r->rel,++r->nbrel) = gc_upto(av, F);
2767 83825 : }
2768 :
2769 : /* Let p=C^2+c
2770 : * Solve h = (C+x)*(C+a)-p = 0 [mod l]
2771 : * h= -c+x*(C+a)+C*a = 0 [mod l]
2772 : * x = (c-C*a)/(C+a) [mod l]
2773 : * h = -c+C*(x+a)+a*x */
2774 : GEN
2775 30261 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
2776 : {
2777 30261 : pari_sp ltop = avma;
2778 30261 : long i, j, th, n = lg(pi)-1, rel = 1, ab = labs(a), ae;
2779 30261 : GEN sieve = zero_zv(2*ab+2)+1+ab;
2780 30261 : GEN L = cgetg(1+2*ab+2, t_VEC);
2781 30261 : pari_sp av = avma;
2782 30261 : GEN z, h = addis(C,a);
2783 30261 : if ((z = Z_issmooth_fact(h, prmax)))
2784 : {
2785 2170 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
2786 2170 : av = avma;
2787 : }
2788 14021056 : for (i=1; i<=n; i++)
2789 : {
2790 13990795 : ulong li = pi[i], s = sz[i], al = smodss(a,li);
2791 13990795 : ulong iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
2792 : long u;
2793 13990795 : if (!iv) continue;
2794 13673338 : u = Fl_add(Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv, li),ab%li,li)-ab;
2795 50391033 : for(j = u; j<=ab; j+=li) sieve[j] += s;
2796 : }
2797 30261 : if (a)
2798 : {
2799 30198 : long e = expi(mulis(C,a));
2800 30198 : th = e - expu(e) - 1;
2801 63 : } else th = -1;
2802 30261 : ae = a>=0 ? ab-1: ab;
2803 15539713 : for (j = 1-ab; j <= ae; j++)
2804 15509452 : if (sieve[j]>=th)
2805 : {
2806 109053 : GEN h = absi(addis(subii(mulis(C,a+j),c), a*j));
2807 109053 : if ((z = Z_issmooth_fact(h, prmax)))
2808 : {
2809 106757 : gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
2810 106757 : av = avma;
2811 2296 : } else set_avma(av);
2812 : }
2813 : /* j = a */
2814 30261 : if (sieve[a]>=th)
2815 : {
2816 448 : GEN h = absi(addiu(subii(mulis(C,2*a),c), a*a));
2817 448 : if ((z = Z_issmooth_fact(h, prmax)))
2818 364 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
2819 : }
2820 30261 : setlg(L, rel); return gc_GEN(ltop, L);
2821 : }
2822 :
2823 : static long
2824 63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
2825 : {
2826 : struct pari_mt pt;
2827 63 : long i, j, nb = 0;
2828 63 : GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
2829 : mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
2830 63 : long running, pending = 0;
2831 63 : GEN W = zerovec(r->nbgen);
2832 63 : mt_queue_start_lim(&pt, worker, r->nbgen);
2833 30459 : for (i = 0; (running = (i < r->nbgen)) || pending; i++)
2834 : {
2835 : GEN done;
2836 : long idx;
2837 30396 : mt_queue_submit(&pt, i, running ? mkvec(stoi(tr(i))): NULL);
2838 30396 : done = mt_queue_get(&pt, &idx, &pending);
2839 30396 : if (!done || lg(done)==1) continue;
2840 27636 : gel(W, idx+1) = done;
2841 27636 : nb += lg(done)-1;
2842 27636 : if (DEBUGLEVEL && (i&127)==0)
2843 0 : err_printf("%ld%% ",100*nb/r->nbmax);
2844 : }
2845 63 : mt_queue_end(&pt);
2846 26362 : for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
2847 : {
2848 : long ll, m;
2849 26299 : GEN L = gel(W,j);
2850 26299 : if (isintzero(L)) continue;
2851 23681 : ll = lg(L);
2852 109669 : for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
2853 : {
2854 85988 : GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
2855 85988 : if (v[1] == 1)
2856 2163 : addifsmooth1(r, h, v[2], v[3]);
2857 : else
2858 83825 : addifsmooth2(r, h, v[2], v[3]);
2859 : }
2860 : }
2861 63 : return j;
2862 : }
2863 :
2864 : static GEN
2865 837 : ECP_psi(GEN x, GEN y)
2866 : {
2867 837 : long prec = realprec(x);
2868 837 : GEN lx = glog(x, prec), ly = glog(y, prec);
2869 837 : GEN u = gdiv(lx, ly);
2870 837 : return gpow(u, gneg(u),prec);
2871 : }
2872 :
2873 : struct computeG
2874 : {
2875 : GEN C;
2876 : long bnd, nbi;
2877 : };
2878 :
2879 : static GEN
2880 837 : _computeG(void *E, GEN gen)
2881 : {
2882 837 : struct computeG * d = (struct computeG *) E;
2883 837 : GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
2884 837 : return gsub(gmul(gsqr(gen),ps),gmulgs(gaddgs(gen,d->nbi),3));
2885 : }
2886 :
2887 : static long
2888 63 : compute_nbgen(GEN C, long bnd, long nbi)
2889 : {
2890 : struct computeG d;
2891 63 : d.C = shifti(C, 1);
2892 63 : d.bnd = bnd;
2893 63 : d.nbi = nbi;
2894 63 : return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
2895 : }
2896 :
2897 : static GEN
2898 1714 : _psi(void*E, GEN y)
2899 : {
2900 1714 : GEN lx = (GEN) E;
2901 1714 : long prec = realprec(lx);
2902 1714 : GEN ly = glog(y, prec);
2903 1714 : GEN u = gdiv(lx, ly);
2904 1714 : return gsub(gdiv(y, ly), gpow(u, u, prec));
2905 : }
2906 :
2907 : static GEN
2908 63 : opt_param(GEN x, long prec)
2909 : {
2910 63 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
2911 : }
2912 :
2913 : static GEN
2914 63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
2915 : {
2916 63 : pari_sp av = avma;
2917 63 : long lM = lg(M)-1, nbcol = lM;
2918 63 : long tbs = maxss(1, expu(nbg/expi(m)));
2919 : for (;;)
2920 42 : {
2921 105 : GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
2922 : GEN tab;
2923 105 : long i, f=0;
2924 105 : long l = lg(K), lm = lgefint(m);
2925 105 : GEN idx = diviiexact(subiu(p,1),m), g;
2926 : pari_timer ti;
2927 105 : if (DEBUGLEVEL) timer_start(&ti);
2928 210 : for(i=1; i<l; i++)
2929 210 : if (signe(gel(K,i)))
2930 105 : break;
2931 105 : g = Fp_pow(utoi(i), idx, p);
2932 105 : tab = Fp_pow_init(g, p, tbs, p);
2933 105 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
2934 121520 : for(i=1; i<l; i++)
2935 : {
2936 121415 : GEN k = gel(K,i);
2937 121415 : GEN j = i<=prmax ? utoi(i): addis(C,tr(i-(prmax+1)));
2938 121415 : if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
2939 82369 : gel(K,i) = cgetineg(lm);
2940 : else
2941 39046 : f++;
2942 : }
2943 105 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
2944 105 : if(f > (nbg>>1)) return gc_upto(av, K);
2945 10024 : for(i=1; i<=nbcol; i++)
2946 : {
2947 9982 : long a = 1+random_Fl(lM);
2948 9982 : swap(gel(M,a),gel(M,i));
2949 : }
2950 42 : if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
2951 : }
2952 : }
2953 :
2954 : static GEN
2955 126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
2956 : {
2957 126 : pari_sp av=avma;
2958 126 : GEN aa = gen_1;
2959 126 : long AV = 0;
2960 : for(;;)
2961 0 : {
2962 126 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
2963 126 : GEN F = gel(A,1), E = gel(A,2);
2964 126 : GEN Ao = gen_0;
2965 126 : long i, l = lg(F);
2966 807 : for(i=1; i<l; i++)
2967 : {
2968 681 : GEN Ki = gel(K,F[i]);
2969 681 : if (signe(Ki)<0) break;
2970 681 : Ao = addii(Ao, mulis(Ki, E[i]));
2971 : }
2972 126 : if (i==l) return Fp_divu(Ao, AV, m);
2973 0 : aa = gc_INT(av, aa);
2974 : }
2975 : }
2976 :
2977 : static GEN
2978 63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
2979 : {
2980 63 : pari_sp av = avma, av2;
2981 63 : long i, j, nbi, nbr = 0, nbrow, nbg;
2982 : GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
2983 : pari_timer ti;
2984 : struct Fp_log_rel r;
2985 63 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
2986 63 : ulong bnd = 4*bnds;
2987 63 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
2988 :
2989 63 : p_1 = subiu(p,1);
2990 63 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
2991 0 : m = diviiexact(p_1, Z_ppo(p_1, m));
2992 63 : pr = primes_upto_zv(bnd);
2993 63 : nbi = lg(pr)-1;
2994 63 : C = sqrtremi(p, &c);
2995 63 : av2 = avma;
2996 12796 : for (i = 1; i <= nbi; ++i)
2997 : {
2998 12733 : ulong lp = pr[i];
2999 26894 : while (lp <= bnd)
3000 : {
3001 14161 : nbr++;
3002 14161 : lp *= pr[i];
3003 : }
3004 : }
3005 63 : pi = cgetg(nbr+1,t_VECSMALL);
3006 63 : Ci = cgetg(nbr+1,t_VECSMALL);
3007 63 : ci = cgetg(nbr+1,t_VECSMALL);
3008 63 : sz = cgetg(nbr+1,t_VECSMALL);
3009 12796 : for (i = 1, j = 1; i <= nbi; ++i)
3010 : {
3011 12733 : ulong lp = pr[i], sp = expu(2*lp-1);
3012 26894 : while (lp <= bnd)
3013 : {
3014 14161 : pi[j] = lp;
3015 14161 : Ci[j] = umodiu(C, lp);
3016 14161 : ci[j] = umodiu(c, lp);
3017 14161 : sz[j] = sp;
3018 14161 : lp *= pr[i];
3019 14161 : j++;
3020 : }
3021 : }
3022 63 : r.nbrel = 0;
3023 63 : r.nbgen = compute_nbgen(C, bnd, nbi);
3024 63 : r.nbmax = 2*(nbi+r.nbgen);
3025 63 : r.rel = cgetg(r.nbmax+1,t_VEC);
3026 63 : r.prmax = pr[nbi];
3027 63 : if (DEBUGLEVEL)
3028 : {
3029 0 : err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
3030 0 : timer_start(&ti);
3031 : }
3032 63 : nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
3033 63 : nbrow = r.prmax + nbg;
3034 63 : if (DEBUGLEVEL)
3035 : {
3036 0 : err_printf("\n");
3037 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
3038 : }
3039 63 : setlg(r.rel,r.nbrel+1);
3040 63 : r.rel = gc_GEN(av2, r.rel);
3041 63 : K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
3042 63 : if (DEBUGLEVEL) timer_start(&ti);
3043 63 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
3044 63 : if (DEBUGLEVEL) timer_printf(&ti," log element");
3045 63 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
3046 63 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
3047 63 : d = gcdii(Ao,Bo);
3048 63 : l = Fp_div(diviiexact(Ao, d), diviiexact(Bo, d), m);
3049 63 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
3050 63 : return gc_INT(av, l);
3051 : }
3052 :
3053 : static int
3054 5039619 : Fp_log_use_index(long e, long p)
3055 : {
3056 5039619 : return (e >= 27 && 20*(p+6)<=e*e);
3057 : }
3058 :
3059 : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
3060 : static GEN
3061 8819103 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
3062 : {
3063 8819103 : pari_sp av = avma;
3064 8819103 : GEN p = (GEN)E;
3065 : /* assume a reduced mod p, p not necessarily prime */
3066 8819103 : if (equali1(a)) return gen_0;
3067 : /* p > 2 */
3068 5589725 : if (equalii(subiu(p,1), a)) /* -1 */
3069 : {
3070 : pari_sp av2;
3071 : GEN t;
3072 1382881 : ord = get_arith_Z(ord);
3073 1382881 : if (mpodd(ord)) retgc_const(av, cgetg(1, t_VEC)); /* no solution */
3074 1382867 : t = shifti(ord,-1); /* only possible solution */
3075 1382867 : av2 = avma;
3076 1382867 : if (!equalii(Fp_pow(g, t, p), a)) retgc_const(av, cgetg(1, t_VEC));
3077 1382839 : set_avma(av2); return gc_INT(av, t);
3078 : }
3079 4206844 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
3080 63 : return Fp_log_index(a, g, ord, p);
3081 4206781 : return gc_NULL(av); /* not easy */
3082 : }
3083 :
3084 : GEN
3085 4259498 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
3086 : {
3087 4259498 : GEN v = get_arith_ZZM(ord);
3088 4259470 : GEN F = gmael(v,2,1);
3089 4259470 : long lF = lg(F)-1, lmax;
3090 4259470 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
3091 4259442 : lmax = expi(gel(F,lF));
3092 4259442 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
3093 91 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
3094 4259442 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
3095 : }
3096 :
3097 : /* assume !(p & HIGHMASK) */
3098 : static ulong
3099 132376 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
3100 : {
3101 132376 : ulong i, h=1;
3102 364436 : for (i = 0; i < ord; i++, h = (h * g) % p)
3103 364436 : if (a==h) return i;
3104 0 : return ~0UL;
3105 : }
3106 :
3107 : static ulong
3108 29661 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
3109 : {
3110 29661 : ulong i, h=1;
3111 73901 : for (i = 0; i < ord; i++, h = Fl_mul_pre(h, g, p, pi))
3112 73901 : if (a==h) return i;
3113 0 : return ~0UL;
3114 : }
3115 :
3116 : static ulong
3117 0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
3118 : {
3119 0 : pari_sp av = avma;
3120 0 : GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
3121 0 : return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
3122 : }
3123 :
3124 : /* allow pi = 0 */
3125 : ulong
3126 30073 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
3127 : {
3128 30073 : if (!pi) return Fl_log(a, g, ord, p);
3129 29661 : if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
3130 0 : return Fl_log_Fp(a, g, ord, p);
3131 : }
3132 :
3133 : ulong
3134 132376 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
3135 : {
3136 132376 : if (ord <= 200)
3137 0 : return (p&HIGHMASK)? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
3138 132376 : : Fl_log_naive(a, g, ord, p);
3139 0 : return Fl_log_Fp(a, g, ord, p);
3140 : }
3141 :
3142 : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
3143 : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
3144 : static GEN
3145 126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
3146 : {
3147 126 : long l = lg(P) - 1, e = E[l];
3148 126 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
3149 : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
3150 :
3151 126 : if (l == 1) {
3152 98 : hpe = h;
3153 98 : gpe = g;
3154 : } else {
3155 28 : hpe = modii(h, pe);
3156 28 : gpe = modii(g, pe);
3157 : }
3158 126 : if (e == 1) {
3159 42 : hp = hpe;
3160 42 : gp = gpe;
3161 : } else {
3162 84 : hp = remii(hpe, p);
3163 84 : gp = remii(gpe, p);
3164 : }
3165 126 : if (hp == gen_0 || gp == gen_0) return NULL;
3166 105 : if (absequaliu(p, 2))
3167 : {
3168 35 : GEN n = int2n(e);
3169 35 : ogpe = Zp_order(gpe, gen_2, e, n);
3170 35 : a = Fp_log(hpe, gpe, ogpe, n);
3171 35 : if (typ(a) != t_INT) return NULL;
3172 : }
3173 : else
3174 : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
3175 : is trivial */
3176 : /* [order(gp), factor(order(gp))] */
3177 70 : GEN v = Fp_factored_order(gp, subiu(p,1), p);
3178 70 : GEN ogp = gel(v,1);
3179 70 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
3180 70 : a = Fp_log(hp, gp, v, p);
3181 70 : if (typ(a) != t_INT) return NULL;
3182 70 : if (e == 1) ogpe = ogp;
3183 : else
3184 : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
3185 : /* use p-adic log: O(log p + e) mul*/
3186 : long vpogpe, vpohpe;
3187 :
3188 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
3189 28 : gpe = Fp_pow(gpe, ogp, pe);
3190 : /* g,h = 1 mod p; compute b s.t. h = g^b */
3191 :
3192 : /* v_p(order g mod pe) */
3193 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
3194 : /* v_p(order h mod pe) */
3195 28 : vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
3196 28 : if (vpohpe > vpogpe) return NULL;
3197 :
3198 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
3199 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
3200 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
3201 28 : a = addii(a, mulii(ogp, padic_to_Q(b)));
3202 : }
3203 : }
3204 : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
3205 91 : if (l == 1) return a;
3206 :
3207 28 : N = diviiexact(N, pe); /* make N coprime to p */
3208 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
3209 28 : g = Fp_pow(g, modii(ogpe, phi), N);
3210 28 : setlg(P, l); /* remove last element */
3211 28 : setlg(E, l);
3212 28 : b = znlog_rec(h, g, N, P, E, PHI);
3213 28 : if (!b) return NULL;
3214 28 : return addmulii(a, b, ogpe);
3215 : }
3216 :
3217 : static GEN
3218 98 : get_PHI(GEN P, GEN E)
3219 : {
3220 98 : long i, l = lg(P);
3221 98 : GEN PHI = cgetg(l, t_VEC);
3222 98 : gel(PHI,1) = gen_1;
3223 126 : for (i=1; i<l-1; i++)
3224 : {
3225 28 : GEN t, p = gel(P,i);
3226 28 : long e = E[i];
3227 28 : t = mulii(powiu(p, e-1), subiu(p,1));
3228 28 : if (i > 1) t = mulii(t, gel(PHI,i));
3229 28 : gel(PHI,i+1) = t;
3230 : }
3231 98 : return PHI;
3232 : }
3233 :
3234 : GEN
3235 238 : znlog(GEN h, GEN g, GEN o)
3236 : {
3237 238 : pari_sp av = avma;
3238 : GEN N, fa, P, E, x;
3239 238 : switch (typ(g))
3240 : {
3241 28 : case t_PADIC:
3242 : {
3243 28 : GEN p = padic_p(g);
3244 28 : long v = valp(g);
3245 28 : if (v < 0) pari_err_DIM("znlog");
3246 28 : if (v > 0) {
3247 0 : long k = gvaluation(h, p);
3248 0 : if (k % v) return cgetg(1,t_VEC);
3249 0 : k /= v;
3250 0 : if (!gequal(h, gpowgs(g,k))) retgc_const(av, cgetg(1, t_VEC));
3251 0 : return gc_stoi(av, k);
3252 : }
3253 28 : N = padic_pd(g);
3254 28 : g = Rg_to_Fp(g, N);
3255 28 : break;
3256 : }
3257 203 : case t_INTMOD:
3258 203 : N = gel(g,1);
3259 203 : g = gel(g,2); break;
3260 7 : default: pari_err_TYPE("znlog", g);
3261 : return NULL; /* LCOV_EXCL_LINE */
3262 : }
3263 231 : if (equali1(N)) { set_avma(av); return gen_0; }
3264 231 : h = Rg_to_Fp(h, N);
3265 224 : if (o) return gc_upto(av, Fp_log(h, g, o, N));
3266 98 : fa = Z_factor(N);
3267 98 : P = gel(fa,1);
3268 98 : E = vec_to_vecsmall(gel(fa,2));
3269 98 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
3270 98 : if (!x) retgc_const(av, cgetg(1, t_VEC));
3271 63 : return gc_INT(av, x);
3272 : }
3273 :
3274 : GEN
3275 173539 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
3276 : {
3277 173539 : if (lgefint(p)==3)
3278 : {
3279 172921 : long nn = itos_or_0(n);
3280 172921 : if (nn)
3281 : {
3282 172921 : ulong pp = p[2];
3283 : ulong uz;
3284 172921 : ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
3285 172900 : if (r==ULONG_MAX) return NULL;
3286 172858 : if (zeta) *zeta = utoi(uz);
3287 172858 : return utoi(r);
3288 : }
3289 : }
3290 618 : a = modii(a,p);
3291 618 : if (!signe(a))
3292 : {
3293 0 : if (zeta) *zeta = gen_1;
3294 0 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
3295 0 : return gen_0;
3296 : }
3297 618 : if (absequaliu(n,2))
3298 : {
3299 418 : if (zeta) *zeta = subiu(p,1);
3300 418 : return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
3301 : }
3302 200 : return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
3303 : }
3304 :
3305 : /*********************************************************************/
3306 : /** FACTORIAL **/
3307 : /*********************************************************************/
3308 : GEN
3309 92879 : mulu_interval_step(ulong a, ulong b, ulong step)
3310 : {
3311 92879 : pari_sp av = avma;
3312 : ulong k, l, N, n;
3313 : long lx;
3314 : GEN x;
3315 :
3316 92879 : if (!a) return gen_0;
3317 92879 : if (step == 1) return mulu_interval(a, b);
3318 92879 : n = 1 + (b-a) / step;
3319 92879 : b -= (b-a) % step;
3320 92879 : if (n < 61)
3321 : {
3322 91495 : if (n == 1) return utoipos(a);
3323 70410 : x = muluu(a,a+step); if (n == 2) return x;
3324 548362 : for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
3325 55221 : return gc_INT(av, x);
3326 : }
3327 : /* step | b-a */
3328 1384 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3329 1384 : N = b + a;
3330 1384 : for (k = a;; k += step)
3331 : {
3332 227455 : l = N - k; if (l <= k) break;
3333 226071 : gel(x,lx++) = muluu(k,l);
3334 : }
3335 1384 : if (l == k) gel(x,lx++) = utoipos(k);
3336 1384 : setlg(x, lx);
3337 1384 : return gc_INT(av, ZV_prod(x));
3338 : }
3339 : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
3340 : * first is slower ... ] */
3341 : GEN
3342 166976 : mulu_interval(ulong a, ulong b)
3343 : {
3344 166976 : pari_sp av = avma;
3345 : ulong k, l, N, n;
3346 : long lx;
3347 : GEN x;
3348 :
3349 166976 : if (!a) return gen_0;
3350 166976 : n = b - a + 1;
3351 166976 : if (n < 61)
3352 : {
3353 166241 : if (n == 1) return utoipos(a);
3354 97914 : x = muluu(a,a+1); if (n == 2) return x;
3355 375701 : for (k=a+2; k<b; k++) x = mului(k,x);
3356 : /* avoid k <= b: broken if b = ULONG_MAX */
3357 80463 : return gc_INT(av, mului(b,x));
3358 : }
3359 735 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3360 735 : N = b + a;
3361 735 : for (k = a;; k++)
3362 : {
3363 29974 : l = N - k; if (l <= k) break;
3364 29239 : gel(x,lx++) = muluu(k,l);
3365 : }
3366 735 : if (l == k) gel(x,lx++) = utoipos(k);
3367 735 : setlg(x, lx);
3368 735 : return gc_INT(av, ZV_prod(x));
3369 : }
3370 : GEN
3371 595 : muls_interval(long a, long b)
3372 : {
3373 595 : pari_sp av = avma;
3374 595 : long lx, k, l, N, n = b - a + 1;
3375 : GEN x;
3376 :
3377 595 : if (a <= 0 && b >= 0) return gen_0;
3378 322 : if (n < 61)
3379 : {
3380 322 : x = stoi(a);
3381 518 : for (k=a+1; k<=b; k++) x = mulsi(k,x);
3382 322 : return gc_INT(av, x);
3383 : }
3384 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3385 0 : N = b + a;
3386 0 : for (k = a;; k++)
3387 : {
3388 0 : l = N - k; if (l <= k) break;
3389 0 : gel(x,lx++) = mulss(k,l);
3390 : }
3391 0 : if (l == k) gel(x,lx++) = stoi(k);
3392 0 : setlg(x, lx);
3393 0 : return gc_INT(av, ZV_prod(x));
3394 : }
3395 :
3396 : GEN
3397 105 : mpprimorial(long n)
3398 : {
3399 105 : pari_sp av = avma;
3400 105 : if (n <= 12) switch(n)
3401 : {
3402 14 : case 0: case 1: return gen_1;
3403 7 : case 2: return gen_2;
3404 14 : case 3: case 4: return utoipos(6);
3405 14 : case 5: case 6: return utoipos(30);
3406 28 : case 7: case 8: case 9: case 10: return utoipos(210);
3407 14 : case 11: case 12: return utoipos(2310);
3408 7 : default: pari_err_DOMAIN("primorial", "argument","<",gen_0,stoi(n));
3409 : }
3410 7 : return gc_INT(av, zv_prod_Z(primes_upto_zv(n)));
3411 : }
3412 :
3413 : GEN
3414 501083 : mpfact(long n)
3415 : {
3416 501083 : pari_sp av = avma;
3417 : GEN a, v;
3418 : long k;
3419 501083 : if (n <= 12) switch(n)
3420 : {
3421 431559 : case 0: case 1: return gen_1;
3422 25171 : case 2: return gen_2;
3423 3556 : case 3: return utoipos(6);
3424 4145 : case 4: return utoipos(24);
3425 2887 : case 5: return utoipos(120);
3426 2563 : case 6: return utoipos(720);
3427 2448 : case 7: return utoipos(5040);
3428 2451 : case 8: return utoipos(40320);
3429 2458 : case 9: return utoipos(362880);
3430 2715 : case 10:return utoipos(3628800);
3431 1409 : case 11:return utoipos(39916800);
3432 591 : case 12:return utoipos(479001600);
3433 0 : default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
3434 : }
3435 19130 : v = cgetg(expu(n) + 2, t_VEC);
3436 19130 : for (k = 1;; k++)
3437 89071 : {
3438 108201 : long m = n >> (k-1), l;
3439 108201 : if (m <= 2) break;
3440 89071 : l = (1 + (n >> k)) | 1;
3441 : /* product of odd numbers in ]n / 2^k, n / 2^(k-1)] */
3442 89071 : a = mulu_interval_step(l, m, 2);
3443 89071 : gel(v,k) = k == 1? a: powiu(a, k);
3444 : }
3445 89071 : a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
3446 19130 : a = shifti(a, factorial_lval(n, 2));
3447 19130 : return gc_INT(av, a);
3448 : }
3449 :
3450 : ulong
3451 57255 : factorial_Fl(long n, ulong p)
3452 : {
3453 : long k;
3454 : ulong v;
3455 57255 : if (p <= (ulong)n) return 0;
3456 57255 : v = Fl_powu(2, factorial_lval(n, 2), p);
3457 57255 : for (k = 1;; k++)
3458 143736 : {
3459 200991 : long m = n >> (k-1), l, i;
3460 200991 : ulong a = 1;
3461 200991 : if (m <= 2) break;
3462 143736 : l = (1 + (n >> k)) | 1;
3463 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3464 787514 : for (i=l; i<=m; i+=2)
3465 643778 : a = Fl_mul(a, i, p);
3466 143736 : v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
3467 : }
3468 57255 : return v;
3469 : }
3470 :
3471 : GEN
3472 186 : factorial_Fp(long n, GEN p)
3473 : {
3474 186 : pari_sp av = avma;
3475 : long k;
3476 186 : GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
3477 186 : for (k = 1;; k++)
3478 456 : {
3479 642 : long m = n >> (k-1), l, i;
3480 642 : GEN a = gen_1;
3481 642 : if (m <= 2) break;
3482 456 : l = (1 + (n >> k)) | 1;
3483 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3484 1690 : for (i=l; i<=m; i+=2)
3485 1234 : a = Fp_mulu(a, i, p);
3486 456 : v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
3487 456 : v = gc_INT(av, v);
3488 : }
3489 186 : return v;
3490 : }
3491 :
3492 : /*******************************************************************/
3493 : /** LUCAS & FIBONACCI **/
3494 : /*******************************************************************/
3495 : static void
3496 56 : lucas(ulong n, GEN *a, GEN *b)
3497 : {
3498 : GEN z, t, zt;
3499 56 : if (!n) { *a = gen_2; *b = gen_1; return; }
3500 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
3501 49 : switch(n & 3) {
3502 14 : case 0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
3503 14 : case 1: *a = subiu(zt,1); *b = addiu(sqri(t),2); break;
3504 7 : case 2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
3505 14 : case 3: *a = addiu(zt,1); *b = subiu(sqri(t),2);
3506 : }
3507 : }
3508 :
3509 : GEN
3510 7 : fibo(long n)
3511 : {
3512 7 : pari_sp av = avma;
3513 : GEN a, b;
3514 7 : if (!n) return gen_0;
3515 7 : lucas((ulong)(labs(n)-1), &a, &b);
3516 7 : a = diviuexact(addii(shifti(a,1),b), 5);
3517 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
3518 7 : return gc_INT(av, a);
3519 : }
3520 :
3521 : /*******************************************************************/
3522 : /* CONTINUED FRACTIONS */
3523 : /*******************************************************************/
3524 : static GEN
3525 3137064 : icopy_lg(GEN x, long l)
3526 : {
3527 3137064 : long lx = lgefint(x);
3528 : GEN y;
3529 :
3530 3137064 : if (lx >= l) return icopy(x);
3531 49 : y = cgeti(l); affii(x, y); return y;
3532 : }
3533 :
3534 : /* continued fraction of a/b. If y != NULL, stop when partial quotients
3535 : * differ from y */
3536 : static GEN
3537 3137414 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
3538 : {
3539 : GEN z, c;
3540 3137414 : ulong i, l, ly = lgefint(b);
3541 :
3542 : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
3543 3137414 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
3544 3137414 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
3545 3137414 : if (l > LGBITS) l = LGBITS;
3546 :
3547 3137414 : z = cgetg(l,t_VEC);
3548 3137414 : l--;
3549 3137414 : if (y) {
3550 350 : pari_sp av = avma;
3551 350 : if (l >= (ulong)lg(y)) l = lg(y)-1;
3552 25209 : for (i = 1; i <= l; i++)
3553 : {
3554 24985 : GEN q = gel(y,i);
3555 24985 : gel(z,i) = q;
3556 24985 : c = b; if (!gequal1(q)) c = mulii(q, b);
3557 24985 : c = subii(a, c);
3558 24985 : if (signe(c) < 0)
3559 : { /* partial quotient too large */
3560 96 : c = addii(c, b);
3561 96 : if (signe(c) >= 0) i++; /* by 1 */
3562 96 : break;
3563 : }
3564 24889 : if (cmpii(c, b) >= 0)
3565 : { /* partial quotient too small */
3566 30 : c = subii(c, b);
3567 30 : if (cmpii(c, b) < 0) {
3568 : /* by 1. If next quotient is 1 in y, add 1 */
3569 12 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
3570 12 : i++;
3571 : }
3572 30 : break;
3573 : }
3574 24859 : if ((i & 0xff) == 0) (void)gc_all(av, 2, &b, &c);
3575 24859 : a = b; b = c;
3576 : }
3577 : } else {
3578 3137064 : a = icopy_lg(a, ly);
3579 3137064 : b = icopy(b);
3580 24524443 : for (i = 1; i <= l; i++)
3581 : {
3582 24524125 : gel(z,i) = truedvmdii(a,b,&c);
3583 24524125 : if (c == gen_0) { i++; break; }
3584 21387379 : affii(c, a); cgiv(c); c = a;
3585 21387379 : a = b; b = c;
3586 : }
3587 : }
3588 3137414 : i--;
3589 3137414 : if (i > 1 && gequal1(gel(z,i)))
3590 : {
3591 101 : cgiv(gel(z,i)); --i;
3592 101 : gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
3593 : }
3594 3137414 : setlg(z,i+1); return z;
3595 : }
3596 :
3597 : static GEN
3598 0 : sersfcont(GEN a, GEN b, long k)
3599 : {
3600 0 : long i, l = typ(a) == t_POL? lg(a): 3;
3601 : GEN y, c;
3602 0 : if (lg(b) > l) l = lg(b);
3603 0 : if (k > 0 && l > k+1) l = k+1;
3604 0 : y = cgetg(l,t_VEC);
3605 0 : for (i=1; i<l; i++)
3606 : {
3607 0 : gel(y,i) = poldivrem(a,b,&c);
3608 0 : if (gequal0(c)) { i++; break; }
3609 0 : a = b; b = c;
3610 : }
3611 0 : setlg(y, i); return y;
3612 : }
3613 : static GEN
3614 7 : quadsfcontbound(GEN a, long k)
3615 : {
3616 7 : pari_sp av = avma;
3617 7 : long i, l = k+1;
3618 7 : GEN y = cgetg(l,t_VEC);
3619 147 : for (i=1; i<l; i++)
3620 : {
3621 140 : GEN c = gfloor(a);
3622 140 : gel(y,i) = c;
3623 140 : a = ginv(gsub(a,c));
3624 : }
3625 7 : return gc_GEN(av, y);
3626 : }
3627 :
3628 : static int
3629 21 : quad_isreduced(GEN x)
3630 : {
3631 21 : GEN c = conj_i(x);
3632 21 : return gcmp(x, gen_1) > 0 && gcmp(c,gen_0) < 0 && gcmp(c,gen_m1) > 0;
3633 : }
3634 :
3635 : static GEN
3636 14 : quadsfcont(GEN a)
3637 : {
3638 14 : pari_sp av = avma;
3639 14 : GEN a0 = NULL, V, W;
3640 14 : long i, l = 16;
3641 14 : V = cgetg(l+1, t_VEC);
3642 14 : for (i = 1;;)
3643 7 : {
3644 : GEN c;
3645 21 : if (quad_isreduced(a))
3646 14 : break;
3647 7 : c = gfloor(a);
3648 7 : gel(V,i++) = c;
3649 7 : a = ginv(gsub(a, c));
3650 7 : if (i==l+1)
3651 : {
3652 0 : l *= 2; V = vec_lengthen(V, l);
3653 : }
3654 : }
3655 14 : setlg(V, i);
3656 14 : l = 16; a0 = a;
3657 14 : W = cgetg(l+1, t_VEC);
3658 14 : for (i = 1;; i++)
3659 63 : {
3660 77 : GEN c = gfloor(a);
3661 77 : gel(W,i) = c;
3662 77 : a = ginv(gsub(a, c));
3663 77 : if (gequal(a, a0))
3664 14 : break;
3665 63 : if (i==l)
3666 : {
3667 0 : l *= 2; W = vec_lengthen(W, l);
3668 : }
3669 : }
3670 14 : setlg(W, i+1); return gc_GEN(av, mkvec2(V,W));
3671 : }
3672 :
3673 : GEN
3674 3142426 : gboundcf(GEN x, long k)
3675 : {
3676 : pari_sp av;
3677 3142426 : long tx = typ(x), e;
3678 : GEN y, a, b, c;
3679 :
3680 3142426 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
3681 3142419 : if (is_scalar_t(tx))
3682 : {
3683 3142405 : if (gequal0(x)) return mkvec(gen_0);
3684 3142286 : switch(tx)
3685 : {
3686 5194 : case t_INT: return mkveccopy(x);
3687 357 : case t_REAL:
3688 357 : av = avma;
3689 357 : c = mantissa_real(x,&e);
3690 357 : if (e < 0) pari_err_PREC("gboundcf");
3691 350 : y = int2n(e);
3692 350 : a = Qsfcont(c,y, NULL, k);
3693 350 : b = addsi(signe(x), c);
3694 350 : return gc_GEN(av, Qsfcont(b,y, a, k));
3695 :
3696 3136714 : case t_FRAC:
3697 3136714 : av = avma;
3698 3136714 : return gc_upto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
3699 21 : case t_QUAD:
3700 21 : if (signe(quad_disc(x)) <= 0) pari_err_DOMAIN("contfrac","x.disc","<",gen_0,x);
3701 21 : return k ? quadsfcontbound(x, k): quadsfcont(x);
3702 : }
3703 0 : pari_err_TYPE("gboundcf",x);
3704 : }
3705 :
3706 14 : switch(tx)
3707 : {
3708 14 : case t_QFB:
3709 14 : if (signe(qfb_disc(x)) <= 0) pari_err_DOMAIN("contfrac","x.disc","<",gen_0,x);
3710 14 : return k ? qfr_boundcf(x,k): qfr_cf(x);
3711 0 : case t_POL: return mkveccopy(x);
3712 0 : case t_SER:
3713 0 : av = avma;
3714 0 : return gc_upto(av, gboundcf(ser2rfrac_i(x), k));
3715 0 : case t_RFRAC:
3716 0 : av = avma;
3717 0 : return gc_GEN(av, sersfcont(gel(x,1), gel(x,2), k));
3718 : }
3719 0 : pari_err_TYPE("gboundcf",x);
3720 : return NULL; /* LCOV_EXCL_LINE */
3721 : }
3722 :
3723 : static GEN
3724 14 : sfcont2(GEN b, GEN x, long k)
3725 : {
3726 14 : pari_sp av = avma;
3727 14 : long lb = lg(b), tx = typ(x), i;
3728 : GEN y,p1;
3729 :
3730 14 : if (k)
3731 : {
3732 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
3733 0 : lb = k+1;
3734 : }
3735 7 : y = cgetg(lb,t_VEC);
3736 7 : if (lb==1) return y;
3737 7 : if (is_scalar_t(tx))
3738 : {
3739 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
3740 : }
3741 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
3742 :
3743 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
3744 7 : for (i = 1;;)
3745 : {
3746 35 : if (tx == t_REAL)
3747 : {
3748 35 : long e = expo(x);
3749 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
3750 35 : gel(y,i) = floorr(x);
3751 35 : p1 = subri(x, gel(y,i));
3752 : }
3753 : else
3754 : {
3755 0 : gel(y,i) = gfloor(x);
3756 0 : p1 = gsub(x, gel(y,i));
3757 : }
3758 35 : if (++i >= lb) break;
3759 28 : if (gequal0(p1)) break;
3760 28 : x = gdiv(gel(b,i),p1);
3761 : }
3762 7 : setlg(y,i);
3763 7 : return gc_GEN(av,y);
3764 : }
3765 :
3766 : GEN
3767 126 : gcf(GEN x) { return gboundcf(x,0); }
3768 : GEN
3769 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
3770 : GEN
3771 84 : contfrac0(GEN x, GEN b, long nmax)
3772 : {
3773 : long tb;
3774 :
3775 84 : if (!b) return gboundcf(x,nmax);
3776 42 : tb = typ(b);
3777 42 : if (tb == t_INT) return gboundcf(x,itos(b));
3778 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
3779 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
3780 14 : return sfcont2(b,x,nmax);
3781 : }
3782 :
3783 : GEN
3784 266 : contfracpnqn(GEN x, long n)
3785 : {
3786 266 : pari_sp av = avma;
3787 266 : long i, lx = lg(x);
3788 : GEN M,A,B, p0,p1, q0,q1;
3789 :
3790 266 : if (lx == 1)
3791 : {
3792 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
3793 21 : if (n >= 0) return cgetg(1,t_MAT);
3794 7 : return matid(2);
3795 : }
3796 238 : switch(typ(x))
3797 : {
3798 196 : case t_VEC: case t_COL: A = x; B = NULL; break;
3799 42 : case t_MAT:
3800 42 : switch(lgcols(x))
3801 : {
3802 0 : case 2: A = row(x,1); B = NULL; break;
3803 35 : case 3: A = row(x,2); B = row(x,1); break;
3804 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
3805 : return NULL; /*LCOV_EXCL_LINE*/
3806 : }
3807 35 : break;
3808 0 : default: pari_err_TYPE("pnqn",x);
3809 : return NULL; /*LCOV_EXCL_LINE*/
3810 : }
3811 231 : p1 = gel(A,1);
3812 231 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
3813 231 : if (n >= 0)
3814 : {
3815 196 : lx = minss(lx, n+2);
3816 196 : if (lx == 2) return gc_GEN(av, mkmat(mkcol2(p1,q1)));
3817 : }
3818 35 : else if (lx == 2)
3819 7 : return gc_GEN(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
3820 : /* lx >= 3 */
3821 119 : p0 = gen_1;
3822 119 : q0 = gen_0; /* p[-1], q[-1] */
3823 119 : M = cgetg(lx, t_MAT);
3824 119 : gel(M,1) = mkcol2(p1,q1);
3825 399 : for (i=2; i<lx; i++)
3826 : {
3827 280 : GEN a = gel(A,i), p2,q2;
3828 280 : if (B) {
3829 84 : GEN b = gel(B,i);
3830 84 : p0 = gmul(b,p0);
3831 84 : q0 = gmul(b,q0);
3832 : }
3833 280 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
3834 280 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
3835 280 : gel(M,i) = mkcol2(p1,q1);
3836 : }
3837 119 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
3838 119 : return gc_GEN(av, M);
3839 : }
3840 : GEN
3841 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
3842 : /* x = [a0, ..., an] from gboundcf, n >= 0;
3843 : * return [[p0, ..., pn], [q0,...,qn]] */
3844 : GEN
3845 894831 : ZV_allpnqn(GEN x)
3846 : {
3847 894831 : long i, lx = lg(x);
3848 894831 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
3849 :
3850 894831 : gel(v,1) = P = cgetg(lx, t_VEC);
3851 894831 : gel(v,2) = Q = cgetg(lx, t_VEC);
3852 894831 : p0 = gen_1; q0 = gen_0;
3853 894831 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
3854 3106250 : for (i=2; i<lx; i++)
3855 : {
3856 2211419 : GEN a = gel(x,i);
3857 2211419 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
3858 2211419 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
3859 : }
3860 894831 : return v;
3861 : }
3862 :
3863 : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
3864 : static GEN
3865 42 : mod_to_frac(GEN x, GEN N, GEN B)
3866 : {
3867 : GEN a, b, A;
3868 42 : if (B) A = divii(shifti(N, -1), B);
3869 : else
3870 : {
3871 14 : A = sqrti(shifti(N, -1));
3872 14 : B = A;
3873 : }
3874 42 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
3875 28 : return equali1(b)? a: mkfrac(a,b);
3876 : }
3877 :
3878 : static GEN
3879 112 : mod_to_rfrac(GEN x, GEN N, long B)
3880 : {
3881 : GEN a, b;
3882 112 : long A, d = degpol(N);
3883 112 : if (B >= 0) A = d-1 - B;
3884 : else
3885 : {
3886 42 : B = d >> 1;
3887 42 : A = odd(d)? B : B-1;
3888 : }
3889 112 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
3890 112 : if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
3891 91 : return gdiv(a,b);
3892 : }
3893 :
3894 : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
3895 : * of the continued fraction of x with b <= k maximal */
3896 : static GEN
3897 7 : bestappr_frac(GEN x, GEN k)
3898 : {
3899 : pari_sp av;
3900 : GEN p0, p1, p, q0, q1, q, a, y;
3901 :
3902 7 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
3903 0 : av = avma; y = x;
3904 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
3905 0 : q1 = gen_0; q0 = gen_1;
3906 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
3907 : for(;;)
3908 : {
3909 0 : x = ginv(x); /* > 1 */
3910 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
3911 0 : if (cmpii(a,k) > 0)
3912 : { /* next partial quotient will overflow limits */
3913 : GEN n, d;
3914 0 : a = divii(subii(k, q1), q0);
3915 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3916 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3917 : /* compare |y-p0/q0|, |y-p1/q1| */
3918 0 : n = gel(y,1);
3919 0 : d = gel(y,2);
3920 0 : if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
3921 : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
3922 0 : { p1 = p0; q1 = q0; }
3923 0 : break;
3924 : }
3925 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3926 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3927 :
3928 0 : if (cmpii(q0,k) > 0) break;
3929 0 : x = gsub(x,a); /* 0 <= x < 1 */
3930 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
3931 :
3932 : }
3933 0 : return gc_upto(av, gdiv(p1,q1));
3934 : }
3935 : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
3936 : * of the continued fraction of x with b <= k maximal */
3937 : static GEN
3938 1426091 : bestappr_real(GEN x, GEN k)
3939 : {
3940 1426091 : pari_sp av = avma;
3941 1426091 : GEN kr, p0, p1, p, q0, q1, q, a, y = x;
3942 :
3943 1426091 : p1 = gen_1; a = p0 = floorr(x);
3944 1426091 : q1 = gen_0; q0 = gen_1;
3945 1426091 : x = subri(x,a); /* 0 <= x < 1 */
3946 1426091 : if (!signe(x)) { cgiv(x); return a; }
3947 1306080 : kr = itor(k, realprec(x));
3948 : for(;;)
3949 9147418 : {
3950 : long d;
3951 10453498 : x = invr(x); /* > 1 */
3952 10453498 : if (cmprr(x,kr) > 0)
3953 : { /* next partial quotient will overflow limits */
3954 1102892 : a = divii(subii(k, q1), q0);
3955 1102892 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3956 1102892 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3957 : /* compare |y-p0/q0|, |y-p1/q1| */
3958 1102892 : if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
3959 : mulir(q0, subri(mulir(q1,y), p1))) < 0)
3960 128047 : { p1 = p0; q1 = q0; }
3961 1102892 : break;
3962 : }
3963 9350606 : d = nbits2prec(expo(x) + 1);
3964 9350606 : if (d > realprec(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
3965 :
3966 9350073 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
3967 9350073 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3968 9350073 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3969 :
3970 9350073 : if (cmpii(q0,k) > 0) break;
3971 9163248 : x = subri(x,a); /* 0 <= x < 1 */
3972 9163248 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
3973 : }
3974 1306080 : if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
3975 1306080 : return gc_GEN(av, equali1(q1)? p1: mkfrac(p1,q1));
3976 : }
3977 :
3978 : /* k t_INT or NULL */
3979 : static GEN
3980 2443909 : bestappr_Q(GEN x, GEN k)
3981 : {
3982 2443909 : long lx, tx = typ(x), i;
3983 : GEN a, y;
3984 :
3985 2443909 : switch(tx)
3986 : {
3987 1260 : case t_INT: return icopy(x);
3988 7 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
3989 1687086 : case t_REAL:
3990 1687086 : if (!signe(x)) return gen_0;
3991 : /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
3992 1426091 : i = bit_prec(x); if (i <= expo(x)) return NULL;
3993 1426091 : return bestappr_real(x, k? k: int2n(i));
3994 :
3995 28 : case t_INTMOD: {
3996 28 : pari_sp av = avma;
3997 28 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
3998 21 : return gc_GEN(av, a);
3999 : }
4000 14 : case t_PADIC: {
4001 14 : pari_sp av = avma;
4002 14 : long v = valp(x);
4003 14 : a = mod_to_frac(padic_u(x), padic_pd(x), k); if (!a) return NULL;
4004 7 : if (v) a = gmul(a, powis(padic_p(x), v));
4005 7 : return gc_GEN(av, a);
4006 : }
4007 :
4008 5474 : case t_COMPLEX: {
4009 5474 : pari_sp av = avma;
4010 5474 : y = cgetg(3, t_COMPLEX);
4011 5474 : gel(y,2) = bestappr(gel(x,2), k);
4012 5474 : gel(y,1) = bestappr(gel(x,1), k);
4013 5474 : if (gequal0(gel(y,2))) return gc_upto(av, gel(y,1));
4014 112 : return y;
4015 : }
4016 0 : case t_SER:
4017 0 : if (ser_isexactzero(x)) return gcopy(x);
4018 : /* fall through */
4019 : case t_POLMOD: case t_POL: case t_RFRAC:
4020 : case t_VEC: case t_COL: case t_MAT:
4021 750040 : y = cgetg_copy(x, &lx);
4022 750054 : for(i = 1; i < lontyp[tx]; i++) y[i] = x[i];
4023 2907079 : for (; i < lx; i++)
4024 : {
4025 2157039 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
4026 2157039 : gel(y,i) = a;
4027 : }
4028 750040 : if (tx == t_POL) return normalizepol(y);
4029 750026 : if (tx == t_SER) return normalizeser(y);
4030 750026 : return y;
4031 : }
4032 0 : pari_err_TYPE("bestappr_Q",x);
4033 : return NULL; /* LCOV_EXCL_LINE */
4034 : }
4035 :
4036 : static GEN
4037 98 : bestappr_ser(GEN x, long B)
4038 : {
4039 98 : long dN, v = valser(x), lx = lg(x);
4040 : GEN t;
4041 98 : x = normalizepol(ser2pol_i(x, lx));
4042 98 : dN = lx-2;
4043 98 : if (v > 0)
4044 : {
4045 21 : x = RgX_shift_shallow(x, v);
4046 21 : dN += v;
4047 : }
4048 77 : else if (v < 0)
4049 : {
4050 14 : if (B >= 0) B = maxss(B+v, 0);
4051 : }
4052 98 : t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
4053 98 : if (!t) return NULL;
4054 77 : if (v < 0)
4055 : {
4056 : GEN a, b;
4057 : long vx;
4058 14 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
4059 : /* t_RFRAC */
4060 14 : vx = varn(x);
4061 14 : a = gel(t,1);
4062 14 : b = gel(t,2);
4063 14 : v -= RgX_valrem(b, &b);
4064 14 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
4065 14 : if (v < 0) b = RgX_shift_shallow(b, -v);
4066 0 : else if (v > 0) {
4067 0 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
4068 0 : a = RgX_shift_shallow(a, v);
4069 : }
4070 14 : t = mkrfraccopy(a, b);
4071 : }
4072 77 : return t;
4073 : }
4074 : static GEN
4075 42 : gc_empty(pari_sp av) { retgc_const(av, cgetg(1, t_VEC)); }
4076 : static GEN
4077 112 : _gc_upto(pari_sp av, GEN x) { return x? gc_upto(av, x): NULL; }
4078 :
4079 : static GEN bestappr_RgX(GEN x, long B);
4080 : /* B >= 0 or < 0 [omit condition on B].
4081 : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b ~ x */
4082 : static GEN
4083 119 : bestappr_RgX(GEN x, long B)
4084 : {
4085 : pari_sp av;
4086 119 : switch(typ(x))
4087 : {
4088 0 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
4089 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
4090 0 : return gcopy(x);
4091 14 : case t_RFRAC:
4092 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
4093 7 : av = avma; return _gc_upto(av, bestappr_ser(rfrac_to_ser_i(x, 2*B+1), B));
4094 14 : case t_POLMOD:
4095 14 : av = avma; return _gc_upto(av, mod_to_rfrac(gel(x,2), gel(x,1), B));
4096 91 : case t_SER:
4097 91 : av = avma; return _gc_upto(av, bestappr_ser(x, B));
4098 0 : case t_VEC: case t_COL: case t_MAT: {
4099 : long i, lx;
4100 0 : GEN y = cgetg_copy(x, &lx);
4101 0 : for (i = 1; i < lx; i++)
4102 : {
4103 0 : GEN t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
4104 0 : gel(y,i) = t;
4105 : }
4106 0 : return y;
4107 : }
4108 : }
4109 0 : pari_err_TYPE("bestappr_RgX",x);
4110 : return NULL; /* LCOV_EXCL_LINE */
4111 : }
4112 :
4113 : /* allow k = NULL: maximal accuracy */
4114 : GEN
4115 286870 : bestappr(GEN x, GEN k)
4116 : {
4117 286870 : pari_sp av = avma;
4118 286870 : if (k) { /* replace by floor(k) */
4119 286149 : switch(typ(k))
4120 : {
4121 214090 : case t_INT:
4122 214090 : break;
4123 72059 : case t_REAL: case t_FRAC:
4124 72059 : k = floor_safe(k); /* left on stack for efficiency */
4125 72059 : if (!signe(k)) k = gen_1;
4126 72059 : break;
4127 0 : default:
4128 0 : pari_err_TYPE("bestappr [bound type]", k);
4129 0 : break;
4130 : }
4131 : }
4132 286870 : x = bestappr_Q(x, k);
4133 286870 : return x? x: gc_empty(av);
4134 : }
4135 : GEN
4136 119 : bestapprPade(GEN x, long B)
4137 : {
4138 119 : pari_sp av = avma;
4139 119 : GEN t = bestappr_RgX(x, B);
4140 119 : return t? t: gc_empty(av);
4141 : }
4142 :
4143 : static GEN
4144 49 : serPade(GEN S, long p, long q)
4145 : {
4146 49 : pari_sp av = avma;
4147 49 : long va, v, t = typ(S);
4148 49 : if (t!=t_SER && t!=t_POL && t!=t_RFRAC) pari_err_TYPE("bestapprPade", S);
4149 49 : va = gvar(S); v = gvaluation(S, pol_x(va));
4150 49 : if (p < 0) pari_err_DOMAIN("bestapprPade", "p", "<", gen_0, stoi(p));
4151 49 : if (q < 0) pari_err_DOMAIN("bestapprPade", "q", "<", gen_0, stoi(q));
4152 49 : if (v == LONG_MAX) return gc_empty(av);
4153 42 : S = gadd(S, zeroser(va, p + q + 1 + v));
4154 42 : return gc_upto(av, bestapprPade(S, q));
4155 : }
4156 :
4157 : GEN
4158 126 : bestapprPade0(GEN x, long p, long q)
4159 : {
4160 77 : return (p >= 0 && q >= 0)? serPade(x, p, q)
4161 203 : : bestapprPade(x, p >= 0? p: q);
4162 : }
|
