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 132320 : 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 132319 : 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 2079 : is_gener_expo(GEN p, GEN L0)
39 : {
40 2079 : GEN L, q = shifti(p,-1);
41 : long i, l;
42 2079 : if (L0) {
43 304 : l = lg(L0);
44 304 : L = cgetg(l, t_VEC);
45 : } else {
46 1775 : L0 = L = odd_prime_divisors(q);
47 1775 : l = lg(L);
48 : }
49 9217 : for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
50 2079 : return L;
51 : }
52 : static GEN
53 204536 : u_is_gener_expo(ulong p, GEN L0)
54 : {
55 204536 : const ulong q = p >> 1;
56 : long i;
57 : GEN L;
58 204536 : if (!L0) L0 = u_odd_prime_divisors(q);
59 204533 : L = cgetg_copy(L0,&i);
60 417954 : while (--i) L[i] = q / uel(L0,i);
61 204532 : return L;
62 : }
63 :
64 : int
65 578933 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
66 : {
67 : long i;
68 578933 : if (krouu(x, p) >= 0) return 0;
69 525777 : for (i=lg(L)-1; i; i--)
70 : {
71 316516 : ulong t = Fl_powu(x, uel(L,i), p);
72 316513 : if (t == p_1 || t == 1) return 0;
73 : }
74 209261 : return 1;
75 : }
76 : /* assume p prime */
77 : ulong
78 493895 : pgener_Fl_local(ulong p, GEN L0)
79 : {
80 493895 : const pari_sp av = avma;
81 493895 : const ulong p_1 = p-1;
82 : long x;
83 : GEN L;
84 493895 : if (p <= 19) switch(p)
85 : { /* quick trivial cases */
86 21 : case 2: return 1;
87 60564 : case 7:
88 60564 : case 17: return 3;
89 228837 : default: return 2;
90 : }
91 204473 : L = u_is_gener_expo(p,L0);
92 204501 : for (x = 2;; x++)
93 571093 : if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
94 : }
95 : ulong
96 201834 : 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 8472 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
102 : {
103 8472 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
104 8472 : if (t >= 0) return 0;
105 16582 : for (i = lg(L)-1; i; i--)
106 : {
107 11904 : GEN t = Fp_pow(x, gel(L,i), p);
108 11904 : if (equalii(t, p_1) || equali1(t)) return 0;
109 : }
110 4678 : return 1;
111 : }
112 :
113 : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
114 : GEN
115 354672 : pgener_Fp_local(GEN p, GEN L0)
116 : {
117 354672 : pari_sp av0 = avma;
118 : GEN x, p_1, L;
119 354672 : if (lgefint(p) == 3)
120 : {
121 : ulong z;
122 352600 : if (p[2] == 2) return gen_1;
123 257778 : if (L0) L0 = ZV_to_nv(L0);
124 257771 : z = pgener_Fl_local(uel(p,2), L0);
125 257819 : set_avma(av0); return utoipos(z);
126 : }
127 2072 : p_1 = subiu(p,1); L = is_gener_expo(p, L0);
128 2074 : x = utoipos(2);
129 4443 : for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
130 2074 : set_avma(av0); return utoipos(uel(x,2));
131 : }
132 :
133 : GEN
134 44016 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
135 :
136 : ulong
137 114214 : pgener_Zl(ulong p)
138 : {
139 114214 : 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 114214 : if (p == 40487) return 10;
142 : #ifndef LONG_IS_64BIT
143 16312 : return pgener_Fl(p);
144 : #else
145 97902 : 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 113687 : pgener_Zp(GEN p)
162 : {
163 113687 : 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 : set_avma(av); return utoipos(uel(x,2));
172 : }
173 : }
174 :
175 : static GEN
176 245 : gener_Zp(GEN q, GEN F)
177 : {
178 245 : GEN p = NULL;
179 245 : long e = 0;
180 245 : 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", "argument","=",F,F);
189 14 : e = itos(gel(E,i));
190 : }
191 14 : if (!p) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
192 : }
193 : else
194 231 : e = Z_isanypower(q, &p);
195 245 : return e > 1? pgener_Zp(p): pgener_Fp(q);
196 : }
197 :
198 : GEN
199 315 : znprimroot(GEN N)
200 : {
201 315 : pari_sp av = avma;
202 : GEN x, n, F;
203 :
204 315 : if ((F = check_arith_non0(N,"znprimroot")))
205 : {
206 14 : F = clean_Z_factor(F);
207 14 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
208 : }
209 308 : N = absi_shallow(N);
210 308 : if (abscmpiu(N, 4) <= 0) { set_avma(av); return mkintmodu(N[2]-1,N[2]); }
211 259 : switch(mod4(N))
212 : {
213 14 : case 0: /* N = 0 mod 4 */
214 14 : pari_err_DOMAIN("znprimroot", "argument","=",N,N);
215 0 : x = NULL; break;
216 21 : case 2: /* N = 2 mod 4 */
217 21 : n = shifti(N,-1); /* becomes odd */
218 21 : x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
219 21 : break;
220 224 : default: /* N odd */
221 224 : x = gener_Zp(N,F);
222 224 : break;
223 : }
224 245 : return gerepilecopy(av, mkintmod(x, N));
225 : }
226 :
227 : /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
228 : GEN
229 0 : rootsof1_Fp(GEN n, GEN p)
230 : {
231 0 : pari_sp av = avma;
232 0 : GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
233 0 : GEN z = pgener_Fp_local(p, L);
234 0 : z = Fp_pow(z, diviiexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
235 0 : return gerepileuptoint(av, z);
236 : }
237 :
238 : GEN
239 217 : rootsof1u_Fp(ulong n, GEN p)
240 : {
241 217 : pari_sp av = avma;
242 217 : GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
243 217 : z = pgener_Fp_local(p, Flv_to_ZV(L));
244 217 : z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
245 217 : return gerepileuptoint(av, z);
246 : }
247 :
248 : ulong
249 32115 : rootsof1_Fl(ulong n, ulong p)
250 : {
251 32115 : pari_sp av = avma;
252 32115 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
253 32115 : ulong z = pgener_Fl_local(p, L);
254 32115 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
255 32115 : return gc_ulong(av,z);
256 : }
257 :
258 : /*********************************************************************/
259 : /** INVERSE TOTIENT FUNCTION **/
260 : /*********************************************************************/
261 : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
262 : * primes. Return factor(N) */
263 : GEN
264 350651 : Z_factor_listP(GEN N, GEN L)
265 : {
266 350651 : long i, k, l = lg(L);
267 350651 : GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
268 1346688 : for (i = k = 1; i < l; i++)
269 : {
270 996037 : GEN p = gel(L,i);
271 996037 : long v = Z_pvalrem(N, p, &N);
272 996037 : if (v)
273 : {
274 792176 : gel(P,k) = p;
275 792176 : gel(E,k) = utoipos(v);
276 792176 : k++;
277 : }
278 : }
279 350651 : setlg(P, k);
280 350651 : setlg(E, k); return mkmat2(P,E);
281 : }
282 :
283 : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
284 : * L is a list of primes containing all prime divisors of n. */
285 : static long
286 621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
287 : {
288 621565 : pari_sp av = avma, av2;
289 : GEN k, D;
290 : long i, v;
291 621565 : if (m && mod2(n))
292 : {
293 270914 : if (!equali1(n)) return 0;
294 69986 : if (px) *px = gen_1;
295 69986 : return 1;
296 : }
297 350651 : D = divisors(Z_factor_listP(shifti(n, -1), L));
298 : /* loop through primes p > m, d = p-1 | n */
299 350651 : av2 = avma;
300 350651 : if (!m)
301 : { /* special case p = 2, d = 1 */
302 69986 : k = n;
303 69986 : for (v = 1;; v++) {
304 69986 : if (istotient_i(k, gen_2, L, px)) {
305 69986 : if (px) *px = shifti(*px, v);
306 69986 : return 1;
307 : }
308 0 : if (mod2(k)) break;
309 0 : k = shifti(k,-1);
310 : }
311 0 : set_avma(av2);
312 : }
313 1099462 : for (i = 1; i < lg(D); ++i)
314 : {
315 1001588 : GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
316 1001588 : if (m && cmpii(d, m) < 0) continue;
317 677782 : p = addiu(d, 1);
318 677782 : if (!isprime(p)) continue;
319 442064 : k = diviiexact(n, d);
320 481593 : for (v = 1;; v++) {
321 : GEN r;
322 481593 : if (istotient_i(k, p, L, px)) {
323 182791 : if (px) *px = mulii(*px, powiu(p, v));
324 182791 : return 1;
325 : }
326 298802 : k = dvmdii(k, p, &r);
327 298802 : if (r != gen_0) break;
328 : }
329 259273 : set_avma(av2);
330 : }
331 97874 : return gc_long(av,0);
332 : }
333 :
334 : /* find x such that phi(x) = n */
335 : long
336 70000 : istotient(GEN n, GEN *px)
337 : {
338 70000 : pari_sp av = avma;
339 70000 : if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
340 70000 : if (signe(n) < 1) return 0;
341 70000 : if (mod2(n))
342 : {
343 14 : if (!equali1(n)) return 0;
344 14 : if (px) *px = gen_1;
345 14 : return 1;
346 : }
347 69986 : if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
348 : {
349 69986 : if (!px) set_avma(av);
350 : else
351 69986 : *px = gerepileuptoint(av, *px);
352 69986 : return 1;
353 : }
354 0 : return gc_long(av,0);
355 : }
356 :
357 : /*********************************************************************/
358 : /** KRONECKER SYMBOL **/
359 : /*********************************************************************/
360 : /* t = 3,5 mod 8 ? (= 2 not a square mod t) */
361 : static int
362 298297947 : ome(long t)
363 : {
364 298297947 : switch(t & 7)
365 : {
366 169433605 : case 3:
367 169433605 : case 5: return 1;
368 128864342 : default: return 0;
369 : }
370 : }
371 : /* t a t_INT, is t = 3,5 mod 8 ? */
372 : static int
373 4098386 : gome(GEN t)
374 4098386 : { return signe(t)? ome( mod2BIL(t) ): 0; }
375 :
376 : /* assume y odd, return kronecker(x,y) * s */
377 : static long
378 212697992 : krouu_s(ulong x, ulong y, long s)
379 : {
380 212697992 : ulong x1 = x, y1 = y, z;
381 955723329 : while (x1)
382 : {
383 743045056 : long r = vals(x1);
384 743006678 : if (r)
385 : {
386 395412975 : if (odd(r) && ome(y1)) s = -s;
387 395431634 : x1 >>= r;
388 : }
389 743025337 : if (x1 & y1 & 2) s = -s;
390 743025337 : z = y1 % x1; y1 = x1; x1 = z;
391 : }
392 212678273 : return (y1 == 1)? s: 0;
393 : }
394 :
395 : long
396 9518524 : kronecker(GEN x, GEN y)
397 : {
398 9518524 : pari_sp av = avma;
399 9518524 : long s = 1, r;
400 : ulong xu;
401 :
402 9518524 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
403 9518524 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
404 9518524 : switch (signe(y))
405 : {
406 63 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
407 0 : case 0: return is_pm1(x);
408 : }
409 9518524 : r = vali(y);
410 9518526 : if (r)
411 : {
412 312749 : if (!mpodd(x)) return gc_long(av,0);
413 310887 : if (odd(r) && gome(x)) s = -s;
414 310887 : y = shifti(y,-r);
415 : }
416 9516663 : x = modii(x,y);
417 10517635 : while (lgefint(x) > 3) /* x < y */
418 : {
419 : GEN z;
420 1001025 : r = vali(x);
421 1001002 : if (r)
422 : {
423 549381 : if (odd(r) && gome(y)) s = -s;
424 549345 : x = shifti(x,-r);
425 : }
426 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
427 999655 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
428 1000008 : z = remii(y,x); y = x; x = z;
429 1001124 : if (gc_needed(av,2))
430 : {
431 0 : if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
432 0 : gerepileall(av, 2, &x, &y);
433 : }
434 : }
435 9516610 : xu = itou(x);
436 9516613 : if (!xu) return is_pm1(y)? s: 0;
437 9465082 : r = vals(xu);
438 9465095 : if (r)
439 : {
440 4705412 : if (odd(r) && gome(y)) s = -s;
441 4705412 : xu >>= r;
442 : }
443 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
444 9465095 : if (xu & mod2BIL(y) & 2) s = -s;
445 9465107 : return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
446 : }
447 :
448 : long
449 35973 : krois(GEN x, long y)
450 : {
451 : ulong yu;
452 35973 : long s = 1;
453 :
454 35973 : if (y <= 0)
455 : {
456 7 : if (y == 0) return is_pm1(x);
457 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
458 : }
459 : else
460 35966 : yu = (ulong)y;
461 35966 : if (!odd(yu))
462 : {
463 : long r;
464 14938 : if (!mpodd(x)) return 0;
465 11088 : r = vals(yu); yu >>= r;
466 11088 : if (odd(r) && gome(x)) s = -s;
467 : }
468 32116 : return krouu_s(umodiu(x, yu), yu, s);
469 : }
470 : /* assume y != 0 */
471 : long
472 27339887 : kroiu(GEN x, ulong y)
473 : {
474 : long r;
475 27339887 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
476 297091 : if (!mpodd(x)) return 0;
477 203131 : r = vals(y); y >>= r;
478 203134 : return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
479 : }
480 :
481 : /* assume y > 0, odd, return s * kronecker(x,y) */
482 : static long
483 163172 : krouodd(ulong x, GEN y, long s)
484 : {
485 : long r;
486 163172 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
487 19820 : if (!x) return 0; /* y != 1 */
488 19820 : r = vals(x);
489 19820 : if (r)
490 : {
491 10220 : if (odd(r) && gome(y)) s = -s;
492 10220 : x >>= r;
493 : }
494 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
495 19820 : if (x & mod2BIL(y) & 2) s = -s;
496 19820 : return krouu_s(umodiu(y,x), x, s);
497 : }
498 :
499 : long
500 154767 : krosi(long x, GEN y)
501 : {
502 154767 : const pari_sp av = avma;
503 154767 : long s = 1, r;
504 154767 : switch (signe(y))
505 : {
506 0 : case -1: y = negi(y); if (x < 0) s = -1; break;
507 0 : case 0: return (x==1 || x==-1);
508 : }
509 154767 : r = vali(y);
510 154767 : if (r)
511 : {
512 16884 : if (!odd(x)) return gc_long(av,0);
513 16884 : if (odd(r) && ome(x)) s = -s;
514 16884 : y = shifti(y,-r);
515 : }
516 154767 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
517 154767 : return gc_long(av, krouodd((ulong)x, y, s));
518 : }
519 :
520 : long
521 8405 : kroui(ulong x, GEN y)
522 : {
523 8405 : const pari_sp av = avma;
524 8405 : long s = 1, r;
525 8405 : switch (signe(y))
526 : {
527 0 : case -1: y = negi(y); break;
528 0 : case 0: return x==1UL;
529 : }
530 8405 : r = vali(y);
531 8405 : if (r)
532 : {
533 0 : if (!odd(x)) return gc_long(av,0);
534 0 : if (odd(r) && ome(x)) s = -s;
535 0 : y = shifti(y,-r);
536 : }
537 8405 : return gc_long(av, krouodd(x, y, s));
538 : }
539 :
540 : long
541 94633498 : kross(long x, long y)
542 : {
543 : ulong yu;
544 94633498 : long s = 1;
545 :
546 94633498 : if (y <= 0)
547 : {
548 68943 : if (y == 0) return (labs(x)==1);
549 68915 : yu = (ulong)-y; if (x < 0) s = -1;
550 : }
551 : else
552 94564555 : yu = (ulong)y;
553 94633470 : if (!odd(yu))
554 : {
555 : long r;
556 24360810 : if (!odd(x)) return 0;
557 17450540 : r = vals(yu); yu >>= r;
558 17450540 : if (odd(r) && ome(x)) s = -s;
559 : }
560 87723200 : x %= (long)yu; if (x < 0) x += yu;
561 87723200 : return krouu_s((ulong)x, yu, s);
562 : }
563 :
564 : long
565 88084008 : krouu(ulong x, ulong y)
566 : {
567 : long r;
568 88084008 : if (odd(y)) return krouu_s(x, y, 1);
569 5707 : if (!odd(x)) return 0;
570 6182 : r = vals(y); y >>= r;
571 6182 : return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
572 : }
573 :
574 : /*********************************************************************/
575 : /** HILBERT SYMBOL **/
576 : /*********************************************************************/
577 : /* x,y are t_INT or t_REAL */
578 : static long
579 7343 : mphilbertoo(GEN x, GEN y)
580 : {
581 7343 : long sx = signe(x), sy = signe(y);
582 7343 : if (!sx || !sy) return 0;
583 7343 : return (sx < 0 && sy < 0)? -1: 1;
584 : }
585 :
586 : long
587 134624 : hilbertii(GEN x, GEN y, GEN p)
588 : {
589 : pari_sp av;
590 : long oddvx, oddvy, z;
591 :
592 134624 : if (!p) return mphilbertoo(x,y);
593 127302 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
594 127302 : if (!signe(x) || !signe(y)) return 0;
595 127281 : av = avma;
596 127281 : oddvx = odd(Z_pvalrem(x,p,&x));
597 127281 : oddvy = odd(Z_pvalrem(y,p,&y));
598 : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
599 127281 : if (absequaliu(p, 2))
600 : {
601 12019 : z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
602 12019 : if (oddvx && gome(y)) z = -z;
603 12019 : if (oddvy && gome(x)) z = -z;
604 : }
605 : else
606 : {
607 115262 : z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
608 115262 : if (oddvx && kronecker(y,p) < 0) z = -z;
609 115262 : if (oddvy && kronecker(x,p) < 0) z = -z;
610 : }
611 127281 : return gc_long(av, z);
612 : }
613 :
614 : static void
615 196 : err_prec(void) { pari_err_PREC("hilbert"); }
616 : static void
617 161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
618 : static void
619 56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
620 :
621 : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
622 : * Return lift(x) provided it's p-adic accuracy is large enough to decide
623 : * hilbert()'s value [ problem at p = 2 ] */
624 : static GEN
625 420 : lift_intmod(GEN x, GEN *pp)
626 : {
627 420 : GEN p = *pp, N = gel(x,1);
628 420 : x = gel(x,2);
629 420 : if (!p)
630 : {
631 266 : *pp = p = N;
632 266 : switch(itos_or_0(p))
633 : {
634 126 : case 2:
635 126 : case 4: err_prec();
636 : }
637 140 : return x;
638 : }
639 154 : if (!signe(p)) err_oo(N);
640 112 : if (absequaliu(p,2))
641 42 : { if (vali(N) <= 2) err_prec(); }
642 : else
643 70 : { if (!dvdii(N,p)) err_p(N,p); }
644 28 : if (!signe(x)) err_prec();
645 21 : return x;
646 : }
647 : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
648 : * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
649 : * to decide hilbert()'s value [ problem at p = 2 ]*/
650 : static GEN
651 210 : lift_padic(GEN x, GEN *pp)
652 : {
653 210 : GEN p = *pp, q = gel(x,2), y = gel(x,4);
654 210 : if (!p) *pp = p = q;
655 147 : else if (!equalii(p,q)) err_p(p, q);
656 105 : if (absequaliu(p,2) && precp(x) <= 2) err_prec();
657 70 : if (!signe(y)) err_prec();
658 70 : return odd(valp(x))? mulii(p,y): y;
659 : }
660 :
661 : long
662 59857 : hilbert(GEN x, GEN y, GEN p)
663 : {
664 59857 : pari_sp av = avma;
665 59857 : long tx = typ(x), ty = typ(y);
666 :
667 59857 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
668 59857 : if (tx == t_REAL)
669 : {
670 77 : if (p && signe(p)) err_oo(p);
671 63 : switch (ty)
672 : {
673 7 : case t_INT:
674 7 : case t_REAL: return mphilbertoo(x,y);
675 0 : case t_FRAC: return mphilbertoo(x,gel(y,1));
676 56 : default: pari_err_TYPE2("hilbert",x,y);
677 : }
678 : }
679 59780 : if (ty == t_REAL)
680 : {
681 14 : if (p && signe(p)) err_oo(p);
682 14 : switch (tx)
683 : {
684 14 : case t_INT:
685 14 : case t_REAL: return mphilbertoo(x,y);
686 0 : case t_FRAC: return mphilbertoo(gel(x,1),y);
687 0 : default: pari_err_TYPE2("hilbert",x,y);
688 : }
689 : }
690 59766 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
691 59563 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
692 :
693 59507 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
694 59444 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
695 :
696 59367 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
697 59367 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
698 :
699 59367 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
700 59367 : if (p && !signe(p)) p = NULL;
701 59367 : return gc_long(av, hilbertii(x,y,p));
702 : }
703 :
704 : /*******************************************************************/
705 : /* SQUARE ROOT MODULO p */
706 : /*******************************************************************/
707 : static void
708 1743913 : checkp(ulong q, ulong p)
709 1743913 : { if (!q) pari_err_PRIME("Fl_nonsquare",utoipos(p)); }
710 : /* p = 1 (mod 4) prime, return the first quadratic nonresidue, a prime */
711 : static ulong
712 10368671 : nonsquare1_Fl(ulong p)
713 : {
714 : forprime_t S;
715 : ulong q;
716 10368671 : if ((p & 7UL) != 1) return 2UL;
717 3745626 : q = p % 3; if (q == 2) return 3UL;
718 1127340 : checkp(q, p);
719 1127412 : q = p % 5; if (q == 2 || q == 3) return 5UL;
720 396420 : checkp(q, p);
721 396420 : q = p % 7; if (q != 4 && q >= 3) return 7UL;
722 138826 : checkp(q, p);
723 : /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
724 138826 : u_forprime_init(&S, 11, 1967);
725 220074 : while ((q = u_forprime_next(&S)))
726 : {
727 220074 : if (krouu(q, p) < 0) return q;
728 81248 : checkp(q, p);
729 : }
730 0 : checkp(0, p);
731 : return 0; /*LCOV_EXCL_LINE*/
732 : }
733 : /* p > 2 a prime */
734 : ulong
735 7853 : nonsquare_Fl(ulong p)
736 7853 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
737 :
738 : ulong
739 155880 : Fl_2gener_pre(ulong p, ulong pi)
740 : {
741 155880 : ulong p1 = p-1;
742 155880 : long e = vals(p1);
743 155879 : if (e == 1) return p1;
744 58337 : return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
745 : }
746 :
747 : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
748 : ulong
749 30272523 : Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
750 : {
751 : long i, e, k;
752 : ulong p1, q, v, w;
753 :
754 30272523 : if (!a) return 0;
755 28930522 : p1 = p - 1; e = vals(p1);
756 28939656 : if (e == 0) /* p = 2 */
757 : {
758 420476 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
759 421973 : return ((a & 1) == 0)? 0: 1;
760 : }
761 28519180 : if (e == 1)
762 : {
763 18207514 : v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
764 18168092 : if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
765 18175483 : p1 = p - v; if (v > p1) v = p1;
766 18175483 : return v;
767 : }
768 10311666 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
769 10311666 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
770 10312036 : if (!p1) return 0;
771 10312036 : v = Fl_mul_pre(a, p1, p, pi);
772 10312124 : w = Fl_mul_pre(v, p1, p, pi);
773 10312140 : if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
774 18574622 : while (w != 1)
775 : { /* a*w = v^2, y primitive 2^e-th root of 1
776 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
777 8264460 : p1 = Fl_sqr_pre(w,p,pi);
778 14746911 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
779 8264507 : if (k == e) return ~0UL;
780 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
781 8262547 : p1 = y;
782 10856836 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
783 8262546 : y = Fl_sqr_pre(p1, p, pi); e = k;
784 8262556 : w = Fl_mul_pre(y, w, p, pi);
785 8262554 : v = Fl_mul_pre(v, p1, p, pi);
786 : }
787 10310162 : p1 = p - v; if (v > p1) v = p1;
788 10310162 : return v;
789 : }
790 :
791 : ulong
792 26769307 : Fl_sqrt(ulong a, ulong p)
793 26769307 : { ulong pi = get_Fl_red(p); return Fl_sqrt_pre_i(a, 0, p, pi); }
794 :
795 : ulong
796 3453453 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
797 3453453 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
798 :
799 : static ulong
800 48606 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
801 : {
802 : ulong x, y, m;
803 48606 : ulong le1 = upowuu(l, e-1);
804 48606 : for (x = 2; ; x++)
805 : {
806 76507 : y = Fl_powu_pre(x, r, p, pi);
807 76508 : if (y==1) continue;
808 59487 : m = Fl_powu_pre(y, le1, p, pi);
809 59486 : if (m != 1) break;
810 : }
811 48606 : *pt_m = m; return y;
812 : }
813 :
814 : /* solve x^l = a , l prime in G of order q.
815 : *
816 : * q = (l^e)*r, e >= 1, (r,l) = 1
817 : * y generates the l-Sylow of G
818 : * m = y^(l^(e-1)) != 1 */
819 : static ulong
820 117063 : Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
821 : {
822 : ulong p1, v, w, z, dl;
823 : ulong u2;
824 117063 : if (a==0) return a;
825 117063 : u2 = Fl_inv(l%r, r);
826 117064 : v = Fl_powu_pre(a, u2, p, pi);
827 117065 : w = Fl_powu_pre(v, l, p, pi);
828 117064 : w = Fl_mul_pre(w, Fl_inv(a, p), p, pi);
829 117047 : if (w==1) return v;
830 47241 : if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
831 67501 : while (w!=1)
832 : {
833 51795 : ulong k = 0;
834 51795 : p1 = w;
835 : do
836 : {
837 76829 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
838 76829 : if (++k == e) return ULONG_MAX;
839 45294 : } while (p1!=1);
840 20260 : dl = Fl_log_pre(z, m, l, p, pi);
841 20260 : dl = Fl_neg(dl, l);
842 20260 : p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
843 20260 : m = Fl_powu_pre(m, dl, p, pi);
844 20260 : e = k;
845 20260 : v = Fl_mul_pre(p1,v,p,pi);
846 20260 : y = Fl_powu_pre(p1,l,p,pi);
847 20260 : w = Fl_mul_pre(y,w,p,pi);
848 : }
849 15706 : return v;
850 : }
851 :
852 : static ulong
853 114713 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
854 : {
855 114713 : ulong r, e = u_lvalrem(p-1, l, &r);
856 114712 : return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
857 : }
858 :
859 : ulong
860 114713 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
861 114713 : { return Fl_sqrtl_i(a, l, p, pi, 0, 0); }
862 :
863 : ulong
864 0 : Fl_sqrtl(ulong a, ulong l, ulong p)
865 0 : { ulong pi = get_Fl_red(p); return Fl_sqrtl_i(a, l, p, pi, 0, 0); }
866 :
867 : ulong
868 85834 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
869 : {
870 85834 : ulong m, q = p-1, z;
871 85834 : ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
872 85834 : if (a==0)
873 : {
874 48139 : if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
875 48132 : if (zetan) *zetan = 1UL;
876 48132 : return 0;
877 : }
878 37695 : if (n==1)
879 : {
880 420 : if (zetan) *zetan = 1;
881 420 : return n < 0? Fl_inv(a,p): a;
882 : }
883 37275 : if (n==2)
884 : {
885 5887 : if (zetan) *zetan = p-1;
886 5887 : return Fl_sqrt_pre_i(a, 0, p, pi);
887 : }
888 31388 : if (a == 1 && !zetan) return a;
889 7969 : m = ugcd(nn, q);
890 7969 : z = 1;
891 7969 : if (m!=1)
892 : {
893 1396 : GEN F = factoru(m);
894 : long i, j, e;
895 : ulong r, zeta, y, l;
896 3118 : for (i = nbrows(F); i; i--)
897 : {
898 1785 : l = ucoeff(F,i,1);
899 1785 : j = ucoeff(F,i,2);
900 1785 : e = u_lvalrem(q,l, &r);
901 1785 : y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
902 1785 : if (zetan)
903 98 : z = Fl_mul_pre(z, Fl_powu_pre(y, upowuu(l,e-j), p, pi), p, pi);
904 1785 : if (a!=1)
905 : do
906 : {
907 2351 : a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
908 2337 : if (a==ULONG_MAX) return ULONG_MAX;
909 2288 : } while (--j);
910 : }
911 : }
912 7906 : if (m != nn)
913 : {
914 6594 : ulong qm = q/m, nm = nn/m;
915 6594 : a = Fl_powu_pre(a, Fl_inv(nm%qm, qm), p, pi);
916 : }
917 7906 : if (n < 0) a = Fl_inv(a, p);
918 7906 : if (zetan) *zetan = z;
919 7906 : return a;
920 : }
921 :
922 : ulong
923 85834 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
924 : {
925 85834 : ulong pi = get_Fl_red(p);
926 85834 : return Fl_sqrtn_pre(a, n, p, pi, zetan);
927 : }
928 :
929 : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
930 : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
931 : * and in average is about 2 or 2.5 times worse. But try both algorithms for
932 : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
933 : *
934 : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
935 : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
936 : * If (a|p)=1, then sqrt(a) is in F_p.
937 : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
938 :
939 : /* compute y^2, y = y[1] + y[2] X */
940 : static GEN
941 449 : sqrt_Cipolla_sqr(void *data, GEN y)
942 : {
943 449 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
944 449 : GEN u2 = sqri(u), v2 = sqri(v);
945 449 : v = subii(sqri(addii(v,u)), addii(u2,v2));
946 449 : u = addii(u2, mulii(v2,n));
947 449 : retmkvec2(modii(u,p), modii(v,p));
948 : }
949 : /* compute (t+X) y^2 */
950 : static GEN
951 23 : sqrt_Cipolla_msqr(void *data, GEN y)
952 : {
953 23 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2);
954 23 : ulong t = gel(data,4)[2];
955 23 : GEN d = addii(u, mului(t,v)), d2 = sqri(d);
956 23 : GEN b = remii(mulii(a,v), p);
957 23 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
958 23 : v = subii(d2, mulii(b,v));
959 23 : retmkvec2(modii(u,p), modii(v,p));
960 : }
961 : /* assume a reduced mod p [ otherwise correct but inefficient ] */
962 : static GEN
963 8 : sqrt_Cipolla(GEN a, GEN p)
964 : {
965 : pari_sp av;
966 : GEN u, v, n, y, pov2;
967 : ulong t;
968 :
969 8 : if (kronecker(a, p) < 0) return NULL;
970 8 : pov2 = shifti(p,-1); /* center to avoid multiplying by huge base*/
971 8 : if (cmpii(a,pov2) > 0) a = subii(a,p);
972 8 : av = avma;
973 41 : for (t=1; ; t++, set_avma(av))
974 : {
975 41 : n = subsi((long)(t*t), a);
976 41 : if (kronecker(n, p) < 0) break;
977 : }
978 :
979 : /* compute (t+X)^((p-1)/2) =: u+vX */
980 8 : u = utoipos(t);
981 8 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
982 : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
983 : /* Now u+vX = (t+X)^((p-1)/2); thus
984 : * (u+vX)(t+X) = sqrt(a) + 0 X
985 : * Whence,
986 : * sqrt(a) = (u+vt)t - v*a
987 : * 0 = (u+vt)
988 : * Thus a square root is v*a */
989 :
990 8 : v = Fp_mul(gel(y, 2), a, p);
991 8 : if (cmpii(v,pov2) > 0) v = subii(p,v);
992 8 : return v;
993 : }
994 :
995 : /* Return NULL if p is found to be composite */
996 : static GEN
997 3190 : Fp_2gener_all(long e, GEN p)
998 : {
999 : GEN y, m;
1000 : long k;
1001 3190 : GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
1002 3190 : if (e==0 && !equaliu(p,2)) return NULL;
1003 3190 : for (k=2; ; k++)
1004 7811 : {
1005 11001 : long i = krosi(k, p);
1006 11001 : if (i >= 0)
1007 : {
1008 7811 : if (i) continue;
1009 0 : return NULL;
1010 : }
1011 3190 : y = m = Fp_pow(utoi(k), q, p);
1012 10928 : for (i=1; i<e; i++)
1013 7738 : if (equali1(m = Fp_sqr(m, p))) break;
1014 3190 : if (i == e) break; /* success */
1015 : }
1016 3190 : return y;
1017 : }
1018 :
1019 : /* Return NULL if p is found to be composite */
1020 : GEN
1021 1120 : Fp_2gener(GEN p)
1022 1120 : { return Fp_2gener_all(vali(subis(p,1)),p); }
1023 :
1024 : /* smallest square root */
1025 : static GEN
1026 32996 : choose_sqrt(GEN v, GEN p)
1027 : {
1028 32996 : pari_sp av = avma;
1029 32996 : GEN q = subii(p,v);
1030 32995 : if (cmpii(v,q) > 0) v = q; else set_avma(av);
1031 32995 : return v;
1032 : }
1033 : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
1034 : GEN
1035 3334227 : Fp_sqrt_i(GEN a, GEN y, GEN p)
1036 : {
1037 3334227 : pari_sp av = avma;
1038 : long i, k, e;
1039 : GEN p1, q, v, w;
1040 :
1041 3334227 : if (lgefint(p) == 3)
1042 : {
1043 3301126 : ulong pp = uel(p,2), u = umodiu(a, pp);
1044 3301157 : if (!u) return gen_0;
1045 3005441 : u = Fl_sqrt(u, pp);
1046 3005501 : return (u == ~0UL)? NULL: utoipos(u);
1047 : }
1048 :
1049 33101 : a = modii(a, p); if (!signe(a)) return gc_const(av, gen_0);
1050 32984 : p1 = subiu(p,1); e = vali(p1);
1051 32998 : if (e <= 2)
1052 : { /* direct formulas more efficient */
1053 : pari_sp av2;
1054 26519 : if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
1055 26519 : if (e == 1)
1056 : {
1057 16011 : q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
1058 16003 : v = Fp_pow(a, q, p);
1059 : }
1060 : else
1061 : { /* Atkin's formula */
1062 10508 : GEN i, a2 = shifti(a,1);
1063 10506 : if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
1064 10510 : q = shifti(p1, -3); /* (p-5)/8 */
1065 10512 : v = Fp_pow(a2, q, p);
1066 10519 : i = Fp_mul(a2, Fp_sqr(v,p), p); /* i^2 = -1 */
1067 10519 : v = Fp_mul(a, Fp_mul(v, subiu(i,1), p), p);
1068 : }
1069 26545 : av2 = avma;
1070 : /* must check equality in case (a/p) = -1 or p not prime */
1071 26545 : e = equalii(Fp_sqr(v,p), a); set_avma(av2);
1072 26546 : return e? gerepileuptoint(av,choose_sqrt(v,p)): NULL;
1073 : }
1074 : /* On average, Cipolla is better than Tonelli/Shanks if and only if
1075 : * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
1076 6479 : if (e*(e-1) > 20 + 8 * expi(p))
1077 : {
1078 8 : v = sqrt_Cipolla(a,p); if (!v) return gc_NULL(av);
1079 8 : return gerepileuptoint(av,v);
1080 : }
1081 6472 : if (!y)
1082 : {
1083 2070 : y = Fp_2gener_all(e, p);
1084 2070 : if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
1085 : }
1086 6472 : q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
1087 6471 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
1088 6472 : v = Fp_mul(a, p1, p);
1089 6473 : w = Fp_mul(v, p1, p);
1090 15684 : while (!equali1(w))
1091 : { /* a*w = v^2, y primitive 2^e-th root of 1
1092 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1093 9212 : p1 = Fp_sqr(w,p);
1094 19896 : for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
1095 9204 : if (k == e) return gc_NULL(av); /* p composite or (a/p) != 1 */
1096 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1097 9204 : p1 = y;
1098 13380 : for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
1099 9204 : y = Fp_sqr(p1, p); e = k;
1100 9211 : w = Fp_mul(y, w, p);
1101 9212 : v = Fp_mul(v, p1, p);
1102 9211 : if (gc_needed(av,1))
1103 : {
1104 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
1105 0 : gerepileall(av,3, &y,&w,&v);
1106 : }
1107 : }
1108 6470 : return gerepileuptoint(av, choose_sqrt(v,p));
1109 : }
1110 :
1111 : GEN
1112 3314854 : Fp_sqrt(GEN a, GEN p)
1113 3314854 : { return Fp_sqrt_i(a, NULL, p); }
1114 :
1115 : /*********************************************************************/
1116 : /** GCD & BEZOUT **/
1117 : /*********************************************************************/
1118 :
1119 : GEN
1120 41157101 : lcmii(GEN x, GEN y)
1121 : {
1122 : pari_sp av;
1123 : GEN a, b;
1124 41157101 : if (!signe(x) || !signe(y)) return gen_0;
1125 41157192 : av = avma; a = gcdii(x,y);
1126 41155654 : if (absequalii(a,y)) { set_avma(av); return absi(x); }
1127 7805305 : if (!equali1(a)) y = diviiexact(y,a);
1128 7805343 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
1129 : }
1130 :
1131 : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
1132 : * set *pd = gcd(x,N) */
1133 : GEN
1134 6587193 : Fp_invgen(GEN x, GEN N, GEN *pd)
1135 : {
1136 : GEN d, d0, e, v;
1137 6587193 : if (lgefint(N) == 3)
1138 : {
1139 6026489 : ulong dd, NN = N[2], xx = umodiu(x,NN);
1140 6027123 : if (!xx) { *pd = N; return gen_0; }
1141 6027123 : xx = Fl_invgen(xx, NN, &dd);
1142 6027718 : *pd = utoi(dd); return utoi(xx);
1143 : }
1144 560704 : *pd = d = bezout(x, N, &v, NULL);
1145 560717 : if (equali1(d)) return v;
1146 : /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
1147 460833 : e = diviiexact(N,d);
1148 460833 : d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
1149 460833 : if (equali1(d0)) return v;
1150 289792 : if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
1151 289792 : return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
1152 : }
1153 :
1154 : /*********************************************************************/
1155 : /** CHINESE REMAINDERS **/
1156 : /*********************************************************************/
1157 :
1158 : /* Chinese Remainder Theorem. x and y must have the same type (integermod,
1159 : * polymod, or polynomial/vector/matrix recursively constructed with these
1160 : * as coefficients). Creates (with the same type) a z in the same residue
1161 : * class as x and the same residue class as y, if it is possible.
1162 : *
1163 : * We also allow (during recursion) two identical objects even if they are
1164 : * not integermod or polymod. For example:
1165 : *
1166 : * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
1167 : * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
1168 : * ? chinese(x, y)
1169 : * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
1170 :
1171 : static GEN
1172 2127668 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
1173 : {
1174 2127668 : GEN z = gassoc_proto(f,x,NULL);
1175 2127661 : if (z == gen_1) retmkintmod(gen_0,gen_1);
1176 2127626 : return z;
1177 : }
1178 :
1179 : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
1180 : * call chinese: makes Mod(0,1) a better "neutral" element */
1181 : static GEN
1182 21 : chinese_intpol(GEN x,GEN y)
1183 : {
1184 21 : pari_sp av = avma;
1185 21 : GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
1186 21 : return gerepileupto(av, chinese(z, y));
1187 : }
1188 :
1189 : GEN
1190 2275 : chinese1(GEN x) { return gen_chinese(x,chinese); }
1191 :
1192 : GEN
1193 5327 : chinese(GEN x, GEN y)
1194 : {
1195 : pari_sp av;
1196 5327 : long tx = typ(x), ty;
1197 : GEN z,p1,p2,d,u,v;
1198 :
1199 5327 : if (!y) return chinese1(x);
1200 5278 : if (gidentical(x,y)) return gcopy(x);
1201 5271 : ty = typ(y);
1202 5271 : if (tx == ty) switch(tx)
1203 : {
1204 3759 : case t_POLMOD:
1205 : {
1206 3759 : GEN A = gel(x,1), B = gel(y,1);
1207 3759 : GEN a = gel(x,2), b = gel(y,2);
1208 3759 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
1209 3759 : if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
1210 3759 : av = avma;
1211 3759 : d = RgX_extgcd(A,B,&u,&v);
1212 3759 : p2 = gsub(b, a);
1213 3759 : if (!gequal0(gmod(p2, d))) break;
1214 3759 : p1 = gdiv(A,d);
1215 3759 : p2 = gadd(a, gmul(gmul(u,p1), p2));
1216 :
1217 3759 : z = cgetg(3, t_POLMOD);
1218 3759 : gel(z,1) = gmul(p1,B);
1219 3759 : gel(z,2) = gmod(p2,gel(z,1));
1220 3759 : return gerepileupto(av, z);
1221 : }
1222 1477 : case t_INTMOD:
1223 : {
1224 1477 : GEN A = gel(x,1), B = gel(y,1);
1225 1477 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
1226 1477 : z = cgetg(3,t_INTMOD);
1227 1477 : Z_chinese_pre(A, B, &C, &U, &d);
1228 1477 : c = Z_chinese_post(a, b, C, U, d);
1229 1477 : if (!c) pari_err_OP("chinese", x,y);
1230 1477 : set_avma((pari_sp)z);
1231 1477 : gel(z,1) = icopy(C);
1232 1477 : gel(z,2) = icopy(c); return z;
1233 : }
1234 7 : case t_POL:
1235 : {
1236 7 : long i, lx = lg(x), ly = lg(y);
1237 7 : if (varn(x) != varn(y)) break;
1238 7 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
1239 7 : z = cgetg(lx, t_POL); z[1] = x[1];
1240 21 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
1241 14 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
1242 7 : return z;
1243 : }
1244 :
1245 7 : case t_VEC: case t_COL: case t_MAT:
1246 : {
1247 : long i, lx;
1248 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
1249 21 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
1250 7 : return z;
1251 : }
1252 : }
1253 21 : if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
1254 7 : if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
1255 0 : pari_err_OP("chinese",x,y);
1256 : return NULL; /* LCOV_EXCL_LINE */
1257 : }
1258 :
1259 : /* init chinese(Mod(.,A), Mod(.,B)) */
1260 : void
1261 266158 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
1262 : {
1263 266158 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
1264 266162 : GEN t = diviiexact(A,d);
1265 266153 : *pU = mulii(u, t);
1266 266146 : *pC = mulii(t, B);
1267 266145 : if (pd) *pd = d;
1268 266145 : }
1269 : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
1270 : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
1271 : * If d not NULL, check whether a = b mod d. */
1272 : GEN
1273 2730476 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
1274 : {
1275 : GEN b_a;
1276 2730476 : if (!signe(a))
1277 : {
1278 785220 : if (d && !dvdii(b, d)) return NULL;
1279 785220 : return Fp_mul(b, U, C);
1280 : }
1281 1945256 : b_a = subii(b,a);
1282 1945256 : if (d && !dvdii(b_a, d)) return NULL;
1283 1945256 : return modii(addii(a, mulii(U, b_a)), C);
1284 : }
1285 : static ulong
1286 2114536 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
1287 : {
1288 2114536 : if (!a) return Fl_mul(b, U, C);
1289 2114536 : return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
1290 : }
1291 :
1292 : GEN
1293 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
1294 : {
1295 2142 : pari_sp av = avma;
1296 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
1297 2142 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
1298 : }
1299 : GEN
1300 262483 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
1301 : {
1302 262483 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
1303 262469 : return Z_chinese_post(a,b, *pC, U, NULL);
1304 : }
1305 :
1306 : /* return lift(chinese(a mod A, b mod B))
1307 : * assume(A,B)=1, a,b,A,B integers and C = A*B */
1308 : GEN
1309 291052 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
1310 : {
1311 291052 : pari_sp av = avma;
1312 291052 : GEN U = mulii(Fp_inv(A,B), A);
1313 291052 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
1314 : }
1315 : ulong
1316 2114521 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
1317 2114521 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
1318 :
1319 : /* chinese1 for coprime moduli in Z */
1320 : static GEN
1321 2173011 : chinese1_coprime_Z_aux(GEN x, GEN y)
1322 : {
1323 2173011 : GEN z = cgetg(3, t_INTMOD);
1324 2173011 : GEN A = gel(x,1), a = gel(x, 2);
1325 2173011 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
1326 2173011 : pari_sp av = avma;
1327 2173011 : GEN U = mulii(Fp_inv(A,B), A);
1328 2173011 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
1329 2173011 : gel(z,1) = C; return z;
1330 : }
1331 : GEN
1332 2125393 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
1333 :
1334 : /*********************************************************************/
1335 : /** MODULAR EXPONENTIATION **/
1336 : /*********************************************************************/
1337 : /* xa ZV or nv */
1338 : GEN
1339 1998471 : ZV_producttree(GEN xa)
1340 : {
1341 1998471 : long n = lg(xa)-1;
1342 1998471 : long m = n==1 ? 1: expu(n-1)+1;
1343 1998470 : GEN T = cgetg(m+1, t_VEC), t;
1344 : long i, j, k;
1345 1998470 : t = cgetg(((n+1)>>1)+1, t_VEC);
1346 1998461 : if (typ(xa)==t_VECSMALL)
1347 : {
1348 2581316 : for (j=1, k=1; k<n; j++, k+=2)
1349 1703417 : gel(t, j) = muluu(xa[k], xa[k+1]);
1350 877899 : if (k==n) gel(t, j) = utoi(xa[k]);
1351 : } else {
1352 2311315 : for (j=1, k=1; k<n; j++, k+=2)
1353 1190781 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
1354 1120534 : if (k==n) gel(t, j) = icopy(gel(xa,k));
1355 : }
1356 1998437 : gel(T,1) = t;
1357 3269065 : for (i=2; i<=m; i++)
1358 : {
1359 1270618 : GEN u = gel(T, i-1);
1360 1270618 : long n = lg(u)-1;
1361 1270618 : t = cgetg(((n+1)>>1)+1, t_VEC);
1362 2830025 : for (j=1, k=1; k<n; j++, k+=2)
1363 1559397 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
1364 1270628 : if (k==n) gel(t, j) = gel(u, k);
1365 1270628 : gel(T, i) = t;
1366 : }
1367 1998447 : return T;
1368 : }
1369 :
1370 : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
1371 : GEN
1372 53043330 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
1373 : {
1374 : long i,j,k;
1375 53043330 : long m = lg(T)-1, n = lg(P)-1;
1376 : GEN t;
1377 53043330 : GEN Tp = cgetg(m+1, t_VEC);
1378 52980002 : gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
1379 111639528 : for (i=m-1; i>=1; i--)
1380 : {
1381 58753897 : GEN u = gel(T, i);
1382 58753897 : GEN v = gel(Tp, i+1);
1383 58753897 : long n = lg(u)-1;
1384 58753897 : t = cgetg(n+1, t_VEC);
1385 146307430 : for (j=1, k=1; k<n; j++, k+=2)
1386 : {
1387 87689976 : gel(t, k) = modii(gel(v, j), gel(u, k));
1388 87720368 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
1389 : }
1390 58617454 : if (k==n) gel(t, k) = gel(v, j);
1391 58617454 : gel(Tp, i) = t;
1392 : }
1393 : {
1394 52885631 : GEN u = gel(T, i+1);
1395 52885631 : GEN v = gel(Tp, i+1);
1396 52885631 : long l = lg(u)-1;
1397 52885631 : if (typ(P)==t_VECSMALL)
1398 : {
1399 50892637 : GEN R = cgetg(n+1, t_VECSMALL);
1400 188783109 : for (j=1, k=1; j<=l; j++, k+=2)
1401 : {
1402 137581366 : uel(R,k) = umodiu(gel(v, j), P[k]);
1403 137851467 : if (k < n)
1404 110527042 : uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
1405 : }
1406 51201743 : return R;
1407 : }
1408 : else
1409 : {
1410 1992994 : GEN R = cgetg(n+1, t_VEC);
1411 5552869 : for (j=1, k=1; j<=l; j++, k+=2)
1412 : {
1413 3554706 : gel(R,k) = modii(gel(v, j), gel(P,k));
1414 3554678 : if (k < n)
1415 2891105 : gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
1416 : }
1417 1998163 : return R;
1418 : }
1419 : }
1420 : }
1421 :
1422 : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
1423 : GEN
1424 34778835 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
1425 : {
1426 34778835 : long m = lg(T)-1, n = lg(A)-1;
1427 : long i,j,k;
1428 34778835 : GEN Tp = cgetg(m+1, t_VEC);
1429 34751463 : GEN M = gel(T, 1);
1430 34751463 : GEN t = cgetg(lg(M), t_VEC);
1431 34709715 : if (typ(P)==t_VECSMALL)
1432 : {
1433 80887346 : for (j=1, k=1; k<n; j++, k+=2)
1434 : {
1435 60411866 : pari_sp av = avma;
1436 60411866 : GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
1437 60191831 : GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
1438 60351224 : gel(t, j) = gerepileuptoint(av, tj);
1439 : }
1440 20475480 : if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
1441 : } else
1442 : {
1443 30318165 : for (j=1, k=1; k<n; j++, k+=2)
1444 : {
1445 16047924 : pari_sp av = avma;
1446 16047924 : GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
1447 16059331 : GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
1448 16085307 : gel(t, j) = gerepileuptoint(av, tj);
1449 : }
1450 14270241 : if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
1451 : }
1452 34736173 : gel(Tp, 1) = t;
1453 66811088 : for (i=2; i<=m; i++)
1454 : {
1455 32060465 : GEN u = gel(T, i-1), M = gel(T, i);
1456 32060465 : GEN t = cgetg(lg(M), t_VEC);
1457 32092324 : GEN v = gel(Tp, i-1);
1458 32092324 : long n = lg(v)-1;
1459 87492007 : for (j=1, k=1; k<n; j++, k+=2)
1460 : {
1461 55417092 : pari_sp av = avma;
1462 55357501 : gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
1463 55417092 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
1464 : }
1465 32074915 : if (k==n) gel(t, j) = gel(v, k);
1466 32074915 : gel(Tp, i) = t;
1467 : }
1468 34750623 : return gmael(Tp,m,1);
1469 : }
1470 :
1471 : static GEN
1472 931544 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1473 : {
1474 931544 : long i, l = lg(gel(vA,1)), n = lg(P);
1475 931544 : GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
1476 30520609 : for (i=1; i < l; i++)
1477 : {
1478 29589191 : pari_sp av = avma;
1479 29589191 : GEN c, A = cgetg(n, typ(P));
1480 : long j;
1481 181956313 : for (j=1; j < n; j++) A[j] = mael(vA,j,i);
1482 29553802 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
1483 29588021 : gel(V,i) = gerepileuptoint(av, c);
1484 : }
1485 931418 : return V;
1486 : }
1487 :
1488 : static GEN
1489 392441 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1490 : {
1491 392441 : long i, j, l, n = lg(P);
1492 392441 : GEN mod = gmael(T, lg(T)-1, 1), V, w;
1493 392441 : w = cgetg(n, t_VECSMALL);
1494 1434443 : for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
1495 392428 : l = vecsmall_max(w);
1496 392433 : V = cgetg(l, t_POL);
1497 392416 : V[1] = mael(vA,1,1);
1498 2591519 : for (i=2; i < l; i++)
1499 : {
1500 2199086 : pari_sp av = avma;
1501 2199086 : GEN c, A = cgetg(n, typ(P));
1502 2198803 : if (typ(P)==t_VECSMALL)
1503 4669270 : for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
1504 : else
1505 3739003 : for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
1506 2198803 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
1507 2199190 : gel(V,i) = gerepileuptoint(av, c);
1508 : }
1509 392433 : return ZX_renormalize(V, l);
1510 : }
1511 :
1512 : static GEN
1513 4614 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1514 : {
1515 4614 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1516 4614 : GEN A = cgetg(n, t_VEC);
1517 4614 : GEN V = cgetg(l, t_COL);
1518 90887 : for (i=1; i < l; i++)
1519 : {
1520 335054 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1521 86273 : gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
1522 : }
1523 4614 : return V;
1524 : }
1525 :
1526 : static GEN
1527 114633 : polint_chinese(GEN worker, GEN mA, GEN P)
1528 : {
1529 114633 : long cnt, pending, n, i, j, l = lg(gel(mA,1));
1530 : struct pari_mt pt;
1531 : GEN done, va, M, A;
1532 : pari_timer ti;
1533 :
1534 114633 : if (l == 1) return cgetg(1, t_MAT);
1535 85581 : cnt = pending = 0;
1536 85581 : n = lg(P);
1537 85581 : A = cgetg(n, t_VEC);
1538 85581 : va = mkvec(A);
1539 85581 : M = cgetg(l, t_MAT);
1540 85581 : if (DEBUGLEVEL>4) timer_start(&ti);
1541 85581 : if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
1542 85581 : mt_queue_start_lim(&pt, worker, l-1);
1543 671529 : for (i=1; i<l || pending; i++)
1544 : {
1545 : long workid;
1546 2421256 : for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
1547 585948 : mt_queue_submit(&pt, i, i<l? va: NULL);
1548 585948 : done = mt_queue_get(&pt, &workid, &pending);
1549 585948 : if (done)
1550 : {
1551 552208 : gel(M,workid) = done;
1552 552208 : if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
1553 : }
1554 : }
1555 85581 : if (DEBUGLEVEL>5) err_printf("\n");
1556 85581 : if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
1557 85581 : mt_queue_end(&pt);
1558 85581 : return M;
1559 : }
1560 :
1561 : GEN
1562 840 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
1563 : {
1564 840 : return nxCV_polint_center_tree(vA, P, T, R, m2);
1565 : }
1566 :
1567 : static GEN
1568 430 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1569 : {
1570 430 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1571 430 : GEN A = cgetg(n, t_VEC);
1572 430 : GEN V = cgetg(l, t_MAT);
1573 4204 : for (i=1; i < l; i++)
1574 : {
1575 15299 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1576 3774 : gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
1577 : }
1578 430 : return V;
1579 : }
1580 :
1581 : static GEN
1582 90 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
1583 : {
1584 90 : GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
1585 90 : return polint_chinese(worker, mA, P);
1586 : }
1587 :
1588 : static GEN
1589 52973 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1590 : {
1591 52973 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1592 52973 : GEN A = cgetg(n, t_VEC);
1593 52973 : GEN V = cgetg(l, t_MAT);
1594 418874 : for (i=1; i < l; i++)
1595 : {
1596 2132700 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1597 365898 : gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
1598 : }
1599 52976 : return V;
1600 : }
1601 :
1602 : GEN
1603 551268 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
1604 : {
1605 551268 : return ncV_polint_center_tree(vA, P, T, R, m2);
1606 : }
1607 :
1608 : static GEN
1609 114543 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
1610 : {
1611 114543 : GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
1612 114543 : return polint_chinese(worker, mA, P);
1613 : }
1614 :
1615 : /* return [A mod P[i], i=1..#P] */
1616 : GEN
1617 0 : Z_ZV_mod(GEN A, GEN P)
1618 : {
1619 0 : pari_sp av = avma;
1620 0 : return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
1621 : }
1622 : /* P a t_VECSMALL */
1623 : GEN
1624 0 : Z_nv_mod(GEN A, GEN P)
1625 : {
1626 0 : pari_sp av = avma;
1627 0 : return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
1628 : }
1629 : /* B a ZX, T = ZV_producttree(P) */
1630 : GEN
1631 1239547 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
1632 : {
1633 : pari_sp av;
1634 1239547 : long i, j, l = lg(B), n = lg(A)-1;
1635 1239547 : GEN V = cgetg(n+1, t_VEC);
1636 6113346 : for (j=1; j <= n; j++)
1637 : {
1638 4873995 : gel(V, j) = cgetg(l, t_VECSMALL);
1639 4873813 : mael(V, j, 1) = B[1]&VARNBITS;
1640 : }
1641 1239351 : av = avma;
1642 12013888 : for (i=2; i < l; i++)
1643 : {
1644 10775780 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
1645 74818998 : for (j=1; j <= n; j++)
1646 64052394 : mael(V, j, i) = v[j];
1647 10766604 : set_avma(av);
1648 : }
1649 6112621 : for (j=1; j <= n; j++)
1650 4874435 : (void) Flx_renormalize(gel(V, j), l);
1651 1238186 : return V;
1652 : }
1653 :
1654 : static GEN
1655 159883 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
1656 :
1657 : GEN
1658 10277 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
1659 : {
1660 10277 : pari_sp av = avma;
1661 10277 : long i, j, l = lg(P), n = lg(xa)-1;
1662 10277 : GEN V = cgetg(n+1, t_VEC);
1663 39626 : for (j=1; j <= n; j++)
1664 : {
1665 29349 : gel(V, j) = cgetg(l, t_POL);
1666 29349 : mael(V, j, 1) = P[1]&VARNBITS;
1667 : }
1668 88474 : for (i=2; i < l; i++)
1669 : {
1670 78198 : GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
1671 321139 : for (j=1; j <= n; j++)
1672 242942 : gmael(V, j, i) = gel(v,j);
1673 : }
1674 39625 : for (j=1; j <= n; j++)
1675 29349 : (void) FlxX_renormalize(gel(V, j), l);
1676 10276 : return gerepilecopy(av, V);
1677 : }
1678 :
1679 : GEN
1680 4047 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
1681 : {
1682 4047 : pari_sp av = avma;
1683 4047 : long i, j, l = lg(C), n = lg(xa)-1;
1684 4047 : GEN V = cgetg(n+1, t_VEC);
1685 16915 : for (j = 1; j <= n; j++)
1686 12870 : gel(V, j) = cgetg(l, t_COL);
1687 85716 : for (i = 1; i < l; i++)
1688 : {
1689 81683 : GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
1690 358530 : for (j = 1; j <= n; j++)
1691 276859 : gmael(V, j, i) = gel(v,j);
1692 : }
1693 4033 : return gerepilecopy(av, V);
1694 : }
1695 :
1696 : GEN
1697 430 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
1698 : {
1699 430 : pari_sp av = avma;
1700 430 : long i, j, l = lg(M), n = lg(xa)-1;
1701 430 : GEN V = cgetg(n+1, t_VEC);
1702 2083 : for (j=1; j <= n; j++)
1703 1653 : gel(V, j) = cgetg(l, t_MAT);
1704 4204 : for (i=1; i < l; i++)
1705 : {
1706 3774 : GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
1707 15299 : for (j=1; j <= n; j++)
1708 11525 : gmael(V, j, i) = gel(v,j);
1709 : }
1710 430 : return gerepilecopy(av, V);
1711 : }
1712 :
1713 : GEN
1714 936200 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
1715 : {
1716 : pari_sp av;
1717 936200 : long i, j, l = lg(B), n = lg(A)-1;
1718 936200 : GEN V = cgetg(n+1, t_VEC);
1719 4637043 : for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
1720 936127 : av = avma;
1721 41199659 : for (i=1; i < l; i++)
1722 : {
1723 40267801 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
1724 226384133 : for (j=1; j <= n; j++) mael(V, j, i) = v[j];
1725 40194472 : set_avma(av);
1726 : }
1727 931858 : return V;
1728 : }
1729 :
1730 : GEN
1731 80586 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
1732 : {
1733 80586 : pari_sp av = avma;
1734 80586 : long i, j, l = lg(M), n = lg(xa)-1;
1735 80586 : GEN V = cgetg(n+1, t_VEC);
1736 315207 : for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_MAT);
1737 1016552 : for (i=1; i < l; i++)
1738 : {
1739 935974 : GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
1740 4637166 : for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
1741 : }
1742 80578 : return gerepilecopy(av, V);
1743 : }
1744 :
1745 : static GEN
1746 1994331 : ZV_sqr(GEN z)
1747 : {
1748 1994331 : long i,l = lg(z);
1749 1994331 : GEN x = cgetg(l, t_VEC);
1750 1994353 : if (typ(z)==t_VECSMALL)
1751 4591041 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
1752 : else
1753 3823221 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
1754 1994273 : return x;
1755 : }
1756 :
1757 : static GEN
1758 10538062 : ZT_sqr(GEN x)
1759 : {
1760 10538062 : if (typ(x) == t_INT) return sqri(x);
1761 13796562 : pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
1762 : }
1763 :
1764 : static GEN
1765 1994323 : ZV_invdivexact(GEN y, GEN x)
1766 : {
1767 1994323 : long i, l = lg(y);
1768 1994323 : GEN z = cgetg(l,t_VEC);
1769 1994344 : if (typ(x)==t_VECSMALL)
1770 4590603 : for (i=1; i<l; i++)
1771 : {
1772 3713014 : pari_sp av = avma;
1773 3713014 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
1774 3713295 : set_avma(av); gel(z,i) = utoi(a);
1775 : }
1776 : else
1777 3823317 : for (i=1; i<l; i++)
1778 2706588 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
1779 1994318 : return z;
1780 : }
1781 :
1782 : /* P t_VECSMALL or t_VEC of t_INT */
1783 : GEN
1784 1994305 : ZV_chinesetree(GEN P, GEN T)
1785 : {
1786 1994305 : GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
1787 1994310 : GEN mod = gmael(T,lg(T)-1,1);
1788 1994310 : return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
1789 : }
1790 :
1791 : static GEN
1792 485179 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
1793 : {
1794 485179 : if (!pt_mod)
1795 11159 : return gerepileupto(av, a);
1796 : else
1797 : {
1798 474020 : GEN mod = gmael(T, lg(T)-1, 1);
1799 474020 : gerepileall(av, 2, &a, &mod);
1800 474021 : *pt_mod = mod;
1801 474021 : return a;
1802 : }
1803 : }
1804 :
1805 : GEN
1806 160311 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1807 : {
1808 160311 : pari_sp av = avma;
1809 160311 : GEN T = ZV_producttree(P);
1810 160311 : GEN R = ZV_chinesetree(P, T);
1811 160311 : GEN a = ZV_chinese_tree(A, P, T, R);
1812 160311 : GEN mod = gmael(T, lg(T)-1, 1);
1813 160311 : GEN ca = Fp_center(a, mod, shifti(mod,-1));
1814 160311 : return gc_chinese(av, T, ca, pt_mod);
1815 : }
1816 :
1817 : GEN
1818 5088 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
1819 : {
1820 5088 : pari_sp av = avma;
1821 5088 : GEN T = ZV_producttree(P);
1822 5088 : GEN R = ZV_chinesetree(P, T);
1823 5088 : GEN a = ZV_chinese_tree(A, P, T, R);
1824 5088 : return gc_chinese(av, T, a, pt_mod);
1825 : }
1826 :
1827 : GEN
1828 110134 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
1829 : {
1830 110134 : pari_sp av = avma;
1831 110134 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1832 110134 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
1833 110135 : return gerepileupto(av, a);
1834 : }
1835 :
1836 : GEN
1837 196024 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1838 : {
1839 196024 : pari_sp av = avma;
1840 196024 : GEN T = ZV_producttree(P);
1841 196023 : GEN R = ZV_chinesetree(P, T);
1842 196022 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1843 196022 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
1844 196024 : return gc_chinese(av, T, a, pt_mod);
1845 : }
1846 :
1847 : GEN
1848 9123 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1849 : {
1850 9123 : pari_sp av = avma;
1851 9123 : GEN T = ZV_producttree(P);
1852 9123 : GEN R = ZV_chinesetree(P, T);
1853 9123 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1854 9123 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
1855 9123 : return gc_chinese(av, T, a, pt_mod);
1856 : }
1857 :
1858 : GEN
1859 5262 : ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
1860 : {
1861 5262 : pari_sp av = avma;
1862 5262 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1863 5262 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
1864 5262 : return gerepileupto(av, a);
1865 : }
1866 :
1867 : GEN
1868 0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
1869 : {
1870 0 : pari_sp av = avma;
1871 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1872 0 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
1873 0 : return gerepileupto(av, a);
1874 : }
1875 :
1876 : GEN
1877 52975 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
1878 : {
1879 52975 : pari_sp av = avma;
1880 52975 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1881 52973 : GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
1882 52976 : return gerepileupto(av, a);
1883 : }
1884 :
1885 : GEN
1886 114543 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1887 : {
1888 114543 : pari_sp av = avma;
1889 114543 : GEN T = ZV_producttree(P);
1890 114543 : GEN R = ZV_chinesetree(P, T);
1891 114543 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1892 114543 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
1893 114543 : return gc_chinese(av, T, a, pt_mod);
1894 : }
1895 :
1896 : GEN
1897 0 : nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
1898 : {
1899 0 : pari_sp av = avma;
1900 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1901 0 : GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
1902 0 : return gerepileupto(av, a);
1903 : }
1904 :
1905 : GEN
1906 0 : nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1907 : {
1908 0 : pari_sp av = avma;
1909 0 : GEN T = ZV_producttree(P);
1910 0 : GEN R = ZV_chinesetree(P, T);
1911 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1912 0 : GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
1913 0 : return gc_chinese(av, T, a, pt_mod);
1914 : }
1915 :
1916 : GEN
1917 430 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
1918 : {
1919 430 : pari_sp av = avma;
1920 430 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1921 430 : GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
1922 430 : return gerepileupto(av, a);
1923 : }
1924 :
1925 : GEN
1926 90 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1927 : {
1928 90 : pari_sp av = avma;
1929 90 : GEN T = ZV_producttree(P);
1930 90 : GEN R = ZV_chinesetree(P, T);
1931 90 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1932 90 : GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
1933 90 : return gc_chinese(av, T, a, pt_mod);
1934 : }
1935 :
1936 : /**********************************************************************
1937 : ** Powering over (Z/NZ)^*, small N **
1938 : **********************************************************************/
1939 :
1940 : /* 2^n mod p; assume n > 1 */
1941 : static ulong
1942 11576477 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
1943 : {
1944 11576477 : ulong y = 2;
1945 11576477 : int j = 1+bfffo(n);
1946 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
1947 11576477 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
1948 390056122 : for (; j; n<<=1,j--)
1949 : {
1950 378480014 : y = Fl_sqr_pre(y,p,pi);
1951 378449766 : if (n & HIGHBIT) y = Fl_double(y, p);
1952 : }
1953 11576108 : return y;
1954 : }
1955 :
1956 : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
1957 : static ulong
1958 1452803 : Fl_2powu(ulong n, ulong p)
1959 : {
1960 1452803 : ulong y = 2;
1961 1452803 : int j = 1+bfffo(n);
1962 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
1963 1452803 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
1964 5597052 : for (; j; n<<=1,j--)
1965 : {
1966 4144240 : y = (y*y) % p;
1967 4144240 : if (n & HIGHBIT) y = Fl_double(y, p);
1968 : }
1969 1452812 : return y;
1970 : }
1971 :
1972 : ulong
1973 52218133 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
1974 : {
1975 : ulong y, z, n;
1976 52218133 : if (n0 <= 1)
1977 : { /* frequent special cases */
1978 5063779 : if (n0 == 1) return x;
1979 1346211 : if (n0 == 0) return 1;
1980 : }
1981 47154293 : if (x <= 2)
1982 : {
1983 12124379 : if (x == 2) return Fl_2powu_pre(n0, p, pi);
1984 546823 : return x; /* 0 or 1 */
1985 : }
1986 35029914 : y = 1; z = x; n = n0;
1987 : for(;;)
1988 : {
1989 448381005 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
1990 448556944 : n>>=1; if (!n) return y;
1991 413556304 : z = Fl_sqr_pre(z,p,pi);
1992 : }
1993 : }
1994 :
1995 : ulong
1996 147114994 : Fl_powu(ulong x, ulong n0, ulong p)
1997 : {
1998 : ulong y, z, n;
1999 147114994 : if (n0 <= 2)
2000 : { /* frequent special cases */
2001 92252795 : if (n0 == 2) return Fl_sqr(x,p);
2002 23395317 : if (n0 == 1) return x;
2003 500151 : if (n0 == 0) return 1;
2004 : }
2005 54820297 : if (x <= 1) return x; /* 0 or 1 */
2006 54605532 : if (p & HIGHMASK)
2007 6317281 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
2008 48288251 : if (x == 2) return Fl_2powu(n0, p);
2009 46835447 : y = 1; z = x; n = n0;
2010 : for(;;)
2011 : {
2012 253193491 : if (n&1) y = (y*z) % p;
2013 253193491 : n>>=1; if (!n) return y;
2014 206358044 : z = (z*z) % p;
2015 : }
2016 : }
2017 :
2018 : /* Reduce data dependency to maximize internal parallelism */
2019 : GEN
2020 11330899 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
2021 : {
2022 : long i, k;
2023 11330899 : GEN powers = cgetg(n + 2, t_VECSMALL);
2024 11318671 : powers[1] = 1; if (n == 0) return powers;
2025 11318671 : powers[2] = x;
2026 282845272 : for (i = 3, k=2; i <= n; i+=2, k++)
2027 : {
2028 271510444 : powers[i] = Fl_sqr_pre(powers[k], p, pi);
2029 271511209 : powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
2030 : }
2031 11334828 : if (i==n+1)
2032 9799232 : powers[i] = Fl_sqr_pre(powers[k], p, pi);
2033 11335632 : return powers;
2034 : }
2035 :
2036 : GEN
2037 24981 : Fl_powers(ulong x, long n, ulong p)
2038 : {
2039 24981 : return Fl_powers_pre(x, n, p, get_Fl_red(p));
2040 : }
2041 :
2042 : /**********************************************************************
2043 : ** Powering over (Z/NZ)^*, large N **
2044 : **********************************************************************/
2045 :
2046 : static GEN
2047 4708642 : Fp_dblsqr(GEN x, GEN N)
2048 : {
2049 4708642 : GEN z = shifti(Fp_sqr(x, N), 1);
2050 4587246 : return cmpii(z, N) >= 0? subii(z, N): z;
2051 : }
2052 :
2053 : typedef struct muldata {
2054 : GEN (*sqr)(void * E, GEN x);
2055 : GEN (*mul)(void * E, GEN x, GEN y);
2056 : GEN (*mul2)(void * E, GEN x);
2057 : } muldata;
2058 :
2059 : /* modified Barrett reduction with one fold */
2060 : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
2061 :
2062 : static GEN
2063 8611 : Fp_invmBarrett(GEN p, long s)
2064 : {
2065 8611 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
2066 8611 : return mkvec2(Q,R);
2067 : }
2068 :
2069 : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
2070 : * a = r (mod N) */
2071 : static GEN
2072 4698940 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
2073 : {
2074 4698940 : pari_sp av = avma;
2075 4698940 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
2076 4698940 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
2077 4698940 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
2078 4698940 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
2079 4698940 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
2080 4698940 : GEN r = subii(A, mulii(q, N));
2081 4698940 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
2082 4698940 : if (signe(sr)<0) return gerepileuptoint(av, r);
2083 2787452 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
2084 2787452 : if (signe(sr)<0) return gerepileuptoint(av, r);
2085 107472 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
2086 107472 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
2087 : }
2088 :
2089 : /* Montgomery reduction */
2090 :
2091 : INLINE ulong
2092 807181 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
2093 :
2094 : struct montred
2095 : {
2096 : GEN N;
2097 : ulong inv;
2098 : };
2099 :
2100 : /* Montgomery reduction */
2101 : static GEN
2102 49771465 : _sqr_montred(void * E, GEN x)
2103 : {
2104 49771465 : struct montred * D = (struct montred *) E;
2105 49771465 : return red_montgomery(sqri(x), D->N, D->inv);
2106 : }
2107 :
2108 : /* Montgomery reduction */
2109 : static GEN
2110 4656997 : _mul_montred(void * E, GEN x, GEN y)
2111 : {
2112 4656997 : struct montred * D = (struct montred *) E;
2113 4656997 : return red_montgomery(mulii(x, y), D->N, D->inv);
2114 : }
2115 :
2116 : static GEN
2117 7528967 : _mul2_montred(void * E, GEN x)
2118 : {
2119 7528967 : struct montred * D = (struct montred *) E;
2120 7528967 : GEN z = shifti(_sqr_montred(E, x), 1);
2121 7526257 : long l = lgefint(D->N);
2122 7969358 : while (lgefint(z) > l) z = subii(z, D->N);
2123 7526832 : return z;
2124 : }
2125 :
2126 : static GEN
2127 19737303 : _sqr_remii(void* N, GEN x)
2128 19737303 : { return remii(sqri(x), (GEN) N); }
2129 :
2130 : static GEN
2131 1432600 : _mul_remii(void* N, GEN x, GEN y)
2132 1432600 : { return remii(mulii(x, y), (GEN) N); }
2133 :
2134 : static GEN
2135 3292437 : _mul2_remii(void* N, GEN x)
2136 3292437 : { return Fp_dblsqr(x, (GEN) N); }
2137 :
2138 : struct redbarrett
2139 : {
2140 : GEN iM, N;
2141 : long s;
2142 : };
2143 :
2144 : static GEN
2145 4243321 : _sqr_remiibar(void *E, GEN x)
2146 : {
2147 4243321 : struct redbarrett * D = (struct redbarrett *) E;
2148 4243321 : return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
2149 : }
2150 :
2151 : static GEN
2152 455619 : _mul_remiibar(void *E, GEN x, GEN y)
2153 : {
2154 455619 : struct redbarrett * D = (struct redbarrett *) E;
2155 455619 : return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
2156 : }
2157 :
2158 : static GEN
2159 1413630 : _mul2_remiibar(void *E, GEN x)
2160 : {
2161 1413630 : struct redbarrett * D = (struct redbarrett *) E;
2162 1413630 : return Fp_dblsqr(x, D->N);
2163 : }
2164 :
2165 : static long
2166 1017199 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
2167 : {
2168 1017199 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
2169 : {
2170 8611 : struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
2171 8611 : D->sqr = &_sqr_remiibar;
2172 8611 : D->mul = &_mul_remiibar;
2173 8611 : D->mul2 = &_mul2_remiibar;
2174 8611 : E->N = N;
2175 8611 : E->s = 1+(expi(N)>>1);
2176 8611 : E->iM = Fp_invmBarrett(N, E->s);
2177 8611 : *pt_E = (void*) E;
2178 8611 : return 0;
2179 : }
2180 1008588 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
2181 : {
2182 807190 : struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
2183 807190 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
2184 807185 : D->sqr = &_sqr_montred;
2185 807185 : D->mul = &_mul_montred;
2186 807185 : D->mul2 = &_mul2_montred;
2187 807185 : E->N = N;
2188 807185 : E->inv = init_montdata(N);
2189 807175 : *pt_E = (void*) E;
2190 807175 : return 1;
2191 : }
2192 : else
2193 : {
2194 201393 : D->sqr = &_sqr_remii;
2195 201393 : D->mul = &_mul_remii;
2196 201393 : D->mul2 = &_mul2_remii;
2197 201393 : *pt_E = (void*) N;
2198 201393 : return 0;
2199 : }
2200 : }
2201 :
2202 : GEN
2203 3055308 : Fp_powu(GEN A, ulong k, GEN N)
2204 : {
2205 3055308 : long lN = lgefint(N);
2206 : int base_is_2, use_montgomery;
2207 : muldata D;
2208 : void *E;
2209 : pari_sp av;
2210 :
2211 3055308 : if (lN == 3) {
2212 1432839 : ulong n = uel(N,2);
2213 1432839 : return utoi( Fl_powu(umodiu(A, n), k, n) );
2214 : }
2215 1622469 : if (k <= 2)
2216 : { /* frequent special cases */
2217 814860 : if (k == 2) return Fp_sqr(A,N);
2218 275650 : if (k == 1) return A;
2219 0 : if (k == 0) return gen_1;
2220 : }
2221 807608 : av = avma; A = modii(A,N);
2222 807609 : base_is_2 = 0;
2223 807609 : if (lgefint(A) == 3) switch(A[2])
2224 : {
2225 770 : case 1: set_avma(av); return gen_1;
2226 38967 : case 2: base_is_2 = 1; break;
2227 : }
2228 :
2229 : /* TODO: Move this out of here and use for general modular computations */
2230 806839 : use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
2231 806839 : if (base_is_2)
2232 38967 : A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
2233 : else
2234 767872 : A = gen_powu_i(A, k, E, D.sqr, D.mul);
2235 806839 : if (use_montgomery)
2236 : {
2237 686756 : A = red_montgomery(A, N, ((struct montred *) E)->inv);
2238 686756 : if (cmpii(A, N) >= 0) A = subii(A,N);
2239 : }
2240 806839 : return gerepileuptoint(av, A);
2241 : }
2242 :
2243 : GEN
2244 22428 : Fp_pows(GEN A, long k, GEN N)
2245 : {
2246 22428 : if (lgefint(N) == 3) {
2247 7939 : ulong n = N[2];
2248 7939 : ulong a = umodiu(A, n);
2249 7939 : if (k < 0) {
2250 133 : a = Fl_inv(a, n);
2251 133 : k = -k;
2252 : }
2253 7939 : return utoi( Fl_powu(a, (ulong)k, n) );
2254 : }
2255 14489 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
2256 14489 : return Fp_powu(A, (ulong)k, N);
2257 : }
2258 :
2259 : /* A^K mod N */
2260 : GEN
2261 34430762 : Fp_pow(GEN A, GEN K, GEN N)
2262 : {
2263 : pari_sp av;
2264 34430762 : long s, lN = lgefint(N), sA, sy;
2265 : int base_is_2, use_montgomery;
2266 : GEN y;
2267 : muldata D;
2268 : void *E;
2269 :
2270 34430762 : s = signe(K);
2271 34430762 : if (!s) return dvdii(A,N)? gen_0: gen_1;
2272 33401821 : if (lN == 3 && lgefint(K) == 3)
2273 : {
2274 32508094 : ulong n = N[2], a = umodiu(A, n);
2275 32508451 : if (s < 0) a = Fl_inv(a, n);
2276 32508520 : if (a <= 1) return utoi(a); /* 0 or 1 */
2277 29547888 : return utoi(Fl_powu(a, uel(K,2), n));
2278 : }
2279 :
2280 893727 : av = avma;
2281 893727 : if (s < 0) y = Fp_inv(A,N);
2282 : else
2283 : {
2284 891770 : y = modii(A,N);
2285 892008 : if (!signe(y)) { set_avma(av); return gen_0; }
2286 : }
2287 893966 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
2288 :
2289 210536 : base_is_2 = 0;
2290 210536 : sy = abscmpii(y, shifti(N,-1)) > 0;
2291 210561 : if (sy) y = subii(N,y);
2292 210565 : sA = sy && mod2(K);
2293 210563 : if (lgefint(y) == 3) switch(y[2])
2294 : {
2295 199 : case 1: set_avma(av); return sA ? subis(N,1): gen_1;
2296 139823 : case 2: base_is_2 = 1; break;
2297 : }
2298 :
2299 : /* TODO: Move this out of here and use for general modular computations */
2300 210364 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
2301 210348 : if (base_is_2)
2302 139833 : y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
2303 : else
2304 70515 : y = gen_pow_i(y, K, E, D.sqr, D.mul);
2305 210376 : if (use_montgomery)
2306 : {
2307 120442 : y = red_montgomery(y, N, ((struct montred *) E)->inv);
2308 120443 : if (cmpii(y,N) >= 0) y = subii(y,N);
2309 : }
2310 210376 : if (sA) y = subii(N, y);
2311 210376 : return gerepileuptoint(av,y);
2312 : }
2313 :
2314 : static GEN
2315 14030630 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
2316 : static GEN
2317 8126876 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
2318 : static GEN
2319 56163 : _Fp_one(void *E) { (void) E; return gen_1; }
2320 :
2321 : GEN
2322 77 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
2323 77 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
2324 :
2325 : GEN
2326 54096 : Fp_pow_table(GEN R, GEN n, GEN p)
2327 54096 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
2328 :
2329 : GEN
2330 4523 : Fp_powers(GEN x, long n, GEN p)
2331 : {
2332 4523 : if (lgefint(p) == 3)
2333 2456 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
2334 2067 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
2335 : }
2336 :
2337 : GEN
2338 434 : FpV_prod(GEN V, GEN p) { return gen_product(V, (void *)p, &_Fp_mul); }
2339 :
2340 : static GEN
2341 27110731 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
2342 : static GEN
2343 147 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
2344 :
2345 : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
2346 : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
2347 : equalii,equali1,Fp_easylog};
2348 :
2349 : static GEN
2350 761553 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
2351 : static GEN
2352 853646 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
2353 : static GEN
2354 767753 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
2355 : static GEN
2356 497784 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
2357 : static GEN
2358 33090 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
2359 : static int
2360 218619 : _Fp_equal0(GEN x) { return signe(x)==0; }
2361 : static GEN
2362 28761 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
2363 :
2364 : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
2365 : _Fp_inv,_Fp_equal0,_Fp_s};
2366 :
2367 7077 : const struct bb_field *get_Fp_field(void **E, GEN p)
2368 7077 : { *E = (void*)p; return &Fp_field; }
2369 :
2370 : /*********************************************************************/
2371 : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
2372 : /*********************************************************************/
2373 : ulong
2374 15078 : Fl_order(ulong a, ulong o, ulong p)
2375 : {
2376 15078 : pari_sp av = avma;
2377 : GEN m, P, E;
2378 : long i;
2379 15078 : if (a==1) return 1;
2380 9954 : if (!o) o = p-1;
2381 9954 : m = factoru(o);
2382 9954 : P = gel(m,1);
2383 9954 : E = gel(m,2);
2384 25018 : for (i = lg(P)-1; i; i--)
2385 : {
2386 15064 : ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
2387 15064 : if (y == 1) o = t;
2388 16779 : else for (j = 1; j < e; j++)
2389 : {
2390 5257 : y = Fl_powu(y, l, p);
2391 5257 : if (y == 1) { o = t * upowuu(l, j); break; }
2392 : }
2393 : }
2394 9954 : return gc_ulong(av, o);
2395 : }
2396 :
2397 : /*Find the exact order of a assuming a^o==1*/
2398 : GEN
2399 10438 : Fp_order(GEN a, GEN o, GEN p) {
2400 10438 : if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
2401 : {
2402 21 : ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
2403 21 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
2404 : }
2405 10417 : return gen_order(a, o, (void*)p, &Fp_star);
2406 : }
2407 : GEN
2408 70 : Fp_factored_order(GEN a, GEN o, GEN p)
2409 70 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
2410 :
2411 : /* return order of a mod p^e, e > 0, pe = p^e */
2412 : static GEN
2413 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
2414 : {
2415 : GEN ap, op;
2416 70 : if (absequaliu(p, 2))
2417 : {
2418 56 : if (e == 1) return gen_1;
2419 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
2420 49 : if (mod4(a) == 1) op = gen_1; else { op = gen_2; a = Fp_sqr(a, pe); }
2421 : } else {
2422 14 : ap = (e == 1)? a: remii(a,p);
2423 14 : op = Fp_order(ap, subiu(p,1), p);
2424 14 : if (e == 1) return op;
2425 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
2426 : }
2427 49 : if (equali1(a)) return op;
2428 7 : return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
2429 : }
2430 :
2431 : GEN
2432 63 : znorder(GEN x, GEN o)
2433 : {
2434 63 : pari_sp av = avma;
2435 : GEN b, a;
2436 :
2437 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
2438 56 : b = gel(x,1); a = gel(x,2);
2439 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
2440 49 : if (!o)
2441 : {
2442 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
2443 35 : long i, l = lg(P);
2444 35 : o = gen_1;
2445 70 : for (i = 1; i < l; i++)
2446 : {
2447 35 : GEN p = gel(P,i);
2448 35 : long e = itos(gel(E,i));
2449 :
2450 35 : if (l == 2)
2451 35 : o = Zp_order(a, p, e, b);
2452 : else {
2453 0 : GEN pe = powiu(p,e);
2454 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
2455 : }
2456 : }
2457 35 : return gerepileuptoint(av, o);
2458 : }
2459 14 : return Fp_order(a, o, b);
2460 : }
2461 :
2462 : /*********************************************************************/
2463 : /** DISCRETE LOGARITHM in (Z/nZ)* **/
2464 : /*********************************************************************/
2465 : static GEN
2466 58614 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
2467 : {
2468 58614 : pari_sp av = avma;
2469 : GEN h1, h2, F, G;
2470 58614 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
2471 35205 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
2472 : {
2473 126 : GEN M = cgetg(3, t_MAT);
2474 126 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
2475 126 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
2476 126 : return gerepileupto(av, M);
2477 : }
2478 35079 : return gc_NULL(av);
2479 : }
2480 :
2481 : static GEN
2482 58614 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
2483 : {
2484 : GEN rel;
2485 58614 : do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
2486 58614 : while (!rel);
2487 126 : return rel;
2488 : }
2489 :
2490 : struct Fp_log_rel
2491 : {
2492 : GEN rel;
2493 : ulong prmax;
2494 : long nbrel, nbmax, nbgen;
2495 : };
2496 :
2497 : /* add u^e */
2498 : static void
2499 3157 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
2500 : {
2501 3157 : pari_sp av = avma;
2502 3157 : long off = r->prmax+1;
2503 3157 : GEN F = cgetg(3, t_MAT);
2504 3157 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
2505 3157 : gel(F,2) = vecsmall_append(gel(z,2), e);
2506 3157 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
2507 3157 : }
2508 :
2509 : /* add u^-1 v^-1 */
2510 : static void
2511 104083 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
2512 : {
2513 104083 : pari_sp av = avma;
2514 104083 : long off = r->prmax+1;
2515 104083 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
2516 104083 : GEN F = cgetg(3, t_MAT);
2517 104083 : gel(F,1) = vecsmall_concat(gel(z,1), P);
2518 104083 : gel(F,2) = vecsmall_concat(gel(z,2), E);
2519 104083 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
2520 104083 : }
2521 :
2522 : /*
2523 : Let p=C^2+c
2524 : Solve h = (C+x)*(C+a)-p = 0 [mod l]
2525 : h= -c+x*(C+a)+C*a = 0 [mod l]
2526 : x = (c-C*a)/(C+a) [mod l]
2527 : h = -c+C*(x+a)+a*x
2528 : */
2529 :
2530 : GEN
2531 40827 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
2532 : {
2533 40827 : pari_sp ltop = avma;
2534 40827 : long i, j, th, n = lg(pi)-1, rel = 1;
2535 40827 : GEN sieve = zero_zv(a+2)+1;
2536 40850 : GEN L = cgetg(1+a+2, t_VEC);
2537 40837 : pari_sp av = avma;
2538 40837 : GEN z, h = addis(C,a);
2539 40817 : if ((z = Z_issmooth_fact(h, prmax)))
2540 : {
2541 3007 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
2542 3016 : av = avma;
2543 : }
2544 16835364 : for (i=1; i<=n; i++)
2545 : {
2546 16809850 : ulong li = pi[i], s = sz[i], al = a % li;
2547 16809850 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
2548 17154451 : if (!iv) continue;
2549 16730061 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
2550 77389417 : for(j = u; j<=a; j+=li) sieve[j] += s;
2551 : }
2552 35046 : if (a)
2553 : {
2554 40737 : long e = expi(mulis(C,a));
2555 40798 : th = e - expu(e) - 1;
2556 54 : } else th = -1;
2557 28023325 : for (j=0; j<a; j++)
2558 27981761 : if (sieve[j]>=th)
2559 : {
2560 119453 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
2561 119334 : if ((z = Z_issmooth_fact(h, prmax)))
2562 : {
2563 109798 : gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
2564 109993 : av = avma;
2565 9431 : } else set_avma(av);
2566 : }
2567 : /* j = a */
2568 41564 : if (sieve[a]>=th)
2569 : {
2570 476 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
2571 476 : if ((z = Z_issmooth_fact(h, prmax)))
2572 385 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
2573 : }
2574 41564 : setlg(L, rel); return gerepilecopy(ltop, L);
2575 : }
2576 :
2577 : static long
2578 63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
2579 : {
2580 : struct pari_mt pt;
2581 63 : long i, j, nb = 0;
2582 63 : GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
2583 : mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
2584 63 : long running, pending = 0;
2585 63 : GEN W = zerovec(r->nbgen);
2586 63 : mt_queue_start_lim(&pt, worker, r->nbgen);
2587 41229 : for (i = 0; (running = (i < r->nbgen)) || pending; i++)
2588 : {
2589 : GEN done;
2590 : long idx;
2591 41166 : mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
2592 41166 : done = mt_queue_get(&pt, &idx, &pending);
2593 41166 : if (!done || lg(done)==1) continue;
2594 35917 : gel(W, idx+1) = done;
2595 35917 : nb += lg(done)-1;
2596 35917 : if (DEBUGLEVEL && (i&127)==0)
2597 0 : err_printf("%ld%% ",100*nb/r->nbmax);
2598 : }
2599 63 : mt_queue_end(&pt);
2600 39550 : for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
2601 : {
2602 : long ll, m;
2603 39487 : GEN L = gel(W,j);
2604 39487 : if (isintzero(L)) continue;
2605 34531 : ll = lg(L);
2606 141771 : for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
2607 : {
2608 107240 : GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
2609 107240 : if (v[1] == 1)
2610 3157 : addifsmooth1(r, h, v[2], v[3]);
2611 : else
2612 104083 : addifsmooth2(r, h, v[2], v[3]);
2613 : }
2614 : }
2615 63 : return j;
2616 : }
2617 :
2618 : static GEN
2619 735 : ECP_psi(GEN x, GEN y)
2620 : {
2621 735 : long prec = realprec(x);
2622 735 : GEN lx = glog(x, prec), ly = glog(y, prec);
2623 735 : GEN u = gdiv(lx, ly);
2624 735 : return gpow(u, gneg(u),prec);
2625 : }
2626 :
2627 : struct computeG
2628 : {
2629 : GEN C;
2630 : long bnd, nbi;
2631 : };
2632 :
2633 : static GEN
2634 735 : _computeG(void *E, GEN gen)
2635 : {
2636 735 : struct computeG * d = (struct computeG *) E;
2637 735 : GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
2638 735 : return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
2639 : }
2640 :
2641 : static long
2642 63 : compute_nbgen(GEN C, long bnd, long nbi)
2643 : {
2644 : struct computeG d;
2645 63 : d.C = shifti(C, 1);
2646 63 : d.bnd = bnd;
2647 63 : d.nbi = nbi;
2648 63 : return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
2649 : }
2650 :
2651 : static GEN
2652 1714 : _psi(void*E, GEN y)
2653 : {
2654 1714 : GEN lx = (GEN) E;
2655 1714 : long prec = realprec(lx);
2656 1714 : GEN ly = glog(y, prec);
2657 1714 : GEN u = gdiv(lx, ly);
2658 1714 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
2659 : }
2660 :
2661 : static GEN
2662 63 : opt_param(GEN x, long prec)
2663 : {
2664 63 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
2665 : }
2666 :
2667 : static GEN
2668 63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
2669 : {
2670 63 : pari_sp av = avma;
2671 63 : long lM = lg(M)-1, nbcol = lM;
2672 63 : long tbs = maxss(1, expu(nbg/expi(m)));
2673 : for (;;)
2674 14 : {
2675 77 : GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
2676 : GEN tab;
2677 77 : long i, f=0;
2678 77 : long l = lg(K), lm = lgefint(m);
2679 77 : GEN idx = diviiexact(subiu(p,1),m), g;
2680 : pari_timer ti;
2681 77 : if (DEBUGLEVEL) timer_start(&ti);
2682 154 : for(i=1; i<l; i++)
2683 154 : if (signe(gel(K,i)))
2684 77 : break;
2685 77 : g = Fp_pow(utoi(i), idx, p);
2686 77 : tab = Fp_pow_init(g, p, tbs, p);
2687 77 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
2688 128464 : for(i=1; i<l; i++)
2689 : {
2690 128387 : GEN k = gel(K,i);
2691 128387 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
2692 128387 : if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
2693 76391 : gel(K,i) = cgetineg(lm);
2694 : else
2695 51996 : f++;
2696 : }
2697 77 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
2698 77 : if(f > (nbg>>1)) return gerepileupto(av, K);
2699 4585 : for(i=1; i<=nbcol; i++)
2700 : {
2701 4571 : long a = 1+random_Fl(lM);
2702 4571 : swap(gel(M,a),gel(M,i));
2703 : }
2704 14 : if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
2705 : }
2706 : }
2707 :
2708 : static GEN
2709 126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
2710 : {
2711 126 : pari_sp av=avma;
2712 126 : GEN aa = gen_1;
2713 126 : long AV = 0;
2714 : for(;;)
2715 0 : {
2716 126 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
2717 126 : GEN F = gel(A,1), E = gel(A,2);
2718 126 : GEN Ao = gen_0;
2719 126 : long i, l = lg(F);
2720 962 : for(i=1; i<l; i++)
2721 : {
2722 836 : GEN Ki = gel(K,F[i]);
2723 836 : if (signe(Ki)<0) break;
2724 836 : Ao = addii(Ao, mulis(Ki, E[i]));
2725 : }
2726 126 : if (i==l) return Fp_divu(Ao, AV, m);
2727 0 : aa = gerepileuptoint(av, aa);
2728 : }
2729 : }
2730 :
2731 : static GEN
2732 63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
2733 : {
2734 63 : pari_sp av = avma, av2;
2735 63 : long i, j, nbi, nbr = 0, nbrow, nbg;
2736 : GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
2737 : pari_timer ti;
2738 : struct Fp_log_rel r;
2739 63 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
2740 63 : ulong bnd = 4*bnds;
2741 63 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
2742 :
2743 63 : p_1 = subiu(p,1);
2744 63 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
2745 0 : m = diviiexact(p_1, Z_ppo(p_1, m));
2746 63 : pr = primes_upto_zv(bnd);
2747 63 : nbi = lg(pr)-1;
2748 63 : C = sqrtremi(p, &c);
2749 63 : av2 = avma;
2750 12796 : for (i = 1; i <= nbi; ++i)
2751 : {
2752 12733 : ulong lp = pr[i];
2753 26894 : while (lp <= bnd)
2754 : {
2755 14161 : nbr++;
2756 14161 : lp *= pr[i];
2757 : }
2758 : }
2759 63 : pi = cgetg(nbr+1,t_VECSMALL);
2760 63 : Ci = cgetg(nbr+1,t_VECSMALL);
2761 63 : ci = cgetg(nbr+1,t_VECSMALL);
2762 63 : sz = cgetg(nbr+1,t_VECSMALL);
2763 12796 : for (i = 1, j = 1; i <= nbi; ++i)
2764 : {
2765 12733 : ulong lp = pr[i], sp = expu(2*lp-1);
2766 26894 : while (lp <= bnd)
2767 : {
2768 14161 : pi[j] = lp;
2769 14161 : Ci[j] = umodiu(C, lp);
2770 14161 : ci[j] = umodiu(c, lp);
2771 14161 : sz[j] = sp;
2772 14161 : lp *= pr[i];
2773 14161 : j++;
2774 : }
2775 : }
2776 63 : r.nbrel = 0;
2777 63 : r.nbgen = compute_nbgen(C, bnd, nbi);
2778 63 : r.nbmax = 2*(nbi+r.nbgen);
2779 63 : r.rel = cgetg(r.nbmax+1,t_VEC);
2780 63 : r.prmax = pr[nbi];
2781 63 : if (DEBUGLEVEL)
2782 : {
2783 0 : err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
2784 0 : timer_start(&ti);
2785 : }
2786 63 : nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
2787 63 : nbrow = r.prmax + nbg;
2788 63 : if (DEBUGLEVEL)
2789 : {
2790 0 : err_printf("\n");
2791 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
2792 : }
2793 63 : setlg(r.rel,r.nbrel+1);
2794 63 : r.rel = gerepilecopy(av2, r.rel);
2795 63 : K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
2796 63 : if (DEBUGLEVEL) timer_start(&ti);
2797 63 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
2798 63 : if (DEBUGLEVEL) timer_printf(&ti," log element");
2799 63 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
2800 63 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
2801 63 : d = gcdii(Ao,Bo);
2802 63 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
2803 63 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
2804 63 : return gerepileuptoint(av, l);
2805 : }
2806 :
2807 : static int
2808 4472902 : Fp_log_use_index(long e, long p)
2809 : {
2810 4472902 : return (e >= 27 && 20*(p+6)<=e*e);
2811 : }
2812 :
2813 : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
2814 : static GEN
2815 8258051 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
2816 : {
2817 8258051 : pari_sp av = avma;
2818 8258051 : GEN p = (GEN)E;
2819 : /* assume a reduced mod p, p not necessarily prime */
2820 8258051 : if (equali1(a)) return gen_0;
2821 : /* p > 2 */
2822 5264538 : if (equalii(subiu(p,1), a)) /* -1 */
2823 : {
2824 : pari_sp av2;
2825 : GEN t;
2826 1292351 : ord = get_arith_Z(ord);
2827 1292351 : if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
2828 1292337 : t = shifti(ord,-1); /* only possible solution */
2829 1292337 : av2 = avma;
2830 1292337 : if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
2831 1292309 : set_avma(av2); return gerepileuptoint(av, t);
2832 : }
2833 3972187 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
2834 63 : return Fp_log_index(a, g, ord, p);
2835 3972124 : return gc_NULL(av); /* not easy */
2836 : }
2837 :
2838 : GEN
2839 3803343 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
2840 : {
2841 3803343 : GEN v = get_arith_ZZM(ord);
2842 3803317 : GEN F = gmael(v,2,1);
2843 3803317 : long lF = lg(F)-1, lmax;
2844 3803317 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
2845 3803289 : lmax = expi(gel(F,lF));
2846 3803291 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
2847 91 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
2848 3803290 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
2849 : }
2850 :
2851 : static ulong
2852 125870 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
2853 : {
2854 125870 : ulong i, h=1;
2855 349360 : for(i=0; i<ord; i++, h = Fl_mul(h, g, p))
2856 349361 : if(a==h) return i;
2857 0 : return ~0UL;
2858 : }
2859 :
2860 : static ulong
2861 20260 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
2862 : {
2863 20260 : ulong i, h=1;
2864 51007 : for(i=0; i<ord; i++, h = Fl_mul_pre(h, g, p, pi))
2865 51007 : if(a==h) return i;
2866 0 : return ~0UL;
2867 : }
2868 :
2869 : static ulong
2870 0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
2871 : {
2872 0 : pari_sp av = avma;
2873 0 : GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
2874 0 : return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
2875 : }
2876 :
2877 : ulong
2878 20260 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
2879 : {
2880 20260 : if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
2881 0 : return Fl_log_Fp(a, g, ord, p);
2882 : }
2883 :
2884 : ulong
2885 125870 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
2886 : {
2887 125870 : if (ord <= 200)
2888 0 : return (p&HIGHMASK) ? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
2889 125870 : : Fl_log_naive(a, g, ord, p);
2890 0 : return Fl_log_Fp(a, g, ord, p);
2891 : }
2892 :
2893 : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
2894 : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
2895 : static GEN
2896 126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
2897 : {
2898 126 : long l = lg(P) - 1, e = E[l];
2899 126 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
2900 : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
2901 :
2902 126 : if (l == 1) {
2903 98 : hpe = h;
2904 98 : gpe = g;
2905 : } else {
2906 28 : hpe = modii(h, pe);
2907 28 : gpe = modii(g, pe);
2908 : }
2909 126 : if (e == 1) {
2910 42 : hp = hpe;
2911 42 : gp = gpe;
2912 : } else {
2913 84 : hp = remii(hpe, p);
2914 84 : gp = remii(gpe, p);
2915 : }
2916 126 : if (hp == gen_0 || gp == gen_0) return NULL;
2917 105 : if (absequaliu(p, 2))
2918 : {
2919 35 : GEN n = int2n(e);
2920 35 : ogpe = Zp_order(gpe, gen_2, e, n);
2921 35 : a = Fp_log(hpe, gpe, ogpe, n);
2922 35 : if (typ(a) != t_INT) return NULL;
2923 : }
2924 : else
2925 : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
2926 : is trivial */
2927 : /* [order(gp), factor(order(gp))] */
2928 70 : GEN v = Fp_factored_order(gp, subiu(p,1), p);
2929 70 : GEN ogp = gel(v,1);
2930 70 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
2931 70 : a = Fp_log(hp, gp, v, p);
2932 70 : if (typ(a) != t_INT) return NULL;
2933 70 : if (e == 1) ogpe = ogp;
2934 : else
2935 : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
2936 : /* use p-adic log: O(log p + e) mul*/
2937 : long vpogpe, vpohpe;
2938 :
2939 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
2940 28 : gpe = Fp_pow(gpe, ogp, pe);
2941 : /* g,h = 1 mod p; compute b s.t. h = g^b */
2942 :
2943 : /* v_p(order g mod pe) */
2944 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
2945 : /* v_p(order h mod pe) */
2946 28 : vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
2947 28 : if (vpohpe > vpogpe) return NULL;
2948 :
2949 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
2950 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
2951 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
2952 28 : a = addii(a, mulii(ogp, padic_to_Q(b)));
2953 : }
2954 : }
2955 : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
2956 91 : if (l == 1) return a;
2957 :
2958 28 : N = diviiexact(N, pe); /* make N coprime to p */
2959 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
2960 28 : g = Fp_pow(g, modii(ogpe, phi), N);
2961 28 : setlg(P, l); /* remove last element */
2962 28 : setlg(E, l);
2963 28 : b = znlog_rec(h, g, N, P, E, PHI);
2964 28 : if (!b) return NULL;
2965 28 : return addmulii(a, b, ogpe);
2966 : }
2967 :
2968 : static GEN
2969 98 : get_PHI(GEN P, GEN E)
2970 : {
2971 98 : long i, l = lg(P);
2972 98 : GEN PHI = cgetg(l, t_VEC);
2973 98 : gel(PHI,1) = gen_1;
2974 126 : for (i=1; i<l-1; i++)
2975 : {
2976 28 : GEN t, p = gel(P,i);
2977 28 : long e = E[i];
2978 28 : t = mulii(powiu(p, e-1), subiu(p,1));
2979 28 : if (i > 1) t = mulii(t, gel(PHI,i));
2980 28 : gel(PHI,i+1) = t;
2981 : }
2982 98 : return PHI;
2983 : }
2984 :
2985 : GEN
2986 238 : znlog(GEN h, GEN g, GEN o)
2987 : {
2988 238 : pari_sp av = avma;
2989 : GEN N, fa, P, E, x;
2990 238 : switch (typ(g))
2991 : {
2992 28 : case t_PADIC:
2993 : {
2994 28 : GEN p = gel(g,2);
2995 28 : long v = valp(g);
2996 28 : if (v < 0) pari_err_DIM("znlog");
2997 28 : if (v > 0) {
2998 0 : long k = gvaluation(h, p);
2999 0 : if (k % v) return cgetg(1,t_VEC);
3000 0 : k /= v;
3001 0 : if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
3002 0 : set_avma(av); return stoi(k);
3003 : }
3004 28 : N = gel(g,3);
3005 28 : g = Rg_to_Fp(g, N);
3006 28 : break;
3007 : }
3008 203 : case t_INTMOD:
3009 203 : N = gel(g,1);
3010 203 : g = gel(g,2); break;
3011 7 : default: pari_err_TYPE("znlog", g);
3012 : return NULL; /* LCOV_EXCL_LINE */
3013 : }
3014 231 : if (equali1(N)) { set_avma(av); return gen_0; }
3015 231 : h = Rg_to_Fp(h, N);
3016 224 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
3017 98 : fa = Z_factor(N);
3018 98 : P = gel(fa,1);
3019 98 : E = vec_to_vecsmall(gel(fa,2));
3020 98 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
3021 98 : if (!x) { set_avma(av); return cgetg(1,t_VEC); }
3022 63 : return gerepileuptoint(av, x);
3023 : }
3024 :
3025 : GEN
3026 62790 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
3027 : {
3028 62790 : if (lgefint(p)==3)
3029 : {
3030 62370 : long nn = itos_or_0(n);
3031 62370 : if (nn)
3032 : {
3033 62370 : ulong pp = p[2];
3034 : ulong uz;
3035 62370 : ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
3036 62349 : if (r==ULONG_MAX) return NULL;
3037 62293 : if (zeta) *zeta = utoi(uz);
3038 62293 : return utoi(r);
3039 : }
3040 : }
3041 420 : a = modii(a,p);
3042 420 : if (!signe(a))
3043 : {
3044 0 : if (zeta) *zeta = gen_1;
3045 0 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
3046 0 : return gen_0;
3047 : }
3048 420 : if (absequaliu(n,2))
3049 : {
3050 224 : if (zeta) *zeta = subiu(p,1);
3051 224 : return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
3052 : }
3053 196 : return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
3054 : }
3055 :
3056 : /*********************************************************************/
3057 : /** FACTORIAL **/
3058 : /*********************************************************************/
3059 : GEN
3060 75484 : mulu_interval_step(ulong a, ulong b, ulong step)
3061 : {
3062 75484 : pari_sp av = avma;
3063 : ulong k, l, N, n;
3064 : long lx;
3065 : GEN x;
3066 :
3067 75484 : if (!a) return gen_0;
3068 75484 : if (step == 1) return mulu_interval(a, b);
3069 75484 : n = 1 + (b-a) / step;
3070 75484 : b -= (b-a) % step;
3071 75484 : if (n < 61)
3072 : {
3073 74035 : if (n == 1) return utoipos(a);
3074 57645 : x = muluu(a,a+step); if (n == 2) return x;
3075 497071 : for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
3076 46029 : return gerepileuptoint(av, x);
3077 : }
3078 : /* step | b-a */
3079 1449 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3080 1449 : N = b + a;
3081 1449 : for (k = a;; k += step)
3082 : {
3083 230195 : l = N - k; if (l <= k) break;
3084 228746 : gel(x,lx++) = muluu(k,l);
3085 : }
3086 1449 : if (l == k) gel(x,lx++) = utoipos(k);
3087 1449 : setlg(x, lx);
3088 1449 : return gerepileuptoint(av, ZV_prod(x));
3089 : }
3090 : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
3091 : * first is slower ... ] */
3092 : GEN
3093 182397 : mulu_interval(ulong a, ulong b)
3094 : {
3095 182397 : pari_sp av = avma;
3096 : ulong k, l, N, n;
3097 : long lx;
3098 : GEN x;
3099 :
3100 182397 : if (!a) return gen_0;
3101 182397 : n = b - a + 1;
3102 182397 : if (n < 61)
3103 : {
3104 182383 : if (n == 1) return utoipos(a);
3105 126103 : x = muluu(a,a+1); if (n == 2) return x;
3106 371697 : for (k=a+2; k<b; k++) x = mului(k,x);
3107 : /* avoid k <= b: broken if b = ULONG_MAX */
3108 92111 : return gerepileuptoint(av, mului(b,x));
3109 : }
3110 14 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3111 14 : N = b + a;
3112 14 : for (k = a;; k++)
3113 : {
3114 7007 : l = N - k; if (l <= k) break;
3115 6993 : gel(x,lx++) = muluu(k,l);
3116 : }
3117 14 : if (l == k) gel(x,lx++) = utoipos(k);
3118 14 : setlg(x, lx);
3119 14 : return gerepileuptoint(av, ZV_prod(x));
3120 : }
3121 : GEN
3122 588 : muls_interval(long a, long b)
3123 : {
3124 588 : pari_sp av = avma;
3125 588 : long lx, k, l, N, n = b - a + 1;
3126 : GEN x;
3127 :
3128 588 : if (a <= 0 && b >= 0) return gen_0;
3129 315 : if (n < 61)
3130 : {
3131 315 : x = stoi(a);
3132 553 : for (k=a+1; k<=b; k++) x = mulsi(k,x);
3133 315 : return gerepileuptoint(av, x);
3134 : }
3135 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3136 0 : N = b + a;
3137 0 : for (k = a;; k++)
3138 : {
3139 0 : l = N - k; if (l <= k) break;
3140 0 : gel(x,lx++) = mulss(k,l);
3141 : }
3142 0 : if (l == k) gel(x,lx++) = stoi(k);
3143 0 : setlg(x, lx);
3144 0 : return gerepileuptoint(av, ZV_prod(x));
3145 : }
3146 :
3147 : GEN
3148 495888 : mpfact(long n)
3149 : {
3150 495888 : pari_sp av = avma;
3151 : GEN a, v;
3152 : long k;
3153 495888 : if (n <= 12) switch(n)
3154 : {
3155 440002 : case 0: case 1: return gen_1;
3156 34554 : case 2: return gen_2;
3157 1304 : case 3: return utoipos(6);
3158 2133 : case 4: return utoipos(24);
3159 780 : case 5: return utoipos(120);
3160 486 : case 6: return utoipos(720);
3161 387 : case 7: return utoipos(5040);
3162 388 : case 8: return utoipos(40320);
3163 416 : case 9: return utoipos(362880);
3164 661 : case 10:return utoipos(3628800);
3165 249 : case 11:return utoipos(39916800);
3166 236 : case 12:return utoipos(479001600);
3167 0 : default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
3168 : }
3169 14292 : v = cgetg(expu(n) + 2, t_VEC);
3170 14292 : for (k = 1;; k++)
3171 71697 : {
3172 85989 : long m = n >> (k-1), l;
3173 85989 : if (m <= 2) break;
3174 71697 : l = (1 + (n >> k)) | 1;
3175 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3176 71697 : a = mulu_interval_step(l, m, 2);
3177 71696 : gel(v,k) = k == 1? a: powiu(a, k);
3178 : }
3179 71697 : a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
3180 14292 : a = shifti(a, factorial_lval(n, 2));
3181 14292 : return gerepileuptoint(av, a);
3182 : }
3183 :
3184 : ulong
3185 49897 : factorial_Fl(long n, ulong p)
3186 : {
3187 : long k;
3188 : ulong v;
3189 49897 : if (p <= (ulong)n) return 0;
3190 49897 : v = Fl_powu(2, factorial_lval(n, 2), p);
3191 49966 : for (k = 1;; k++)
3192 128676 : {
3193 178642 : long m = n >> (k-1), l, i;
3194 178642 : ulong a = 1;
3195 178642 : if (m <= 2) break;
3196 128712 : l = (1 + (n >> k)) | 1;
3197 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3198 748448 : for (i=l; i<=m; i+=2)
3199 619761 : a = Fl_mul(a, i, p);
3200 128687 : v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
3201 : }
3202 49930 : return v;
3203 : }
3204 :
3205 : GEN
3206 60 : factorial_Fp(long n, GEN p)
3207 : {
3208 60 : pari_sp av = avma;
3209 : long k;
3210 60 : GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
3211 60 : for (k = 1;; k++)
3212 134 : {
3213 194 : long m = n >> (k-1), l, i;
3214 194 : GEN a = gen_1;
3215 194 : if (m <= 2) break;
3216 134 : l = (1 + (n >> k)) | 1;
3217 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3218 402 : for (i=l; i<=m; i+=2)
3219 268 : a = Fp_mulu(a, i, p);
3220 134 : v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
3221 134 : v = gerepileuptoint(av, v);
3222 : }
3223 60 : return v;
3224 : }
3225 :
3226 : /*******************************************************************/
3227 : /** LUCAS & FIBONACCI **/
3228 : /*******************************************************************/
3229 : static void
3230 56 : lucas(ulong n, GEN *a, GEN *b)
3231 : {
3232 : GEN z, t, zt;
3233 56 : if (!n) { *a = gen_2; *b = gen_1; return; }
3234 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
3235 49 : switch(n & 3) {
3236 14 : case 0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
3237 14 : case 1: *a = subiu(zt,1); *b = addiu(sqri(t),2); break;
3238 7 : case 2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
3239 14 : case 3: *a = addiu(zt,1); *b = subiu(sqri(t),2);
3240 : }
3241 49 : }
3242 :
3243 : GEN
3244 7 : fibo(long n)
3245 : {
3246 7 : pari_sp av = avma;
3247 : GEN a, b;
3248 7 : if (!n) return gen_0;
3249 7 : lucas((ulong)(labs(n)-1), &a, &b);
3250 7 : a = diviuexact(addii(shifti(a,1),b), 5);
3251 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
3252 7 : return gerepileuptoint(av, a);
3253 : }
3254 :
3255 : /*******************************************************************/
3256 : /* CONTINUED FRACTIONS */
3257 : /*******************************************************************/
3258 : static GEN
3259 2830727 : icopy_lg(GEN x, long l)
3260 : {
3261 2830727 : long lx = lgefint(x);
3262 : GEN y;
3263 :
3264 2830727 : if (lx >= l) return icopy(x);
3265 35 : y = cgeti(l); affii(x, y); return y;
3266 : }
3267 :
3268 : /* continued fraction of a/b. If y != NULL, stop when partial quotients
3269 : * differ from y */
3270 : static GEN
3271 2831060 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
3272 : {
3273 : GEN z, c;
3274 2831060 : ulong i, l, ly = lgefint(b);
3275 :
3276 : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
3277 2831060 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
3278 2831060 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
3279 2831060 : if (l > LGBITS) l = LGBITS;
3280 :
3281 2831060 : z = cgetg(l,t_VEC);
3282 2831060 : l--;
3283 2831060 : if (y) {
3284 333 : pari_sp av = avma;
3285 333 : if (l >= (ulong)lg(y)) l = lg(y)-1;
3286 27537 : for (i = 1; i <= l; i++)
3287 : {
3288 27327 : GEN q = gel(y,i);
3289 27327 : gel(z,i) = q;
3290 27327 : c = b; if (!gequal1(q)) c = mulii(q, b);
3291 27327 : c = subii(a, c);
3292 27327 : if (signe(c) < 0)
3293 : { /* partial quotient too large */
3294 117 : c = addii(c, b);
3295 117 : if (signe(c) >= 0) i++; /* by 1 */
3296 117 : break;
3297 : }
3298 27210 : if (cmpii(c, b) >= 0)
3299 : { /* partial quotient too small */
3300 6 : c = subii(c, b);
3301 6 : if (cmpii(c, b) < 0) {
3302 : /* by 1. If next quotient is 1 in y, add 1 */
3303 6 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
3304 6 : i++;
3305 : }
3306 6 : break;
3307 : }
3308 27204 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
3309 27204 : a = b; b = c;
3310 : }
3311 : } else {
3312 2830727 : a = icopy_lg(a, ly);
3313 2830727 : b = icopy(b);
3314 23450113 : for (i = 1; i <= l; i++)
3315 : {
3316 23449804 : gel(z,i) = truedvmdii(a,b,&c);
3317 23449804 : if (c == gen_0) { i++; break; }
3318 20619386 : affii(c, a); cgiv(c); c = a;
3319 20619386 : a = b; b = c;
3320 : }
3321 : }
3322 2831060 : i--;
3323 2831060 : if (i > 1 && gequal1(gel(z,i)))
3324 : {
3325 100 : cgiv(gel(z,i)); --i;
3326 100 : gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
3327 : }
3328 2831060 : setlg(z,i+1); return z;
3329 : }
3330 :
3331 : static GEN
3332 0 : sersfcont(GEN a, GEN b, long k)
3333 : {
3334 0 : long i, l = typ(a) == t_POL? lg(a): 3;
3335 : GEN y, c;
3336 0 : if (lg(b) > l) l = lg(b);
3337 0 : if (k > 0 && l > k+1) l = k+1;
3338 0 : y = cgetg(l,t_VEC);
3339 0 : for (i=1; i<l; i++)
3340 : {
3341 0 : gel(y,i) = poldivrem(a,b,&c);
3342 0 : if (gequal0(c)) { i++; break; }
3343 0 : a = b; b = c;
3344 : }
3345 0 : setlg(y, i); return y;
3346 : }
3347 :
3348 : GEN
3349 2831742 : gboundcf(GEN x, long k)
3350 : {
3351 : pari_sp av;
3352 2831742 : long tx = typ(x), e;
3353 : GEN y, a, b, c;
3354 :
3355 2831742 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
3356 2831735 : if (is_scalar_t(tx))
3357 : {
3358 2831735 : if (gequal0(x)) return mkvec(gen_0);
3359 2831630 : switch(tx)
3360 : {
3361 896 : case t_INT: return mkveccopy(x);
3362 340 : case t_REAL:
3363 340 : av = avma;
3364 340 : c = mantissa_real(x,&e);
3365 340 : if (e < 0) pari_err_PREC("gboundcf");
3366 333 : y = int2n(e);
3367 333 : a = Qsfcont(c,y, NULL, k);
3368 333 : b = addsi(signe(x), c);
3369 333 : return gerepilecopy(av, Qsfcont(b,y, a, k));
3370 :
3371 2830394 : case t_FRAC:
3372 2830394 : av = avma;
3373 2830394 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
3374 : }
3375 0 : pari_err_TYPE("gboundcf",x);
3376 : }
3377 :
3378 0 : switch(tx)
3379 : {
3380 0 : case t_POL: return mkveccopy(x);
3381 0 : case t_SER:
3382 0 : av = avma;
3383 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
3384 0 : case t_RFRAC:
3385 0 : av = avma;
3386 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
3387 : }
3388 0 : pari_err_TYPE("gboundcf",x);
3389 : return NULL; /* LCOV_EXCL_LINE */
3390 : }
3391 :
3392 : static GEN
3393 14 : sfcont2(GEN b, GEN x, long k)
3394 : {
3395 14 : pari_sp av = avma;
3396 14 : long lb = lg(b), tx = typ(x), i;
3397 : GEN y,p1;
3398 :
3399 14 : if (k)
3400 : {
3401 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
3402 0 : lb = k+1;
3403 : }
3404 7 : y = cgetg(lb,t_VEC);
3405 7 : if (lb==1) return y;
3406 7 : if (is_scalar_t(tx))
3407 : {
3408 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
3409 : }
3410 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
3411 :
3412 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
3413 7 : for (i = 1;;)
3414 : {
3415 35 : if (tx == t_REAL)
3416 : {
3417 35 : long e = expo(x);
3418 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
3419 35 : gel(y,i) = floorr(x);
3420 35 : p1 = subri(x, gel(y,i));
3421 : }
3422 : else
3423 : {
3424 0 : gel(y,i) = gfloor(x);
3425 0 : p1 = gsub(x, gel(y,i));
3426 : }
3427 35 : if (++i >= lb) break;
3428 28 : if (gequal0(p1)) break;
3429 28 : x = gdiv(gel(b,i),p1);
3430 : }
3431 7 : setlg(y,i);
3432 7 : return gerepilecopy(av,y);
3433 : }
3434 :
3435 : GEN
3436 112 : gcf(GEN x) { return gboundcf(x,0); }
3437 : GEN
3438 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
3439 : GEN
3440 49 : contfrac0(GEN x, GEN b, long nmax)
3441 : {
3442 : long tb;
3443 :
3444 49 : if (!b) return gboundcf(x,nmax);
3445 28 : tb = typ(b);
3446 28 : if (tb == t_INT) return gboundcf(x,itos(b));
3447 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
3448 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
3449 14 : return sfcont2(b,x,nmax);
3450 : }
3451 :
3452 : GEN
3453 252 : contfracpnqn(GEN x, long n)
3454 : {
3455 252 : pari_sp av = avma;
3456 252 : long i, lx = lg(x);
3457 : GEN M,A,B, p0,p1, q0,q1;
3458 :
3459 252 : if (lx == 1)
3460 : {
3461 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
3462 21 : if (n >= 0) return cgetg(1,t_MAT);
3463 7 : return matid(2);
3464 : }
3465 224 : switch(typ(x))
3466 : {
3467 182 : case t_VEC: case t_COL: A = x; B = NULL; break;
3468 42 : case t_MAT:
3469 42 : switch(lgcols(x))
3470 : {
3471 0 : case 2: A = row(x,1); B = NULL; break;
3472 35 : case 3: A = row(x,2); B = row(x,1); break;
3473 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
3474 : return NULL; /*LCOV_EXCL_LINE*/
3475 : }
3476 35 : break;
3477 0 : default: pari_err_TYPE("pnqn",x);
3478 : return NULL; /*LCOV_EXCL_LINE*/
3479 : }
3480 217 : p1 = gel(A,1);
3481 217 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
3482 217 : if (n >= 0)
3483 : {
3484 182 : lx = minss(lx, n+2);
3485 182 : if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
3486 : }
3487 35 : else if (lx == 2)
3488 7 : return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
3489 : /* lx >= 3 */
3490 119 : p0 = gen_1;
3491 119 : q0 = gen_0; /* p[-1], q[-1] */
3492 119 : M = cgetg(lx, t_MAT);
3493 119 : gel(M,1) = mkcol2(p1,q1);
3494 399 : for (i=2; i<lx; i++)
3495 : {
3496 280 : GEN a = gel(A,i), p2,q2;
3497 280 : if (B) {
3498 84 : GEN b = gel(B,i);
3499 84 : p0 = gmul(b,p0);
3500 84 : q0 = gmul(b,q0);
3501 : }
3502 280 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
3503 280 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
3504 280 : gel(M,i) = mkcol2(p1,q1);
3505 : }
3506 119 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
3507 119 : return gerepilecopy(av, M);
3508 : }
3509 : GEN
3510 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
3511 : /* x = [a0, ..., an] from gboundcf, n >= 0;
3512 : * return [[p0, ..., pn], [q0,...,qn]] */
3513 : GEN
3514 609308 : ZV_allpnqn(GEN x)
3515 : {
3516 609308 : long i, lx = lg(x);
3517 609308 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
3518 :
3519 609308 : gel(v,1) = P = cgetg(lx, t_VEC);
3520 609308 : gel(v,2) = Q = cgetg(lx, t_VEC);
3521 609308 : p0 = gen_1; q0 = gen_0;
3522 609308 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
3523 2092209 : for (i=2; i<lx; i++)
3524 : {
3525 1482901 : GEN a = gel(x,i);
3526 1482901 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
3527 1482901 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
3528 : }
3529 609308 : return v;
3530 : }
3531 :
3532 : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
3533 : static GEN
3534 42 : mod_to_frac(GEN x, GEN N, GEN B)
3535 : {
3536 : GEN a, b, A;
3537 42 : if (B) A = divii(shifti(N, -1), B);
3538 : else
3539 : {
3540 14 : A = sqrti(shifti(N, -1));
3541 14 : B = A;
3542 : }
3543 42 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
3544 28 : return equali1(b)? a: mkfrac(a,b);
3545 : }
3546 :
3547 : static GEN
3548 70 : mod_to_rfrac(GEN x, GEN N, long B)
3549 : {
3550 : GEN a, b;
3551 70 : long A, d = degpol(N);
3552 70 : if (B >= 0) A = d-1 - B;
3553 : else
3554 : {
3555 42 : B = d >> 1;
3556 42 : A = odd(d)? B : B-1;
3557 : }
3558 70 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
3559 70 : if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
3560 56 : return gdiv(a,b);
3561 : }
3562 :
3563 : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
3564 : * of the continued fraction of x with b <= k maximal */
3565 : static GEN
3566 7 : bestappr_frac(GEN x, GEN k)
3567 : {
3568 : pari_sp av;
3569 : GEN p0, p1, p, q0, q1, q, a, y;
3570 :
3571 7 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
3572 0 : av = avma; y = x;
3573 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
3574 0 : q1 = gen_0; q0 = gen_1;
3575 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
3576 : for(;;)
3577 : {
3578 0 : x = ginv(x); /* > 1 */
3579 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
3580 0 : if (cmpii(a,k) > 0)
3581 : { /* next partial quotient will overflow limits */
3582 : GEN n, d;
3583 0 : a = divii(subii(k, q1), q0);
3584 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3585 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3586 : /* compare |y-p0/q0|, |y-p1/q1| */
3587 0 : n = gel(y,1);
3588 0 : d = gel(y,2);
3589 0 : if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
3590 : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
3591 0 : { p1 = p0; q1 = q0; }
3592 0 : break;
3593 : }
3594 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3595 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3596 :
3597 0 : if (cmpii(q0,k) > 0) break;
3598 0 : x = gsub(x,a); /* 0 <= x < 1 */
3599 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
3600 :
3601 : }
3602 0 : return gerepileupto(av, gdiv(p1,q1));
3603 : }
3604 : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
3605 : * of the continued fraction of x with b <= k maximal */
3606 : static GEN
3607 1428430 : bestappr_real(GEN x, GEN k)
3608 : {
3609 1428430 : pari_sp av = avma;
3610 1428430 : GEN kr, p0, p1, p, q0, q1, q, a, y = x;
3611 :
3612 1428430 : p1 = gen_1; a = p0 = floorr(x);
3613 1428333 : q1 = gen_0; q0 = gen_1;
3614 1428333 : x = subri(x,a); /* 0 <= x < 1 */
3615 1428325 : if (!signe(x)) { cgiv(x); return a; }
3616 1324976 : kr = itor(k, realprec(x));
3617 : for(;;)
3618 1471032 : {
3619 : long d;
3620 2796033 : x = invr(x); /* > 1 */
3621 2795915 : if (cmprr(x,kr) > 0)
3622 : { /* next partial quotient will overflow limits */
3623 1312453 : a = divii(subii(k, q1), q0);
3624 1312371 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3625 1312377 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3626 : /* compare |y-p0/q0|, |y-p1/q1| */
3627 1312314 : if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
3628 : mulir(q0, subri(mulir(q1,y), p1))) < 0)
3629 124576 : { p1 = p0; q1 = q0; }
3630 1312457 : break;
3631 : }
3632 1483583 : d = nbits2prec(expo(x) + 1);
3633 1483583 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
3634 :
3635 1483388 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
3636 1483351 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3637 1483371 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3638 :
3639 1483373 : if (cmpii(q0,k) > 0) break;
3640 1477460 : x = subri(x,a); /* 0 <= x < 1 */
3641 1477463 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
3642 : }
3643 1324999 : if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
3644 1324999 : return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
3645 : }
3646 :
3647 : /* k t_INT or NULL */
3648 : static GEN
3649 2377894 : bestappr_Q(GEN x, GEN k)
3650 : {
3651 2377894 : long lx, tx = typ(x), i;
3652 : GEN a, y;
3653 :
3654 2377894 : switch(tx)
3655 : {
3656 77 : case t_INT: return icopy(x);
3657 7 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
3658 1669001 : case t_REAL:
3659 1669001 : if (!signe(x)) return gen_0;
3660 : /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
3661 1428415 : i = bit_prec(x); if (i <= expo(x)) return NULL;
3662 1428430 : return bestappr_real(x, k? k: int2n(i));
3663 :
3664 28 : case t_INTMOD: {
3665 28 : pari_sp av = avma;
3666 28 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
3667 21 : return gerepilecopy(av, a);
3668 : }
3669 14 : case t_PADIC: {
3670 14 : pari_sp av = avma;
3671 14 : long v = valp(x);
3672 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
3673 7 : if (v) a = gmul(a, powis(gel(x,2), v));
3674 7 : return gerepilecopy(av, a);
3675 : }
3676 :
3677 322 : case t_COMPLEX: {
3678 322 : pari_sp av = avma;
3679 322 : y = cgetg(3, t_COMPLEX);
3680 322 : gel(y,2) = bestappr(gel(x,2), k);
3681 322 : gel(y,1) = bestappr(gel(x,1), k);
3682 322 : if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
3683 91 : return y;
3684 : }
3685 0 : case t_SER:
3686 0 : if (ser_isexactzero(x)) return gcopy(x);
3687 : /* fall through */
3688 : case t_POLMOD: case t_POL: case t_RFRAC:
3689 : case t_VEC: case t_COL: case t_MAT:
3690 708445 : y = cgetg_copy(x, &lx);
3691 708480 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
3692 3020145 : for (; i<lx; i++)
3693 : {
3694 2311678 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
3695 2311665 : gel(y,i) = a;
3696 : }
3697 708467 : if (tx == t_POL) return normalizepol(y);
3698 708453 : if (tx == t_SER) return normalizeser(y);
3699 708453 : return y;
3700 : }
3701 0 : pari_err_TYPE("bestappr_Q",x);
3702 : return NULL; /* LCOV_EXCL_LINE */
3703 : }
3704 :
3705 : static GEN
3706 56 : bestappr_ser(GEN x, long B)
3707 : {
3708 56 : long dN, v = valp(x), lx = lg(x);
3709 : GEN t;
3710 56 : x = normalizepol(ser2pol_i(x, lx));
3711 56 : dN = lx-2;
3712 56 : if (v > 0)
3713 : {
3714 14 : x = RgX_shift_shallow(x, v);
3715 14 : dN += v;
3716 : }
3717 42 : else if (v < 0)
3718 : {
3719 7 : if (B >= 0) B = maxss(B+v, 0);
3720 : }
3721 56 : t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
3722 56 : if (!t) return NULL;
3723 42 : if (v < 0)
3724 : {
3725 : GEN a, b;
3726 : long vx;
3727 7 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
3728 : /* t_RFRAC */
3729 7 : vx = varn(x);
3730 7 : a = gel(t,1);
3731 7 : b = gel(t,2);
3732 7 : v -= RgX_valrem(b, &b);
3733 7 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
3734 7 : if (v < 0) b = RgX_shift(b, -v);
3735 0 : else if (v > 0) {
3736 0 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
3737 0 : a = RgX_shift(a, v);
3738 : }
3739 7 : t = mkrfraccopy(a, b);
3740 : }
3741 42 : return t;
3742 : }
3743 : static GEN bestappr_RgX(GEN x, long B);
3744 : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
3745 : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
3746 : static GEN
3747 77 : bestappr_RgX(GEN x, long B)
3748 : {
3749 77 : long i, lx, tx = typ(x);
3750 : GEN y, t;
3751 77 : switch(tx)
3752 : {
3753 0 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
3754 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
3755 0 : return gcopy(x);
3756 :
3757 14 : case t_RFRAC: {
3758 14 : pari_sp av = avma;
3759 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
3760 7 : x = rfrac_to_ser_i(x, 2*B+1);
3761 7 : t = bestappr_ser(x, B); if (!t) return NULL;
3762 0 : return gerepileupto(av, t);
3763 : }
3764 14 : case t_POLMOD: {
3765 14 : pari_sp av = avma;
3766 14 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
3767 14 : return gerepileupto(av, t);
3768 : }
3769 49 : case t_SER: {
3770 49 : pari_sp av = avma;
3771 49 : t = bestappr_ser(x, B); if (!t) return NULL;
3772 42 : return gerepileupto(av, t);
3773 : }
3774 :
3775 0 : case t_VEC: case t_COL: case t_MAT:
3776 0 : y = cgetg_copy(x, &lx);
3777 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
3778 0 : for (; i<lx; i++)
3779 : {
3780 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
3781 0 : gel(y,i) = t;
3782 : }
3783 0 : return y;
3784 : }
3785 0 : pari_err_TYPE("bestappr_RgX",x);
3786 : return NULL; /* LCOV_EXCL_LINE */
3787 : }
3788 :
3789 : /* allow k = NULL: maximal accuracy */
3790 : GEN
3791 66215 : bestappr(GEN x, GEN k)
3792 : {
3793 66215 : pari_sp av = avma;
3794 66215 : if (k) { /* replace by floor(k) */
3795 65928 : switch(typ(k))
3796 : {
3797 4424 : case t_INT:
3798 4424 : break;
3799 61504 : case t_REAL: case t_FRAC:
3800 61504 : k = floor_safe(k); /* left on stack for efficiency */
3801 61507 : if (!signe(k)) k = gen_1;
3802 61507 : break;
3803 0 : default:
3804 0 : pari_err_TYPE("bestappr [bound type]", k);
3805 0 : break;
3806 : }
3807 287 : }
3808 66218 : x = bestappr_Q(x, k);
3809 66218 : if (!x) { set_avma(av); return cgetg(1,t_VEC); }
3810 66203 : return x;
3811 : }
3812 : GEN
3813 77 : bestapprPade(GEN x, long B)
3814 : {
3815 77 : pari_sp av = avma;
3816 77 : GEN t = bestappr_RgX(x, B);
3817 77 : if (!t) { set_avma(av); return cgetg(1,t_VEC); }
3818 63 : return t;
3819 : }
|
