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 464816 : 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 464816 : 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 4895 : is_gener_expo(GEN p, GEN L0)
39 : {
40 4895 : GEN L, q = shifti(p,-1);
41 : long i, l;
42 4895 : if (L0) {
43 3120 : l = lg(L0);
44 3120 : L = cgetg(l, t_VEC);
45 : } else {
46 1775 : L0 = L = odd_prime_divisors(q);
47 1775 : l = lg(L);
48 : }
49 14171 : for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
50 4895 : return L;
51 : }
52 : static GEN
53 532669 : u_is_gener_expo(ulong p, GEN L0)
54 : {
55 532669 : const ulong q = p >> 1;
56 : long i;
57 : GEN L;
58 532669 : if (!L0) L0 = u_odd_prime_divisors(q);
59 532671 : L = cgetg_copy(L0,&i);
60 1154949 : while (--i) L[i] = q / uel(L0,i);
61 532670 : return L;
62 : }
63 :
64 : int
65 1629492 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
66 : {
67 : long i;
68 1629492 : if (krouu(x, p) >= 0) return 0;
69 1379928 : for (i=lg(L)-1; i; i--)
70 : {
71 842462 : ulong t = Fl_powu(x, uel(L,i), p);
72 842460 : if (t == p_1 || t == 1) return 0;
73 : }
74 537466 : return 1;
75 : }
76 : /* assume p prime */
77 : ulong
78 1054026 : pgener_Fl_local(ulong p, GEN L0)
79 : {
80 1054026 : const pari_sp av = avma;
81 1054026 : const ulong p_1 = p-1;
82 : long x;
83 : GEN L;
84 1054026 : if (p <= 19) switch(p)
85 : { /* quick trivial cases */
86 28 : case 2: return 1;
87 116176 : case 7:
88 116176 : case 17: return 3;
89 405199 : default: return 2;
90 : }
91 532623 : L = u_is_gener_expo(p,L0);
92 1621675 : for (x = 2;; x++)
93 1621675 : if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
94 : }
95 : ulong
96 579292 : 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 13755 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
102 : {
103 13755 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
104 13755 : if (t >= 0) return 0;
105 21797 : for (i = lg(L)-1; i; i--)
106 : {
107 14302 : GEN t = Fp_pow(x, gel(L,i), p);
108 14302 : if (equalii(t, p_1) || equali1(t)) return 0;
109 : }
110 7495 : return 1;
111 : }
112 :
113 : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
114 : GEN
115 357761 : pgener_Fp_local(GEN p, GEN L0)
116 : {
117 357761 : pari_sp av0 = avma;
118 : GEN x, p_1, L;
119 357761 : if (lgefint(p) == 3)
120 : {
121 : ulong z;
122 352872 : if (p[2] == 2) return gen_1;
123 257990 : if (L0) L0 = ZV_to_nv(L0);
124 257985 : z = pgener_Fl_local(uel(p,2), L0);
125 258013 : return gc_utoipos(av0, z);
126 : }
127 4889 : p_1 = subiu(p,1); L = is_gener_expo(p, L0);
128 4890 : x = utoipos(2);
129 9899 : for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
130 4890 : return gc_utoipos(av0, uel(x,2));
131 : }
132 :
133 : GEN
134 44134 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
135 :
136 : ulong
137 203274 : pgener_Zl(ulong p)
138 : {
139 203274 : 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 203274 : if (p == 40487) return 10;
142 : #ifndef LONG_IS_64BIT
143 29460 : return pgener_Fl(p);
144 : #else
145 173814 : 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 102 : 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 168539 : pgener_Zp(GEN p)
162 : {
163 168539 : if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
164 : else
165 : {
166 5 : const pari_sp av = avma;
167 5 : GEN p_1 = subiu(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
168 5 : GEN x = utoipos(2);
169 12 : for (;; x[2]++)
170 17 : if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
171 5 : return gc_utoipos(av, uel(x,2));
172 : }
173 : }
174 :
175 : static GEN
176 259 : gener_Zp(GEN q, GEN F)
177 : {
178 259 : GEN p = NULL;
179 259 : long e = 0;
180 259 : if (F)
181 : {
182 14 : GEN P = gel(F,1), E = gel(F,2);
183 14 : long i, l = lg(P);
184 42 : for (i = 1; i < l; i++)
185 : {
186 28 : p = gel(P,i);
187 28 : if (absequaliu(p, 2)) continue;
188 14 : if (i < l-1) pari_err_DOMAIN("znprimroot", "n","=",F,F);
189 14 : e = itos(gel(E,i));
190 : }
191 14 : if (!p) pari_err_DOMAIN("znprimroot", "n","=",F,F);
192 : }
193 : else
194 245 : e = Z_isanypower(q, &p);
195 259 : if (!BPSW_psp(e? p: q)) pari_err_DOMAIN("znprimroot", "n","=", q,q);
196 245 : return e > 1? pgener_Zp(p): pgener_Fp(q);
197 : }
198 :
199 : GEN
200 329 : znprimroot(GEN N)
201 : {
202 329 : pari_sp av = avma;
203 : GEN x, n, F;
204 :
205 329 : if ((F = check_arith_non0(N,"znprimroot")))
206 : {
207 14 : F = clean_Z_factor(F);
208 14 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
209 : }
210 322 : N = absi_shallow(N);
211 322 : if (abscmpiu(N, 4) <= 0) { set_avma(av); return mkintmodu(N[2]-1,N[2]); }
212 273 : switch(mod4(N))
213 : {
214 14 : case 0: /* N = 0 mod 4 */
215 14 : pari_err_DOMAIN("znprimroot", "n","=",N,N);
216 0 : x = NULL; break;
217 28 : case 2: /* N = 2 mod 4 */
218 28 : n = shifti(N,-1); /* becomes odd */
219 28 : x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
220 21 : break;
221 231 : default: /* N odd */
222 231 : x = gener_Zp(N,F);
223 224 : break;
224 : }
225 245 : return gerepilecopy(av, mkintmod(x, N));
226 : }
227 :
228 : /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
229 : GEN
230 0 : rootsof1_Fp(GEN n, GEN p)
231 : {
232 0 : pari_sp av = avma;
233 0 : GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
234 0 : GEN z = pgener_Fp_local(p, L);
235 0 : z = Fp_pow(z, diviiexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
236 0 : return gerepileuptoint(av, z);
237 : }
238 :
239 : GEN
240 3033 : rootsof1u_Fp(ulong n, GEN p)
241 : {
242 3033 : pari_sp av = avma;
243 3033 : GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
244 3033 : z = pgener_Fp_local(p, Flv_to_ZV(L));
245 3033 : z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
246 3033 : return gerepileuptoint(av, z);
247 : }
248 :
249 : ulong
250 214579 : rootsof1_Fl(ulong n, ulong p)
251 : {
252 214579 : pari_sp av = avma;
253 214579 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
254 214579 : ulong z = pgener_Fl_local(p, L);
255 214579 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
256 214579 : return gc_ulong(av,z);
257 : }
258 :
259 : /*********************************************************************/
260 : /** INVERSE TOTIENT FUNCTION **/
261 : /*********************************************************************/
262 : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
263 : * primes. Return factor(N) */
264 : GEN
265 350651 : Z_factor_listP(GEN N, GEN L)
266 : {
267 350651 : long i, k, l = lg(L);
268 350651 : GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
269 1346688 : for (i = k = 1; i < l; i++)
270 : {
271 996037 : GEN p = gel(L,i);
272 996037 : long v = Z_pvalrem(N, p, &N);
273 996037 : if (v)
274 : {
275 792176 : gel(P,k) = p;
276 792176 : gel(E,k) = utoipos(v);
277 792176 : k++;
278 : }
279 : }
280 350651 : setlg(P, k);
281 350651 : setlg(E, k); return mkmat2(P,E);
282 : }
283 :
284 : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
285 : * L is a list of primes containing all prime divisors of n. */
286 : static long
287 621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
288 : {
289 621565 : pari_sp av = avma, av2;
290 : GEN k, D;
291 : long i, v;
292 621565 : if (m && mod2(n))
293 : {
294 270914 : if (!equali1(n)) return 0;
295 69986 : if (px) *px = gen_1;
296 69986 : return 1;
297 : }
298 350651 : D = divisors(Z_factor_listP(shifti(n, -1), L));
299 : /* loop through primes p > m, d = p-1 | n */
300 350651 : av2 = avma;
301 350651 : if (!m)
302 : { /* special case p = 2, d = 1 */
303 69986 : k = n;
304 69986 : for (v = 1;; v++) {
305 69986 : if (istotient_i(k, gen_2, L, px)) {
306 69986 : if (px) *px = shifti(*px, v);
307 69986 : return 1;
308 : }
309 0 : if (mod2(k)) break;
310 0 : k = shifti(k,-1);
311 : }
312 0 : set_avma(av2);
313 : }
314 1099462 : for (i = 1; i < lg(D); ++i)
315 : {
316 1001588 : GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
317 1001588 : if (m && cmpii(d, m) < 0) continue;
318 677782 : p = addiu(d, 1);
319 677782 : if (!isprime(p)) continue;
320 442064 : k = diviiexact(n, d);
321 481593 : for (v = 1;; v++) {
322 : GEN r;
323 481593 : if (istotient_i(k, p, L, px)) {
324 182791 : if (px) *px = mulii(*px, powiu(p, v));
325 182791 : return 1;
326 : }
327 298802 : k = dvmdii(k, p, &r);
328 298802 : if (r != gen_0) break;
329 : }
330 259273 : set_avma(av2);
331 : }
332 97874 : return gc_long(av,0);
333 : }
334 :
335 : /* find x such that phi(x) = n */
336 : long
337 70000 : istotient(GEN n, GEN *px)
338 : {
339 70000 : pari_sp av = avma;
340 70000 : if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
341 70000 : if (signe(n) < 1) return 0;
342 70000 : if (mod2(n))
343 : {
344 14 : if (!equali1(n)) return 0;
345 14 : if (px) *px = gen_1;
346 14 : return 1;
347 : }
348 69986 : if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
349 : {
350 69986 : if (!px) set_avma(av);
351 : else
352 69986 : *px = gerepileuptoint(av, *px);
353 69986 : return 1;
354 : }
355 0 : return gc_long(av,0);
356 : }
357 :
358 : /*********************************************************************/
359 : /** KRONECKER SYMBOL **/
360 : /*********************************************************************/
361 : /* t = 3,5 mod 8 ? (= 2 not a square mod t) */
362 : static int
363 312100188 : ome(long t)
364 : {
365 312100188 : switch(t & 7)
366 : {
367 176905919 : case 3:
368 176905919 : case 5: return 1;
369 135194269 : default: return 0;
370 : }
371 : }
372 : /* t a t_INT, is t = 3,5 mod 8 ? */
373 : static int
374 5492364 : gome(GEN t)
375 5492364 : { return signe(t)? ome( mod2BIL(t) ): 0; }
376 :
377 : /* assume y odd, return kronecker(x,y) * s */
378 : static long
379 220107216 : krouu_s(ulong x, ulong y, long s)
380 : {
381 220107216 : ulong x1 = x, y1 = y, z;
382 999401012 : while (x1)
383 : {
384 779329853 : long r = vals(x1);
385 779401933 : if (r)
386 : {
387 414825280 : if (odd(r) && ome(y1)) s = -s;
388 414717143 : x1 >>= r;
389 : }
390 779293796 : if (x1 & y1 & 2) s = -s;
391 779293796 : z = y1 % x1; y1 = x1; x1 = z;
392 : }
393 220071159 : return (y1 == 1)? s: 0;
394 : }
395 :
396 : long
397 11612034 : kronecker(GEN x, GEN y)
398 : {
399 11612034 : pari_sp av = avma;
400 11612034 : long s = 1, r;
401 : ulong xu;
402 :
403 11612034 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
404 11612034 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
405 11612034 : switch (signe(y))
406 : {
407 63 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
408 0 : case 0: return is_pm1(x);
409 : }
410 11612034 : r = vali(y);
411 11612034 : if (r)
412 : {
413 1347351 : if (!mpodd(x)) return gc_long(av,0);
414 322649 : if (odd(r) && gome(x)) s = -s;
415 322649 : y = shifti(y,-r);
416 : }
417 10587331 : x = modii(x,y);
418 12969113 : while (lgefint(x) > 3) /* x < y */
419 : {
420 : GEN z;
421 2381831 : r = vali(x);
422 2381844 : if (r)
423 : {
424 1301516 : if (odd(r) && gome(y)) s = -s;
425 1301465 : x = shifti(x,-r);
426 : }
427 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
428 2380393 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
429 2380648 : z = remii(y,x); y = x; x = z;
430 2381825 : if (gc_needed(av,2))
431 : {
432 0 : if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
433 0 : gerepileall(av, 2, &x, &y);
434 : }
435 : }
436 10587282 : xu = itou(x);
437 10587286 : if (!xu) return is_pm1(y)? s: 0;
438 10521333 : r = vals(xu);
439 10521339 : if (r)
440 : {
441 5600235 : if (odd(r) && gome(y)) s = -s;
442 5600235 : xu >>= r;
443 : }
444 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
445 10521339 : if (xu & mod2BIL(y) & 2) s = -s;
446 10521346 : return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
447 : }
448 :
449 : long
450 35973 : krois(GEN x, long y)
451 : {
452 : ulong yu;
453 35973 : long s = 1;
454 :
455 35973 : if (y <= 0)
456 : {
457 7 : if (y == 0) return is_pm1(x);
458 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
459 : }
460 : else
461 35966 : yu = (ulong)y;
462 35966 : if (!odd(yu))
463 : {
464 : long r;
465 14938 : if (!mpodd(x)) return 0;
466 11088 : r = vals(yu); yu >>= r;
467 11088 : if (odd(r) && gome(x)) s = -s;
468 : }
469 32116 : return krouu_s(umodiu(x, yu), yu, s);
470 : }
471 : /* assume y != 0 */
472 : long
473 27375932 : kroiu(GEN x, ulong y)
474 : {
475 : long r;
476 27375932 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
477 299211 : if (!mpodd(x)) return 0;
478 204606 : r = vals(y); y >>= r;
479 204614 : return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
480 : }
481 :
482 : /* assume y > 0, odd, return s * kronecker(x,y) */
483 : static long
484 178605 : krouodd(ulong x, GEN y, long s)
485 : {
486 : long r;
487 178605 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
488 28776 : if (!x) return 0; /* y != 1 */
489 28776 : r = vals(x);
490 28776 : if (r)
491 : {
492 15679 : if (odd(r) && gome(y)) s = -s;
493 15679 : x >>= r;
494 : }
495 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
496 28776 : if (x & mod2BIL(y) & 2) s = -s;
497 28776 : return krouu_s(umodiu(y,x), x, s);
498 : }
499 :
500 : long
501 162697 : krosi(long x, GEN y)
502 : {
503 162697 : const pari_sp av = avma;
504 162697 : long s = 1, r;
505 162697 : switch (signe(y))
506 : {
507 0 : case -1: y = negi(y); if (x < 0) s = -1; break;
508 0 : case 0: return (x==1 || x==-1);
509 : }
510 162697 : r = vali(y);
511 162697 : if (r)
512 : {
513 16884 : if (!odd(x)) return gc_long(av,0);
514 16884 : if (odd(r) && ome(x)) s = -s;
515 16884 : y = shifti(y,-r);
516 : }
517 162697 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
518 162697 : return gc_long(av, krouodd((ulong)x, y, s));
519 : }
520 :
521 : long
522 15908 : kroui(ulong x, GEN y)
523 : {
524 15908 : const pari_sp av = avma;
525 15908 : long s = 1, r;
526 15908 : switch (signe(y))
527 : {
528 0 : case -1: y = negi(y); break;
529 0 : case 0: return x==1UL;
530 : }
531 15908 : r = vali(y);
532 15908 : if (r)
533 : {
534 0 : if (!odd(x)) return gc_long(av,0);
535 0 : if (odd(r) && ome(x)) s = -s;
536 0 : y = shifti(y,-r);
537 : }
538 15908 : return gc_long(av, krouodd(x, y, s));
539 : }
540 :
541 : long
542 95816617 : kross(long x, long y)
543 : {
544 : ulong yu;
545 95816617 : long s = 1;
546 :
547 95816617 : if (y <= 0)
548 : {
549 68943 : if (y == 0) return (labs(x)==1);
550 68915 : yu = (ulong)-y; if (x < 0) s = -1;
551 : }
552 : else
553 95747674 : yu = (ulong)y;
554 95816589 : if (!odd(yu))
555 : {
556 : long r;
557 23352911 : if (!odd(x)) return 0;
558 16438914 : r = vals(yu); yu >>= r;
559 16438914 : if (odd(r) && ome(x)) s = -s;
560 : }
561 88902592 : x %= (long)yu; if (x < 0) x += yu;
562 88902592 : return krouu_s((ulong)x, yu, s);
563 : }
564 :
565 : long
566 93217874 : krouu(ulong x, ulong y)
567 : {
568 : long r;
569 93217874 : if (odd(y)) return krouu_s(x, y, 1);
570 5520 : if (!odd(x)) return 0;
571 6217 : r = vals(y); y >>= r;
572 6217 : return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
573 : }
574 :
575 : /*********************************************************************/
576 : /** HILBERT SYMBOL **/
577 : /*********************************************************************/
578 : /* x,y are t_INT or t_REAL */
579 : static long
580 7343 : mphilbertoo(GEN x, GEN y)
581 : {
582 7343 : long sx = signe(x), sy = signe(y);
583 7343 : if (!sx || !sy) return 0;
584 7343 : return (sx < 0 && sy < 0)? -1: 1;
585 : }
586 :
587 : long
588 135628 : hilbertii(GEN x, GEN y, GEN p)
589 : {
590 : pari_sp av;
591 : long oddvx, oddvy, z;
592 :
593 135628 : if (!p) return mphilbertoo(x,y);
594 128306 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
595 128306 : if (!signe(x) || !signe(y)) return 0;
596 128285 : av = avma;
597 128285 : oddvx = odd(Z_pvalrem(x,p,&x));
598 128285 : oddvy = odd(Z_pvalrem(y,p,&y));
599 : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
600 128285 : if (absequaliu(p, 2))
601 : {
602 12117 : z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
603 12117 : if (oddvx && gome(y)) z = -z;
604 12117 : if (oddvy && gome(x)) z = -z;
605 : }
606 : else
607 : {
608 116168 : z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
609 116168 : if (oddvx && kronecker(y,p) < 0) z = -z;
610 116168 : if (oddvy && kronecker(x,p) < 0) z = -z;
611 : }
612 128285 : return gc_long(av, z);
613 : }
614 :
615 : static void
616 196 : err_prec(void) { pari_err_PREC("hilbert"); }
617 : static void
618 161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
619 : static void
620 56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
621 :
622 : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
623 : * Return lift(x) provided it's p-adic accuracy is large enough to decide
624 : * hilbert()'s value [ problem at p = 2 ] */
625 : static GEN
626 420 : lift_intmod(GEN x, GEN *pp)
627 : {
628 420 : GEN p = *pp, N = gel(x,1);
629 420 : x = gel(x,2);
630 420 : if (!p)
631 : {
632 266 : *pp = p = N;
633 266 : switch(itos_or_0(p))
634 : {
635 126 : case 2:
636 126 : case 4: err_prec();
637 : }
638 140 : return x;
639 : }
640 154 : if (!signe(p)) err_oo(N);
641 112 : if (absequaliu(p,2))
642 42 : { if (vali(N) <= 2) err_prec(); }
643 : else
644 70 : { if (!dvdii(N,p)) err_p(N,p); }
645 28 : if (!signe(x)) err_prec();
646 21 : return x;
647 : }
648 : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
649 : * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
650 : * to decide hilbert()'s value [ problem at p = 2 ]*/
651 : static GEN
652 210 : lift_padic(GEN x, GEN *pp)
653 : {
654 210 : GEN p = *pp, q = gel(x,2), y = gel(x,4);
655 210 : if (!p) *pp = p = q;
656 147 : else if (!equalii(p,q)) err_p(p, q);
657 105 : if (absequaliu(p,2) && precp(x) <= 2) err_prec();
658 70 : if (!signe(y)) err_prec();
659 70 : return odd(valp(x))? mulii(p,y): y;
660 : }
661 :
662 : long
663 60665 : hilbert(GEN x, GEN y, GEN p)
664 : {
665 60665 : pari_sp av = avma;
666 60665 : long tx = typ(x), ty = typ(y);
667 :
668 60665 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
669 60665 : if (tx == t_REAL)
670 : {
671 77 : if (p && signe(p)) err_oo(p);
672 63 : switch (ty)
673 : {
674 7 : case t_INT:
675 7 : case t_REAL: return mphilbertoo(x,y);
676 0 : case t_FRAC: return mphilbertoo(x,gel(y,1));
677 56 : default: pari_err_TYPE2("hilbert",x,y);
678 : }
679 : }
680 60588 : if (ty == t_REAL)
681 : {
682 14 : if (p && signe(p)) err_oo(p);
683 14 : switch (tx)
684 : {
685 14 : case t_INT:
686 14 : case t_REAL: return mphilbertoo(x,y);
687 0 : case t_FRAC: return mphilbertoo(gel(x,1),y);
688 0 : default: pari_err_TYPE2("hilbert",x,y);
689 : }
690 : }
691 60574 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
692 60371 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
693 :
694 60315 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
695 60252 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
696 :
697 60175 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
698 60175 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
699 :
700 60175 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
701 60175 : if (p && !signe(p)) p = NULL;
702 60175 : return gc_long(av, hilbertii(x,y,p));
703 : }
704 :
705 : /*******************************************************************/
706 : /* SQUARE ROOT MODULO p */
707 : /*******************************************************************/
708 : static void
709 2078507 : checkp(ulong q, ulong p)
710 2078507 : { if (!q) pari_err_PRIME("Fl_nonsquare",utoipos(p)); }
711 : /* p = 1 (mod 4) prime, return the first quadratic nonresidue, a prime */
712 : static ulong
713 11078356 : nonsquare1_Fl(ulong p)
714 : {
715 : forprime_t S;
716 : ulong q;
717 11078356 : if ((p & 7UL) != 1) return 2UL;
718 4104123 : q = p % 3; if (q == 2) return 3UL;
719 1303148 : checkp(q, p);
720 1309710 : q = p % 5; if (q == 2 || q == 3) return 5UL;
721 482051 : checkp(q, p);
722 482051 : q = p % 7; if (q != 4 && q >= 3) return 7UL;
723 176525 : checkp(q, p);
724 : /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
725 176547 : u_forprime_init(&S, 11, 1967);
726 286988 : while ((q = u_forprime_next(&S)))
727 : {
728 286986 : if (krouu(q, p) < 0) return q;
729 110435 : checkp(q, p);
730 : }
731 0 : checkp(0, p);
732 : return 0; /*LCOV_EXCL_LINE*/
733 : }
734 : /* p > 2 a prime */
735 : ulong
736 7893 : nonsquare_Fl(ulong p)
737 7893 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
738 :
739 : /* allow pi = 0 */
740 : ulong
741 173626 : Fl_2gener_pre(ulong p, ulong pi)
742 : {
743 173626 : ulong p1 = p-1;
744 173626 : long e = vals(p1);
745 173614 : if (e == 1) return p1;
746 64698 : return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
747 : }
748 :
749 : ulong
750 67793 : Fl_2gener_pre_i(ulong ns, ulong p, ulong pi)
751 : {
752 67793 : ulong p1 = p-1;
753 67793 : long e = vals(p1);
754 67793 : if (e == 1) return p1;
755 26020 : return Fl_powu_pre(ns, p1 >> e, p, pi);
756 : }
757 :
758 : static ulong
759 11495909 : Fl_sqrt_i(ulong a, ulong y, ulong p)
760 : {
761 : long i, e, k;
762 : ulong p1, q, v, w;
763 :
764 11495909 : if (!a) return 0;
765 10243682 : p1 = p - 1; e = vals(p1);
766 10243450 : if (e == 0) /* p = 2 */
767 : {
768 421559 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
769 422606 : return ((a & 1) == 0)? 0: 1;
770 : }
771 9821891 : if (e == 1)
772 : {
773 4728353 : v = Fl_powu(a, (p+1) >> 2, p);
774 4728542 : if (Fl_sqr(v, p) != a) return ~0UL;
775 4723699 : p1 = p - v; if (v > p1) v = p1;
776 4723699 : return v;
777 : }
778 5093538 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
779 5093538 : p1 = Fl_powu(a, q >> 1, p); /* a ^ [(q-1)/2] */
780 5093632 : if (!p1) return 0;
781 5093632 : v = Fl_mul(a, p1, p);
782 5093646 : w = Fl_mul(v, p1, p);
783 5093738 : if (!y) y = Fl_powu(nonsquare1_Fl(p), q, p);
784 8663633 : while (w != 1)
785 : { /* a*w = v^2, y primitive 2^e-th root of 1
786 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
787 3571826 : p1 = Fl_sqr(w, p);
788 5881970 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr(p1, p);
789 3571833 : if (k == e) return ~0UL;
790 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
791 3569756 : p1 = y;
792 4688866 : for (i=1; i < e-k; i++) p1 = Fl_sqr(p1, p);
793 3569756 : y = Fl_sqr(p1, p); e = k;
794 3569781 : w = Fl_mul(y, w, p);
795 3569815 : v = Fl_mul(v, p1, p);
796 : }
797 5091807 : p1 = p - v; if (v > p1) v = p1;
798 5091807 : return v;
799 : }
800 :
801 : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. Allow pi = 0 */
802 : ulong
803 32290597 : Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
804 : {
805 : long i, e, k;
806 : ulong p1, q, v, w;
807 :
808 32290597 : if (!pi) return Fl_sqrt_i(a, y, p);
809 20794801 : if (!a) return 0;
810 20670902 : p1 = p - 1; e = vals(p1);
811 20678218 : if (e == 0) /* p = 2 */
812 : {
813 0 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
814 0 : return ((a & 1) == 0)? 0: 1;
815 : }
816 20689568 : if (e == 1)
817 : {
818 14716161 : v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
819 14672229 : if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
820 14686448 : p1 = p - v; if (v > p1) v = p1;
821 14686448 : return v;
822 : }
823 5973407 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
824 5973407 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
825 5968696 : if (!p1) return 0;
826 5968696 : v = Fl_mul_pre(a, p1, p, pi);
827 5969926 : w = Fl_mul_pre(v, p1, p, pi);
828 5969577 : if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
829 11362986 : while (w != 1)
830 : { /* a*w = v^2, y primitive 2^e-th root of 1
831 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
832 5393988 : p1 = Fl_sqr_pre(w,p,pi);
833 10137092 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
834 5395855 : if (k == e) return ~0UL;
835 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
836 5395763 : p1 = y;
837 7068512 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
838 5395769 : y = Fl_sqr_pre(p1, p, pi); e = k;
839 5396146 : w = Fl_mul_pre(y, w, p, pi);
840 5395690 : v = Fl_mul_pre(v, p1, p, pi);
841 : }
842 5968998 : p1 = p - v; if (v > p1) v = p1;
843 5968998 : return v;
844 : }
845 :
846 : ulong
847 11511734 : Fl_sqrt(ulong a, ulong p)
848 11511734 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0; return Fl_sqrt_pre_i(a, 0, p, pi); }
849 :
850 : ulong
851 20622190 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
852 20622190 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
853 :
854 : /* allow pi = 0 */
855 : static ulong
856 51035 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
857 : {
858 51035 : ulong x, y, m, le1 = upowuu(l, e-1);
859 51034 : for (x = 2; ; x++)
860 : {
861 79000 : y = Fl_powu_pre(x, r, p, pi);
862 78998 : if (y==1) continue;
863 62341 : m = Fl_powu_pre(y, le1, p, pi);
864 62346 : if (m != 1) break;
865 : }
866 51037 : *pt_m = m; return y;
867 : }
868 :
869 : /* solve x^l = a , l prime in G of order q.
870 : *
871 : * q = (l^e)*r, e >= 1, (r,l) = 1
872 : * y generates the l-Sylow of G
873 : * m = y^(l^(e-1)) != 1 */
874 : static ulong
875 124189 : Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
876 : {
877 : ulong u2, p1, v, w, z, dl;
878 124189 : if (a==0) return a;
879 124183 : u2 = Fl_inv(l%r, r);
880 124184 : v = Fl_powu_pre(a, u2, p, pi);
881 124182 : w = Fl_powu_pre(v, l, p, pi);
882 124177 : w = pi? Fl_mul_pre(w, Fl_inv(a, p), p, pi): Fl_div(w, a, p);
883 124166 : if (w==1) return v;
884 49684 : if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
885 72279 : while (w!=1)
886 : {
887 54923 : ulong k = 0;
888 54923 : p1 = w;
889 : do
890 : {
891 81422 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
892 81422 : if (++k == e) return ULONG_MAX;
893 49092 : } while (p1!=1);
894 22593 : dl = Fl_log_pre(z, m, l, p, pi);
895 22593 : dl = Fl_neg(dl, l);
896 22593 : p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
897 22593 : m = Fl_powu_pre(m, dl, p, pi);
898 22593 : e = k;
899 22593 : v = pi? Fl_mul_pre(p1,v,p,pi): Fl_mul(p1,v,p);
900 22593 : y = Fl_powu_pre(p1,l,p,pi);
901 22593 : w = pi? Fl_mul_pre(y,w,p,pi): Fl_mul(y,w,p);
902 : }
903 17356 : return v;
904 : }
905 :
906 : /* allow pi = 0 */
907 : static ulong
908 121855 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
909 : {
910 121855 : ulong r, e = u_lvalrem(p-1, l, &r);
911 121856 : return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
912 : }
913 : /* allow pi = 0 */
914 : ulong
915 121854 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
916 121854 : { return Fl_sqrtl_i(a, l, p, pi, 0, 0); }
917 :
918 : ulong
919 0 : Fl_sqrtl(ulong a, ulong l, ulong p)
920 0 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
921 0 : return Fl_sqrtl_i(a, l, p, pi, 0, 0); }
922 :
923 : /* allow pi = 0 */
924 : ulong
925 85685 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
926 : {
927 85685 : ulong m, q = p-1, z;
928 85685 : ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
929 85685 : if (a==0)
930 : {
931 48153 : if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
932 48146 : if (zetan) *zetan = 1UL;
933 48146 : return 0;
934 : }
935 37532 : if (n==1)
936 : {
937 420 : if (zetan) *zetan = 1;
938 420 : return n < 0? Fl_inv(a,p): a;
939 : }
940 37112 : if (n==2)
941 : {
942 5710 : if (zetan) *zetan = p-1;
943 5710 : return Fl_sqrt_pre_i(a, 0, p, pi);
944 : }
945 31402 : if (a == 1 && !zetan) return a;
946 7973 : m = ugcd(nn, q);
947 7973 : z = 1;
948 7973 : if (m!=1)
949 : {
950 1400 : GEN F = factoru(m);
951 : long i, j, e;
952 : ulong r, zeta, y, l;
953 3108 : for (i = nbrows(F); i; i--)
954 : {
955 1771 : l = ucoeff(F,i,1);
956 1771 : j = ucoeff(F,i,2);
957 1771 : e = u_lvalrem(q,l, &r);
958 1771 : y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
959 1771 : if (zetan)
960 : {
961 112 : ulong Y = Fl_powu_pre(y, upowuu(l,e-j), p, pi);
962 112 : z = pi? Fl_mul_pre(z, Y, p, pi): Fl_mul(z, Y, p);
963 : }
964 1771 : if (a!=1)
965 : do
966 : {
967 2331 : a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
968 2317 : if (a==ULONG_MAX) return ULONG_MAX;
969 2268 : } while (--j);
970 : }
971 : }
972 7910 : if (m != nn)
973 : {
974 6594 : ulong qm = q/m, nm = (nn/m) % qm;
975 6594 : a = Fl_powu_pre(a, Fl_inv(nm, qm), p, pi);
976 : }
977 7910 : if (n < 0) a = Fl_inv(a, p);
978 7910 : if (zetan) *zetan = z;
979 7910 : return a;
980 : }
981 :
982 : ulong
983 85685 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
984 : {
985 85685 : ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
986 85685 : return Fl_sqrtn_pre(a, n, p, pi, zetan);
987 : }
988 :
989 : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
990 : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
991 : * and in average is about 2 or 2.5 times worse. But try both algorithms for
992 : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
993 : *
994 : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
995 : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
996 : * If (a|p)=1, then sqrt(a) is in F_p.
997 : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
998 :
999 : /* compute y^2, y = y[1] + y[2] X */
1000 : static GEN
1001 449 : sqrt_Cipolla_sqr(void *data, GEN y)
1002 : {
1003 449 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
1004 449 : GEN u2 = sqri(u), v2 = sqri(v);
1005 449 : v = subii(sqri(addii(v,u)), addii(u2,v2));
1006 449 : u = addii(u2, mulii(v2,n));
1007 449 : retmkvec2(modii(u,p), modii(v,p));
1008 : }
1009 : /* compute (t+X) y^2 */
1010 : static GEN
1011 23 : sqrt_Cipolla_msqr(void *data, GEN y)
1012 : {
1013 23 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2);
1014 23 : ulong t = gel(data,4)[2];
1015 23 : GEN d = addii(u, mului(t,v)), d2 = sqri(d);
1016 23 : GEN b = remii(mulii(a,v), p);
1017 23 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
1018 23 : v = subii(d2, mulii(b,v));
1019 23 : retmkvec2(modii(u,p), modii(v,p));
1020 : }
1021 : /* assume a reduced mod p [ otherwise correct but inefficient ] */
1022 : static GEN
1023 8 : sqrt_Cipolla(GEN a, GEN p)
1024 : {
1025 : pari_sp av;
1026 : GEN u, v, n, y, pov2;
1027 : ulong t;
1028 :
1029 8 : if (kronecker(a, p) < 0) return NULL;
1030 8 : pov2 = shifti(p,-1); /* center to avoid multiplying by huge base*/
1031 8 : if (cmpii(a,pov2) > 0) a = subii(a,p);
1032 8 : av = avma;
1033 41 : for (t=1; ; t++, set_avma(av))
1034 : {
1035 41 : n = subsi((long)(t*t), a);
1036 41 : if (kronecker(n, p) < 0) break;
1037 : }
1038 :
1039 : /* compute (t+X)^((p-1)/2) =: u+vX */
1040 8 : u = utoipos(t);
1041 8 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
1042 : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
1043 : /* Now u+vX = (t+X)^((p-1)/2); thus
1044 : * (u+vX)(t+X) = sqrt(a) + 0 X
1045 : * Whence,
1046 : * sqrt(a) = (u+vt)t - v*a
1047 : * 0 = (u+vt)
1048 : * Thus a square root is v*a */
1049 :
1050 8 : v = Fp_mul(gel(y, 2), a, p);
1051 8 : if (cmpii(v,pov2) > 0) v = subii(p,v);
1052 8 : return v;
1053 : }
1054 :
1055 : /* Return NULL if p is found to be composite */
1056 : static GEN
1057 6097 : Fp_2gener_all(long e, GEN p)
1058 : {
1059 : GEN y, m;
1060 : long k;
1061 6097 : GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
1062 6097 : if (e==0 && !equaliu(p,2)) return NULL;
1063 6097 : for (k=2; ; k++)
1064 12336 : {
1065 18433 : long i = krosi(k, p);
1066 18433 : if (i >= 0)
1067 : {
1068 12336 : if (i) continue;
1069 0 : return NULL;
1070 : }
1071 6097 : y = m = Fp_pow(utoi(k), q, p);
1072 17986 : for (i=1; i<e; i++)
1073 11889 : if (equali1(m = Fp_sqr(m, p))) break;
1074 6097 : if (i == e) break; /* success */
1075 : }
1076 6097 : return y;
1077 : }
1078 :
1079 : /* Return NULL if p is found to be composite */
1080 : GEN
1081 3185 : Fp_2gener(GEN p)
1082 3185 : { return Fp_2gener_all(vali(subis(p,1)),p); }
1083 :
1084 : GEN
1085 19931 : Fp_2gener_i(GEN ns, GEN p)
1086 : {
1087 19931 : GEN p1 = subiu(p,1);
1088 19931 : long e = vali(p1);
1089 19931 : if (e == 1) return p1;
1090 18546 : return Fp_pow(ns, shifti(p1,-e), p);
1091 : }
1092 :
1093 : /* smallest square root */
1094 : static GEN
1095 78327 : choose_sqrt(GEN v, GEN p)
1096 : {
1097 78327 : pari_sp av = avma;
1098 78327 : GEN q = subii(p,v);
1099 78323 : if (cmpii(v,q) > 0) v = q; else set_avma(av);
1100 78325 : return v;
1101 : }
1102 : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
1103 : GEN
1104 3382295 : Fp_sqrt_i(GEN a, GEN y, GEN p)
1105 : {
1106 3382295 : pari_sp av = avma;
1107 : long i, k, e;
1108 : GEN p1, q, v, w;
1109 :
1110 3382295 : if (lgefint(p) == 3)
1111 : {
1112 3303454 : ulong pp = uel(p,2), u = umodiu(a, pp);
1113 3303487 : if (!u) return gen_0;
1114 3006980 : u = Fl_sqrt(u, pp);
1115 3007022 : return (u == ~0UL)? NULL: utoipos(u);
1116 : }
1117 :
1118 78841 : a = modii(a, p); if (!signe(a)) return gc_const(av, gen_0);
1119 78720 : p1 = subiu(p,1); e = vali(p1);
1120 78738 : if (e <= 2)
1121 : { /* direct formulas more efficient */
1122 : pari_sp av2;
1123 38392 : if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
1124 38394 : if (e == 1)
1125 : {
1126 19638 : q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
1127 19627 : v = Fp_pow(a, q, p);
1128 : }
1129 : else
1130 : { /* Atkin's formula */
1131 18756 : GEN i, a2 = shifti(a,1);
1132 18750 : if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
1133 18755 : q = shifti(p1, -3); /* (p-5)/8 */
1134 18757 : v = Fp_pow(a2, q, p);
1135 18760 : i = Fp_mul(a2, Fp_sqr(v,p), p); /* i^2 = -1 */
1136 18760 : v = Fp_mul(a, Fp_mul(v, subiu(i,1), p), p);
1137 : }
1138 38418 : av2 = avma;
1139 : /* must check equality in case (a/p) = -1 or p not prime */
1140 38418 : e = equalii(Fp_sqr(v,p), a); set_avma(av2);
1141 38418 : return e? gerepileuptoint(av,choose_sqrt(v,p)): NULL;
1142 : }
1143 : /* On average, Cipolla is better than Tonelli/Shanks if and only if
1144 : * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
1145 40346 : if (e*(e-1) > 20 + 8 * expi(p))
1146 : {
1147 8 : v = sqrt_Cipolla(a,p); if (!v) return gc_NULL(av);
1148 8 : return gerepileuptoint(av,v);
1149 : }
1150 40339 : if (!y)
1151 : {
1152 2912 : y = Fp_2gener_all(e, p);
1153 2912 : if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
1154 : }
1155 40339 : q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
1156 40338 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
1157 40342 : v = Fp_mul(a, p1, p);
1158 40343 : w = Fp_mul(v, p1, p);
1159 96095 : while (!equali1(w))
1160 : { /* a*w = v^2, y primitive 2^e-th root of 1
1161 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1162 55796 : p1 = Fp_sqr(w,p);
1163 115050 : for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
1164 55791 : if (k == e) return gc_NULL(av); /* p composite or (a/p) != 1 */
1165 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1166 55748 : p1 = y;
1167 79619 : for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
1168 55748 : y = Fp_sqr(p1, p); e = k;
1169 55754 : w = Fp_mul(y, w, p);
1170 55752 : v = Fp_mul(v, p1, p);
1171 55753 : if (gc_needed(av,1))
1172 : {
1173 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
1174 0 : gerepileall(av,3, &y,&w,&v);
1175 : }
1176 : }
1177 40298 : return gerepileuptoint(av, choose_sqrt(v,p));
1178 : }
1179 :
1180 : GEN
1181 3323902 : Fp_sqrt(GEN a, GEN p)
1182 3323902 : { return Fp_sqrt_i(a, NULL, p); }
1183 :
1184 : /*********************************************************************/
1185 : /** GCD & BEZOUT **/
1186 : /*********************************************************************/
1187 :
1188 : GEN
1189 43641036 : lcmii(GEN x, GEN y)
1190 : {
1191 : pari_sp av;
1192 : GEN a, b;
1193 43641036 : if (!signe(x) || !signe(y)) return gen_0;
1194 43641116 : av = avma; a = gcdii(x,y);
1195 43639449 : if (absequalii(a,y)) { set_avma(av); return absi(x); }
1196 8639931 : if (!equali1(a)) y = diviiexact(y,a);
1197 8639946 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
1198 : }
1199 :
1200 : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
1201 : * set *pd = gcd(x,N) */
1202 : GEN
1203 6610419 : Fp_invgen(GEN x, GEN N, GEN *pd)
1204 : {
1205 : GEN d, d0, e, v;
1206 6610419 : if (lgefint(N) == 3)
1207 : {
1208 5829200 : ulong dd, NN = N[2], xx = umodiu(x,NN);
1209 5829760 : if (!xx) { *pd = N; return gen_0; }
1210 5829760 : xx = Fl_invgen(xx, NN, &dd);
1211 5830378 : *pd = utoi(dd); return utoi(xx);
1212 : }
1213 781219 : *pd = d = bezout(x, N, &v, NULL);
1214 781232 : if (equali1(d)) return v;
1215 : /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
1216 691263 : e = diviiexact(N,d);
1217 691263 : d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
1218 691263 : if (equali1(d0)) return v;
1219 527976 : if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
1220 527976 : return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
1221 : }
1222 :
1223 : /*********************************************************************/
1224 : /** CHINESE REMAINDERS **/
1225 : /*********************************************************************/
1226 :
1227 : /* Chinese Remainder Theorem. x and y must have the same type (integermod,
1228 : * polymod, or polynomial/vector/matrix recursively constructed with these
1229 : * as coefficients). Creates (with the same type) a z in the same residue
1230 : * class as x and the same residue class as y, if it is possible.
1231 : *
1232 : * We also allow (during recursion) two identical objects even if they are
1233 : * not integermod or polymod. For example:
1234 : *
1235 : * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
1236 : * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
1237 : * ? chinese(x, y)
1238 : * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
1239 :
1240 : static GEN
1241 2401623 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
1242 : {
1243 2401623 : GEN z = gassoc_proto(f,x,NULL);
1244 2401616 : if (z == gen_1) retmkintmod(gen_0,gen_1);
1245 2401581 : return z;
1246 : }
1247 :
1248 : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
1249 : * call chinese: makes Mod(0,1) a better "neutral" element */
1250 : static GEN
1251 21 : chinese_intpol(GEN x,GEN y)
1252 : {
1253 21 : pari_sp av = avma;
1254 21 : GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
1255 21 : return gerepileupto(av, chinese(z, y));
1256 : }
1257 :
1258 : GEN
1259 2345 : chinese1(GEN x) { return gen_chinese(x,chinese); }
1260 :
1261 : GEN
1262 5432 : chinese(GEN x, GEN y)
1263 : {
1264 : pari_sp av;
1265 5432 : long tx = typ(x), ty;
1266 : GEN z,p1,p2,d,u,v;
1267 :
1268 5432 : if (!y) return chinese1(x);
1269 5383 : if (gidentical(x,y)) return gcopy(x);
1270 5376 : ty = typ(y);
1271 5376 : if (tx == ty) switch(tx)
1272 : {
1273 3829 : case t_POLMOD:
1274 : {
1275 3829 : GEN A = gel(x,1), B = gel(y,1);
1276 3829 : GEN a = gel(x,2), b = gel(y,2);
1277 3829 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
1278 3829 : if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
1279 3829 : av = avma;
1280 3829 : d = RgX_extgcd(A,B,&u,&v);
1281 3829 : p2 = gsub(b, a);
1282 3829 : if (!gequal0(gmod(p2, d))) break;
1283 3829 : p1 = gdiv(A,d);
1284 3829 : p2 = gadd(a, gmul(gmul(u,p1), p2));
1285 :
1286 3829 : z = cgetg(3, t_POLMOD);
1287 3829 : gel(z,1) = gmul(p1,B);
1288 3829 : gel(z,2) = gmod(p2,gel(z,1));
1289 3829 : return gerepileupto(av, z);
1290 : }
1291 1505 : case t_INTMOD:
1292 : {
1293 1505 : GEN A = gel(x,1), B = gel(y,1);
1294 1505 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
1295 1505 : z = cgetg(3,t_INTMOD);
1296 1505 : Z_chinese_pre(A, B, &C, &U, &d);
1297 1505 : c = Z_chinese_post(a, b, C, U, d);
1298 1505 : if (!c) pari_err_OP("chinese", x,y);
1299 1505 : set_avma((pari_sp)z);
1300 1505 : gel(z,1) = icopy(C);
1301 1505 : gel(z,2) = icopy(c); return z;
1302 : }
1303 14 : case t_POL:
1304 : {
1305 14 : long i, lx = lg(x), ly = lg(y);
1306 14 : if (varn(x) != varn(y)) break;
1307 14 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
1308 14 : z = cgetg(lx, t_POL); z[1] = x[1];
1309 42 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
1310 14 : if (i < lx)
1311 : {
1312 14 : GEN _0 = Rg_get_0(y);
1313 28 : for ( ; i<lx; i++) gel(z,i) = chinese(gel(x,i),_0);
1314 : }
1315 14 : return z;
1316 : }
1317 7 : case t_VEC: case t_COL: case t_MAT:
1318 : {
1319 : long i, lx;
1320 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
1321 21 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
1322 7 : return z;
1323 : }
1324 : }
1325 21 : if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
1326 7 : if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
1327 0 : pari_err_OP("chinese",x,y);
1328 : return NULL; /* LCOV_EXCL_LINE */
1329 : }
1330 :
1331 : /* init chinese(Mod(.,A), Mod(.,B)) */
1332 : void
1333 267283 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
1334 : {
1335 267283 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
1336 267285 : GEN t = diviiexact(A,d);
1337 267262 : *pU = mulii(u, t);
1338 267265 : *pC = mulii(t, B);
1339 267257 : if (pd) *pd = d;
1340 267257 : }
1341 : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
1342 : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
1343 : * If d not NULL, check whether a = b mod d. */
1344 : GEN
1345 2981548 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
1346 : {
1347 : GEN b_a;
1348 2981548 : if (!signe(a))
1349 : {
1350 788655 : if (d && !dvdii(b, d)) return NULL;
1351 788655 : return Fp_mul(b, U, C);
1352 : }
1353 2192893 : b_a = subii(b,a);
1354 2192893 : if (d && !dvdii(b_a, d)) return NULL;
1355 2192893 : return modii(addii(a, mulii(U, b_a)), C);
1356 : }
1357 : static ulong
1358 2178065 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
1359 : {
1360 2178065 : if (!a) return Fl_mul(b, U, C);
1361 2178065 : return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
1362 : }
1363 :
1364 : GEN
1365 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
1366 : {
1367 2142 : pari_sp av = avma;
1368 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
1369 2142 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
1370 : }
1371 : GEN
1372 263579 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
1373 : {
1374 263579 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
1375 263554 : return Z_chinese_post(a,b, *pC, U, NULL);
1376 : }
1377 :
1378 : /* return lift(chinese(a mod A, b mod B))
1379 : * assume(A,B)=1, a,b,A,B integers and C = A*B */
1380 : GEN
1381 529236 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
1382 : {
1383 529236 : pari_sp av = avma;
1384 529236 : GEN U = mulii(Fp_inv(A,B), A);
1385 529236 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
1386 : }
1387 : ulong
1388 2178045 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
1389 2178045 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
1390 :
1391 : /* chinese1 for coprime moduli in Z */
1392 : static GEN
1393 2184786 : chinese1_coprime_Z_aux(GEN x, GEN y)
1394 : {
1395 2184786 : GEN z = cgetg(3, t_INTMOD);
1396 2184786 : GEN A = gel(x,1), a = gel(x, 2);
1397 2184786 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
1398 2184786 : pari_sp av = avma;
1399 2184786 : GEN U = mulii(Fp_inv(A,B), A);
1400 2184786 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
1401 2184786 : gel(z,1) = C; return z;
1402 : }
1403 : GEN
1404 2399278 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
1405 :
1406 : /*********************************************************************/
1407 : /** MODULAR EXPONENTIATION **/
1408 : /*********************************************************************/
1409 : /* xa ZV or nv */
1410 : GEN
1411 2170985 : ZV_producttree(GEN xa)
1412 : {
1413 2170985 : long n = lg(xa)-1;
1414 2170985 : long m = n==1 ? 1: expu(n-1)+1;
1415 2170986 : GEN T = cgetg(m+1, t_VEC), t;
1416 : long i, j, k;
1417 2170991 : t = cgetg(((n+1)>>1)+1, t_VEC);
1418 2170983 : if (typ(xa)==t_VECSMALL)
1419 : {
1420 2849274 : for (j=1, k=1; k<n; j++, k+=2)
1421 1805005 : gel(t, j) = muluu(xa[k], xa[k+1]);
1422 1044269 : if (k==n) gel(t, j) = utoi(xa[k]);
1423 : } else {
1424 2323895 : for (j=1, k=1; k<n; j++, k+=2)
1425 1197198 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
1426 1126697 : if (k==n) gel(t, j) = icopy(gel(xa,k));
1427 : }
1428 2170965 : gel(T,1) = t;
1429 3451700 : for (i=2; i<=m; i++)
1430 : {
1431 1280725 : GEN u = gel(T, i-1);
1432 1280725 : long n = lg(u)-1;
1433 1280725 : t = cgetg(((n+1)>>1)+1, t_VEC);
1434 2857149 : for (j=1, k=1; k<n; j++, k+=2)
1435 1576414 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
1436 1280735 : if (k==n) gel(t, j) = gel(u, k);
1437 1280735 : gel(T, i) = t;
1438 : }
1439 2170975 : return T;
1440 : }
1441 :
1442 : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
1443 : GEN
1444 53532992 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
1445 : {
1446 : long i,j,k;
1447 53532992 : long m = lg(T)-1, n = lg(P)-1;
1448 : GEN t;
1449 53532992 : GEN Tp = cgetg(m+1, t_VEC);
1450 53485673 : gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
1451 112802011 : for (i=m-1; i>=1; i--)
1452 : {
1453 59398998 : GEN u = gel(T, i);
1454 59398998 : GEN v = gel(Tp, i+1);
1455 59398998 : long n = lg(u)-1;
1456 59398998 : t = cgetg(n+1, t_VEC);
1457 147930176 : for (j=1, k=1; k<n; j++, k+=2)
1458 : {
1459 88655007 : gel(t, k) = modii(gel(v, j), gel(u, k));
1460 88672986 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
1461 : }
1462 59275169 : if (k==n) gel(t, k) = gel(v, j);
1463 59275169 : gel(Tp, i) = t;
1464 : }
1465 : {
1466 53403013 : GEN u = gel(T, i+1);
1467 53403013 : GEN v = gel(Tp, i+1);
1468 53403013 : long l = lg(u)-1;
1469 53403013 : if (typ(P)==t_VECSMALL)
1470 : {
1471 51237348 : GEN R = cgetg(n+1, t_VECSMALL);
1472 190788161 : for (j=1, k=1; j<=l; j++, k+=2)
1473 : {
1474 139178585 : uel(R,k) = umodiu(gel(v, j), P[k]);
1475 139523503 : if (k < n)
1476 111973078 : uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
1477 : }
1478 51609576 : return R;
1479 : }
1480 : else
1481 : {
1482 2165665 : GEN R = cgetg(n+1, t_VEC);
1483 5914667 : for (j=1, k=1; j<=l; j++, k+=2)
1484 : {
1485 3744081 : gel(R,k) = modii(gel(v, j), gel(P,k));
1486 3744077 : if (k < n)
1487 2999007 : gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
1488 : }
1489 2170586 : return R;
1490 : }
1491 : }
1492 : }
1493 :
1494 : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
1495 : GEN
1496 35292624 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
1497 : {
1498 35292624 : long m = lg(T)-1, n = lg(A)-1;
1499 : long i,j,k;
1500 35292624 : GEN Tp = cgetg(m+1, t_VEC);
1501 35270293 : GEN M = gel(T, 1);
1502 35270293 : GEN t = cgetg(lg(M), t_VEC);
1503 35241126 : if (typ(P)==t_VECSMALL)
1504 : {
1505 82120195 : for (j=1, k=1; k<n; j++, k+=2)
1506 : {
1507 61163012 : pari_sp av = avma;
1508 61163012 : GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
1509 60982336 : GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
1510 61071895 : gel(t, j) = gerepileuptoint(av, tj);
1511 : }
1512 20957183 : if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
1513 : } else
1514 : {
1515 30386276 : for (j=1, k=1; k<n; j++, k+=2)
1516 : {
1517 16079264 : pari_sp av = avma;
1518 16079264 : GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
1519 16077761 : GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
1520 16103639 : gel(t, j) = gerepileuptoint(av, tj);
1521 : }
1522 14307012 : if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
1523 : }
1524 35260194 : gel(Tp, 1) = t;
1525 67481218 : for (i=2; i<=m; i++)
1526 : {
1527 32207580 : GEN u = gel(T, i-1), M = gel(T, i);
1528 32207580 : GEN t = cgetg(lg(M), t_VEC);
1529 32256413 : GEN v = gel(Tp, i-1);
1530 32256413 : long n = lg(v)-1;
1531 88168843 : for (j=1, k=1; k<n; j++, k+=2)
1532 : {
1533 55947819 : pari_sp av = avma;
1534 55890080 : gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
1535 55947819 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
1536 : }
1537 32221024 : if (k==n) gel(t, j) = gel(v, k);
1538 32221024 : gel(Tp, i) = t;
1539 : }
1540 35273638 : return gmael(Tp,m,1);
1541 : }
1542 :
1543 : static GEN
1544 937784 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1545 : {
1546 937784 : long i, l = lg(gel(vA,1)), n = lg(P);
1547 937784 : GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
1548 30564388 : for (i=1; i < l; i++)
1549 : {
1550 29626712 : pari_sp av = avma;
1551 29626712 : GEN c, A = cgetg(n, typ(P));
1552 : long j;
1553 183053771 : for (j=1; j < n; j++) A[j] = mael(vA,j,i);
1554 29596419 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
1555 29616502 : gel(V,i) = gerepileuptoint(av, c);
1556 : }
1557 937676 : return V;
1558 : }
1559 :
1560 : static GEN
1561 553476 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1562 : {
1563 553476 : long i, j, l, n = lg(P);
1564 553476 : GEN mod = gmael(T, lg(T)-1, 1), V, w;
1565 553476 : w = cgetg(n, t_VECSMALL);
1566 1847779 : for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
1567 553473 : l = vecsmall_max(w);
1568 553474 : V = cgetg(l, t_POL);
1569 553454 : V[1] = mael(vA,1,1);
1570 3225304 : for (i=2; i < l; i++)
1571 : {
1572 2671837 : pari_sp av = avma;
1573 2671837 : GEN c, A = cgetg(n, typ(P));
1574 2671569 : if (typ(P)==t_VECSMALL)
1575 5973746 : for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
1576 : else
1577 3809957 : for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
1578 2671569 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
1579 2671814 : gel(V,i) = gerepileuptoint(av, c);
1580 : }
1581 553467 : return ZX_renormalize(V, l);
1582 : }
1583 :
1584 : static GEN
1585 4614 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1586 : {
1587 4614 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1588 4614 : GEN A = cgetg(n, t_VEC);
1589 4614 : GEN V = cgetg(l, t_COL);
1590 90903 : for (i=1; i < l; i++)
1591 : {
1592 335111 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1593 86289 : gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
1594 : }
1595 4614 : return V;
1596 : }
1597 :
1598 : static GEN
1599 114988 : polint_chinese(GEN worker, GEN mA, GEN P)
1600 : {
1601 114988 : long cnt, pending, n, i, j, l = lg(gel(mA,1));
1602 : struct pari_mt pt;
1603 : GEN done, va, M, A;
1604 : pari_timer ti;
1605 :
1606 114988 : if (l == 1) return cgetg(1, t_MAT);
1607 85936 : cnt = pending = 0;
1608 85936 : n = lg(P);
1609 85936 : A = cgetg(n, t_VEC);
1610 85936 : va = mkvec(A);
1611 85936 : M = cgetg(l, t_MAT);
1612 85936 : if (DEBUGLEVEL>4) timer_start(&ti);
1613 85936 : if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
1614 85936 : mt_queue_start_lim(&pt, worker, l-1);
1615 676057 : for (i=1; i<l || pending; i++)
1616 : {
1617 : long workid;
1618 2469587 : for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
1619 590121 : mt_queue_submit(&pt, i, i<l? va: NULL);
1620 590121 : done = mt_queue_get(&pt, &workid, &pending);
1621 590121 : if (done)
1622 : {
1623 555808 : gel(M,workid) = done;
1624 555808 : if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
1625 : }
1626 : }
1627 85936 : if (DEBUGLEVEL>5) err_printf("\n");
1628 85936 : if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
1629 85936 : mt_queue_end(&pt);
1630 85936 : return M;
1631 : }
1632 :
1633 : GEN
1634 840 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
1635 : {
1636 840 : return nxCV_polint_center_tree(vA, P, T, R, m2);
1637 : }
1638 :
1639 : static GEN
1640 430 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1641 : {
1642 430 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1643 430 : GEN A = cgetg(n, t_VEC);
1644 430 : GEN V = cgetg(l, t_MAT);
1645 4204 : for (i=1; i < l; i++)
1646 : {
1647 15299 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1648 3774 : gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
1649 : }
1650 430 : return V;
1651 : }
1652 :
1653 : static GEN
1654 90 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
1655 : {
1656 90 : GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
1657 90 : return polint_chinese(worker, mA, P);
1658 : }
1659 :
1660 : static GEN
1661 53115 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1662 : {
1663 53115 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1664 53115 : GEN A = cgetg(n, t_VEC);
1665 53114 : GEN V = cgetg(l, t_MAT);
1666 420230 : for (i=1; i < l; i++)
1667 : {
1668 2137025 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1669 367115 : gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
1670 : }
1671 53115 : return V;
1672 : }
1673 :
1674 : GEN
1675 554898 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
1676 : {
1677 554898 : return ncV_polint_center_tree(vA, P, T, R, m2);
1678 : }
1679 :
1680 : static GEN
1681 114898 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
1682 : {
1683 114898 : GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
1684 114898 : return polint_chinese(worker, mA, P);
1685 : }
1686 :
1687 : /* return [A mod P[i], i=1..#P] */
1688 : GEN
1689 0 : Z_ZV_mod(GEN A, GEN P)
1690 : {
1691 0 : pari_sp av = avma;
1692 0 : return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
1693 : }
1694 : /* P a t_VECSMALL */
1695 : GEN
1696 0 : Z_nv_mod(GEN A, GEN P)
1697 : {
1698 0 : pari_sp av = avma;
1699 0 : return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
1700 : }
1701 : /* B a ZX, T = ZV_producttree(P) */
1702 : GEN
1703 1263349 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
1704 : {
1705 : pari_sp av;
1706 1263349 : long i, j, l = lg(B), n = lg(A)-1;
1707 1263349 : GEN V = cgetg(n+1, t_VEC);
1708 6210893 : for (j=1; j <= n; j++)
1709 : {
1710 4947633 : gel(V, j) = cgetg(l, t_VECSMALL);
1711 4947560 : mael(V, j, 1) = B[1]&VARNBITS;
1712 : }
1713 1263260 : av = avma;
1714 12136546 : for (i=2; i < l; i++)
1715 : {
1716 10873695 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
1717 75290334 : for (j=1; j <= n; j++)
1718 64421866 : mael(V, j, i) = v[j];
1719 10868468 : set_avma(av);
1720 : }
1721 6211100 : for (j=1; j <= n; j++)
1722 4948200 : (void) Flx_renormalize(gel(V, j), l);
1723 1262900 : return V;
1724 : }
1725 :
1726 : static GEN
1727 169745 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
1728 :
1729 : GEN
1730 11939 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
1731 : {
1732 11939 : pari_sp av = avma;
1733 11939 : long i, j, l = lg(P), n = lg(xa)-1;
1734 11939 : GEN V = cgetg(n+1, t_VEC);
1735 46965 : for (j=1; j <= n; j++)
1736 : {
1737 35026 : gel(V, j) = cgetg(l, t_POL);
1738 35026 : mael(V, j, 1) = P[1]&VARNBITS;
1739 : }
1740 99692 : for (i=2; i < l; i++)
1741 : {
1742 87753 : GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
1743 362381 : for (j=1; j <= n; j++)
1744 274628 : gmael(V, j, i) = gel(v,j);
1745 : }
1746 46965 : for (j=1; j <= n; j++)
1747 35026 : (void) FlxX_renormalize(gel(V, j), l);
1748 11939 : return gerepilecopy(av, V);
1749 : }
1750 :
1751 : GEN
1752 4045 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
1753 : {
1754 4045 : pari_sp av = avma;
1755 4045 : long i, j, l = lg(C), n = lg(xa)-1;
1756 4045 : GEN V = cgetg(n+1, t_VEC);
1757 16915 : for (j = 1; j <= n; j++)
1758 12870 : gel(V, j) = cgetg(l, t_COL);
1759 86037 : for (i = 1; i < l; i++)
1760 : {
1761 81993 : GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
1762 360579 : for (j = 1; j <= n; j++)
1763 278587 : gmael(V, j, i) = gel(v,j);
1764 : }
1765 4044 : return gerepilecopy(av, V);
1766 : }
1767 :
1768 : GEN
1769 430 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
1770 : {
1771 430 : pari_sp av = avma;
1772 430 : long i, j, l = lg(M), n = lg(xa)-1;
1773 430 : GEN V = cgetg(n+1, t_VEC);
1774 2083 : for (j=1; j <= n; j++)
1775 1653 : gel(V, j) = cgetg(l, t_MAT);
1776 4204 : for (i=1; i < l; i++)
1777 : {
1778 3774 : GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
1779 15299 : for (j=1; j <= n; j++)
1780 11525 : gmael(V, j, i) = gel(v,j);
1781 : }
1782 430 : return gerepilecopy(av, V);
1783 : }
1784 :
1785 : GEN
1786 941093 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
1787 : {
1788 : pari_sp av;
1789 941093 : long i, j, l = lg(B), n = lg(A)-1;
1790 941093 : GEN V = cgetg(n+1, t_VEC);
1791 4667977 : for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
1792 941028 : av = avma;
1793 41408892 : for (i=1; i < l; i++)
1794 : {
1795 40476190 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
1796 228740005 : for (j=1; j <= n; j++) mael(V, j, i) = v[j];
1797 40440222 : set_avma(av);
1798 : }
1799 932702 : return V;
1800 : }
1801 :
1802 : GEN
1803 80998 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
1804 : {
1805 80998 : pari_sp av = avma;
1806 80998 : long i, j, l = lg(M), n = lg(xa)-1;
1807 80998 : GEN V = cgetg(n+1, t_VEC);
1808 317049 : for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_MAT);
1809 1021856 : for (i=1; i < l; i++)
1810 : {
1811 940866 : GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
1812 4667441 : for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
1813 : }
1814 80990 : return gerepilecopy(av, V);
1815 : }
1816 :
1817 : static GEN
1818 2166733 : ZV_sqr(GEN z)
1819 : {
1820 2166733 : long i,l = lg(z);
1821 2166733 : GEN x = cgetg(l, t_VEC);
1822 2166753 : if (typ(z)==t_VECSMALL)
1823 5041388 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
1824 : else
1825 3842661 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
1826 2166709 : return x;
1827 : }
1828 :
1829 : static GEN
1830 11101186 : ZT_sqr(GEN x)
1831 : {
1832 11101186 : if (typ(x) == t_INT) return sqri(x);
1833 14542189 : pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
1834 : }
1835 :
1836 : static GEN
1837 2166730 : ZV_invdivexact(GEN y, GEN x)
1838 : {
1839 2166730 : long i, l = lg(y);
1840 2166730 : GEN z = cgetg(l,t_VEC);
1841 2166766 : if (typ(x)==t_VECSMALL)
1842 5041117 : for (i=1; i<l; i++)
1843 : {
1844 3997236 : pari_sp av = avma;
1845 3997236 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
1846 3997493 : set_avma(av); gel(z,i) = utoi(a);
1847 : }
1848 : else
1849 3842774 : for (i=1; i<l; i++)
1850 2719903 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
1851 2166752 : return z;
1852 : }
1853 :
1854 : /* P t_VECSMALL or t_VEC of t_INT */
1855 : GEN
1856 2166741 : ZV_chinesetree(GEN P, GEN T)
1857 : {
1858 2166741 : GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
1859 2166729 : GEN mod = gmael(T,lg(T)-1,1);
1860 2166729 : return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
1861 : }
1862 :
1863 : static GEN
1864 646491 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
1865 : {
1866 646491 : if (!pt_mod)
1867 12209 : return gerepileupto(av, a);
1868 : else
1869 : {
1870 634282 : GEN mod = gmael(T, lg(T)-1, 1);
1871 634282 : gerepileall(av, 2, &a, &mod);
1872 634281 : *pt_mod = mod;
1873 634281 : return a;
1874 : }
1875 : }
1876 :
1877 : GEN
1878 162125 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1879 : {
1880 162125 : pari_sp av = avma;
1881 162125 : GEN T = ZV_producttree(P);
1882 162126 : GEN R = ZV_chinesetree(P, T);
1883 162126 : GEN a = ZV_chinese_tree(A, P, T, R);
1884 162126 : GEN mod = gmael(T, lg(T)-1, 1);
1885 162126 : GEN ca = Fp_center(a, mod, shifti(mod,-1));
1886 162126 : return gc_chinese(av, T, ca, pt_mod);
1887 : }
1888 :
1889 : GEN
1890 5095 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
1891 : {
1892 5095 : pari_sp av = avma;
1893 5095 : GEN T = ZV_producttree(P);
1894 5095 : GEN R = ZV_chinesetree(P, T);
1895 5095 : GEN a = ZV_chinese_tree(A, P, T, R);
1896 5095 : return gc_chinese(av, T, a, pt_mod);
1897 : }
1898 :
1899 : GEN
1900 113129 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
1901 : {
1902 113129 : pari_sp av = avma;
1903 113129 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1904 113129 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
1905 113129 : return gerepileupto(av, a);
1906 : }
1907 :
1908 : GEN
1909 354046 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1910 : {
1911 354046 : pari_sp av = avma;
1912 354046 : GEN T = ZV_producttree(P);
1913 354046 : GEN R = ZV_chinesetree(P, T);
1914 354045 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1915 354044 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
1916 354045 : return gc_chinese(av, T, a, pt_mod);
1917 : }
1918 :
1919 : GEN
1920 10237 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1921 : {
1922 10237 : pari_sp av = avma;
1923 10237 : GEN T = ZV_producttree(P);
1924 10237 : GEN R = ZV_chinesetree(P, T);
1925 10237 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1926 10237 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
1927 10237 : return gc_chinese(av, T, a, pt_mod);
1928 : }
1929 :
1930 : GEN
1931 5530 : ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
1932 : {
1933 5530 : pari_sp av = avma;
1934 5530 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1935 5530 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
1936 5530 : return gerepileupto(av, a);
1937 : }
1938 :
1939 : GEN
1940 0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
1941 : {
1942 0 : pari_sp av = avma;
1943 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1944 0 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
1945 0 : return gerepileupto(av, a);
1946 : }
1947 :
1948 : GEN
1949 53115 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
1950 : {
1951 53115 : pari_sp av = avma;
1952 53115 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1953 53115 : GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
1954 53115 : return gerepileupto(av, a);
1955 : }
1956 :
1957 : GEN
1958 114898 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1959 : {
1960 114898 : pari_sp av = avma;
1961 114898 : GEN T = ZV_producttree(P);
1962 114898 : GEN R = ZV_chinesetree(P, T);
1963 114898 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1964 114898 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
1965 114898 : return gc_chinese(av, T, a, pt_mod);
1966 : }
1967 :
1968 : GEN
1969 0 : nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
1970 : {
1971 0 : pari_sp av = avma;
1972 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1973 0 : GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
1974 0 : return gerepileupto(av, a);
1975 : }
1976 :
1977 : GEN
1978 0 : nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1979 : {
1980 0 : pari_sp av = avma;
1981 0 : GEN T = ZV_producttree(P);
1982 0 : GEN R = ZV_chinesetree(P, T);
1983 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1984 0 : GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
1985 0 : return gc_chinese(av, T, a, pt_mod);
1986 : }
1987 :
1988 : GEN
1989 430 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
1990 : {
1991 430 : pari_sp av = avma;
1992 430 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
1993 430 : GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
1994 430 : return gerepileupto(av, a);
1995 : }
1996 :
1997 : GEN
1998 90 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1999 : {
2000 90 : pari_sp av = avma;
2001 90 : GEN T = ZV_producttree(P);
2002 90 : GEN R = ZV_chinesetree(P, T);
2003 90 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2004 90 : GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
2005 90 : return gc_chinese(av, T, a, pt_mod);
2006 : }
2007 :
2008 : /**********************************************************************
2009 : ** Powering over (Z/NZ)^*, small N **
2010 : **********************************************************************/
2011 :
2012 : /* 2^n mod p; assume n > 1 */
2013 : static ulong
2014 10406985 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
2015 : {
2016 10406985 : ulong y = 2;
2017 10406985 : int j = 1+bfffo(n);
2018 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
2019 10406985 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
2020 444967317 : for (; j; n<<=1,j--)
2021 : {
2022 434611911 : y = Fl_sqr_pre(y,p,pi);
2023 434675646 : if (n & HIGHBIT) y = Fl_double(y, p);
2024 : }
2025 10355406 : return y;
2026 : }
2027 :
2028 : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
2029 : static ulong
2030 4172176 : Fl_2powu(ulong n, ulong p)
2031 : {
2032 4172176 : ulong y = 2;
2033 4172176 : int j = 1+bfffo(n);
2034 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
2035 4172176 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
2036 23878999 : for (; j; n<<=1,j--)
2037 : {
2038 19706820 : y = (y*y) % p;
2039 19706820 : if (n & HIGHBIT) y = Fl_double(y, p);
2040 : }
2041 4172179 : return y;
2042 : }
2043 :
2044 : /* allow pi = 0 */
2045 : ulong
2046 92138713 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
2047 : {
2048 : ulong y, z, n;
2049 92138713 : if (!pi) return Fl_powu(x, n0, p);
2050 89806538 : if (n0 <= 1)
2051 : { /* frequent special cases */
2052 5520693 : if (n0 == 1) return x;
2053 81176 : if (n0 == 0) return 1;
2054 : }
2055 84285847 : if (x <= 2)
2056 : {
2057 10450806 : if (x == 2) return Fl_2powu_pre(n0, p, pi);
2058 42689 : return x; /* 0 or 1 */
2059 : }
2060 73835041 : y = 1; z = x; n = n0;
2061 : for(;;)
2062 : {
2063 529918699 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
2064 529808295 : n>>=1; if (!n) return y;
2065 455835890 : z = Fl_sqr_pre(z,p,pi);
2066 : }
2067 : }
2068 :
2069 : ulong
2070 144946162 : Fl_powu(ulong x, ulong n0, ulong p)
2071 : {
2072 : ulong y, z, n;
2073 144946162 : if (n0 <= 2)
2074 : { /* frequent special cases */
2075 65326515 : if (n0 == 2) return Fl_sqr(x,p);
2076 31910467 : if (n0 == 1) return x;
2077 1883736 : if (n0 == 0) return 1;
2078 : }
2079 79588034 : if (x <= 1) return x; /* 0 or 1 */
2080 79167748 : if (p & HIGHMASK)
2081 7881972 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
2082 71285776 : if (x == 2) return Fl_2powu(n0, p);
2083 67113591 : y = 1; z = x; n = n0;
2084 : for(;;)
2085 : {
2086 354463454 : if (n&1) y = (y*z) % p;
2087 354463454 : n>>=1; if (!n) return y;
2088 287349863 : z = (z*z) % p;
2089 : }
2090 : }
2091 :
2092 : /* Reduce data dependency to maximize internal parallelism; allow pi = 0 */
2093 : GEN
2094 12768563 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
2095 : {
2096 : long i, k;
2097 12768563 : GEN z = cgetg(n + 2, t_VECSMALL);
2098 12757698 : z[1] = 1; if (n == 0) return z;
2099 12757698 : z[2] = x;
2100 12757698 : if (pi)
2101 : {
2102 90048449 : for (i = 3, k=2; i <= n; i+=2, k++)
2103 : {
2104 77494034 : z[i] = Fl_sqr_pre(z[k], p, pi);
2105 77499431 : z[i+1] = Fl_mul_pre(z[k], z[k+1], p, pi);
2106 : }
2107 12554415 : if (i==n+1) z[i] = Fl_sqr_pre(z[k], p, pi);
2108 : }
2109 212378 : else if (p & HIGHMASK)
2110 : {
2111 0 : for (i = 3, k=2; i <= n; i+=2, k++)
2112 : {
2113 0 : z[i] = Fl_sqr(z[k], p);
2114 0 : z[i+1] = Fl_mul(z[k], z[k+1], p);
2115 : }
2116 0 : if (i==n+1) z[i] = Fl_sqr(z[k], p);
2117 : }
2118 : else
2119 400517421 : for (i = 2; i <= n; i++) z[i+1] = (z[i] * x) % p;
2120 12768124 : return z;
2121 : }
2122 :
2123 : GEN
2124 295824 : Fl_powers(ulong x, long n, ulong p)
2125 : {
2126 295824 : return Fl_powers_pre(x, n, p, (p & HIGHMASK)? get_Fl_red(p): 0);
2127 : }
2128 :
2129 : /**********************************************************************
2130 : ** Powering over (Z/NZ)^*, large N **
2131 : **********************************************************************/
2132 :
2133 : static GEN
2134 4565239 : Fp_dblsqr(GEN x, GEN N)
2135 : {
2136 4565239 : GEN z = shifti(Fp_sqr(x, N), 1);
2137 4466924 : return cmpii(z, N) >= 0? subii(z, N): z;
2138 : }
2139 :
2140 : typedef struct muldata {
2141 : GEN (*sqr)(void * E, GEN x);
2142 : GEN (*mul)(void * E, GEN x, GEN y);
2143 : GEN (*mul2)(void * E, GEN x);
2144 : } muldata;
2145 :
2146 : /* modified Barrett reduction with one fold */
2147 : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
2148 :
2149 : static GEN
2150 12951 : Fp_invmBarrett(GEN p, long s)
2151 : {
2152 12951 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
2153 12951 : return mkvec2(Q,R);
2154 : }
2155 :
2156 : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
2157 : * a = r (mod N) */
2158 : static GEN
2159 6401240 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
2160 : {
2161 6401240 : pari_sp av = avma;
2162 6401240 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
2163 6401240 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
2164 6401240 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
2165 6401240 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
2166 6401240 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
2167 6401240 : GEN r = subii(A, mulii(q, N));
2168 6401240 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
2169 6401240 : if (signe(sr)<0) return gerepileuptoint(av, r);
2170 3343618 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
2171 3343618 : if (signe(sr)<0) return gerepileuptoint(av, r);
2172 105378 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
2173 105378 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
2174 : }
2175 :
2176 : /* Montgomery reduction */
2177 :
2178 : INLINE ulong
2179 864342 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
2180 :
2181 : struct montred
2182 : {
2183 : GEN N;
2184 : ulong inv;
2185 : };
2186 :
2187 : /* Montgomery reduction */
2188 : static GEN
2189 61131223 : _sqr_montred(void * E, GEN x)
2190 : {
2191 61131223 : struct montred * D = (struct montred *) E;
2192 61131223 : return red_montgomery(sqri(x), D->N, D->inv);
2193 : }
2194 :
2195 : /* Montgomery reduction */
2196 : static GEN
2197 7371143 : _mul_montred(void * E, GEN x, GEN y)
2198 : {
2199 7371143 : struct montred * D = (struct montred *) E;
2200 7371143 : return red_montgomery(mulii(x, y), D->N, D->inv);
2201 : }
2202 :
2203 : static GEN
2204 6863019 : _mul2_montred(void * E, GEN x)
2205 : {
2206 6863019 : struct montred * D = (struct montred *) E;
2207 6863019 : GEN z = shifti(_sqr_montred(E, x), 1);
2208 6862567 : long l = lgefint(D->N);
2209 7289776 : while (lgefint(z) > l) z = subii(z, D->N);
2210 6862808 : return z;
2211 : }
2212 :
2213 : static GEN
2214 20454433 : _sqr_remii(void* N, GEN x)
2215 20454433 : { return remii(sqri(x), (GEN) N); }
2216 :
2217 : static GEN
2218 1598651 : _mul_remii(void* N, GEN x, GEN y)
2219 1598651 : { return remii(mulii(x, y), (GEN) N); }
2220 :
2221 : static GEN
2222 3221120 : _mul2_remii(void* N, GEN x)
2223 3221120 : { return Fp_dblsqr(x, (GEN) N); }
2224 :
2225 : struct redbarrett
2226 : {
2227 : GEN iM, N;
2228 : long s;
2229 : };
2230 :
2231 : static GEN
2232 5626446 : _sqr_remiibar(void *E, GEN x)
2233 : {
2234 5626446 : struct redbarrett * D = (struct redbarrett *) E;
2235 5626446 : return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
2236 : }
2237 :
2238 : static GEN
2239 774794 : _mul_remiibar(void *E, GEN x, GEN y)
2240 : {
2241 774794 : struct redbarrett * D = (struct redbarrett *) E;
2242 774794 : return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
2243 : }
2244 :
2245 : static GEN
2246 1342530 : _mul2_remiibar(void *E, GEN x)
2247 : {
2248 1342530 : struct redbarrett * D = (struct redbarrett *) E;
2249 1342530 : return Fp_dblsqr(x, D->N);
2250 : }
2251 :
2252 : static long
2253 1084957 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
2254 : {
2255 1084957 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
2256 : {
2257 12951 : struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
2258 12951 : D->sqr = &_sqr_remiibar;
2259 12951 : D->mul = &_mul_remiibar;
2260 12951 : D->mul2 = &_mul2_remiibar;
2261 12951 : E->N = N;
2262 12951 : E->s = 1+(expi(N)>>1);
2263 12951 : E->iM = Fp_invmBarrett(N, E->s);
2264 12951 : *pt_E = (void*) E;
2265 12951 : return 0;
2266 : }
2267 1072006 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
2268 : {
2269 864341 : struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
2270 864341 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
2271 864345 : D->sqr = &_sqr_montred;
2272 864345 : D->mul = &_mul_montred;
2273 864345 : D->mul2 = &_mul2_montred;
2274 864345 : E->N = N;
2275 864345 : E->inv = init_montdata(N);
2276 864335 : *pt_E = (void*) E;
2277 864335 : return 1;
2278 : }
2279 : else
2280 : {
2281 207660 : D->sqr = &_sqr_remii;
2282 207660 : D->mul = &_mul_remii;
2283 207660 : D->mul2 = &_mul2_remii;
2284 207660 : *pt_E = (void*) N;
2285 207660 : return 0;
2286 : }
2287 : }
2288 :
2289 : GEN
2290 3107374 : Fp_powu(GEN A, ulong k, GEN N)
2291 : {
2292 3107374 : long lN = lgefint(N);
2293 : int base_is_2, use_montgomery;
2294 : muldata D;
2295 : void *E;
2296 : pari_sp av;
2297 :
2298 3107374 : if (lN == 3) {
2299 1443355 : ulong n = uel(N,2);
2300 1443355 : return utoi( Fl_powu(umodiu(A, n), k, n) );
2301 : }
2302 1664019 : if (k <= 2)
2303 : { /* frequent special cases */
2304 837425 : if (k == 2) return Fp_sqr(A,N);
2305 279860 : if (k == 1) return A;
2306 0 : if (k == 0) return gen_1;
2307 : }
2308 826593 : av = avma; A = modii(A,N);
2309 826594 : base_is_2 = 0;
2310 826594 : if (lgefint(A) == 3) switch(A[2])
2311 : {
2312 832 : case 1: set_avma(av); return gen_1;
2313 42003 : case 2: base_is_2 = 1; break;
2314 : }
2315 :
2316 : /* TODO: Move this out of here and use for general modular computations */
2317 825762 : use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
2318 825762 : if (base_is_2)
2319 42003 : A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
2320 : else
2321 783759 : A = gen_powu_i(A, k, E, D.sqr, D.mul);
2322 825762 : if (use_montgomery)
2323 : {
2324 704844 : A = red_montgomery(A, N, ((struct montred *) E)->inv);
2325 704844 : if (cmpii(A, N) >= 0) A = subii(A,N);
2326 : }
2327 825762 : return gerepileuptoint(av, A);
2328 : }
2329 :
2330 : GEN
2331 29410 : Fp_pows(GEN A, long k, GEN N)
2332 : {
2333 29410 : if (lgefint(N) == 3) {
2334 14921 : ulong n = N[2];
2335 14921 : ulong a = umodiu(A, n);
2336 14921 : if (k < 0) {
2337 133 : a = Fl_inv(a, n);
2338 133 : k = -k;
2339 : }
2340 14921 : return utoi( Fl_powu(a, (ulong)k, n) );
2341 : }
2342 14489 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
2343 14489 : return Fp_powu(A, (ulong)k, N);
2344 : }
2345 :
2346 : /* A^K mod N */
2347 : GEN
2348 35964868 : Fp_pow(GEN A, GEN K, GEN N)
2349 : {
2350 : pari_sp av;
2351 35964868 : long s, lN = lgefint(N), sA, sy;
2352 : int base_is_2, use_montgomery;
2353 : GEN y;
2354 : muldata D;
2355 : void *E;
2356 :
2357 35964868 : s = signe(K);
2358 35964868 : if (!s) return dvdii(A,N)? gen_0: gen_1;
2359 34933592 : if (lN == 3 && lgefint(K) == 3)
2360 : {
2361 33979222 : ulong n = N[2], a = umodiu(A, n);
2362 33979484 : if (s < 0) a = Fl_inv(a, n);
2363 33979548 : if (a <= 1) return utoi(a); /* 0 or 1 */
2364 30481069 : return utoi(Fl_powu(a, uel(K,2), n));
2365 : }
2366 :
2367 954370 : av = avma;
2368 954370 : if (s < 0) y = Fp_inv(A,N);
2369 : else
2370 : {
2371 952415 : y = modii(A,N);
2372 952616 : if (!signe(y)) { set_avma(av); return gen_0; }
2373 : }
2374 954571 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
2375 :
2376 259365 : base_is_2 = 0;
2377 259365 : sy = abscmpii(y, shifti(N,-1)) > 0;
2378 259387 : if (sy) y = subii(N,y);
2379 259398 : sA = sy && mod2(K);
2380 259398 : if (lgefint(y) == 3) switch(y[2])
2381 : {
2382 200 : case 1: set_avma(av); return sA ? subis(N,1): gen_1;
2383 132923 : case 2: base_is_2 = 1; break;
2384 : }
2385 :
2386 : /* TODO: Move this out of here and use for general modular computations */
2387 259198 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
2388 259185 : if (base_is_2)
2389 132923 : y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
2390 : else
2391 126262 : y = gen_pow_i(y, K, E, D.sqr, D.mul);
2392 259205 : if (use_montgomery)
2393 : {
2394 159509 : y = red_montgomery(y, N, ((struct montred *) E)->inv);
2395 159508 : if (cmpii(y,N) >= 0) y = subii(y,N);
2396 : }
2397 259203 : if (sA) y = subii(N, y);
2398 259203 : return gerepileuptoint(av,y);
2399 : }
2400 :
2401 : static GEN
2402 14207608 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
2403 : static GEN
2404 8134365 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
2405 : static GEN
2406 57564 : _Fp_one(void *E) { (void) E; return gen_1; }
2407 :
2408 : GEN
2409 77 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
2410 77 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
2411 :
2412 : GEN
2413 54096 : Fp_pow_table(GEN R, GEN n, GEN p)
2414 54096 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
2415 :
2416 : GEN
2417 5924 : Fp_powers(GEN x, long n, GEN p)
2418 : {
2419 5924 : if (lgefint(p) == 3)
2420 2456 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
2421 3468 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
2422 : }
2423 :
2424 : GEN
2425 497 : FpV_prod(GEN V, GEN p) { return gen_product(V, (void *)p, &_Fp_mul); }
2426 :
2427 : static GEN
2428 27743391 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
2429 : static GEN
2430 154 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
2431 :
2432 : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
2433 : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
2434 : equalii,equali1,Fp_easylog};
2435 :
2436 : static GEN
2437 782379 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
2438 : static GEN
2439 936403 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
2440 : static GEN
2441 844296 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
2442 : static GEN
2443 513167 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
2444 : static GEN
2445 33113 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
2446 : static int
2447 222261 : _Fp_equal0(GEN x) { return signe(x)==0; }
2448 : static GEN
2449 28726 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
2450 :
2451 : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
2452 : _Fp_inv,_Fp_equal0,_Fp_s};
2453 :
2454 7093 : const struct bb_field *get_Fp_field(void **E, GEN p)
2455 7093 : { *E = (void*)p; return &Fp_field; }
2456 :
2457 : /*********************************************************************/
2458 : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
2459 : /*********************************************************************/
2460 : ulong
2461 481619 : Fl_order(ulong a, ulong o, ulong p)
2462 : {
2463 481619 : pari_sp av = avma;
2464 : GEN m, P, E;
2465 : long i;
2466 481619 : if (a==1) return 1;
2467 438733 : if (!o) o = p-1;
2468 438733 : m = factoru(o);
2469 438733 : P = gel(m,1);
2470 438733 : E = gel(m,2);
2471 1250931 : for (i = lg(P)-1; i; i--)
2472 : {
2473 812198 : ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
2474 812198 : if (y == 1) o = t;
2475 772762 : else for (j = 1; j < e; j++)
2476 : {
2477 385554 : y = Fl_powu(y, l, p);
2478 385554 : if (y == 1) { o = t * upowuu(l, j); break; }
2479 : }
2480 : }
2481 438733 : return gc_ulong(av, o);
2482 : }
2483 :
2484 : /*Find the exact order of a assuming a^o==1*/
2485 : GEN
2486 10656 : Fp_order(GEN a, GEN o, GEN p) {
2487 10656 : if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
2488 : {
2489 21 : ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
2490 21 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
2491 : }
2492 10635 : return gen_order(a, o, (void*)p, &Fp_star);
2493 : }
2494 : GEN
2495 70 : Fp_factored_order(GEN a, GEN o, GEN p)
2496 70 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
2497 :
2498 : /* return order of a mod p^e, e > 0, pe = p^e */
2499 : static GEN
2500 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
2501 : {
2502 : GEN ap, op;
2503 70 : if (absequaliu(p, 2))
2504 : {
2505 56 : if (e == 1) return gen_1;
2506 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
2507 49 : if (mod4(a) == 1) op = gen_1; else { op = gen_2; a = Fp_sqr(a, pe); }
2508 : } else {
2509 14 : ap = (e == 1)? a: remii(a,p);
2510 14 : op = Fp_order(ap, subiu(p,1), p);
2511 14 : if (e == 1) return op;
2512 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
2513 : }
2514 49 : if (equali1(a)) return op;
2515 7 : return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
2516 : }
2517 :
2518 : GEN
2519 63 : znorder(GEN x, GEN o)
2520 : {
2521 63 : pari_sp av = avma;
2522 : GEN b, a;
2523 :
2524 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
2525 56 : b = gel(x,1); a = gel(x,2);
2526 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
2527 49 : if (!o)
2528 : {
2529 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
2530 35 : long i, l = lg(P);
2531 35 : o = gen_1;
2532 70 : for (i = 1; i < l; i++)
2533 : {
2534 35 : GEN p = gel(P,i);
2535 35 : long e = itos(gel(E,i));
2536 :
2537 35 : if (l == 2)
2538 35 : o = Zp_order(a, p, e, b);
2539 : else {
2540 0 : GEN pe = powiu(p,e);
2541 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
2542 : }
2543 : }
2544 35 : return gerepileuptoint(av, o);
2545 : }
2546 14 : return Fp_order(a, o, b);
2547 : }
2548 :
2549 : /*********************************************************************/
2550 : /** DISCRETE LOGARITHM in (Z/nZ)* **/
2551 : /*********************************************************************/
2552 : static GEN
2553 58305 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
2554 : {
2555 58305 : pari_sp av = avma;
2556 : GEN h1, h2, F, G;
2557 58305 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
2558 35013 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
2559 : {
2560 126 : GEN M = cgetg(3, t_MAT);
2561 126 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
2562 126 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
2563 126 : return gerepileupto(av, M);
2564 : }
2565 34887 : return gc_NULL(av);
2566 : }
2567 :
2568 : static GEN
2569 58305 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
2570 : {
2571 : GEN rel;
2572 58305 : do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
2573 58305 : while (!rel);
2574 126 : return rel;
2575 : }
2576 :
2577 : struct Fp_log_rel
2578 : {
2579 : GEN rel;
2580 : ulong prmax;
2581 : long nbrel, nbmax, nbgen;
2582 : };
2583 :
2584 : /* add u^e */
2585 : static void
2586 3157 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
2587 : {
2588 3157 : pari_sp av = avma;
2589 3157 : long off = r->prmax+1;
2590 3157 : GEN F = cgetg(3, t_MAT);
2591 3157 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
2592 3157 : gel(F,2) = vecsmall_append(gel(z,2), e);
2593 3157 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
2594 3157 : }
2595 :
2596 : /* add u^-1 v^-1 */
2597 : static void
2598 104083 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
2599 : {
2600 104083 : pari_sp av = avma;
2601 104083 : long off = r->prmax+1;
2602 104083 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
2603 104083 : GEN F = cgetg(3, t_MAT);
2604 104083 : gel(F,1) = vecsmall_concat(gel(z,1), P);
2605 104083 : gel(F,2) = vecsmall_concat(gel(z,2), E);
2606 104083 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
2607 104083 : }
2608 :
2609 : /*
2610 : Let p=C^2+c
2611 : Solve h = (C+x)*(C+a)-p = 0 [mod l]
2612 : h= -c+x*(C+a)+C*a = 0 [mod l]
2613 : x = (c-C*a)/(C+a) [mod l]
2614 : h = -c+C*(x+a)+a*x
2615 : */
2616 :
2617 : GEN
2618 40831 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
2619 : {
2620 40831 : pari_sp ltop = avma;
2621 40831 : long i, j, th, n = lg(pi)-1, rel = 1;
2622 40831 : GEN sieve = zero_zv(a+2)+1;
2623 40866 : GEN L = cgetg(1+a+2, t_VEC);
2624 40853 : pari_sp av = avma;
2625 40853 : GEN z, h = addis(C,a);
2626 40818 : if ((z = Z_issmooth_fact(h, prmax)))
2627 : {
2628 3006 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
2629 3014 : av = avma;
2630 : }
2631 16818295 : for (i=1; i<=n; i++)
2632 : {
2633 16802230 : ulong li = pi[i], s = sz[i], al = a % li;
2634 16802230 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
2635 17102966 : if (!iv) continue;
2636 16679079 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
2637 77374746 : for(j = u; j<=a; j+=li) sieve[j] += s;
2638 : }
2639 35046 : if (a)
2640 : {
2641 40740 : long e = expi(mulis(C,a));
2642 40777 : th = e - expu(e) - 1;
2643 54 : } else th = -1;
2644 27986456 : for (j=0; j<a; j++)
2645 27944026 : if (sieve[j]>=th)
2646 : {
2647 119437 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
2648 119343 : if ((z = Z_issmooth_fact(h, prmax)))
2649 : {
2650 109745 : gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
2651 110003 : av = avma;
2652 9423 : } else set_avma(av);
2653 : }
2654 : /* j = a */
2655 42430 : if (sieve[a]>=th)
2656 : {
2657 476 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
2658 476 : if ((z = Z_issmooth_fact(h, prmax)))
2659 385 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
2660 : }
2661 42430 : setlg(L, rel); return gerepilecopy(ltop, L);
2662 : }
2663 :
2664 : static long
2665 63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
2666 : {
2667 : struct pari_mt pt;
2668 63 : long i, j, nb = 0;
2669 63 : GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
2670 : mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
2671 63 : long running, pending = 0;
2672 63 : GEN W = zerovec(r->nbgen);
2673 63 : mt_queue_start_lim(&pt, worker, r->nbgen);
2674 41229 : for (i = 0; (running = (i < r->nbgen)) || pending; i++)
2675 : {
2676 : GEN done;
2677 : long idx;
2678 41166 : mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
2679 41166 : done = mt_queue_get(&pt, &idx, &pending);
2680 41166 : if (!done || lg(done)==1) continue;
2681 35917 : gel(W, idx+1) = done;
2682 35917 : nb += lg(done)-1;
2683 35917 : if (DEBUGLEVEL && (i&127)==0)
2684 0 : err_printf("%ld%% ",100*nb/r->nbmax);
2685 : }
2686 63 : mt_queue_end(&pt);
2687 39550 : for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
2688 : {
2689 : long ll, m;
2690 39487 : GEN L = gel(W,j);
2691 39487 : if (isintzero(L)) continue;
2692 34531 : ll = lg(L);
2693 141771 : for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
2694 : {
2695 107240 : GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
2696 107240 : if (v[1] == 1)
2697 3157 : addifsmooth1(r, h, v[2], v[3]);
2698 : else
2699 104083 : addifsmooth2(r, h, v[2], v[3]);
2700 : }
2701 : }
2702 63 : return j;
2703 : }
2704 :
2705 : static GEN
2706 735 : ECP_psi(GEN x, GEN y)
2707 : {
2708 735 : long prec = realprec(x);
2709 735 : GEN lx = glog(x, prec), ly = glog(y, prec);
2710 735 : GEN u = gdiv(lx, ly);
2711 735 : return gpow(u, gneg(u),prec);
2712 : }
2713 :
2714 : struct computeG
2715 : {
2716 : GEN C;
2717 : long bnd, nbi;
2718 : };
2719 :
2720 : static GEN
2721 735 : _computeG(void *E, GEN gen)
2722 : {
2723 735 : struct computeG * d = (struct computeG *) E;
2724 735 : GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
2725 735 : return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
2726 : }
2727 :
2728 : static long
2729 63 : compute_nbgen(GEN C, long bnd, long nbi)
2730 : {
2731 : struct computeG d;
2732 63 : d.C = shifti(C, 1);
2733 63 : d.bnd = bnd;
2734 63 : d.nbi = nbi;
2735 63 : return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
2736 : }
2737 :
2738 : static GEN
2739 1714 : _psi(void*E, GEN y)
2740 : {
2741 1714 : GEN lx = (GEN) E;
2742 1714 : long prec = realprec(lx);
2743 1714 : GEN ly = glog(y, prec);
2744 1714 : GEN u = gdiv(lx, ly);
2745 1714 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
2746 : }
2747 :
2748 : static GEN
2749 63 : opt_param(GEN x, long prec)
2750 : {
2751 63 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
2752 : }
2753 :
2754 : static GEN
2755 63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
2756 : {
2757 63 : pari_sp av = avma;
2758 63 : long lM = lg(M)-1, nbcol = lM;
2759 63 : long tbs = maxss(1, expu(nbg/expi(m)));
2760 : for (;;)
2761 14 : {
2762 77 : GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
2763 : GEN tab;
2764 77 : long i, f=0;
2765 77 : long l = lg(K), lm = lgefint(m);
2766 77 : GEN idx = diviiexact(subiu(p,1),m), g;
2767 : pari_timer ti;
2768 77 : if (DEBUGLEVEL) timer_start(&ti);
2769 154 : for(i=1; i<l; i++)
2770 154 : if (signe(gel(K,i)))
2771 77 : break;
2772 77 : g = Fp_pow(utoi(i), idx, p);
2773 77 : tab = Fp_pow_init(g, p, tbs, p);
2774 77 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
2775 128464 : for(i=1; i<l; i++)
2776 : {
2777 128387 : GEN k = gel(K,i);
2778 128387 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
2779 128387 : if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
2780 76391 : gel(K,i) = cgetineg(lm);
2781 : else
2782 51996 : f++;
2783 : }
2784 77 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
2785 77 : if(f > (nbg>>1)) return gerepileupto(av, K);
2786 4585 : for(i=1; i<=nbcol; i++)
2787 : {
2788 4571 : long a = 1+random_Fl(lM);
2789 4571 : swap(gel(M,a),gel(M,i));
2790 : }
2791 14 : if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
2792 : }
2793 : }
2794 :
2795 : static GEN
2796 126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
2797 : {
2798 126 : pari_sp av=avma;
2799 126 : GEN aa = gen_1;
2800 126 : long AV = 0;
2801 : for(;;)
2802 0 : {
2803 126 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
2804 126 : GEN F = gel(A,1), E = gel(A,2);
2805 126 : GEN Ao = gen_0;
2806 126 : long i, l = lg(F);
2807 953 : for(i=1; i<l; i++)
2808 : {
2809 827 : GEN Ki = gel(K,F[i]);
2810 827 : if (signe(Ki)<0) break;
2811 827 : Ao = addii(Ao, mulis(Ki, E[i]));
2812 : }
2813 126 : if (i==l) return Fp_divu(Ao, AV, m);
2814 0 : aa = gerepileuptoint(av, aa);
2815 : }
2816 : }
2817 :
2818 : static GEN
2819 63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
2820 : {
2821 63 : pari_sp av = avma, av2;
2822 63 : long i, j, nbi, nbr = 0, nbrow, nbg;
2823 : GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
2824 : pari_timer ti;
2825 : struct Fp_log_rel r;
2826 63 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
2827 63 : ulong bnd = 4*bnds;
2828 63 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
2829 :
2830 63 : p_1 = subiu(p,1);
2831 63 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
2832 0 : m = diviiexact(p_1, Z_ppo(p_1, m));
2833 63 : pr = primes_upto_zv(bnd);
2834 63 : nbi = lg(pr)-1;
2835 63 : C = sqrtremi(p, &c);
2836 63 : av2 = avma;
2837 12796 : for (i = 1; i <= nbi; ++i)
2838 : {
2839 12733 : ulong lp = pr[i];
2840 26894 : while (lp <= bnd)
2841 : {
2842 14161 : nbr++;
2843 14161 : lp *= pr[i];
2844 : }
2845 : }
2846 63 : pi = cgetg(nbr+1,t_VECSMALL);
2847 63 : Ci = cgetg(nbr+1,t_VECSMALL);
2848 63 : ci = cgetg(nbr+1,t_VECSMALL);
2849 63 : sz = cgetg(nbr+1,t_VECSMALL);
2850 12796 : for (i = 1, j = 1; i <= nbi; ++i)
2851 : {
2852 12733 : ulong lp = pr[i], sp = expu(2*lp-1);
2853 26894 : while (lp <= bnd)
2854 : {
2855 14161 : pi[j] = lp;
2856 14161 : Ci[j] = umodiu(C, lp);
2857 14161 : ci[j] = umodiu(c, lp);
2858 14161 : sz[j] = sp;
2859 14161 : lp *= pr[i];
2860 14161 : j++;
2861 : }
2862 : }
2863 63 : r.nbrel = 0;
2864 63 : r.nbgen = compute_nbgen(C, bnd, nbi);
2865 63 : r.nbmax = 2*(nbi+r.nbgen);
2866 63 : r.rel = cgetg(r.nbmax+1,t_VEC);
2867 63 : r.prmax = pr[nbi];
2868 63 : if (DEBUGLEVEL)
2869 : {
2870 0 : err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
2871 0 : timer_start(&ti);
2872 : }
2873 63 : nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
2874 63 : nbrow = r.prmax + nbg;
2875 63 : if (DEBUGLEVEL)
2876 : {
2877 0 : err_printf("\n");
2878 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
2879 : }
2880 63 : setlg(r.rel,r.nbrel+1);
2881 63 : r.rel = gerepilecopy(av2, r.rel);
2882 63 : K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
2883 63 : if (DEBUGLEVEL) timer_start(&ti);
2884 63 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
2885 63 : if (DEBUGLEVEL) timer_printf(&ti," log element");
2886 63 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
2887 63 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
2888 63 : d = gcdii(Ao,Bo);
2889 63 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
2890 63 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
2891 63 : return gerepileuptoint(av, l);
2892 : }
2893 :
2894 : static int
2895 4642731 : Fp_log_use_index(long e, long p)
2896 : {
2897 4642731 : return (e >= 27 && 20*(p+6)<=e*e);
2898 : }
2899 :
2900 : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
2901 : static GEN
2902 8442668 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
2903 : {
2904 8442668 : pari_sp av = avma;
2905 8442668 : GEN p = (GEN)E;
2906 : /* assume a reduced mod p, p not necessarily prime */
2907 8442668 : if (equali1(a)) return gen_0;
2908 : /* p > 2 */
2909 5421247 : if (equalii(subiu(p,1), a)) /* -1 */
2910 : {
2911 : pari_sp av2;
2912 : GEN t;
2913 1321379 : ord = get_arith_Z(ord);
2914 1321379 : if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
2915 1321365 : t = shifti(ord,-1); /* only possible solution */
2916 1321364 : av2 = avma;
2917 1321364 : if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
2918 1321335 : set_avma(av2); return gerepileuptoint(av, t);
2919 : }
2920 4099870 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
2921 63 : return Fp_log_index(a, g, ord, p);
2922 4099807 : return gc_NULL(av); /* not easy */
2923 : }
2924 :
2925 : GEN
2926 3913070 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
2927 : {
2928 3913070 : GEN v = get_arith_ZZM(ord);
2929 3913048 : GEN F = gmael(v,2,1);
2930 3913048 : long lF = lg(F)-1, lmax;
2931 3913048 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
2932 3913020 : lmax = expi(gel(F,lF));
2933 3913022 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
2934 91 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
2935 3913018 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
2936 : }
2937 :
2938 : /* assume !(p & HIGHMASK) */
2939 : static ulong
2940 130283 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
2941 : {
2942 130283 : ulong i, h=1;
2943 357162 : for (i = 0; i < ord; i++, h = (h * g) % p)
2944 357162 : if (a==h) return i;
2945 0 : return ~0UL;
2946 : }
2947 :
2948 : static ulong
2949 22215 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
2950 : {
2951 22215 : ulong i, h=1;
2952 56012 : for (i = 0; i < ord; i++, h = Fl_mul_pre(h, g, p, pi))
2953 56012 : if (a==h) return i;
2954 0 : return ~0UL;
2955 : }
2956 :
2957 : static ulong
2958 0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
2959 : {
2960 0 : pari_sp av = avma;
2961 0 : GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
2962 0 : return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
2963 : }
2964 :
2965 : /* allow pi = 0 */
2966 : ulong
2967 22593 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
2968 : {
2969 22593 : if (!pi) return Fl_log(a, g, ord, p);
2970 22215 : if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
2971 0 : return Fl_log_Fp(a, g, ord, p);
2972 : }
2973 :
2974 : ulong
2975 130283 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
2976 : {
2977 130283 : if (ord <= 200)
2978 0 : return (p&HIGHMASK)? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
2979 130283 : : Fl_log_naive(a, g, ord, p);
2980 0 : return Fl_log_Fp(a, g, ord, p);
2981 : }
2982 :
2983 : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
2984 : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
2985 : static GEN
2986 126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
2987 : {
2988 126 : long l = lg(P) - 1, e = E[l];
2989 126 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
2990 : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
2991 :
2992 126 : if (l == 1) {
2993 98 : hpe = h;
2994 98 : gpe = g;
2995 : } else {
2996 28 : hpe = modii(h, pe);
2997 28 : gpe = modii(g, pe);
2998 : }
2999 126 : if (e == 1) {
3000 42 : hp = hpe;
3001 42 : gp = gpe;
3002 : } else {
3003 84 : hp = remii(hpe, p);
3004 84 : gp = remii(gpe, p);
3005 : }
3006 126 : if (hp == gen_0 || gp == gen_0) return NULL;
3007 105 : if (absequaliu(p, 2))
3008 : {
3009 35 : GEN n = int2n(e);
3010 35 : ogpe = Zp_order(gpe, gen_2, e, n);
3011 35 : a = Fp_log(hpe, gpe, ogpe, n);
3012 35 : if (typ(a) != t_INT) return NULL;
3013 : }
3014 : else
3015 : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
3016 : is trivial */
3017 : /* [order(gp), factor(order(gp))] */
3018 70 : GEN v = Fp_factored_order(gp, subiu(p,1), p);
3019 70 : GEN ogp = gel(v,1);
3020 70 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
3021 70 : a = Fp_log(hp, gp, v, p);
3022 70 : if (typ(a) != t_INT) return NULL;
3023 70 : if (e == 1) ogpe = ogp;
3024 : else
3025 : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
3026 : /* use p-adic log: O(log p + e) mul*/
3027 : long vpogpe, vpohpe;
3028 :
3029 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
3030 28 : gpe = Fp_pow(gpe, ogp, pe);
3031 : /* g,h = 1 mod p; compute b s.t. h = g^b */
3032 :
3033 : /* v_p(order g mod pe) */
3034 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
3035 : /* v_p(order h mod pe) */
3036 28 : vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
3037 28 : if (vpohpe > vpogpe) return NULL;
3038 :
3039 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
3040 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
3041 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
3042 28 : a = addii(a, mulii(ogp, padic_to_Q(b)));
3043 : }
3044 : }
3045 : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
3046 91 : if (l == 1) return a;
3047 :
3048 28 : N = diviiexact(N, pe); /* make N coprime to p */
3049 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
3050 28 : g = Fp_pow(g, modii(ogpe, phi), N);
3051 28 : setlg(P, l); /* remove last element */
3052 28 : setlg(E, l);
3053 28 : b = znlog_rec(h, g, N, P, E, PHI);
3054 28 : if (!b) return NULL;
3055 28 : return addmulii(a, b, ogpe);
3056 : }
3057 :
3058 : static GEN
3059 98 : get_PHI(GEN P, GEN E)
3060 : {
3061 98 : long i, l = lg(P);
3062 98 : GEN PHI = cgetg(l, t_VEC);
3063 98 : gel(PHI,1) = gen_1;
3064 126 : for (i=1; i<l-1; i++)
3065 : {
3066 28 : GEN t, p = gel(P,i);
3067 28 : long e = E[i];
3068 28 : t = mulii(powiu(p, e-1), subiu(p,1));
3069 28 : if (i > 1) t = mulii(t, gel(PHI,i));
3070 28 : gel(PHI,i+1) = t;
3071 : }
3072 98 : return PHI;
3073 : }
3074 :
3075 : GEN
3076 238 : znlog(GEN h, GEN g, GEN o)
3077 : {
3078 238 : pari_sp av = avma;
3079 : GEN N, fa, P, E, x;
3080 238 : switch (typ(g))
3081 : {
3082 28 : case t_PADIC:
3083 : {
3084 28 : GEN p = gel(g,2);
3085 28 : long v = valp(g);
3086 28 : if (v < 0) pari_err_DIM("znlog");
3087 28 : if (v > 0) {
3088 0 : long k = gvaluation(h, p);
3089 0 : if (k % v) return cgetg(1,t_VEC);
3090 0 : k /= v;
3091 0 : if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
3092 0 : return gc_stoi(av, k);
3093 : }
3094 28 : N = gel(g,3);
3095 28 : g = Rg_to_Fp(g, N);
3096 28 : break;
3097 : }
3098 203 : case t_INTMOD:
3099 203 : N = gel(g,1);
3100 203 : g = gel(g,2); break;
3101 7 : default: pari_err_TYPE("znlog", g);
3102 : return NULL; /* LCOV_EXCL_LINE */
3103 : }
3104 231 : if (equali1(N)) { set_avma(av); return gen_0; }
3105 231 : h = Rg_to_Fp(h, N);
3106 224 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
3107 98 : fa = Z_factor(N);
3108 98 : P = gel(fa,1);
3109 98 : E = vec_to_vecsmall(gel(fa,2));
3110 98 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
3111 98 : if (!x) { set_avma(av); return cgetg(1,t_VEC); }
3112 63 : return gerepileuptoint(av, x);
3113 : }
3114 :
3115 : GEN
3116 62832 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
3117 : {
3118 62832 : if (lgefint(p)==3)
3119 : {
3120 62398 : long nn = itos_or_0(n);
3121 62398 : if (nn)
3122 : {
3123 62398 : ulong pp = p[2];
3124 : ulong uz;
3125 62398 : ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
3126 62377 : if (r==ULONG_MAX) return NULL;
3127 62321 : if (zeta) *zeta = utoi(uz);
3128 62321 : return utoi(r);
3129 : }
3130 : }
3131 434 : a = modii(a,p);
3132 434 : if (!signe(a))
3133 : {
3134 0 : if (zeta) *zeta = gen_1;
3135 0 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
3136 0 : return gen_0;
3137 : }
3138 434 : if (absequaliu(n,2))
3139 : {
3140 238 : if (zeta) *zeta = subiu(p,1);
3141 238 : return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
3142 : }
3143 196 : return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
3144 : }
3145 :
3146 : /*********************************************************************/
3147 : /** FACTORIAL **/
3148 : /*********************************************************************/
3149 : GEN
3150 84157 : mulu_interval_step(ulong a, ulong b, ulong step)
3151 : {
3152 84157 : pari_sp av = avma;
3153 : ulong k, l, N, n;
3154 : long lx;
3155 : GEN x;
3156 :
3157 84157 : if (!a) return gen_0;
3158 84157 : if (step == 1) return mulu_interval(a, b);
3159 84157 : n = 1 + (b-a) / step;
3160 84157 : b -= (b-a) % step;
3161 84157 : if (n < 61)
3162 : {
3163 82708 : if (n == 1) return utoipos(a);
3164 63539 : x = muluu(a,a+step); if (n == 2) return x;
3165 511199 : for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
3166 49656 : return gerepileuptoint(av, x);
3167 : }
3168 : /* step | b-a */
3169 1449 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3170 1449 : N = b + a;
3171 1449 : for (k = a;; k += step)
3172 : {
3173 230241 : l = N - k; if (l <= k) break;
3174 228792 : gel(x,lx++) = muluu(k,l);
3175 : }
3176 1449 : if (l == k) gel(x,lx++) = utoipos(k);
3177 1449 : setlg(x, lx);
3178 1449 : return gerepileuptoint(av, ZV_prod(x));
3179 : }
3180 : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
3181 : * first is slower ... ] */
3182 : GEN
3183 156497 : mulu_interval(ulong a, ulong b)
3184 : {
3185 156497 : pari_sp av = avma;
3186 : ulong k, l, N, n;
3187 : long lx;
3188 : GEN x;
3189 :
3190 156497 : if (!a) return gen_0;
3191 156497 : n = b - a + 1;
3192 156497 : if (n < 61)
3193 : {
3194 156497 : if (n == 1) return utoipos(a);
3195 106174 : x = muluu(a,a+1); if (n == 2) return x;
3196 371718 : for (k=a+2; k<b; k++) x = mului(k,x);
3197 : /* avoid k <= b: broken if b = ULONG_MAX */
3198 92132 : return gerepileuptoint(av, mului(b,x));
3199 : }
3200 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3201 0 : N = b + a;
3202 0 : for (k = a;; k++)
3203 : {
3204 0 : l = N - k; if (l <= k) break;
3205 0 : gel(x,lx++) = muluu(k,l);
3206 : }
3207 0 : if (l == k) gel(x,lx++) = utoipos(k);
3208 0 : setlg(x, lx);
3209 0 : return gerepileuptoint(av, ZV_prod(x));
3210 : }
3211 : GEN
3212 588 : muls_interval(long a, long b)
3213 : {
3214 588 : pari_sp av = avma;
3215 588 : long lx, k, l, N, n = b - a + 1;
3216 : GEN x;
3217 :
3218 588 : if (a <= 0 && b >= 0) return gen_0;
3219 315 : if (n < 61)
3220 : {
3221 315 : x = stoi(a);
3222 553 : for (k=a+1; k<=b; k++) x = mulsi(k,x);
3223 315 : return gerepileuptoint(av, x);
3224 : }
3225 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3226 0 : N = b + a;
3227 0 : for (k = a;; k++)
3228 : {
3229 0 : l = N - k; if (l <= k) break;
3230 0 : gel(x,lx++) = mulss(k,l);
3231 : }
3232 0 : if (l == k) gel(x,lx++) = stoi(k);
3233 0 : setlg(x, lx);
3234 0 : return gerepileuptoint(av, ZV_prod(x));
3235 : }
3236 :
3237 : GEN
3238 525376 : mpfact(long n)
3239 : {
3240 525376 : pari_sp av = avma;
3241 : GEN a, v;
3242 : long k;
3243 525376 : if (n <= 12) switch(n)
3244 : {
3245 446373 : case 0: case 1: return gen_1;
3246 36916 : case 2: return gen_2;
3247 3390 : case 3: return utoipos(6);
3248 4233 : case 4: return utoipos(24);
3249 2859 : case 5: return utoipos(120);
3250 2516 : case 6: return utoipos(720);
3251 2410 : case 7: return utoipos(5040);
3252 2411 : case 8: return utoipos(40320);
3253 2432 : case 9: return utoipos(362880);
3254 2684 : case 10:return utoipos(3628800);
3255 1425 : case 11:return utoipos(39916800);
3256 579 : case 12:return utoipos(479001600);
3257 0 : default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
3258 : }
3259 17148 : v = cgetg(expu(n) + 2, t_VEC);
3260 17148 : for (k = 1;; k++)
3261 80363 : {
3262 97511 : long m = n >> (k-1), l;
3263 97511 : if (m <= 2) break;
3264 80363 : l = (1 + (n >> k)) | 1;
3265 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3266 80363 : a = mulu_interval_step(l, m, 2);
3267 80363 : gel(v,k) = k == 1? a: powiu(a, k);
3268 : }
3269 80363 : a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
3270 17148 : a = shifti(a, factorial_lval(n, 2));
3271 17148 : return gerepileuptoint(av, a);
3272 : }
3273 :
3274 : ulong
3275 56752 : factorial_Fl(long n, ulong p)
3276 : {
3277 : long k;
3278 : ulong v;
3279 56752 : if (p <= (ulong)n) return 0;
3280 56752 : v = Fl_powu(2, factorial_lval(n, 2), p);
3281 56796 : for (k = 1;; k++)
3282 142389 : {
3283 199185 : long m = n >> (k-1), l, i;
3284 199185 : ulong a = 1;
3285 199185 : if (m <= 2) break;
3286 142408 : l = (1 + (n >> k)) | 1;
3287 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3288 779701 : for (i=l; i<=m; i+=2)
3289 637306 : a = Fl_mul(a, i, p);
3290 142395 : v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
3291 : }
3292 56777 : return v;
3293 : }
3294 :
3295 : GEN
3296 158 : factorial_Fp(long n, GEN p)
3297 : {
3298 158 : pari_sp av = avma;
3299 : long k;
3300 158 : GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
3301 158 : for (k = 1;; k++)
3302 344 : {
3303 502 : long m = n >> (k-1), l, i;
3304 502 : GEN a = gen_1;
3305 502 : if (m <= 2) break;
3306 344 : l = (1 + (n >> k)) | 1;
3307 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3308 850 : for (i=l; i<=m; i+=2)
3309 506 : a = Fp_mulu(a, i, p);
3310 344 : v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
3311 344 : v = gerepileuptoint(av, v);
3312 : }
3313 158 : return v;
3314 : }
3315 :
3316 : /*******************************************************************/
3317 : /** LUCAS & FIBONACCI **/
3318 : /*******************************************************************/
3319 : static void
3320 56 : lucas(ulong n, GEN *a, GEN *b)
3321 : {
3322 : GEN z, t, zt;
3323 56 : if (!n) { *a = gen_2; *b = gen_1; return; }
3324 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
3325 49 : switch(n & 3) {
3326 14 : case 0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
3327 14 : case 1: *a = subiu(zt,1); *b = addiu(sqri(t),2); break;
3328 7 : case 2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
3329 14 : case 3: *a = addiu(zt,1); *b = subiu(sqri(t),2);
3330 : }
3331 49 : }
3332 :
3333 : GEN
3334 7 : fibo(long n)
3335 : {
3336 7 : pari_sp av = avma;
3337 : GEN a, b;
3338 7 : if (!n) return gen_0;
3339 7 : lucas((ulong)(labs(n)-1), &a, &b);
3340 7 : a = diviuexact(addii(shifti(a,1),b), 5);
3341 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
3342 7 : return gerepileuptoint(av, a);
3343 : }
3344 :
3345 : /*******************************************************************/
3346 : /* CONTINUED FRACTIONS */
3347 : /*******************************************************************/
3348 : static GEN
3349 2830748 : icopy_lg(GEN x, long l)
3350 : {
3351 2830748 : long lx = lgefint(x);
3352 : GEN y;
3353 :
3354 2830748 : if (lx >= l) return icopy(x);
3355 49 : y = cgeti(l); affii(x, y); return y;
3356 : }
3357 :
3358 : /* continued fraction of a/b. If y != NULL, stop when partial quotients
3359 : * differ from y */
3360 : static GEN
3361 2831102 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
3362 : {
3363 : GEN z, c;
3364 2831102 : ulong i, l, ly = lgefint(b);
3365 :
3366 : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
3367 2831102 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
3368 2831102 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
3369 2831102 : if (l > LGBITS) l = LGBITS;
3370 :
3371 2831102 : z = cgetg(l,t_VEC);
3372 2831102 : l--;
3373 2831102 : if (y) {
3374 354 : pari_sp av = avma;
3375 354 : if (l >= (ulong)lg(y)) l = lg(y)-1;
3376 28713 : for (i = 1; i <= l; i++)
3377 : {
3378 28482 : GEN q = gel(y,i);
3379 28482 : gel(z,i) = q;
3380 28482 : c = b; if (!gequal1(q)) c = mulii(q, b);
3381 28482 : c = subii(a, c);
3382 28482 : if (signe(c) < 0)
3383 : { /* partial quotient too large */
3384 117 : c = addii(c, b);
3385 117 : if (signe(c) >= 0) i++; /* by 1 */
3386 117 : break;
3387 : }
3388 28365 : if (cmpii(c, b) >= 0)
3389 : { /* partial quotient too small */
3390 6 : c = subii(c, b);
3391 6 : if (cmpii(c, b) < 0) {
3392 : /* by 1. If next quotient is 1 in y, add 1 */
3393 6 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
3394 6 : i++;
3395 : }
3396 6 : break;
3397 : }
3398 28359 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
3399 28359 : a = b; b = c;
3400 : }
3401 : } else {
3402 2830748 : a = icopy_lg(a, ly);
3403 2830748 : b = icopy(b);
3404 23451289 : for (i = 1; i <= l; i++)
3405 : {
3406 23450959 : gel(z,i) = truedvmdii(a,b,&c);
3407 23450959 : if (c == gen_0) { i++; break; }
3408 20620541 : affii(c, a); cgiv(c); c = a;
3409 20620541 : a = b; b = c;
3410 : }
3411 : }
3412 2831102 : i--;
3413 2831102 : if (i > 1 && gequal1(gel(z,i)))
3414 : {
3415 100 : cgiv(gel(z,i)); --i;
3416 100 : gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
3417 : }
3418 2831102 : setlg(z,i+1); return z;
3419 : }
3420 :
3421 : static GEN
3422 0 : sersfcont(GEN a, GEN b, long k)
3423 : {
3424 0 : long i, l = typ(a) == t_POL? lg(a): 3;
3425 : GEN y, c;
3426 0 : if (lg(b) > l) l = lg(b);
3427 0 : if (k > 0 && l > k+1) l = k+1;
3428 0 : y = cgetg(l,t_VEC);
3429 0 : for (i=1; i<l; i++)
3430 : {
3431 0 : gel(y,i) = poldivrem(a,b,&c);
3432 0 : if (gequal0(c)) { i++; break; }
3433 0 : a = b; b = c;
3434 : }
3435 0 : setlg(y, i); return y;
3436 : }
3437 :
3438 : GEN
3439 2831763 : gboundcf(GEN x, long k)
3440 : {
3441 : pari_sp av;
3442 2831763 : long tx = typ(x), e;
3443 : GEN y, a, b, c;
3444 :
3445 2831763 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
3446 2831756 : if (is_scalar_t(tx))
3447 : {
3448 2831756 : if (gequal0(x)) return mkvec(gen_0);
3449 2831651 : switch(tx)
3450 : {
3451 896 : case t_INT: return mkveccopy(x);
3452 361 : case t_REAL:
3453 361 : av = avma;
3454 361 : c = mantissa_real(x,&e);
3455 361 : if (e < 0) pari_err_PREC("gboundcf");
3456 354 : y = int2n(e);
3457 354 : a = Qsfcont(c,y, NULL, k);
3458 354 : b = addsi(signe(x), c);
3459 354 : return gerepilecopy(av, Qsfcont(b,y, a, k));
3460 :
3461 2830394 : case t_FRAC:
3462 2830394 : av = avma;
3463 2830394 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
3464 : }
3465 0 : pari_err_TYPE("gboundcf",x);
3466 : }
3467 :
3468 0 : switch(tx)
3469 : {
3470 0 : case t_POL: return mkveccopy(x);
3471 0 : case t_SER:
3472 0 : av = avma;
3473 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
3474 0 : case t_RFRAC:
3475 0 : av = avma;
3476 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
3477 : }
3478 0 : pari_err_TYPE("gboundcf",x);
3479 : return NULL; /* LCOV_EXCL_LINE */
3480 : }
3481 :
3482 : static GEN
3483 14 : sfcont2(GEN b, GEN x, long k)
3484 : {
3485 14 : pari_sp av = avma;
3486 14 : long lb = lg(b), tx = typ(x), i;
3487 : GEN y,p1;
3488 :
3489 14 : if (k)
3490 : {
3491 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
3492 0 : lb = k+1;
3493 : }
3494 7 : y = cgetg(lb,t_VEC);
3495 7 : if (lb==1) return y;
3496 7 : if (is_scalar_t(tx))
3497 : {
3498 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
3499 : }
3500 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
3501 :
3502 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
3503 7 : for (i = 1;;)
3504 : {
3505 35 : if (tx == t_REAL)
3506 : {
3507 35 : long e = expo(x);
3508 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
3509 35 : gel(y,i) = floorr(x);
3510 35 : p1 = subri(x, gel(y,i));
3511 : }
3512 : else
3513 : {
3514 0 : gel(y,i) = gfloor(x);
3515 0 : p1 = gsub(x, gel(y,i));
3516 : }
3517 35 : if (++i >= lb) break;
3518 28 : if (gequal0(p1)) break;
3519 28 : x = gdiv(gel(b,i),p1);
3520 : }
3521 7 : setlg(y,i);
3522 7 : return gerepilecopy(av,y);
3523 : }
3524 :
3525 : GEN
3526 112 : gcf(GEN x) { return gboundcf(x,0); }
3527 : GEN
3528 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
3529 : GEN
3530 49 : contfrac0(GEN x, GEN b, long nmax)
3531 : {
3532 : long tb;
3533 :
3534 49 : if (!b) return gboundcf(x,nmax);
3535 28 : tb = typ(b);
3536 28 : if (tb == t_INT) return gboundcf(x,itos(b));
3537 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
3538 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
3539 14 : return sfcont2(b,x,nmax);
3540 : }
3541 :
3542 : GEN
3543 252 : contfracpnqn(GEN x, long n)
3544 : {
3545 252 : pari_sp av = avma;
3546 252 : long i, lx = lg(x);
3547 : GEN M,A,B, p0,p1, q0,q1;
3548 :
3549 252 : if (lx == 1)
3550 : {
3551 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
3552 21 : if (n >= 0) return cgetg(1,t_MAT);
3553 7 : return matid(2);
3554 : }
3555 224 : switch(typ(x))
3556 : {
3557 182 : case t_VEC: case t_COL: A = x; B = NULL; break;
3558 42 : case t_MAT:
3559 42 : switch(lgcols(x))
3560 : {
3561 0 : case 2: A = row(x,1); B = NULL; break;
3562 35 : case 3: A = row(x,2); B = row(x,1); break;
3563 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
3564 : return NULL; /*LCOV_EXCL_LINE*/
3565 : }
3566 35 : break;
3567 0 : default: pari_err_TYPE("pnqn",x);
3568 : return NULL; /*LCOV_EXCL_LINE*/
3569 : }
3570 217 : p1 = gel(A,1);
3571 217 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
3572 217 : if (n >= 0)
3573 : {
3574 182 : lx = minss(lx, n+2);
3575 182 : if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
3576 : }
3577 35 : else if (lx == 2)
3578 7 : return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
3579 : /* lx >= 3 */
3580 119 : p0 = gen_1;
3581 119 : q0 = gen_0; /* p[-1], q[-1] */
3582 119 : M = cgetg(lx, t_MAT);
3583 119 : gel(M,1) = mkcol2(p1,q1);
3584 399 : for (i=2; i<lx; i++)
3585 : {
3586 280 : GEN a = gel(A,i), p2,q2;
3587 280 : if (B) {
3588 84 : GEN b = gel(B,i);
3589 84 : p0 = gmul(b,p0);
3590 84 : q0 = gmul(b,q0);
3591 : }
3592 280 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
3593 280 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
3594 280 : gel(M,i) = mkcol2(p1,q1);
3595 : }
3596 119 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
3597 119 : return gerepilecopy(av, M);
3598 : }
3599 : GEN
3600 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
3601 : /* x = [a0, ..., an] from gboundcf, n >= 0;
3602 : * return [[p0, ..., pn], [q0,...,qn]] */
3603 : GEN
3604 609308 : ZV_allpnqn(GEN x)
3605 : {
3606 609308 : long i, lx = lg(x);
3607 609308 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
3608 :
3609 609308 : gel(v,1) = P = cgetg(lx, t_VEC);
3610 609308 : gel(v,2) = Q = cgetg(lx, t_VEC);
3611 609308 : p0 = gen_1; q0 = gen_0;
3612 609308 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
3613 2092209 : for (i=2; i<lx; i++)
3614 : {
3615 1482901 : GEN a = gel(x,i);
3616 1482901 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
3617 1482901 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
3618 : }
3619 609308 : return v;
3620 : }
3621 :
3622 : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
3623 : static GEN
3624 42 : mod_to_frac(GEN x, GEN N, GEN B)
3625 : {
3626 : GEN a, b, A;
3627 42 : if (B) A = divii(shifti(N, -1), B);
3628 : else
3629 : {
3630 14 : A = sqrti(shifti(N, -1));
3631 14 : B = A;
3632 : }
3633 42 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
3634 28 : return equali1(b)? a: mkfrac(a,b);
3635 : }
3636 :
3637 : static GEN
3638 70 : mod_to_rfrac(GEN x, GEN N, long B)
3639 : {
3640 : GEN a, b;
3641 70 : long A, d = degpol(N);
3642 70 : if (B >= 0) A = d-1 - B;
3643 : else
3644 : {
3645 42 : B = d >> 1;
3646 42 : A = odd(d)? B : B-1;
3647 : }
3648 70 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
3649 70 : if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
3650 56 : return gdiv(a,b);
3651 : }
3652 :
3653 : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
3654 : * of the continued fraction of x with b <= k maximal */
3655 : static GEN
3656 7 : bestappr_frac(GEN x, GEN k)
3657 : {
3658 : pari_sp av;
3659 : GEN p0, p1, p, q0, q1, q, a, y;
3660 :
3661 7 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
3662 0 : av = avma; y = x;
3663 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
3664 0 : q1 = gen_0; q0 = gen_1;
3665 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
3666 : for(;;)
3667 : {
3668 0 : x = ginv(x); /* > 1 */
3669 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
3670 0 : if (cmpii(a,k) > 0)
3671 : { /* next partial quotient will overflow limits */
3672 : GEN n, d;
3673 0 : a = divii(subii(k, q1), q0);
3674 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3675 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3676 : /* compare |y-p0/q0|, |y-p1/q1| */
3677 0 : n = gel(y,1);
3678 0 : d = gel(y,2);
3679 0 : if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
3680 : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
3681 0 : { p1 = p0; q1 = q0; }
3682 0 : break;
3683 : }
3684 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3685 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3686 :
3687 0 : if (cmpii(q0,k) > 0) break;
3688 0 : x = gsub(x,a); /* 0 <= x < 1 */
3689 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
3690 :
3691 : }
3692 0 : return gerepileupto(av, gdiv(p1,q1));
3693 : }
3694 : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
3695 : * of the continued fraction of x with b <= k maximal */
3696 : static GEN
3697 1522524 : bestappr_real(GEN x, GEN k)
3698 : {
3699 1522524 : pari_sp av = avma;
3700 1522524 : GEN kr, p0, p1, p, q0, q1, q, a, y = x;
3701 :
3702 1522524 : p1 = gen_1; a = p0 = floorr(x);
3703 1522429 : q1 = gen_0; q0 = gen_1;
3704 1522429 : x = subri(x,a); /* 0 <= x < 1 */
3705 1522431 : if (!signe(x)) { cgiv(x); return a; }
3706 1416551 : kr = itor(k, realprec(x));
3707 : for(;;)
3708 1687472 : {
3709 : long d;
3710 3104040 : x = invr(x); /* > 1 */
3711 3103933 : if (cmprr(x,kr) > 0)
3712 : { /* next partial quotient will overflow limits */
3713 1403615 : a = divii(subii(k, q1), q0);
3714 1403551 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3715 1403601 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3716 : /* compare |y-p0/q0|, |y-p1/q1| */
3717 1403555 : if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
3718 : mulir(q0, subri(mulir(q1,y), p1))) < 0)
3719 124346 : { p1 = p0; q1 = q0; }
3720 1403641 : break;
3721 : }
3722 1700406 : d = nbits2prec(expo(x) + 1);
3723 1700408 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
3724 :
3725 1700212 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
3726 1700174 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3727 1700205 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3728 :
3729 1700192 : if (cmpii(q0,k) > 0) break;
3730 1694176 : x = subri(x,a); /* 0 <= x < 1 */
3731 1694169 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
3732 : }
3733 1416555 : if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
3734 1416555 : return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
3735 : }
3736 :
3737 : /* k t_INT or NULL */
3738 : static GEN
3739 2510959 : bestappr_Q(GEN x, GEN k)
3740 : {
3741 2510959 : long lx, tx = typ(x), i;
3742 : GEN a, y;
3743 :
3744 2510959 : switch(tx)
3745 : {
3746 136 : case t_INT: return icopy(x);
3747 7 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
3748 1773762 : case t_REAL:
3749 1773762 : if (!signe(x)) return gen_0;
3750 : /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
3751 1522502 : i = bit_prec(x); if (i <= expo(x)) return NULL;
3752 1522522 : return bestappr_real(x, k? k: int2n(i));
3753 :
3754 28 : case t_INTMOD: {
3755 28 : pari_sp av = avma;
3756 28 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
3757 21 : return gerepilecopy(av, a);
3758 : }
3759 14 : case t_PADIC: {
3760 14 : pari_sp av = avma;
3761 14 : long v = valp(x);
3762 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
3763 7 : if (v) a = gmul(a, powis(gel(x,2), v));
3764 7 : return gerepilecopy(av, a);
3765 : }
3766 :
3767 312 : case t_COMPLEX: {
3768 312 : pari_sp av = avma;
3769 312 : y = cgetg(3, t_COMPLEX);
3770 312 : gel(y,2) = bestappr(gel(x,2), k);
3771 312 : gel(y,1) = bestappr(gel(x,1), k);
3772 312 : if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
3773 91 : return y;
3774 : }
3775 0 : case t_SER:
3776 0 : if (ser_isexactzero(x)) return gcopy(x);
3777 : /* fall through */
3778 : case t_POLMOD: case t_POL: case t_RFRAC:
3779 : case t_VEC: case t_COL: case t_MAT:
3780 736700 : y = cgetg_copy(x, &lx);
3781 736745 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
3782 3180968 : for (; i<lx; i++)
3783 : {
3784 2444238 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
3785 2444223 : gel(y,i) = a;
3786 : }
3787 736730 : if (tx == t_POL) return normalizepol(y);
3788 736716 : if (tx == t_SER) return normalizeser(y);
3789 736716 : return y;
3790 : }
3791 0 : pari_err_TYPE("bestappr_Q",x);
3792 : return NULL; /* LCOV_EXCL_LINE */
3793 : }
3794 :
3795 : static GEN
3796 56 : bestappr_ser(GEN x, long B)
3797 : {
3798 56 : long dN, v = valser(x), lx = lg(x);
3799 : GEN t;
3800 56 : x = normalizepol(ser2pol_i(x, lx));
3801 56 : dN = lx-2;
3802 56 : if (v > 0)
3803 : {
3804 14 : x = RgX_shift_shallow(x, v);
3805 14 : dN += v;
3806 : }
3807 42 : else if (v < 0)
3808 : {
3809 7 : if (B >= 0) B = maxss(B+v, 0);
3810 : }
3811 56 : t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
3812 56 : if (!t) return NULL;
3813 42 : if (v < 0)
3814 : {
3815 : GEN a, b;
3816 : long vx;
3817 7 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
3818 : /* t_RFRAC */
3819 7 : vx = varn(x);
3820 7 : a = gel(t,1);
3821 7 : b = gel(t,2);
3822 7 : v -= RgX_valrem(b, &b);
3823 7 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
3824 7 : if (v < 0) b = RgX_shift(b, -v);
3825 0 : else if (v > 0) {
3826 0 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
3827 0 : a = RgX_shift(a, v);
3828 : }
3829 7 : t = mkrfraccopy(a, b);
3830 : }
3831 42 : return t;
3832 : }
3833 : static GEN bestappr_RgX(GEN x, long B);
3834 : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
3835 : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
3836 : static GEN
3837 77 : bestappr_RgX(GEN x, long B)
3838 : {
3839 77 : long i, lx, tx = typ(x);
3840 : GEN y, t;
3841 77 : switch(tx)
3842 : {
3843 0 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
3844 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
3845 0 : return gcopy(x);
3846 :
3847 14 : case t_RFRAC: {
3848 14 : pari_sp av = avma;
3849 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
3850 7 : x = rfrac_to_ser_i(x, 2*B+1);
3851 7 : t = bestappr_ser(x, B); if (!t) return NULL;
3852 0 : return gerepileupto(av, t);
3853 : }
3854 14 : case t_POLMOD: {
3855 14 : pari_sp av = avma;
3856 14 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
3857 14 : return gerepileupto(av, t);
3858 : }
3859 49 : case t_SER: {
3860 49 : pari_sp av = avma;
3861 49 : t = bestappr_ser(x, B); if (!t) return NULL;
3862 42 : return gerepileupto(av, t);
3863 : }
3864 :
3865 0 : case t_VEC: case t_COL: case t_MAT:
3866 0 : y = cgetg_copy(x, &lx);
3867 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
3868 0 : for (; i<lx; i++)
3869 : {
3870 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
3871 0 : gel(y,i) = t;
3872 : }
3873 0 : return y;
3874 : }
3875 0 : pari_err_TYPE("bestappr_RgX",x);
3876 : return NULL; /* LCOV_EXCL_LINE */
3877 : }
3878 :
3879 : /* allow k = NULL: maximal accuracy */
3880 : GEN
3881 66711 : bestappr(GEN x, GEN k)
3882 : {
3883 66711 : pari_sp av = avma;
3884 66711 : if (k) { /* replace by floor(k) */
3885 66389 : switch(typ(k))
3886 : {
3887 4404 : case t_INT:
3888 4404 : break;
3889 61985 : case t_REAL: case t_FRAC:
3890 61985 : k = floor_safe(k); /* left on stack for efficiency */
3891 61986 : if (!signe(k)) k = gen_1;
3892 61986 : break;
3893 0 : default:
3894 0 : pari_err_TYPE("bestappr [bound type]", k);
3895 0 : break;
3896 : }
3897 322 : }
3898 66712 : x = bestappr_Q(x, k);
3899 66712 : if (!x) { set_avma(av); return cgetg(1,t_VEC); }
3900 66697 : return x;
3901 : }
3902 : GEN
3903 77 : bestapprPade(GEN x, long B)
3904 : {
3905 77 : pari_sp av = avma;
3906 77 : GEN t = bestappr_RgX(x, B);
3907 77 : if (!t) { set_avma(av); return cgetg(1,t_VEC); }
3908 63 : return t;
3909 : }
|
