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 465056 : 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 465056 : 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 532842 : u_is_gener_expo(ulong p, GEN L0)
54 : {
55 532842 : const ulong q = p >> 1;
56 : long i;
57 : GEN L;
58 532842 : if (!L0) L0 = u_odd_prime_divisors(q);
59 532842 : L = cgetg_copy(L0,&i);
60 1155274 : while (--i) L[i] = q / uel(L0,i);
61 532842 : return L;
62 : }
63 :
64 : int
65 1629724 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
66 : {
67 : long i;
68 1629724 : if (krouu(x, p) >= 0) return 0;
69 1380256 : for (i=lg(L)-1; i; i--)
70 : {
71 842625 : ulong t = Fl_powu(x, uel(L,i), p);
72 842625 : if (t == p_1 || t == 1) return 0;
73 : }
74 537631 : return 1;
75 : }
76 : /* assume p prime */
77 : ulong
78 1054306 : pgener_Fl_local(ulong p, GEN L0)
79 : {
80 1054306 : const pari_sp av = avma;
81 1054306 : const ulong p_1 = p-1;
82 : long x;
83 : GEN L;
84 1054306 : if (p <= 19) switch(p)
85 : { /* quick trivial cases */
86 28 : case 2: return 1;
87 116168 : case 7:
88 116168 : case 17: return 3;
89 405321 : default: return 2;
90 : }
91 532789 : L = u_is_gener_expo(p,L0);
92 1622042 : for (x = 2;; x++)
93 1622042 : if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
94 : }
95 : ulong
96 579423 : 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 13878 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
102 : {
103 13878 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
104 13878 : if (t >= 0) return 0;
105 21913 : for (i = lg(L)-1; i; i--)
106 : {
107 14396 : GEN t = Fp_pow(x, gel(L,i), p);
108 14396 : if (equalii(t, p_1) || equali1(t)) return 0;
109 : }
110 7517 : return 1;
111 : }
112 :
113 : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
114 : GEN
115 357748 : pgener_Fp_local(GEN p, GEN L0)
116 : {
117 357748 : pari_sp av0 = avma;
118 : GEN x, p_1, L;
119 357748 : if (lgefint(p) == 3)
120 : {
121 : ulong z;
122 352858 : if (p[2] == 2) return gen_1;
123 257939 : if (L0) L0 = ZV_to_nv(L0);
124 257944 : z = pgener_Fl_local(uel(p,2), L0);
125 257981 : return gc_utoipos(av0, z);
126 : }
127 4890 : 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 44114 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
135 :
136 : ulong
137 203316 : pgener_Zl(ulong p)
138 : {
139 203316 : 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 203316 : if (p == 40487) return 10;
142 : #ifndef LONG_IS_64BIT
143 29466 : return pgener_Fl(p);
144 : #else
145 173850 : 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 168567 : pgener_Zp(GEN p)
162 : {
163 168567 : 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 214770 : rootsof1_Fl(ulong n, ulong p)
251 : {
252 214770 : pari_sp av = avma;
253 214770 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
254 214768 : ulong z = pgener_Fl_local(p, L);
255 214768 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
256 214768 : 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 317849826 : ome(long t)
364 : {
365 317849826 : switch(t & 7)
366 : {
367 180155326 : case 3:
368 180155326 : case 5: return 1;
369 137694500 : default: return 0;
370 : }
371 : }
372 : /* t a t_INT, is t = 3,5 mod 8 ? */
373 : static int
374 5593818 : gome(GEN t)
375 5593818 : { return signe(t)? ome( mod2BIL(t) ): 0; }
376 :
377 : /* assume y odd, return kronecker(x,y) * s */
378 : static long
379 224140351 : krouu_s(ulong x, ulong y, long s)
380 : {
381 224140351 : ulong x1 = x, y1 = y, z;
382 1018915201 : while (x1)
383 : {
384 794765391 : long r = vals(x1);
385 794879253 : if (r)
386 : {
387 423165240 : if (odd(r) && ome(y1)) s = -s;
388 423060837 : x1 >>= r;
389 : }
390 794774850 : if (x1 & y1 & 2) s = -s;
391 794774850 : z = y1 % x1; y1 = x1; x1 = z;
392 : }
393 224149810 : return (y1 == 1)? s: 0;
394 : }
395 :
396 : long
397 11963856 : kronecker(GEN x, GEN y)
398 : {
399 11963856 : pari_sp av = avma;
400 11963856 : long s = 1, r;
401 : ulong xu;
402 :
403 11963856 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
404 11963856 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
405 11963856 : switch (signe(y))
406 : {
407 63 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
408 133 : case 0: return is_pm1(x);
409 : }
410 11963723 : r = vali(y);
411 11963711 : if (r)
412 : {
413 1348898 : if (!mpodd(x)) return gc_long(av,0);
414 321711 : if (odd(r) && gome(x)) s = -s;
415 321711 : y = shifti(y,-r);
416 : }
417 10936524 : x = modii(x,y);
418 13320310 : while (lgefint(x) > 3) /* x < y */
419 : {
420 : GEN z;
421 2383818 : r = vali(x);
422 2383892 : if (r)
423 : {
424 1302371 : if (odd(r) && gome(y)) s = -s;
425 1302372 : x = shifti(x,-r);
426 : }
427 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
428 2382524 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
429 2382791 : z = remii(y,x); y = x; x = z;
430 2383811 : 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 10936492 : xu = itou(x);
437 10936488 : if (!xu) return is_pm1(y)? s: 0;
438 10838186 : r = vals(xu);
439 10838184 : if (r)
440 : {
441 5755974 : if (odd(r) && gome(y)) s = -s;
442 5755974 : xu >>= r;
443 : }
444 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
445 10838184 : if (xu & mod2BIL(y) & 2) s = -s;
446 10838186 : return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
447 : }
448 :
449 : long
450 31542 : krois(GEN x, long y)
451 : {
452 : ulong yu;
453 31542 : long s = 1;
454 :
455 31542 : if (y <= 0)
456 : {
457 28 : if (y == 0) return is_pm1(x);
458 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
459 : }
460 : else
461 31514 : yu = (ulong)y;
462 31514 : if (!odd(yu))
463 : {
464 : long r;
465 15008 : if (!mpodd(x)) return 0;
466 11158 : r = vals(yu); yu >>= r;
467 11158 : if (odd(r) && gome(x)) s = -s;
468 : }
469 27664 : return krouu_s(umodiu(x, yu), yu, s);
470 : }
471 : /* assume y != 0 */
472 : long
473 27332590 : kroiu(GEN x, ulong y)
474 : {
475 : long r;
476 27332590 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
477 299061 : if (!mpodd(x)) return 0;
478 204499 : r = vals(y); y >>= r;
479 204502 : 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 179183 : krouodd(ulong x, GEN y, long s)
485 : {
486 : long r;
487 179183 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
488 29241 : if (!x) return 0; /* y != 1 */
489 29241 : r = vals(x);
490 29241 : if (r)
491 : {
492 15230 : if (odd(r) && gome(y)) s = -s;
493 15230 : x >>= r;
494 : }
495 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
496 29241 : if (x & mod2BIL(y) & 2) s = -s;
497 29241 : return krouu_s(umodiu(y,x), x, s);
498 : }
499 :
500 : long
501 144257 : krosi(long x, GEN y)
502 : {
503 144257 : const pari_sp av = avma;
504 144257 : long s = 1, r;
505 144257 : 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 144257 : r = vali(y);
511 144257 : 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 144257 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
518 144257 : return gc_long(av, krouodd((ulong)x, y, s));
519 : }
520 :
521 : long
522 34926 : kroui(ulong x, GEN y)
523 : {
524 34926 : const pari_sp av = avma;
525 34926 : long s = 1, r;
526 34926 : switch (signe(y))
527 : {
528 0 : case -1: y = negi(y); break;
529 0 : case 0: return x==1UL;
530 : }
531 34926 : r = vali(y);
532 34926 : 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 34926 : return gc_long(av, krouodd(x, y, s));
539 : }
540 :
541 : long
542 96276415 : kross(long x, long y)
543 : {
544 : ulong yu;
545 96276415 : long s = 1;
546 :
547 96276415 : 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 96207472 : yu = (ulong)y;
554 96276387 : if (!odd(yu))
555 : {
556 : long r;
557 23414183 : if (!odd(x)) return 0;
558 16495643 : r = vals(yu); yu >>= r;
559 16495643 : if (odd(r) && ome(x)) s = -s;
560 : }
561 89357847 : x %= (long)yu; if (x < 0) x += yu;
562 89357847 : return krouu_s((ulong)x, yu, s);
563 : }
564 :
565 : long
566 96528223 : krouu(ulong x, ulong y)
567 : {
568 : long r;
569 96528223 : if (odd(y)) return krouu_s(x, y, 1);
570 5055 : if (!odd(x)) return 0;
571 5632 : r = vals(y); y >>= r;
572 5632 : return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
573 : }
574 :
575 : /*********************************************************************/
576 : /** HILBERT SYMBOL **/
577 : /*********************************************************************/
578 : /* x,y are t_INT or t_REAL */
579 : static long
580 7329 : mphilbertoo(GEN x, GEN y)
581 : {
582 7329 : long sx = signe(x), sy = signe(y);
583 7329 : if (!sx || !sy) return 0;
584 7329 : return (sx < 0 && sy < 0)? -1: 1;
585 : }
586 :
587 : long
588 135331 : hilbertii(GEN x, GEN y, GEN p)
589 : {
590 : pari_sp av;
591 : long oddvx, oddvy, z;
592 :
593 135331 : if (!p) return mphilbertoo(x,y);
594 128023 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
595 128023 : if (!signe(x) || !signe(y)) return 0;
596 128002 : av = avma;
597 128002 : oddvx = odd(Z_pvalrem(x,p,&x));
598 128002 : oddvy = odd(Z_pvalrem(y,p,&y));
599 : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
600 128002 : if (absequaliu(p, 2))
601 : {
602 12103 : z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
603 12103 : if (oddvx && gome(y)) z = -z;
604 12103 : if (oddvy && gome(x)) z = -z;
605 : }
606 : else
607 : {
608 115899 : z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
609 115899 : if (oddvx && kronecker(y,p) < 0) z = -z;
610 115899 : if (oddvy && kronecker(x,p) < 0) z = -z;
611 : }
612 128002 : 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 60403 : hilbert(GEN x, GEN y, GEN p)
664 : {
665 60403 : pari_sp av = avma;
666 60403 : long tx = typ(x), ty = typ(y);
667 :
668 60403 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
669 60403 : 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 60326 : 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 60312 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
692 60109 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
693 :
694 60053 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
695 59990 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
696 :
697 59913 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
698 59913 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
699 :
700 59913 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
701 59913 : if (p && !signe(p)) p = NULL;
702 59913 : return gc_long(av, hilbertii(x,y,p));
703 : }
704 :
705 : /*******************************************************************/
706 : /* SQUARE ROOT MODULO p */
707 : /*******************************************************************/
708 : static void
709 2087861 : checkp(ulong q, ulong p)
710 2087861 : { 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 11228975 : nonsquare1_Fl(ulong p)
714 : {
715 : forprime_t S;
716 : ulong q;
717 11228975 : if ((p & 7UL) != 1) return 2UL;
718 4172179 : q = p % 3; if (q == 2) return 3UL;
719 1309419 : checkp(q, p);
720 1314996 : q = p % 5; if (q == 2 || q == 3) return 5UL;
721 482165 : checkp(q, p);
722 482159 : q = p % 7; if (q != 4 && q >= 3) return 7UL;
723 178087 : checkp(q, p);
724 : /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
725 178121 : u_forprime_init(&S, 11, 1967);
726 290904 : while ((q = u_forprime_next(&S)))
727 : {
728 290902 : if (krouu(q, p) < 0) return q;
729 112778 : checkp(q, p);
730 : }
731 0 : checkp(0, p);
732 : return 0; /*LCOV_EXCL_LINE*/
733 : }
734 : /* p > 2 a prime */
735 : ulong
736 7895 : nonsquare_Fl(ulong p)
737 7895 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
738 :
739 : /* allow pi = 0 */
740 : ulong
741 174870 : Fl_2gener_pre(ulong p, ulong pi)
742 : {
743 174870 : ulong p1 = p-1;
744 174870 : long e = vals(p1);
745 174859 : if (e == 1) return p1;
746 65051 : return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
747 : }
748 :
749 : ulong
750 67856 : Fl_2gener_pre_i(ulong ns, ulong p, ulong pi)
751 : {
752 67856 : ulong p1 = p-1;
753 67856 : long e = vals(p1);
754 67856 : if (e == 1) return p1;
755 26020 : return Fl_powu_pre(ns, p1 >> e, p, pi);
756 : }
757 :
758 : static ulong
759 11742712 : Fl_sqrt_i(ulong a, ulong y, ulong p)
760 : {
761 : long i, e, k;
762 : ulong p1, q, v, w;
763 :
764 11742712 : if (!a) return 0;
765 10484853 : p1 = p - 1; e = vals(p1);
766 10484825 : if (e == 0) /* p = 2 */
767 : {
768 539339 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
769 540614 : return ((a & 1) == 0)? 0: 1;
770 : }
771 9945486 : if (e == 1)
772 : {
773 4800722 : v = Fl_powu(a, (p+1) >> 2, p);
774 4800942 : if (Fl_sqr(v, p) != a) return ~0UL;
775 4796068 : p1 = p - v; if (v > p1) v = p1;
776 4796068 : return v;
777 : }
778 5144764 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
779 5144764 : p1 = Fl_powu(a, q >> 1, p); /* a ^ [(q-1)/2] */
780 5144867 : if (!p1) return 0;
781 5144867 : v = Fl_mul(a, p1, p);
782 5144850 : w = Fl_mul(v, p1, p);
783 5144882 : if (!y) y = Fl_powu(nonsquare1_Fl(p), q, p);
784 8782185 : 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 3639126 : p1 = Fl_sqr(w, p);
788 5975066 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr(p1, p);
789 3639104 : if (k == e) return ~0UL;
790 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
791 3637027 : p1 = y;
792 4774392 : for (i=1; i < e-k; i++) p1 = Fl_sqr(p1, p);
793 3637027 : y = Fl_sqr(p1, p); e = k;
794 3637085 : w = Fl_mul(y, w, p);
795 3637115 : v = Fl_mul(v, p1, p);
796 : }
797 5143059 : p1 = p - v; if (v > p1) v = p1;
798 5143059 : return v;
799 : }
800 :
801 : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. Allow pi = 0 */
802 : ulong
803 32902430 : 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 32902430 : if (!pi) return Fl_sqrt_i(a, y, p);
809 21159891 : if (!a) return 0;
810 21037034 : p1 = p - 1; e = vals(p1);
811 21041883 : 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 21053501 : if (e == 1)
817 : {
818 14983451 : v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
819 14962865 : if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
820 14966006 : p1 = p - v; if (v > p1) v = p1;
821 14966006 : return v;
822 : }
823 6070050 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
824 6070050 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
825 6067377 : if (!p1) return 0;
826 6067377 : v = Fl_mul_pre(a, p1, p, pi);
827 6068767 : w = Fl_mul_pre(v, p1, p, pi);
828 6068857 : if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
829 11555935 : 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 5488777 : p1 = Fl_sqr_pre(w,p,pi);
833 10275005 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
834 5489816 : if (k == e) return ~0UL;
835 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
836 5489724 : p1 = y;
837 7195860 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
838 5489780 : y = Fl_sqr_pre(p1, p, pi); e = k;
839 5490409 : w = Fl_mul_pre(y, w, p, pi);
840 5490609 : v = Fl_mul_pre(v, p1, p, pi);
841 : }
842 6067158 : p1 = p - v; if (v > p1) v = p1;
843 6067158 : return v;
844 : }
845 :
846 : ulong
847 11757550 : Fl_sqrt(ulong a, ulong p)
848 11757550 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0; return Fl_sqrt_pre_i(a, 0, p, pi); }
849 :
850 : ulong
851 20986685 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
852 20986685 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
853 :
854 : /* allow pi = 0 */
855 : static ulong
856 140557 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
857 : {
858 140557 : ulong x, y, m, le1 = upowuu(l, e-1);
859 140558 : for (x = 2; ; x++)
860 : {
861 171355 : y = Fl_powu_pre(x, r, p, pi);
862 171354 : if (y==1) continue;
863 154392 : m = Fl_powu_pre(y, le1, p, pi);
864 154392 : if (m != 1) break;
865 : }
866 140557 : *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 223784 : 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 223784 : if (a==0) return a;
879 223778 : u2 = Fl_inv(l%r, r);
880 223779 : v = Fl_powu_pre(a, u2, p, pi);
881 223776 : w = Fl_powu_pre(v, l, p, pi);
882 223778 : w = pi? Fl_mul_pre(w, Fl_inv(a, p), p, pi): Fl_div(w, a, p);
883 223765 : if (w==1) return v;
884 138476 : if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
885 163147 : while (w!=1)
886 : {
887 143725 : ulong k = 0;
888 143725 : p1 = w;
889 : do
890 : {
891 187121 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
892 187122 : if (++k == e) return ULONG_MAX;
893 68067 : } while (p1!=1);
894 24671 : dl = Fl_log_pre(z, m, l, p, pi);
895 24671 : dl = Fl_neg(dl, l);
896 24671 : p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
897 24671 : m = Fl_powu_pre(m, dl, p, pi);
898 24671 : e = k;
899 24671 : v = pi? Fl_mul_pre(p1,v,p,pi): Fl_mul(p1,v,p);
900 24671 : y = Fl_powu_pre(p1,l,p,pi);
901 24671 : w = pi? Fl_mul_pre(y,w,p,pi): Fl_mul(y,w,p);
902 : }
903 19422 : return v;
904 : }
905 :
906 : /* allow pi = 0 */
907 : static ulong
908 221417 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
909 : {
910 221417 : ulong r, e = u_lvalrem(p-1, l, &r);
911 221417 : return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
912 : }
913 : /* allow pi = 0 */
914 : ulong
915 221415 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
916 221415 : { 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 142100 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
926 : {
927 142100 : ulong m, q = p-1, z;
928 142100 : ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
929 142100 : if (a==0)
930 : {
931 82467 : if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
932 82460 : if (zetan) *zetan = 1UL;
933 82460 : return 0;
934 : }
935 59633 : if (n==1)
936 : {
937 420 : if (zetan) *zetan = 1;
938 420 : return n < 0? Fl_inv(a,p): a;
939 : }
940 59213 : if (n==2)
941 : {
942 6482 : if (zetan) *zetan = p-1;
943 6482 : return Fl_sqrt_pre_i(a, 0, p, pi);
944 : }
945 52731 : if (a == 1 && !zetan) return a;
946 26219 : m = ugcd(nn, q);
947 26219 : z = 1;
948 26219 : if (m!=1)
949 : {
950 2076 : GEN F = factoru(m);
951 : long i, j, e;
952 : ulong r, zeta, y, l;
953 4535 : for (i = nbrows(F); i; i--)
954 : {
955 2522 : l = ucoeff(F,i,1);
956 2522 : j = ucoeff(F,i,2);
957 2522 : e = u_lvalrem(q,l, &r);
958 2522 : y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
959 2522 : if (zetan)
960 : {
961 770 : ulong Y = Fl_powu_pre(y, upowuu(l,e-j), p, pi);
962 770 : z = pi? Fl_mul_pre(z, Y, p, pi): Fl_mul(z, Y, p);
963 : }
964 2522 : if (a!=1)
965 : do
966 : {
967 2367 : a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
968 2353 : if (a==ULONG_MAX) return ULONG_MAX;
969 2304 : } while (--j);
970 : }
971 : }
972 26156 : if (m != nn)
973 : {
974 24164 : ulong qm = q/m, nm = (nn/m) % qm;
975 24164 : a = Fl_powu_pre(a, Fl_inv(nm, qm), p, pi);
976 : }
977 26156 : if (n < 0) a = Fl_inv(a, p);
978 26156 : if (zetan) *zetan = z;
979 26156 : return a;
980 : }
981 :
982 : ulong
983 142100 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
984 : {
985 142100 : ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
986 142100 : 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 0 : sqrt_Cipolla_sqr(void *data, GEN y)
1002 : {
1003 0 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
1004 0 : GEN u2 = sqri(u), v2 = sqri(v);
1005 0 : v = subii(sqri(addii(v,u)), addii(u2,v2));
1006 0 : u = addii(u2, mulii(v2,n));
1007 0 : retmkvec2(modii(u,p), modii(v,p));
1008 : }
1009 : /* compute (t+X) y^2 */
1010 : static GEN
1011 0 : sqrt_Cipolla_msqr(void *data, GEN y)
1012 : {
1013 0 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2);
1014 0 : ulong t = gel(data,4)[2];
1015 0 : GEN d = addii(u, mului(t,v)), d2 = sqri(d);
1016 0 : GEN b = remii(mulii(a,v), p);
1017 0 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
1018 0 : v = subii(d2, mulii(b,v));
1019 0 : retmkvec2(modii(u,p), modii(v,p));
1020 : }
1021 : /* assume a reduced mod p [ otherwise correct but inefficient ] */
1022 : static GEN
1023 0 : sqrt_Cipolla(GEN a, GEN p)
1024 : {
1025 : pari_sp av;
1026 : GEN u, n, y, pov2;
1027 : ulong t;
1028 :
1029 0 : if (kronecker(a, p) < 0) return NULL;
1030 0 : pov2 = shifti(p,-1); /* center to avoid multiplying by huge base*/
1031 0 : if (cmpii(a,pov2) > 0) a = subii(a,p);
1032 0 : av = avma;
1033 0 : for (t=1; ; t++, set_avma(av))
1034 : {
1035 0 : n = subsi((long)(t*t), a);
1036 0 : if (kronecker(n, p) < 0) break;
1037 : }
1038 :
1039 : /* compute (t+X)^((p-1)/2) =: u+vX */
1040 0 : u = utoipos(t);
1041 0 : 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 0 : return Fp_mul(gel(y,2), a, p);
1050 : }
1051 :
1052 : /* Return NULL if p is found to be composite.
1053 : * p odd, q = (p-1)/2^oo is odd */
1054 : static GEN
1055 5904 : Fp_2gener_all(GEN q, GEN p)
1056 : {
1057 : long k;
1058 5904 : for (k = 2;; k++)
1059 11607 : {
1060 17511 : long i = kroui(k, p);
1061 17511 : if (i < 0) return Fp_pow(utoipos(k), q, p);
1062 11607 : if (i == 0) return NULL;
1063 : }
1064 : }
1065 :
1066 : /* Return NULL if p is found to be composite */
1067 : GEN
1068 3192 : Fp_2gener(GEN p)
1069 : {
1070 3192 : GEN q = subiu(p, 1);
1071 3192 : long e = Z_lvalrem(q, 2, &q);
1072 3192 : if (e == 0 && !equaliu(p,2)) return NULL;
1073 3192 : return Fp_2gener_all(q, p);
1074 : }
1075 :
1076 : GEN
1077 19931 : Fp_2gener_i(GEN ns, GEN p)
1078 : {
1079 19931 : GEN q = subiu(p,1);
1080 19931 : long e = vali(q);
1081 19931 : if (e == 1) return q;
1082 18546 : return Fp_pow(ns, shifti(q,-e), p);
1083 : }
1084 :
1085 : static GEN
1086 1465 : nonsquare_Fp(GEN p)
1087 : {
1088 : forprime_t T;
1089 : ulong a;
1090 1465 : if (mod4(p)==3) return gen_m1;
1091 1465 : if (mod8(p)==5) return gen_2;
1092 714 : u_forprime_init(&T, 3, ULONG_MAX);
1093 1375 : while((a = u_forprime_next(&T)))
1094 1375 : if (kroui(a,p) < 0) return utoi(a);
1095 0 : pari_err_PRIME("Fp_sqrt [modulus]",p);
1096 : return NULL; /* LCOV_EXCL_LINE */
1097 : }
1098 :
1099 : static GEN
1100 796 : Fp_rootsof1(ulong l, GEN p)
1101 : {
1102 796 : GEN z, pl = diviuexact(subis(p,1),l);
1103 : ulong a;
1104 : forprime_t T;
1105 796 : u_forprime_init(&T, 3, ULONG_MAX);
1106 1062 : while((a = u_forprime_next(&T)))
1107 : {
1108 1062 : z = Fp_pow(utoi(a), pl, p);
1109 1062 : if (!equali1(z)) return z;
1110 : }
1111 0 : pari_err_PRIME("Fp_sqrt [modulus]",p);
1112 : return NULL; /* LCOV_EXCL_LINE */
1113 : }
1114 :
1115 : static GEN
1116 334 : Fp_gausssum(long D, GEN p)
1117 : {
1118 334 : long i, l = labs(D);
1119 334 : GEN z = Fp_rootsof1(l, p);
1120 334 : GEN s = z, x = z;
1121 3020 : for(i = 2; i < l; i++)
1122 : {
1123 2686 : long k = kross(i,l);
1124 2686 : x = mulii(x, z);
1125 2686 : if (k==1) s = addii(s, x);
1126 1510 : else if (k==-1) s = subii(s, x);
1127 : }
1128 334 : return s;
1129 : }
1130 :
1131 : static GEN
1132 19950 : Fp_sqrts(long a, GEN p)
1133 : {
1134 19950 : long v = vals(a)>>1;
1135 19949 : GEN r = gen_0;
1136 19949 : a >>= v << 1;
1137 19949 : switch(a)
1138 : {
1139 8 : case 1:
1140 8 : r = gen_1;
1141 8 : break;
1142 1121 : case -1:
1143 1121 : if (mod4(p)==1)
1144 1121 : r = Fp_pow(nonsquare_Fp(p), shifti(p,-2),p);
1145 : else
1146 0 : r = NULL;
1147 1121 : break;
1148 140 : case 2:
1149 140 : if (mod8(p)==1)
1150 : {
1151 140 : GEN z = Fp_pow(nonsquare_Fp(p), shifti(p,-3),p);
1152 140 : r = Fp_mul(z,Fp_sub(gen_1,Fp_sqr(z,p),p),p);
1153 0 : } else if (mod8(p)==7)
1154 0 : r = Fp_pow(gen_2, shifti(addiu(p,1),-2),p);
1155 : else
1156 0 : return NULL;
1157 140 : break;
1158 204 : case -2:
1159 204 : if (mod8(p)==1)
1160 : {
1161 204 : GEN z = Fp_pow(nonsquare_Fp(p), shifti(p,-3),p);
1162 204 : r = Fp_mul(z,Fp_add(gen_1,Fp_sqr(z,p),p),p);
1163 0 : } else if (mod8(p)==3)
1164 0 : r = Fp_pow(gen_m2, shifti(addiu(p,1),-2),p);
1165 : else
1166 0 : return NULL;
1167 204 : break;
1168 462 : case -3:
1169 462 : if (umodiu(p,3)==1)
1170 : {
1171 462 : GEN z = Fp_rootsof1(3, p);
1172 462 : r = Fp_sub(z,Fp_sqr(z,p),p);
1173 : }
1174 : else
1175 0 : return NULL;
1176 462 : break;
1177 2212 : case 5: case 13: case 17: case 21: case 29: case 33:
1178 : case -7: case -11: case -15: case -19: case -23:
1179 2212 : if (umodiu(p,labs(a))==1)
1180 334 : r = Fp_gausssum(a,p);
1181 : else
1182 1880 : return gen_0;
1183 334 : break;
1184 15802 : default:
1185 15802 : return gen_0;
1186 : }
1187 2269 : return remii(shifti(r, v), p);
1188 : }
1189 :
1190 : static GEN
1191 78831 : Fp_sqrt_ii(GEN a, GEN y, GEN p)
1192 : {
1193 78831 : pari_sp av = avma;
1194 78831 : GEN q, v, w, p1 = subiu(p,1);
1195 78834 : long i, k, e = vali(p1), as;
1196 :
1197 : /* direct formulas more efficient */
1198 78839 : if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
1199 78840 : if (e == 1)
1200 : {
1201 19702 : q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
1202 19699 : v = Fp_pow(a, q, p);
1203 : /* must check equality in case (a/p) = -1 or p not prime */
1204 19723 : av = avma; e = equalii(Fp_sqr(v,p), a); set_avma(av);
1205 19723 : return e? v: NULL;
1206 : }
1207 59138 : as = itos_or_0(a);
1208 59139 : if (!as) as = itos_or_0(subii(a,p));
1209 59146 : if (as)
1210 : {
1211 19951 : GEN res = Fp_sqrts(as, p);
1212 19950 : if (!res) return gc_NULL(av);
1213 19950 : if (signe(res)) return gerepileupto(av, res);
1214 : }
1215 56876 : if (e == 2)
1216 : { /* Atkin's formula */
1217 17624 : GEN I, a2 = shifti(a,1);
1218 17617 : if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
1219 17621 : q = shifti(p1, -3); /* (p-5)/8 */
1220 17618 : v = Fp_pow(a2, q, p);
1221 17627 : I = Fp_mul(a2, Fp_sqr(v,p), p); /* I^2 = -1 */
1222 17628 : v = Fp_mul(a, Fp_mul(v, subiu(I,1), p), p);
1223 : /* must check equality in case (a/p) = -1 or p not prime */
1224 17627 : av = avma; e = equalii(Fp_sqr(v,p), a); set_avma(av);
1225 17626 : return e? v: NULL;
1226 : }
1227 : /* On average, Cipolla is better than Tonelli/Shanks if and only if
1228 : * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
1229 39252 : if (e*(e-1) > 20 + 8 * expi(p)) return sqrt_Cipolla(a,p);
1230 : /* Tonelli-Shanks */
1231 39252 : av = avma; q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
1232 39251 : if (!y)
1233 : {
1234 2712 : y = Fp_2gener_all(q, p);
1235 2712 : if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
1236 : }
1237 39251 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
1238 39257 : v = Fp_mul(a, p1, p);
1239 39257 : w = Fp_mul(v, p1, p);
1240 93692 : while (!equali1(w))
1241 : { /* a*w = v^2, y primitive 2^e-th root of 1
1242 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1243 54475 : p1 = Fp_sqr(w,p);
1244 112012 : for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
1245 54479 : if (k == e) return NULL; /* p composite or (a/p) != 1 */
1246 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1247 54436 : p1 = y;
1248 77462 : for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
1249 54436 : y = Fp_sqr(p1, p); e = k;
1250 54438 : w = Fp_mul(y, w, p);
1251 54439 : v = Fp_mul(v, p1, p);
1252 54435 : if (gc_needed(av,1))
1253 : {
1254 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
1255 0 : gerepileall(av,3, &y,&w,&v);
1256 : }
1257 : }
1258 39214 : return v;
1259 : }
1260 :
1261 : /* Assume p is prime and return NULL if (a,p) = -1; y = NULL or generator
1262 : * of Fp^* 2-Sylow */
1263 : GEN
1264 4050555 : Fp_sqrt_i(GEN a, GEN y, GEN p)
1265 : {
1266 4050555 : pari_sp av = avma, av2;
1267 : GEN q;
1268 :
1269 4050555 : if (lgefint(p) == 3)
1270 : {
1271 3971615 : ulong pp = uel(p,2), u = umodiu(a, pp);
1272 3971630 : if (!u) return gen_0;
1273 3199089 : u = Fl_sqrt(u, pp);
1274 3199184 : return (u == ~0UL)? NULL: utoipos(u);
1275 : }
1276 78940 : a = modii(a, p); if (!signe(a)) return gen_0;
1277 78814 : a = Fp_sqrt_ii(a, y, p); if (!a) return gc_NULL(av);
1278 : /* smallest square root */
1279 78457 : av2 = avma; q = subii(p, a);
1280 78457 : if (cmpii(a, q) > 0) a = q; else set_avma(av2);
1281 78458 : return gerepileuptoint(av, a);
1282 : }
1283 : GEN
1284 3992178 : Fp_sqrt(GEN a, GEN p) { return Fp_sqrt_i(a, NULL, p); }
1285 :
1286 : /*********************************************************************/
1287 : /** GCD & BEZOUT **/
1288 : /*********************************************************************/
1289 :
1290 : GEN
1291 50355537 : lcmii(GEN x, GEN y)
1292 : {
1293 : pari_sp av;
1294 : GEN a, b;
1295 50355537 : if (!signe(x) || !signe(y)) return gen_0;
1296 50355610 : av = avma; a = gcdii(x,y);
1297 50353787 : if (absequalii(a,y)) { set_avma(av); return absi(x); }
1298 10567265 : if (!equali1(a)) y = diviiexact(y,a);
1299 10567297 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
1300 : }
1301 :
1302 : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
1303 : * set *pd = gcd(x,N) */
1304 : GEN
1305 5726871 : Fp_invgen(GEN x, GEN N, GEN *pd)
1306 : {
1307 : GEN d, d0, e, v;
1308 5726871 : if (lgefint(N) == 3)
1309 : {
1310 4971156 : ulong dd, NN = N[2], xx = umodiu(x,NN);
1311 4971679 : if (!xx) { *pd = N; return gen_0; }
1312 4971679 : xx = Fl_invgen(xx, NN, &dd);
1313 4972315 : *pd = utoi(dd); return utoi(xx);
1314 : }
1315 755715 : *pd = d = bezout(x, N, &v, NULL);
1316 755723 : if (equali1(d)) return v;
1317 : /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
1318 661210 : e = diviiexact(N,d);
1319 661210 : d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
1320 661210 : if (equali1(d0)) return v;
1321 523870 : if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
1322 523870 : return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
1323 : }
1324 :
1325 : /*********************************************************************/
1326 : /** CHINESE REMAINDERS **/
1327 : /*********************************************************************/
1328 :
1329 : /* Chinese Remainder Theorem. x and y must have the same type (integermod,
1330 : * polymod, or polynomial/vector/matrix recursively constructed with these
1331 : * as coefficients). Creates (with the same type) a z in the same residue
1332 : * class as x and the same residue class as y, if it is possible.
1333 : *
1334 : * We also allow (during recursion) two identical objects even if they are
1335 : * not integermod or polymod. For example:
1336 : *
1337 : * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
1338 : * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
1339 : * ? chinese(x, y)
1340 : * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
1341 :
1342 : static GEN
1343 2409938 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
1344 : {
1345 2409938 : GEN z = gassoc_proto(f,x,NULL);
1346 2409932 : if (z == gen_1) retmkintmod(gen_0,gen_1);
1347 2409897 : return z;
1348 : }
1349 :
1350 : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
1351 : * call chinese: makes Mod(0,1) a better "neutral" element */
1352 : static GEN
1353 21 : chinese_intpol(GEN x,GEN y)
1354 : {
1355 21 : pari_sp av = avma;
1356 21 : GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
1357 21 : return gerepileupto(av, chinese(z, y));
1358 : }
1359 :
1360 : GEN
1361 2345 : chinese1(GEN x) { return gen_chinese(x,chinese); }
1362 :
1363 : GEN
1364 5425 : chinese(GEN x, GEN y)
1365 : {
1366 : pari_sp av;
1367 5425 : long tx = typ(x), ty;
1368 : GEN z,p1,p2,d,u,v;
1369 :
1370 5425 : if (!y) return chinese1(x);
1371 5376 : if (gidentical(x,y)) return gcopy(x);
1372 5369 : ty = typ(y);
1373 5369 : if (tx == ty) switch(tx)
1374 : {
1375 3822 : case t_POLMOD:
1376 : {
1377 3822 : GEN A = gel(x,1), B = gel(y,1);
1378 3822 : GEN a = gel(x,2), b = gel(y,2);
1379 3822 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
1380 3822 : if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
1381 3822 : av = avma;
1382 3822 : d = RgX_extgcd(A,B,&u,&v);
1383 3822 : p2 = gsub(b, a);
1384 3822 : if (!gequal0(gmod(p2, d))) break;
1385 3822 : p1 = gdiv(A,d);
1386 3822 : p2 = gadd(a, gmul(gmul(u,p1), p2));
1387 :
1388 3822 : z = cgetg(3, t_POLMOD);
1389 3822 : gel(z,1) = gmul(p1,B);
1390 3822 : gel(z,2) = gmod(p2,gel(z,1));
1391 3822 : return gerepileupto(av, z);
1392 : }
1393 1505 : case t_INTMOD:
1394 : {
1395 1505 : GEN A = gel(x,1), B = gel(y,1);
1396 1505 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
1397 1505 : z = cgetg(3,t_INTMOD);
1398 1505 : Z_chinese_pre(A, B, &C, &U, &d);
1399 1505 : c = Z_chinese_post(a, b, C, U, d);
1400 1505 : if (!c) pari_err_OP("chinese", x,y);
1401 1505 : set_avma((pari_sp)z);
1402 1505 : gel(z,1) = icopy(C);
1403 1505 : gel(z,2) = icopy(c); return z;
1404 : }
1405 14 : case t_POL:
1406 : {
1407 14 : long i, lx = lg(x), ly = lg(y);
1408 14 : if (varn(x) != varn(y)) break;
1409 14 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
1410 14 : z = cgetg(lx, t_POL); z[1] = x[1];
1411 42 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
1412 14 : if (i < lx)
1413 : {
1414 14 : GEN _0 = Rg_get_0(y);
1415 28 : for ( ; i<lx; i++) gel(z,i) = chinese(gel(x,i),_0);
1416 : }
1417 14 : return z;
1418 : }
1419 7 : case t_VEC: case t_COL: case t_MAT:
1420 : {
1421 : long i, lx;
1422 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
1423 21 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
1424 7 : return z;
1425 : }
1426 : }
1427 21 : if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
1428 7 : if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
1429 0 : pari_err_OP("chinese",x,y);
1430 : return NULL; /* LCOV_EXCL_LINE */
1431 : }
1432 :
1433 : /* init chinese(Mod(.,A), Mod(.,B)) */
1434 : void
1435 271319 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
1436 : {
1437 271319 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
1438 271319 : GEN t = diviiexact(A,d);
1439 271306 : *pU = mulii(u, t);
1440 271302 : *pC = mulii(t, B);
1441 271300 : if (pd) *pd = d;
1442 271300 : }
1443 : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
1444 : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
1445 : * If d not NULL, check whether a = b mod d. */
1446 : GEN
1447 2985129 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
1448 : {
1449 : GEN b_a;
1450 2985129 : if (!signe(a))
1451 : {
1452 794745 : if (d && !dvdii(b, d)) return NULL;
1453 794745 : return Fp_mul(b, U, C);
1454 : }
1455 2190384 : b_a = subii(b,a);
1456 2190384 : if (d && !dvdii(b_a, d)) return NULL;
1457 2190384 : return modii(addii(a, mulii(U, b_a)), C);
1458 : }
1459 : static ulong
1460 1519838 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
1461 : {
1462 1519838 : if (!a) return Fl_mul(b, U, C);
1463 1519838 : return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
1464 : }
1465 :
1466 : GEN
1467 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
1468 : {
1469 2142 : pari_sp av = avma;
1470 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
1471 2142 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
1472 : }
1473 : GEN
1474 267614 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
1475 : {
1476 267614 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
1477 267596 : return Z_chinese_post(a,b, *pC, U, NULL);
1478 : }
1479 :
1480 : /* return lift(chinese(a mod A, b mod B))
1481 : * assume(A,B)=1, a,b,A,B integers and C = A*B */
1482 : GEN
1483 525130 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
1484 : {
1485 525130 : pari_sp av = avma;
1486 525130 : GEN U = mulii(Fp_inv(A,B), A);
1487 525130 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
1488 : }
1489 : ulong
1490 1519820 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
1491 1519820 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
1492 :
1493 : /* chinese1 for coprime moduli in Z */
1494 : static GEN
1495 2188432 : chinese1_coprime_Z_aux(GEN x, GEN y)
1496 : {
1497 2188432 : GEN z = cgetg(3, t_INTMOD);
1498 2188432 : GEN A = gel(x,1), a = gel(x, 2);
1499 2188432 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
1500 2188432 : pari_sp av = avma;
1501 2188432 : GEN U = mulii(Fp_inv(A,B), A);
1502 2188432 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
1503 2188432 : gel(z,1) = C; return z;
1504 : }
1505 : GEN
1506 2407593 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
1507 :
1508 : /*********************************************************************/
1509 : /** MODULAR EXPONENTIATION **/
1510 : /*********************************************************************/
1511 : /* xa ZV or nv */
1512 : GEN
1513 2427875 : ZV_producttree(GEN xa)
1514 : {
1515 2427875 : long n = lg(xa)-1;
1516 2427875 : long m = n==1 ? 1: expu(n-1)+1;
1517 2427875 : GEN T = cgetg(m+1, t_VEC), t;
1518 : long i, j, k;
1519 2427874 : t = cgetg(((n+1)>>1)+1, t_VEC);
1520 2427875 : if (typ(xa)==t_VECSMALL)
1521 : {
1522 3197015 : for (j=1, k=1; k<n; j++, k+=2)
1523 2051394 : gel(t, j) = muluu(xa[k], xa[k+1]);
1524 1145621 : if (k==n) gel(t, j) = utoi(xa[k]);
1525 : } else {
1526 2655484 : for (j=1, k=1; k<n; j++, k+=2)
1527 1373225 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
1528 1282259 : if (k==n) gel(t, j) = icopy(gel(xa,k));
1529 : }
1530 2427880 : gel(T,1) = t;
1531 3850388 : for (i=2; i<=m; i++)
1532 : {
1533 1422499 : GEN u = gel(T, i-1);
1534 1422499 : long n = lg(u)-1;
1535 1422499 : t = cgetg(((n+1)>>1)+1, t_VEC);
1536 3192889 : for (j=1, k=1; k<n; j++, k+=2)
1537 1770381 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
1538 1422508 : if (k==n) gel(t, j) = gel(u, k);
1539 1422508 : gel(T, i) = t;
1540 : }
1541 2427889 : return T;
1542 : }
1543 :
1544 : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
1545 : GEN
1546 58133003 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
1547 : {
1548 : long i,j,k;
1549 58133003 : long m = lg(T)-1, n = lg(P)-1;
1550 : GEN t;
1551 58133003 : GEN Tp = cgetg(m+1, t_VEC);
1552 58088788 : gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
1553 122705723 : for (i=m-1; i>=1; i--)
1554 : {
1555 64711656 : GEN u = gel(T, i);
1556 64711656 : GEN v = gel(Tp, i+1);
1557 64711656 : long n = lg(u)-1;
1558 64711656 : t = cgetg(n+1, t_VEC);
1559 159134523 : for (j=1, k=1; k<n; j++, k+=2)
1560 : {
1561 94572679 : gel(t, k) = modii(gel(v, j), gel(u, k));
1562 94584994 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
1563 : }
1564 64561844 : if (k==n) gel(t, k) = gel(v, j);
1565 64561844 : gel(Tp, i) = t;
1566 : }
1567 : {
1568 57994067 : GEN u = gel(T, i+1);
1569 57994067 : GEN v = gel(Tp, i+1);
1570 57994067 : long l = lg(u)-1;
1571 57994067 : if (typ(P)==t_VECSMALL)
1572 : {
1573 55572734 : GEN R = cgetg(n+1, t_VECSMALL);
1574 205105013 : for (j=1, k=1; j<=l; j++, k+=2)
1575 : {
1576 149260617 : uel(R,k) = umodiu(gel(v, j), P[k]);
1577 149514864 : if (k < n)
1578 118339290 : uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
1579 : }
1580 55844396 : return R;
1581 : }
1582 : else
1583 : {
1584 2421333 : GEN R = cgetg(n+1, t_VEC);
1585 6622428 : for (j=1, k=1; j<=l; j++, k+=2)
1586 : {
1587 4194905 : gel(R,k) = modii(gel(v, j), gel(P,k));
1588 4194920 : if (k < n)
1589 3421388 : gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
1590 : }
1591 2427523 : return R;
1592 : }
1593 : }
1594 : }
1595 :
1596 : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
1597 : GEN
1598 39559549 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
1599 : {
1600 39559549 : long m = lg(T)-1, n = lg(A)-1;
1601 : long i,j,k;
1602 39559549 : GEN Tp = cgetg(m+1, t_VEC);
1603 39542923 : GEN M = gel(T, 1);
1604 39542923 : GEN t = cgetg(lg(M), t_VEC);
1605 39507070 : if (typ(P)==t_VECSMALL)
1606 : {
1607 88465756 : for (j=1, k=1; k<n; j++, k+=2)
1608 : {
1609 64985802 : pari_sp av = avma;
1610 64985802 : GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
1611 64779482 : GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
1612 64936789 : gel(t, j) = gerepileuptoint(av, tj);
1613 : }
1614 23479954 : if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
1615 : } else
1616 : {
1617 34185781 : for (j=1, k=1; k<n; j++, k+=2)
1618 : {
1619 18131407 : pari_sp av = avma;
1620 18131407 : GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
1621 18140285 : GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
1622 18166282 : gel(t, j) = gerepileuptoint(av, tj);
1623 : }
1624 16054374 : if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
1625 : }
1626 39524407 : gel(Tp, 1) = t;
1627 74948606 : for (i=2; i<=m; i++)
1628 : {
1629 35417848 : GEN u = gel(T, i-1), M = gel(T, i);
1630 35417848 : GEN t = cgetg(lg(M), t_VEC);
1631 35456700 : GEN v = gel(Tp, i-1);
1632 35456700 : long n = lg(v)-1;
1633 95083643 : for (j=1, k=1; k<n; j++, k+=2)
1634 : {
1635 59659444 : pari_sp av = avma;
1636 59607749 : gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
1637 59659444 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
1638 : }
1639 35424199 : if (k==n) gel(t, j) = gel(v, k);
1640 35424199 : gel(Tp, i) = t;
1641 : }
1642 39530758 : return gmael(Tp,m,1);
1643 : }
1644 :
1645 : static GEN
1646 1464569 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1647 : {
1648 1464569 : long i, l = lg(gel(vA,1)), n = lg(P);
1649 1464569 : GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
1650 35048787 : for (i=1; i < l; i++)
1651 : {
1652 33584305 : pari_sp av = avma;
1653 33584305 : GEN c, A = cgetg(n, typ(P));
1654 : long j;
1655 199406713 : for (j=1; j < n; j++) A[j] = mael(vA,j,i);
1656 33554978 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
1657 33583774 : gel(V,i) = gerepileuptoint(av, c);
1658 : }
1659 1464482 : return V;
1660 : }
1661 :
1662 : static GEN
1663 633295 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1664 : {
1665 633295 : long i, j, l, n = lg(P);
1666 633295 : GEN mod = gmael(T, lg(T)-1, 1), V, w;
1667 633295 : w = cgetg(n, t_VECSMALL);
1668 2173282 : for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
1669 633286 : l = vecsmall_max(w);
1670 633283 : V = cgetg(l, t_POL);
1671 633267 : V[1] = mael(vA,1,1);
1672 3872284 : for (i=2; i < l; i++)
1673 : {
1674 3239005 : pari_sp av = avma;
1675 3239005 : GEN c, A = cgetg(n, typ(P));
1676 3238719 : if (typ(P)==t_VECSMALL)
1677 7769962 : for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
1678 : else
1679 4357071 : for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
1680 3238719 : c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
1681 3239013 : gel(V,i) = gerepileuptoint(av, c);
1682 : }
1683 633279 : return ZX_renormalize(V, l);
1684 : }
1685 :
1686 : static GEN
1687 4614 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1688 : {
1689 4614 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1690 4614 : GEN A = cgetg(n, t_VEC);
1691 4614 : GEN V = cgetg(l, t_COL);
1692 90964 : for (i=1; i < l; i++)
1693 : {
1694 335290 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1695 86350 : gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
1696 : }
1697 4614 : return V;
1698 : }
1699 :
1700 : static GEN
1701 367984 : polint_chinese(GEN worker, GEN mA, GEN P)
1702 : {
1703 367984 : long cnt, pending, n, i, j, l = lg(gel(mA,1));
1704 : struct pari_mt pt;
1705 : GEN done, va, M, A;
1706 : pari_timer ti;
1707 :
1708 367984 : if (l == 1) return cgetg(1, t_MAT);
1709 338932 : cnt = pending = 0;
1710 338932 : n = lg(P);
1711 338932 : A = cgetg(n, t_VEC);
1712 338932 : va = mkvec(A);
1713 338932 : M = cgetg(l, t_MAT);
1714 338932 : if (DEBUGLEVEL>4) timer_start(&ti);
1715 338932 : if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
1716 338932 : mt_queue_start_lim(&pt, worker, l-1);
1717 1319992 : for (i=1; i<l || pending; i++)
1718 : {
1719 : long workid;
1720 3722741 : for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
1721 981060 : mt_queue_submit(&pt, i, i<l? va: NULL);
1722 981060 : done = mt_queue_get(&pt, &workid, &pending);
1723 981060 : if (done)
1724 : {
1725 936000 : gel(M,workid) = done;
1726 936000 : if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
1727 : }
1728 : }
1729 338932 : if (DEBUGLEVEL>5) err_printf("\n");
1730 338932 : if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
1731 338932 : mt_queue_end(&pt);
1732 338932 : return M;
1733 : }
1734 :
1735 : GEN
1736 840 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
1737 : {
1738 840 : return nxCV_polint_center_tree(vA, P, T, R, m2);
1739 : }
1740 :
1741 : static GEN
1742 430 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1743 : {
1744 430 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1745 430 : GEN A = cgetg(n, t_VEC);
1746 430 : GEN V = cgetg(l, t_MAT);
1747 4204 : for (i=1; i < l; i++)
1748 : {
1749 15299 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1750 3774 : gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
1751 : }
1752 430 : return V;
1753 : }
1754 :
1755 : static GEN
1756 90 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
1757 : {
1758 90 : GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
1759 90 : return polint_chinese(worker, mA, P);
1760 : }
1761 :
1762 : static GEN
1763 111717 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
1764 : {
1765 111717 : long i, j, l = lg(gel(vA,1)), n = lg(P);
1766 111717 : GEN A = cgetg(n, t_VEC);
1767 111717 : GEN V = cgetg(l, t_MAT);
1768 625556 : for (i=1; i < l; i++)
1769 : {
1770 3000484 : for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
1771 513838 : gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
1772 : }
1773 111718 : return V;
1774 : }
1775 :
1776 : GEN
1777 935085 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
1778 : {
1779 935085 : return ncV_polint_center_tree(vA, P, T, R, m2);
1780 : }
1781 :
1782 : static GEN
1783 367894 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
1784 : {
1785 367894 : GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
1786 367894 : return polint_chinese(worker, mA, P);
1787 : }
1788 :
1789 : /* return [A mod P[i], i=1..#P] */
1790 : GEN
1791 0 : Z_ZV_mod(GEN A, GEN P)
1792 : {
1793 0 : pari_sp av = avma;
1794 0 : return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
1795 : }
1796 : /* P a t_VECSMALL */
1797 : GEN
1798 0 : Z_nv_mod(GEN A, GEN P)
1799 : {
1800 0 : pari_sp av = avma;
1801 0 : return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
1802 : }
1803 : /* B a ZX, T = ZV_producttree(P) */
1804 : GEN
1805 1461644 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
1806 : {
1807 : pari_sp av;
1808 1461644 : long i, j, l = lg(B), n = lg(A)-1;
1809 1461644 : GEN V = cgetg(n+1, t_VEC);
1810 7164651 : for (j=1; j <= n; j++)
1811 : {
1812 5703140 : gel(V, j) = cgetg(l, t_VECSMALL);
1813 5703054 : mael(V, j, 1) = B[1]&VARNBITS;
1814 : }
1815 1461511 : av = avma;
1816 12772154 : for (i=2; i < l; i++)
1817 : {
1818 11311117 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
1819 78686611 : for (j=1; j <= n; j++)
1820 67382708 : mael(V, j, i) = v[j];
1821 11303903 : set_avma(av);
1822 : }
1823 7164809 : for (j=1; j <= n; j++)
1824 5703634 : (void) Flx_renormalize(gel(V, j), l);
1825 1461175 : return V;
1826 : }
1827 :
1828 : static GEN
1829 299554 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
1830 :
1831 : GEN
1832 31088 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
1833 : {
1834 31088 : pari_sp av = avma;
1835 31088 : long i, j, l = lg(P), n = lg(xa)-1;
1836 31088 : GEN V = cgetg(n+1, t_VEC);
1837 126462 : for (j=1; j <= n; j++)
1838 : {
1839 95374 : gel(V, j) = cgetg(l, t_POL);
1840 95374 : mael(V, j, 1) = P[1]&VARNBITS;
1841 : }
1842 249451 : for (i=2; i < l; i++)
1843 : {
1844 218363 : GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
1845 868229 : for (j=1; j <= n; j++)
1846 649866 : gmael(V, j, i) = gel(v,j);
1847 : }
1848 126462 : for (j=1; j <= n; j++)
1849 95374 : (void) FlxX_renormalize(gel(V, j), l);
1850 31088 : return gerepilecopy(av, V);
1851 : }
1852 :
1853 : GEN
1854 4048 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
1855 : {
1856 4048 : pari_sp av = avma;
1857 4048 : long i, j, l = lg(C), n = lg(xa)-1;
1858 4048 : GEN V = cgetg(n+1, t_VEC);
1859 16948 : for (j = 1; j <= n; j++)
1860 12901 : gel(V, j) = cgetg(l, t_COL);
1861 85229 : for (i = 1; i < l; i++)
1862 : {
1863 81181 : GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
1864 359437 : for (j = 1; j <= n; j++)
1865 278255 : gmael(V, j, i) = gel(v,j);
1866 : }
1867 4048 : return gerepilecopy(av, V);
1868 : }
1869 :
1870 : GEN
1871 430 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
1872 : {
1873 430 : pari_sp av = avma;
1874 430 : long i, j, l = lg(M), n = lg(xa)-1;
1875 430 : GEN V = cgetg(n+1, t_VEC);
1876 2083 : for (j=1; j <= n; j++)
1877 1653 : gel(V, j) = cgetg(l, t_MAT);
1878 4204 : for (i=1; i < l; i++)
1879 : {
1880 3774 : GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
1881 15299 : for (j=1; j <= n; j++)
1882 11525 : gmael(V, j, i) = gel(v,j);
1883 : }
1884 430 : return gerepilecopy(av, V);
1885 : }
1886 :
1887 : GEN
1888 1208147 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
1889 : {
1890 : pari_sp av;
1891 1208147 : long i, j, l = lg(B), n = lg(A)-1;
1892 1208147 : GEN V = cgetg(n+1, t_VEC);
1893 6206281 : for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
1894 1208075 : av = avma;
1895 45584213 : for (i=1; i < l; i++)
1896 : {
1897 44379361 : GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
1898 246548781 : for (j=1; j <= n; j++) mael(V, j, i) = v[j];
1899 44336857 : set_avma(av);
1900 : }
1901 1204852 : return V;
1902 : }
1903 :
1904 : GEN
1905 178606 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
1906 : {
1907 178606 : pari_sp av = avma;
1908 178606 : long i, j, l = lg(M), n = lg(xa)-1;
1909 178606 : GEN V = cgetg(n+1, t_VEC);
1910 893675 : for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_MAT);
1911 1386516 : for (i=1; i < l; i++)
1912 : {
1913 1207916 : GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
1914 6205705 : for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
1915 : }
1916 178600 : return gerepilecopy(av, V);
1917 : }
1918 :
1919 : static GEN
1920 2423668 : ZV_sqr(GEN z)
1921 : {
1922 2423668 : long i,l = lg(z);
1923 2423668 : GEN x = cgetg(l, t_VEC);
1924 2423666 : if (typ(z)==t_VECSMALL)
1925 5657048 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
1926 : else
1927 4357307 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
1928 2423655 : return x;
1929 : }
1930 :
1931 : static GEN
1932 12438961 : ZT_sqr(GEN x)
1933 : {
1934 12438961 : if (typ(x) == t_INT) return sqri(x);
1935 16278675 : pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
1936 : }
1937 :
1938 : static GEN
1939 2423662 : ZV_invdivexact(GEN y, GEN x)
1940 : {
1941 2423662 : long i, l = lg(y);
1942 2423662 : GEN z = cgetg(l,t_VEC);
1943 2423659 : if (typ(x)==t_VECSMALL)
1944 5656772 : for (i=1; i<l; i++)
1945 : {
1946 4511514 : pari_sp av = avma;
1947 4511514 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
1948 4511837 : set_avma(av); gel(z,i) = utoi(a);
1949 : }
1950 : else
1951 4357304 : for (i=1; i<l; i++)
1952 3078908 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
1953 2423654 : return z;
1954 : }
1955 :
1956 : /* P t_VECSMALL or t_VEC of t_INT */
1957 : GEN
1958 2423657 : ZV_chinesetree(GEN P, GEN T)
1959 : {
1960 2423657 : GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
1961 2423657 : GEN mod = gmael(T,lg(T)-1,1);
1962 2423657 : return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
1963 : }
1964 :
1965 : static GEN
1966 928962 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
1967 : {
1968 928962 : if (!pt_mod)
1969 12237 : return gerepileupto(av, a);
1970 : else
1971 : {
1972 916725 : GEN mod = gmael(T, lg(T)-1, 1);
1973 916725 : gerepileall(av, 2, &a, &mod);
1974 916725 : *pt_mod = mod;
1975 916725 : return a;
1976 : }
1977 : }
1978 :
1979 : GEN
1980 157963 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
1981 : {
1982 157963 : pari_sp av = avma;
1983 157963 : GEN T = ZV_producttree(P);
1984 157963 : GEN R = ZV_chinesetree(P, T);
1985 157962 : GEN a = ZV_chinese_tree(A, P, T, R);
1986 157963 : GEN mod = gmael(T, lg(T)-1, 1);
1987 157963 : GEN ca = Fp_center(a, mod, shifti(mod,-1));
1988 157963 : return gc_chinese(av, T, ca, pt_mod);
1989 : }
1990 :
1991 : GEN
1992 5110 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
1993 : {
1994 5110 : pari_sp av = avma;
1995 5110 : GEN T = ZV_producttree(P);
1996 5110 : GEN R = ZV_chinesetree(P, T);
1997 5110 : GEN a = ZV_chinese_tree(A, P, T, R);
1998 5110 : return gc_chinese(av, T, a, pt_mod);
1999 : }
2000 :
2001 : GEN
2002 159223 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2003 : {
2004 159223 : pari_sp av = avma;
2005 159223 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2006 159223 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
2007 159224 : return gerepileupto(av, a);
2008 : }
2009 :
2010 : GEN
2011 387712 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2012 : {
2013 387712 : pari_sp av = avma;
2014 387712 : GEN T = ZV_producttree(P);
2015 387711 : GEN R = ZV_chinesetree(P, T);
2016 387712 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2017 387712 : GEN a = nxV_polint_center_tree(A, P, T, R, m2);
2018 387712 : return gc_chinese(av, T, a, pt_mod);
2019 : }
2020 :
2021 : GEN
2022 10193 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2023 : {
2024 10193 : pari_sp av = avma;
2025 10193 : GEN T = ZV_producttree(P);
2026 10193 : GEN R = ZV_chinesetree(P, T);
2027 10193 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2028 10193 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
2029 10193 : return gc_chinese(av, T, a, pt_mod);
2030 : }
2031 :
2032 : GEN
2033 5457 : ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2034 : {
2035 5457 : pari_sp av = avma;
2036 5457 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2037 5457 : GEN a = ncV_polint_center_tree(A, P, T, R, m2);
2038 5457 : return gerepileupto(av, a);
2039 : }
2040 :
2041 : GEN
2042 0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2043 : {
2044 0 : pari_sp av = avma;
2045 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2046 0 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
2047 0 : return gerepileupto(av, a);
2048 : }
2049 :
2050 : GEN
2051 111717 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
2052 : {
2053 111717 : pari_sp av = avma;
2054 111717 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2055 111717 : GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
2056 111718 : return gerepileupto(av, a);
2057 : }
2058 :
2059 : GEN
2060 367894 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2061 : {
2062 367894 : pari_sp av = avma;
2063 367894 : GEN T = ZV_producttree(P);
2064 367894 : GEN R = ZV_chinesetree(P, T);
2065 367894 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2066 367894 : GEN a = nmV_polint_center_tree(A, P, T, R, m2);
2067 367894 : return gc_chinese(av, T, a, pt_mod);
2068 : }
2069 :
2070 : GEN
2071 0 : nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
2072 : {
2073 0 : pari_sp av = avma;
2074 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2075 0 : GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
2076 0 : return gerepileupto(av, a);
2077 : }
2078 :
2079 : GEN
2080 0 : nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2081 : {
2082 0 : pari_sp av = avma;
2083 0 : GEN T = ZV_producttree(P);
2084 0 : GEN R = ZV_chinesetree(P, T);
2085 0 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2086 0 : GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
2087 0 : return gc_chinese(av, T, a, pt_mod);
2088 : }
2089 :
2090 : GEN
2091 430 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
2092 : {
2093 430 : pari_sp av = avma;
2094 430 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2095 430 : GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
2096 430 : return gerepileupto(av, a);
2097 : }
2098 :
2099 : GEN
2100 90 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2101 : {
2102 90 : pari_sp av = avma;
2103 90 : GEN T = ZV_producttree(P);
2104 90 : GEN R = ZV_chinesetree(P, T);
2105 90 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2106 90 : GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
2107 90 : return gc_chinese(av, T, a, pt_mod);
2108 : }
2109 :
2110 : /**********************************************************************
2111 : ** Powering over (Z/NZ)^*, small N **
2112 : **********************************************************************/
2113 :
2114 : /* 2^n mod p; assume n > 1 */
2115 : static ulong
2116 11968001 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
2117 : {
2118 11968001 : ulong y = 2;
2119 11968001 : int j = 1+bfffo(n);
2120 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
2121 11968001 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
2122 522321349 : for (; j; n<<=1,j--)
2123 : {
2124 510374048 : y = Fl_sqr_pre(y,p,pi);
2125 510345403 : if (n & HIGHBIT) y = Fl_double(y, p);
2126 : }
2127 11947301 : return y;
2128 : }
2129 :
2130 : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
2131 : static ulong
2132 4162443 : Fl_2powu(ulong n, ulong p)
2133 : {
2134 4162443 : ulong y = 2;
2135 4162443 : int j = 1+bfffo(n);
2136 : /* normalize, i.e set highest bit to 1 (we know n != 0) */
2137 4162443 : n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
2138 23256778 : for (; j; n<<=1,j--)
2139 : {
2140 19094311 : y = (y*y) % p;
2141 19094311 : if (n & HIGHBIT) y = Fl_double(y, p);
2142 : }
2143 4162467 : return y;
2144 : }
2145 :
2146 : /* allow pi = 0 */
2147 : ulong
2148 98147159 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
2149 : {
2150 : ulong y, z, n;
2151 98147159 : if (!pi) return Fl_powu(x, n0, p);
2152 95727094 : if (n0 <= 1)
2153 : { /* frequent special cases */
2154 6998757 : if (n0 == 1) return x;
2155 103729 : if (n0 == 0) return 1;
2156 : }
2157 88728324 : if (x <= 2)
2158 : {
2159 12226432 : if (x == 2) return Fl_2powu_pre(n0, p, pi);
2160 257096 : return x; /* 0 or 1 */
2161 : }
2162 76501892 : y = 1; z = x; n = n0;
2163 : for(;;)
2164 : {
2165 540560469 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
2166 540893733 : n>>=1; if (!n) return y;
2167 464176076 : z = Fl_sqr_pre(z,p,pi);
2168 : }
2169 : }
2170 :
2171 : ulong
2172 145795767 : Fl_powu(ulong x, ulong n0, ulong p)
2173 : {
2174 : ulong y, z, n;
2175 145795767 : if (n0 <= 2)
2176 : { /* frequent special cases */
2177 65627039 : if (n0 == 2) return Fl_sqr(x,p);
2178 32068743 : if (n0 == 1) return x;
2179 1905670 : if (n0 == 0) return 1;
2180 : }
2181 80137117 : if (x <= 1) return x; /* 0 or 1 */
2182 79705382 : if (p & HIGHMASK)
2183 7824634 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
2184 71880748 : if (x == 2) return Fl_2powu(n0, p);
2185 67718270 : y = 1; z = x; n = n0;
2186 : for(;;)
2187 : {
2188 358630368 : if (n&1) y = (y*z) % p;
2189 358630368 : n>>=1; if (!n) return y;
2190 290912098 : z = (z*z) % p;
2191 : }
2192 : }
2193 :
2194 : /* Reduce data dependency to maximize internal parallelism; allow pi = 0 */
2195 : GEN
2196 12795123 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
2197 : {
2198 : long i, k;
2199 12795123 : GEN z = cgetg(n + 2, t_VECSMALL);
2200 12785174 : z[1] = 1; if (n == 0) return z;
2201 12785174 : z[2] = x;
2202 12785174 : if (pi)
2203 : {
2204 90087526 : for (i = 3, k=2; i <= n; i+=2, k++)
2205 : {
2206 77502988 : z[i] = Fl_sqr_pre(z[k], p, pi);
2207 77503662 : z[i+1] = Fl_mul_pre(z[k], z[k+1], p, pi);
2208 : }
2209 12584538 : if (i==n+1) z[i] = Fl_sqr_pre(z[k], p, pi);
2210 : }
2211 212231 : else if (p & HIGHMASK)
2212 : {
2213 0 : for (i = 3, k=2; i <= n; i+=2, k++)
2214 : {
2215 0 : z[i] = Fl_sqr(z[k], p);
2216 0 : z[i+1] = Fl_mul(z[k], z[k+1], p);
2217 : }
2218 0 : if (i==n+1) z[i] = Fl_sqr(z[k], p);
2219 : }
2220 : else
2221 400553193 : for (i = 2; i <= n; i++) z[i+1] = (z[i] * x) % p;
2222 12796156 : return z;
2223 : }
2224 :
2225 : GEN
2226 295917 : Fl_powers(ulong x, long n, ulong p)
2227 : {
2228 295917 : return Fl_powers_pre(x, n, p, (p & HIGHMASK)? get_Fl_red(p): 0);
2229 : }
2230 :
2231 : /**********************************************************************
2232 : ** Powering over (Z/NZ)^*, large N **
2233 : **********************************************************************/
2234 : typedef struct muldata {
2235 : GEN (*sqr)(void * E, GEN x);
2236 : GEN (*mul)(void * E, GEN x, GEN y);
2237 : GEN (*mul2)(void * E, GEN x);
2238 : } muldata;
2239 :
2240 : /* modified Barrett reduction with one fold */
2241 : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
2242 :
2243 : static GEN
2244 14989 : Fp_invmBarrett(GEN p, long s)
2245 : {
2246 14989 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
2247 14989 : return mkvec2(Q,R);
2248 : }
2249 :
2250 : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
2251 : * a = r (mod N) */
2252 : static GEN
2253 9203382 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
2254 : {
2255 9203382 : pari_sp av = avma;
2256 9203382 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
2257 9203382 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
2258 9203382 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
2259 9203382 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
2260 9203382 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
2261 9203382 : GEN r = subii(A, mulii(q, N));
2262 9203382 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
2263 9203382 : if (signe(sr)<0) return gerepileuptoint(av, r);
2264 5021321 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
2265 5021321 : if (signe(sr)<0) return gerepileuptoint(av, r);
2266 185579 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
2267 185579 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
2268 : }
2269 :
2270 : /* Montgomery reduction */
2271 :
2272 : INLINE ulong
2273 696029 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
2274 :
2275 : struct montred
2276 : {
2277 : GEN N;
2278 : ulong inv;
2279 : };
2280 :
2281 : /* Montgomery reduction */
2282 : static GEN
2283 69582758 : _sqr_montred(void * E, GEN x)
2284 : {
2285 69582758 : struct montred * D = (struct montred *) E;
2286 69582758 : return red_montgomery(sqri(x), D->N, D->inv);
2287 : }
2288 :
2289 : /* Montgomery reduction */
2290 : static GEN
2291 7194723 : _mul_montred(void * E, GEN x, GEN y)
2292 : {
2293 7194723 : struct montred * D = (struct montred *) E;
2294 7194723 : return red_montgomery(mulii(x, y), D->N, D->inv);
2295 : }
2296 :
2297 : static GEN
2298 11257178 : _mul2_montred(void * E, GEN x)
2299 : {
2300 11257178 : struct montred * D = (struct montred *) E;
2301 11257178 : GEN z = shifti(_sqr_montred(E, x), 1);
2302 11255672 : long l = lgefint(D->N);
2303 11894320 : while (lgefint(z) > l) z = subii(z, D->N);
2304 11256139 : return z;
2305 : }
2306 :
2307 : static GEN
2308 23299752 : _sqr_remii(void* N, GEN x)
2309 23299752 : { return remii(sqri(x), (GEN) N); }
2310 :
2311 : static GEN
2312 1560727 : _mul_remii(void* N, GEN x, GEN y)
2313 1560727 : { return remii(mulii(x, y), (GEN) N); }
2314 :
2315 : static GEN
2316 3227521 : _mul2_remii(void* N, GEN x)
2317 3227521 : { return Fp_double(_sqr_remii(N, x), (GEN)N); }
2318 :
2319 : struct redbarrett
2320 : {
2321 : GEN iM, N;
2322 : long s;
2323 : };
2324 :
2325 : static GEN
2326 8429661 : _sqr_remiibar(void *E, GEN x)
2327 : {
2328 8429661 : struct redbarrett * D = (struct redbarrett *) E;
2329 8429661 : return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
2330 : }
2331 :
2332 : static GEN
2333 773721 : _mul_remiibar(void *E, GEN x, GEN y)
2334 : {
2335 773721 : struct redbarrett * D = (struct redbarrett *) E;
2336 773721 : return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
2337 : }
2338 :
2339 : static GEN
2340 2077445 : _mul2_remiibar(void *E, GEN x)
2341 : {
2342 2077445 : struct redbarrett * D = (struct redbarrett *) E;
2343 2077445 : return Fp_double(_sqr_remiibar(E, x), D->N);
2344 : }
2345 :
2346 : static long
2347 888305 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
2348 : {
2349 888305 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
2350 : {
2351 14989 : struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
2352 14989 : D->sqr = &_sqr_remiibar;
2353 14989 : D->mul = &_mul_remiibar;
2354 14989 : D->mul2 = &_mul2_remiibar;
2355 14989 : E->N = N;
2356 14989 : E->s = 1+(expi(N)>>1);
2357 14989 : E->iM = Fp_invmBarrett(N, E->s);
2358 14989 : *pt_E = (void*) E;
2359 14989 : return 0;
2360 : }
2361 873316 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
2362 : {
2363 696037 : struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
2364 696036 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
2365 696034 : D->sqr = &_sqr_montred;
2366 696034 : D->mul = &_mul_montred;
2367 696034 : D->mul2 = &_mul2_montred;
2368 696034 : E->N = N;
2369 696034 : E->inv = init_montdata(N);
2370 696026 : *pt_E = (void*) E;
2371 696026 : return 1;
2372 : }
2373 : else
2374 : {
2375 177274 : D->sqr = &_sqr_remii;
2376 177274 : D->mul = &_mul_remii;
2377 177274 : D->mul2 = &_mul2_remii;
2378 177274 : *pt_E = (void*) N;
2379 177274 : return 0;
2380 : }
2381 : }
2382 :
2383 : GEN
2384 2870401 : Fp_powu(GEN A, ulong k, GEN N)
2385 : {
2386 2870401 : long lN = lgefint(N);
2387 : int base_is_2, use_montgomery;
2388 : muldata D;
2389 : void *E;
2390 : pari_sp av;
2391 :
2392 2870401 : if (lN == 3) {
2393 1438672 : ulong n = uel(N,2);
2394 1438672 : return utoi( Fl_powu(umodiu(A, n), k, n) );
2395 : }
2396 1431729 : if (k <= 2)
2397 : { /* frequent special cases */
2398 815776 : if (k == 2) return Fp_sqr(A,N);
2399 269529 : if (k == 1) return A;
2400 0 : if (k == 0) return gen_1;
2401 : }
2402 615953 : av = avma; A = modii(A,N);
2403 615953 : base_is_2 = 0;
2404 615953 : if (lgefint(A) == 3) switch(A[2])
2405 : {
2406 832 : case 1: set_avma(av); return gen_1;
2407 42702 : case 2: base_is_2 = 1; break;
2408 : }
2409 :
2410 : /* TODO: Move this out of here and use for general modular computations */
2411 615121 : use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
2412 615122 : if (base_is_2)
2413 42702 : A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
2414 : else
2415 572420 : A = gen_powu_i(A, k, E, D.sqr, D.mul);
2416 615122 : if (use_montgomery)
2417 : {
2418 524119 : A = red_montgomery(A, N, ((struct montred *) E)->inv);
2419 524119 : if (cmpii(A, N) >= 0) A = subii(A,N);
2420 : }
2421 615122 : return gerepileuptoint(av, A);
2422 : }
2423 :
2424 : GEN
2425 29334 : Fp_pows(GEN A, long k, GEN N)
2426 : {
2427 29334 : if (lgefint(N) == 3) {
2428 14876 : ulong n = N[2];
2429 14876 : ulong a = umodiu(A, n);
2430 14876 : if (k < 0) {
2431 133 : a = Fl_inv(a, n);
2432 133 : k = -k;
2433 : }
2434 14876 : return utoi( Fl_powu(a, (ulong)k, n) );
2435 : }
2436 14458 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
2437 14458 : return Fp_powu(A, (ulong)k, N);
2438 : }
2439 :
2440 : /* A^K mod N */
2441 : GEN
2442 35865377 : Fp_pow(GEN A, GEN K, GEN N)
2443 : {
2444 : pari_sp av;
2445 35865377 : long s, lN = lgefint(N), sA, sy;
2446 : int base_is_2, use_montgomery;
2447 : GEN y;
2448 : muldata D;
2449 : void *E;
2450 :
2451 35865377 : s = signe(K);
2452 35865377 : if (!s) return dvdii(A,N)? gen_0: gen_1;
2453 34835341 : if (lN == 3 && lgefint(K) == 3)
2454 : {
2455 34072764 : ulong n = N[2], a = umodiu(A, n);
2456 34072956 : if (s < 0) a = Fl_inv(a, n);
2457 34073032 : if (a <= 1) return utoi(a); /* 0 or 1 */
2458 30559806 : return utoi(Fl_powu(a, uel(K,2), n));
2459 : }
2460 :
2461 762577 : av = avma;
2462 762577 : if (s < 0) y = Fp_inv(A,N);
2463 : else
2464 : {
2465 760654 : y = modii(A,N);
2466 760820 : if (!signe(y)) { set_avma(av); return gen_0; }
2467 : }
2468 762744 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
2469 :
2470 273376 : base_is_2 = 0;
2471 273376 : sy = abscmpii(y, shifti(N,-1)) > 0;
2472 273383 : if (sy) y = subii(N,y);
2473 273395 : sA = sy && mod2(K);
2474 273394 : if (lgefint(y) == 3) switch(y[2])
2475 : {
2476 207 : case 1: set_avma(av); return sA ? subis(N,1): gen_1;
2477 147248 : case 2: base_is_2 = 1; break;
2478 : }
2479 :
2480 : /* TODO: Move this out of here and use for general modular computations */
2481 273187 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
2482 273164 : if (base_is_2)
2483 147246 : y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
2484 : else
2485 125918 : y = gen_pow_i(y, K, E, D.sqr, D.mul);
2486 273200 : if (use_montgomery)
2487 : {
2488 171937 : y = red_montgomery(y, N, ((struct montred *) E)->inv);
2489 171939 : if (cmpii(y,N) >= 0) y = subii(y,N);
2490 : }
2491 273201 : if (sA) y = subii(N, y);
2492 273203 : return gerepileuptoint(av,y);
2493 : }
2494 :
2495 : static GEN
2496 14215874 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
2497 : static GEN
2498 8134365 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
2499 : static GEN
2500 57564 : _Fp_one(void *E) { (void) E; return gen_1; }
2501 :
2502 : GEN
2503 77 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
2504 77 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
2505 :
2506 : GEN
2507 54096 : Fp_pow_table(GEN R, GEN n, GEN p)
2508 54096 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
2509 :
2510 : GEN
2511 5924 : Fp_powers(GEN x, long n, GEN p)
2512 : {
2513 5924 : if (lgefint(p) == 3)
2514 2456 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
2515 3468 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
2516 : }
2517 :
2518 : GEN
2519 497 : FpV_prod(GEN V, GEN p) { return gen_product(V, (void *)p, &_Fp_mul); }
2520 :
2521 : static GEN
2522 27821000 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
2523 : static GEN
2524 168 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
2525 :
2526 : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
2527 : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
2528 : equalii,equali1,Fp_easylog};
2529 :
2530 : static GEN
2531 927374 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
2532 : static GEN
2533 1221969 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
2534 : static GEN
2535 1123231 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
2536 : static GEN
2537 609440 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
2538 : static GEN
2539 36647 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
2540 : static int
2541 265834 : _Fp_equal0(GEN x) { return signe(x)==0; }
2542 : static GEN
2543 33797 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
2544 :
2545 : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
2546 : _Fp_inv,_Fp_equal0,_Fp_s};
2547 :
2548 7870 : const struct bb_field *get_Fp_field(void **E, GEN p)
2549 7870 : { *E = (void*)p; return &Fp_field; }
2550 :
2551 : /*********************************************************************/
2552 : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
2553 : /*********************************************************************/
2554 : ulong
2555 481766 : Fl_order(ulong a, ulong o, ulong p)
2556 : {
2557 481766 : pari_sp av = avma;
2558 : GEN m, P, E;
2559 : long i;
2560 481766 : if (a==1) return 1;
2561 438866 : if (!o) o = p-1;
2562 438866 : m = factoru(o);
2563 438866 : P = gel(m,1);
2564 438866 : E = gel(m,2);
2565 1251386 : for (i = lg(P)-1; i; i--)
2566 : {
2567 812520 : ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
2568 812520 : if (y == 1) o = t;
2569 773077 : else for (j = 1; j < e; j++)
2570 : {
2571 385687 : y = Fl_powu(y, l, p);
2572 385687 : if (y == 1) { o = t * upowuu(l, j); break; }
2573 : }
2574 : }
2575 438866 : return gc_ulong(av, o);
2576 : }
2577 :
2578 : /*Find the exact order of a assuming a^o==1*/
2579 : GEN
2580 76053 : Fp_order(GEN a, GEN o, GEN p) {
2581 76053 : if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
2582 : {
2583 21 : ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
2584 21 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
2585 : }
2586 76032 : return gen_order(a, o, (void*)p, &Fp_star);
2587 : }
2588 : GEN
2589 70 : Fp_factored_order(GEN a, GEN o, GEN p)
2590 70 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
2591 :
2592 : /* return order of a mod p^e, e > 0, pe = p^e */
2593 : static GEN
2594 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
2595 : {
2596 : GEN ap, op;
2597 70 : if (absequaliu(p, 2))
2598 : {
2599 56 : if (e == 1) return gen_1;
2600 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
2601 49 : if (mod4(a) == 1) op = gen_1; else { op = gen_2; a = Fp_sqr(a, pe); }
2602 : } else {
2603 14 : ap = (e == 1)? a: remii(a,p);
2604 14 : op = Fp_order(ap, subiu(p,1), p);
2605 14 : if (e == 1) return op;
2606 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
2607 : }
2608 49 : if (equali1(a)) return op;
2609 7 : return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
2610 : }
2611 :
2612 : GEN
2613 63 : znorder(GEN x, GEN o)
2614 : {
2615 63 : pari_sp av = avma;
2616 : GEN b, a;
2617 :
2618 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
2619 56 : b = gel(x,1); a = gel(x,2);
2620 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
2621 49 : if (!o)
2622 : {
2623 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
2624 35 : long i, l = lg(P);
2625 35 : o = gen_1;
2626 70 : for (i = 1; i < l; i++)
2627 : {
2628 35 : GEN p = gel(P,i);
2629 35 : long e = itos(gel(E,i));
2630 :
2631 35 : if (l == 2)
2632 35 : o = Zp_order(a, p, e, b);
2633 : else {
2634 0 : GEN pe = powiu(p,e);
2635 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
2636 : }
2637 : }
2638 35 : return gerepileuptoint(av, o);
2639 : }
2640 14 : return Fp_order(a, o, b);
2641 : }
2642 :
2643 : /*********************************************************************/
2644 : /** DISCRETE LOGARITHM in (Z/nZ)* **/
2645 : /*********************************************************************/
2646 : static GEN
2647 59181 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
2648 : {
2649 59181 : pari_sp av = avma;
2650 : GEN h1, h2, F, G;
2651 59181 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
2652 35560 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
2653 : {
2654 126 : GEN M = cgetg(3, t_MAT);
2655 126 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
2656 126 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
2657 126 : return gerepileupto(av, M);
2658 : }
2659 35434 : return gc_NULL(av);
2660 : }
2661 :
2662 : static GEN
2663 59181 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
2664 : {
2665 : GEN rel;
2666 59181 : do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
2667 59181 : while (!rel);
2668 126 : return rel;
2669 : }
2670 :
2671 : struct Fp_log_rel
2672 : {
2673 : GEN rel;
2674 : ulong prmax;
2675 : long nbrel, nbmax, nbgen;
2676 : };
2677 :
2678 : /* add u^e */
2679 : static void
2680 3157 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
2681 : {
2682 3157 : pari_sp av = avma;
2683 3157 : long off = r->prmax+1;
2684 3157 : GEN F = cgetg(3, t_MAT);
2685 3157 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
2686 3157 : gel(F,2) = vecsmall_append(gel(z,2), e);
2687 3157 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
2688 3157 : }
2689 :
2690 : /* add u^-1 v^-1 */
2691 : static void
2692 104083 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
2693 : {
2694 104083 : pari_sp av = avma;
2695 104083 : long off = r->prmax+1;
2696 104083 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
2697 104083 : GEN F = cgetg(3, t_MAT);
2698 104083 : gel(F,1) = vecsmall_concat(gel(z,1), P);
2699 104083 : gel(F,2) = vecsmall_concat(gel(z,2), E);
2700 104083 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
2701 104083 : }
2702 :
2703 : /*
2704 : Let p=C^2+c
2705 : Solve h = (C+x)*(C+a)-p = 0 [mod l]
2706 : h= -c+x*(C+a)+C*a = 0 [mod l]
2707 : x = (c-C*a)/(C+a) [mod l]
2708 : h = -c+C*(x+a)+a*x
2709 : */
2710 :
2711 : GEN
2712 40828 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
2713 : {
2714 40828 : pari_sp ltop = avma;
2715 40828 : long i, j, th, n = lg(pi)-1, rel = 1;
2716 40828 : GEN sieve = zero_zv(a+2)+1;
2717 40862 : GEN L = cgetg(1+a+2, t_VEC);
2718 40832 : pari_sp av = avma;
2719 40832 : GEN z, h = addis(C,a);
2720 40822 : if ((z = Z_issmooth_fact(h, prmax)))
2721 : {
2722 3007 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
2723 3006 : av = avma;
2724 : }
2725 16960292 : for (i=1; i<=n; i++)
2726 : {
2727 16935304 : ulong li = pi[i], s = sz[i], al = a % li;
2728 16935304 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
2729 17462969 : if (!iv) continue;
2730 17037998 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
2731 78239604 : for(j = u; j<=a; j+=li) sieve[j] += s;
2732 : }
2733 35046 : if (a)
2734 : {
2735 40756 : long e = expi(mulis(C,a));
2736 40792 : th = e - expu(e) - 1;
2737 54 : } else th = -1;
2738 28025286 : for (j=0; j<a; j++)
2739 27983708 : if (sieve[j]>=th)
2740 : {
2741 119454 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
2742 119413 : if ((z = Z_issmooth_fact(h, prmax)))
2743 : {
2744 109710 : gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
2745 109996 : av = avma;
2746 9426 : } else set_avma(av);
2747 : }
2748 : /* j = a */
2749 41578 : if (sieve[a]>=th)
2750 : {
2751 476 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
2752 476 : if ((z = Z_issmooth_fact(h, prmax)))
2753 385 : gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
2754 : }
2755 41578 : setlg(L, rel); return gerepilecopy(ltop, L);
2756 : }
2757 :
2758 : static long
2759 63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
2760 : {
2761 : struct pari_mt pt;
2762 63 : long i, j, nb = 0;
2763 63 : GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
2764 : mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
2765 63 : long running, pending = 0;
2766 63 : GEN W = zerovec(r->nbgen);
2767 63 : mt_queue_start_lim(&pt, worker, r->nbgen);
2768 41229 : for (i = 0; (running = (i < r->nbgen)) || pending; i++)
2769 : {
2770 : GEN done;
2771 : long idx;
2772 41166 : mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
2773 41166 : done = mt_queue_get(&pt, &idx, &pending);
2774 41166 : if (!done || lg(done)==1) continue;
2775 35917 : gel(W, idx+1) = done;
2776 35917 : nb += lg(done)-1;
2777 35917 : if (DEBUGLEVEL && (i&127)==0)
2778 0 : err_printf("%ld%% ",100*nb/r->nbmax);
2779 : }
2780 63 : mt_queue_end(&pt);
2781 39550 : for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
2782 : {
2783 : long ll, m;
2784 39487 : GEN L = gel(W,j);
2785 39487 : if (isintzero(L)) continue;
2786 34531 : ll = lg(L);
2787 141771 : for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
2788 : {
2789 107240 : GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
2790 107240 : if (v[1] == 1)
2791 3157 : addifsmooth1(r, h, v[2], v[3]);
2792 : else
2793 104083 : addifsmooth2(r, h, v[2], v[3]);
2794 : }
2795 : }
2796 63 : return j;
2797 : }
2798 :
2799 : static GEN
2800 735 : ECP_psi(GEN x, GEN y)
2801 : {
2802 735 : long prec = realprec(x);
2803 735 : GEN lx = glog(x, prec), ly = glog(y, prec);
2804 735 : GEN u = gdiv(lx, ly);
2805 735 : return gpow(u, gneg(u),prec);
2806 : }
2807 :
2808 : struct computeG
2809 : {
2810 : GEN C;
2811 : long bnd, nbi;
2812 : };
2813 :
2814 : static GEN
2815 735 : _computeG(void *E, GEN gen)
2816 : {
2817 735 : struct computeG * d = (struct computeG *) E;
2818 735 : GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
2819 735 : return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
2820 : }
2821 :
2822 : static long
2823 63 : compute_nbgen(GEN C, long bnd, long nbi)
2824 : {
2825 : struct computeG d;
2826 63 : d.C = shifti(C, 1);
2827 63 : d.bnd = bnd;
2828 63 : d.nbi = nbi;
2829 63 : return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
2830 : }
2831 :
2832 : static GEN
2833 1714 : _psi(void*E, GEN y)
2834 : {
2835 1714 : GEN lx = (GEN) E;
2836 1714 : long prec = realprec(lx);
2837 1714 : GEN ly = glog(y, prec);
2838 1714 : GEN u = gdiv(lx, ly);
2839 1714 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
2840 : }
2841 :
2842 : static GEN
2843 63 : opt_param(GEN x, long prec)
2844 : {
2845 63 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
2846 : }
2847 :
2848 : static GEN
2849 63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
2850 : {
2851 63 : pari_sp av = avma;
2852 63 : long lM = lg(M)-1, nbcol = lM;
2853 63 : long tbs = maxss(1, expu(nbg/expi(m)));
2854 : for (;;)
2855 14 : {
2856 77 : GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
2857 : GEN tab;
2858 77 : long i, f=0;
2859 77 : long l = lg(K), lm = lgefint(m);
2860 77 : GEN idx = diviiexact(subiu(p,1),m), g;
2861 : pari_timer ti;
2862 77 : if (DEBUGLEVEL) timer_start(&ti);
2863 154 : for(i=1; i<l; i++)
2864 154 : if (signe(gel(K,i)))
2865 77 : break;
2866 77 : g = Fp_pow(utoi(i), idx, p);
2867 77 : tab = Fp_pow_init(g, p, tbs, p);
2868 77 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
2869 128464 : for(i=1; i<l; i++)
2870 : {
2871 128387 : GEN k = gel(K,i);
2872 128387 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
2873 128387 : if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
2874 76391 : gel(K,i) = cgetineg(lm);
2875 : else
2876 51996 : f++;
2877 : }
2878 77 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
2879 77 : if(f > (nbg>>1)) return gerepileupto(av, K);
2880 4585 : for(i=1; i<=nbcol; i++)
2881 : {
2882 4571 : long a = 1+random_Fl(lM);
2883 4571 : swap(gel(M,a),gel(M,i));
2884 : }
2885 14 : if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
2886 : }
2887 : }
2888 :
2889 : static GEN
2890 126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
2891 : {
2892 126 : pari_sp av=avma;
2893 126 : GEN aa = gen_1;
2894 126 : long AV = 0;
2895 : for(;;)
2896 0 : {
2897 126 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
2898 126 : GEN F = gel(A,1), E = gel(A,2);
2899 126 : GEN Ao = gen_0;
2900 126 : long i, l = lg(F);
2901 959 : for(i=1; i<l; i++)
2902 : {
2903 833 : GEN Ki = gel(K,F[i]);
2904 833 : if (signe(Ki)<0) break;
2905 833 : Ao = addii(Ao, mulis(Ki, E[i]));
2906 : }
2907 126 : if (i==l) return Fp_divu(Ao, AV, m);
2908 0 : aa = gerepileuptoint(av, aa);
2909 : }
2910 : }
2911 :
2912 : static GEN
2913 63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
2914 : {
2915 63 : pari_sp av = avma, av2;
2916 63 : long i, j, nbi, nbr = 0, nbrow, nbg;
2917 : GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
2918 : pari_timer ti;
2919 : struct Fp_log_rel r;
2920 63 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
2921 63 : ulong bnd = 4*bnds;
2922 63 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
2923 :
2924 63 : p_1 = subiu(p,1);
2925 63 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
2926 0 : m = diviiexact(p_1, Z_ppo(p_1, m));
2927 63 : pr = primes_upto_zv(bnd);
2928 63 : nbi = lg(pr)-1;
2929 63 : C = sqrtremi(p, &c);
2930 63 : av2 = avma;
2931 12796 : for (i = 1; i <= nbi; ++i)
2932 : {
2933 12733 : ulong lp = pr[i];
2934 26894 : while (lp <= bnd)
2935 : {
2936 14161 : nbr++;
2937 14161 : lp *= pr[i];
2938 : }
2939 : }
2940 63 : pi = cgetg(nbr+1,t_VECSMALL);
2941 63 : Ci = cgetg(nbr+1,t_VECSMALL);
2942 63 : ci = cgetg(nbr+1,t_VECSMALL);
2943 63 : sz = cgetg(nbr+1,t_VECSMALL);
2944 12796 : for (i = 1, j = 1; i <= nbi; ++i)
2945 : {
2946 12733 : ulong lp = pr[i], sp = expu(2*lp-1);
2947 26894 : while (lp <= bnd)
2948 : {
2949 14161 : pi[j] = lp;
2950 14161 : Ci[j] = umodiu(C, lp);
2951 14161 : ci[j] = umodiu(c, lp);
2952 14161 : sz[j] = sp;
2953 14161 : lp *= pr[i];
2954 14161 : j++;
2955 : }
2956 : }
2957 63 : r.nbrel = 0;
2958 63 : r.nbgen = compute_nbgen(C, bnd, nbi);
2959 63 : r.nbmax = 2*(nbi+r.nbgen);
2960 63 : r.rel = cgetg(r.nbmax+1,t_VEC);
2961 63 : r.prmax = pr[nbi];
2962 63 : if (DEBUGLEVEL)
2963 : {
2964 0 : err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
2965 0 : timer_start(&ti);
2966 : }
2967 63 : nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
2968 63 : nbrow = r.prmax + nbg;
2969 63 : if (DEBUGLEVEL)
2970 : {
2971 0 : err_printf("\n");
2972 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
2973 : }
2974 63 : setlg(r.rel,r.nbrel+1);
2975 63 : r.rel = gerepilecopy(av2, r.rel);
2976 63 : K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
2977 63 : if (DEBUGLEVEL) timer_start(&ti);
2978 63 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
2979 63 : if (DEBUGLEVEL) timer_printf(&ti," log element");
2980 63 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
2981 63 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
2982 63 : d = gcdii(Ao,Bo);
2983 63 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
2984 63 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
2985 63 : return gerepileuptoint(av, l);
2986 : }
2987 :
2988 : static int
2989 4651754 : Fp_log_use_index(long e, long p)
2990 : {
2991 4651754 : return (e >= 27 && 20*(p+6)<=e*e);
2992 : }
2993 :
2994 : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
2995 : static GEN
2996 8446387 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
2997 : {
2998 8446387 : pari_sp av = avma;
2999 8446387 : GEN p = (GEN)E;
3000 : /* assume a reduced mod p, p not necessarily prime */
3001 8446387 : if (equali1(a)) return gen_0;
3002 : /* p > 2 */
3003 5426738 : if (equalii(subiu(p,1), a)) /* -1 */
3004 : {
3005 : pari_sp av2;
3006 : GEN t;
3007 1317847 : ord = get_arith_Z(ord);
3008 1317847 : if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
3009 1317833 : t = shifti(ord,-1); /* only possible solution */
3010 1317833 : av2 = avma;
3011 1317833 : if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
3012 1317806 : set_avma(av2); return gerepileuptoint(av, t);
3013 : }
3014 4108893 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
3015 63 : return Fp_log_index(a, g, ord, p);
3016 4108830 : return gc_NULL(av); /* not easy */
3017 : }
3018 :
3019 : GEN
3020 3917285 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
3021 : {
3022 3917285 : GEN v = get_arith_ZZM(ord);
3023 3917261 : GEN F = gmael(v,2,1);
3024 3917261 : long lF = lg(F)-1, lmax;
3025 3917261 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
3026 3917233 : lmax = expi(gel(F,lF));
3027 3917235 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
3028 91 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
3029 3917220 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
3030 : }
3031 :
3032 : /* assume !(p & HIGHMASK) */
3033 : static ulong
3034 131789 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
3035 : {
3036 131789 : ulong i, h=1;
3037 361303 : for (i = 0; i < ord; i++, h = (h * g) % p)
3038 361303 : if (a==h) return i;
3039 0 : return ~0UL;
3040 : }
3041 :
3042 : static ulong
3043 24272 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
3044 : {
3045 24272 : ulong i, h=1;
3046 62211 : for (i = 0; i < ord; i++, h = Fl_mul_pre(h, g, p, pi))
3047 62211 : if (a==h) return i;
3048 0 : return ~0UL;
3049 : }
3050 :
3051 : static ulong
3052 0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
3053 : {
3054 0 : pari_sp av = avma;
3055 0 : GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
3056 0 : return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
3057 : }
3058 :
3059 : /* allow pi = 0 */
3060 : ulong
3061 24671 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
3062 : {
3063 24671 : if (!pi) return Fl_log(a, g, ord, p);
3064 24272 : if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
3065 0 : return Fl_log_Fp(a, g, ord, p);
3066 : }
3067 :
3068 : ulong
3069 131789 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
3070 : {
3071 131789 : if (ord <= 200)
3072 0 : return (p&HIGHMASK)? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
3073 131789 : : Fl_log_naive(a, g, ord, p);
3074 0 : return Fl_log_Fp(a, g, ord, p);
3075 : }
3076 :
3077 : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
3078 : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
3079 : static GEN
3080 126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
3081 : {
3082 126 : long l = lg(P) - 1, e = E[l];
3083 126 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
3084 : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
3085 :
3086 126 : if (l == 1) {
3087 98 : hpe = h;
3088 98 : gpe = g;
3089 : } else {
3090 28 : hpe = modii(h, pe);
3091 28 : gpe = modii(g, pe);
3092 : }
3093 126 : if (e == 1) {
3094 42 : hp = hpe;
3095 42 : gp = gpe;
3096 : } else {
3097 84 : hp = remii(hpe, p);
3098 84 : gp = remii(gpe, p);
3099 : }
3100 126 : if (hp == gen_0 || gp == gen_0) return NULL;
3101 105 : if (absequaliu(p, 2))
3102 : {
3103 35 : GEN n = int2n(e);
3104 35 : ogpe = Zp_order(gpe, gen_2, e, n);
3105 35 : a = Fp_log(hpe, gpe, ogpe, n);
3106 35 : if (typ(a) != t_INT) return NULL;
3107 : }
3108 : else
3109 : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
3110 : is trivial */
3111 : /* [order(gp), factor(order(gp))] */
3112 70 : GEN v = Fp_factored_order(gp, subiu(p,1), p);
3113 70 : GEN ogp = gel(v,1);
3114 70 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
3115 70 : a = Fp_log(hp, gp, v, p);
3116 70 : if (typ(a) != t_INT) return NULL;
3117 70 : if (e == 1) ogpe = ogp;
3118 : else
3119 : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
3120 : /* use p-adic log: O(log p + e) mul*/
3121 : long vpogpe, vpohpe;
3122 :
3123 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
3124 28 : gpe = Fp_pow(gpe, ogp, pe);
3125 : /* g,h = 1 mod p; compute b s.t. h = g^b */
3126 :
3127 : /* v_p(order g mod pe) */
3128 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
3129 : /* v_p(order h mod pe) */
3130 28 : vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
3131 28 : if (vpohpe > vpogpe) return NULL;
3132 :
3133 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
3134 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
3135 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
3136 28 : a = addii(a, mulii(ogp, padic_to_Q(b)));
3137 : }
3138 : }
3139 : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
3140 91 : if (l == 1) return a;
3141 :
3142 28 : N = diviiexact(N, pe); /* make N coprime to p */
3143 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
3144 28 : g = Fp_pow(g, modii(ogpe, phi), N);
3145 28 : setlg(P, l); /* remove last element */
3146 28 : setlg(E, l);
3147 28 : b = znlog_rec(h, g, N, P, E, PHI);
3148 28 : if (!b) return NULL;
3149 28 : return addmulii(a, b, ogpe);
3150 : }
3151 :
3152 : static GEN
3153 98 : get_PHI(GEN P, GEN E)
3154 : {
3155 98 : long i, l = lg(P);
3156 98 : GEN PHI = cgetg(l, t_VEC);
3157 98 : gel(PHI,1) = gen_1;
3158 126 : for (i=1; i<l-1; i++)
3159 : {
3160 28 : GEN t, p = gel(P,i);
3161 28 : long e = E[i];
3162 28 : t = mulii(powiu(p, e-1), subiu(p,1));
3163 28 : if (i > 1) t = mulii(t, gel(PHI,i));
3164 28 : gel(PHI,i+1) = t;
3165 : }
3166 98 : return PHI;
3167 : }
3168 :
3169 : GEN
3170 238 : znlog(GEN h, GEN g, GEN o)
3171 : {
3172 238 : pari_sp av = avma;
3173 : GEN N, fa, P, E, x;
3174 238 : switch (typ(g))
3175 : {
3176 28 : case t_PADIC:
3177 : {
3178 28 : GEN p = gel(g,2);
3179 28 : long v = valp(g);
3180 28 : if (v < 0) pari_err_DIM("znlog");
3181 28 : if (v > 0) {
3182 0 : long k = gvaluation(h, p);
3183 0 : if (k % v) return cgetg(1,t_VEC);
3184 0 : k /= v;
3185 0 : if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
3186 0 : return gc_stoi(av, k);
3187 : }
3188 28 : N = gel(g,3);
3189 28 : g = Rg_to_Fp(g, N);
3190 28 : break;
3191 : }
3192 203 : case t_INTMOD:
3193 203 : N = gel(g,1);
3194 203 : g = gel(g,2); break;
3195 7 : default: pari_err_TYPE("znlog", g);
3196 : return NULL; /* LCOV_EXCL_LINE */
3197 : }
3198 231 : if (equali1(N)) { set_avma(av); return gen_0; }
3199 231 : h = Rg_to_Fp(h, N);
3200 224 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
3201 98 : fa = Z_factor(N);
3202 98 : P = gel(fa,1);
3203 98 : E = vec_to_vecsmall(gel(fa,2));
3204 98 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
3205 98 : if (!x) { set_avma(av); return cgetg(1,t_VEC); }
3206 63 : return gerepileuptoint(av, x);
3207 : }
3208 :
3209 : GEN
3210 119140 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
3211 : {
3212 119140 : if (lgefint(p)==3)
3213 : {
3214 118706 : long nn = itos_or_0(n);
3215 118706 : if (nn)
3216 : {
3217 118706 : ulong pp = p[2];
3218 : ulong uz;
3219 118706 : ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
3220 118685 : if (r==ULONG_MAX) return NULL;
3221 118629 : if (zeta) *zeta = utoi(uz);
3222 118629 : return utoi(r);
3223 : }
3224 : }
3225 434 : a = modii(a,p);
3226 434 : if (!signe(a))
3227 : {
3228 0 : if (zeta) *zeta = gen_1;
3229 0 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
3230 0 : return gen_0;
3231 : }
3232 434 : if (absequaliu(n,2))
3233 : {
3234 238 : if (zeta) *zeta = subiu(p,1);
3235 238 : return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
3236 : }
3237 196 : return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
3238 : }
3239 :
3240 : /*********************************************************************/
3241 : /** FACTORIAL **/
3242 : /*********************************************************************/
3243 : GEN
3244 90338 : mulu_interval_step(ulong a, ulong b, ulong step)
3245 : {
3246 90338 : pari_sp av = avma;
3247 : ulong k, l, N, n;
3248 : long lx;
3249 : GEN x;
3250 :
3251 90338 : if (!a) return gen_0;
3252 90338 : if (step == 1) return mulu_interval(a, b);
3253 90338 : n = 1 + (b-a) / step;
3254 90338 : b -= (b-a) % step;
3255 90338 : if (n < 61)
3256 : {
3257 88884 : if (n == 1) return utoipos(a);
3258 68397 : x = muluu(a,a+step); if (n == 2) return x;
3259 537327 : for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
3260 53444 : return gerepileuptoint(av, x);
3261 : }
3262 : /* step | b-a */
3263 1454 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3264 1449 : N = b + a;
3265 1449 : for (k = a;; k += step)
3266 : {
3267 230236 : l = N - k; if (l <= k) break;
3268 228787 : gel(x,lx++) = muluu(k,l);
3269 : }
3270 1449 : if (l == k) gel(x,lx++) = utoipos(k);
3271 1449 : setlg(x, lx);
3272 1449 : return gerepileuptoint(av, ZV_prod(x));
3273 : }
3274 : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
3275 : * first is slower ... ] */
3276 : GEN
3277 158205 : mulu_interval(ulong a, ulong b)
3278 : {
3279 158205 : pari_sp av = avma;
3280 : ulong k, l, N, n;
3281 : long lx;
3282 : GEN x;
3283 :
3284 158205 : if (!a) return gen_0;
3285 158205 : n = b - a + 1;
3286 158205 : if (n < 61)
3287 : {
3288 157470 : if (n == 1) return utoipos(a);
3289 107147 : x = muluu(a,a+1); if (n == 2) return x;
3290 402521 : for (k=a+2; k<b; k++) x = mului(k,x);
3291 : /* avoid k <= b: broken if b = ULONG_MAX */
3292 93104 : return gerepileuptoint(av, mului(b,x));
3293 : }
3294 735 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3295 735 : N = b + a;
3296 735 : for (k = a;; k++)
3297 : {
3298 28361 : l = N - k; if (l <= k) break;
3299 27630 : gel(x,lx++) = muluu(k,l);
3300 : }
3301 731 : if (l == k) gel(x,lx++) = utoipos(k);
3302 731 : setlg(x, lx);
3303 732 : return gerepileuptoint(av, ZV_prod(x));
3304 : }
3305 : GEN
3306 588 : muls_interval(long a, long b)
3307 : {
3308 588 : pari_sp av = avma;
3309 588 : long lx, k, l, N, n = b - a + 1;
3310 : GEN x;
3311 :
3312 588 : if (a <= 0 && b >= 0) return gen_0;
3313 315 : if (n < 61)
3314 : {
3315 315 : x = stoi(a);
3316 553 : for (k=a+1; k<=b; k++) x = mulsi(k,x);
3317 315 : return gerepileuptoint(av, x);
3318 : }
3319 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3320 0 : N = b + a;
3321 0 : for (k = a;; k++)
3322 : {
3323 0 : l = N - k; if (l <= k) break;
3324 0 : gel(x,lx++) = mulss(k,l);
3325 : }
3326 0 : if (l == k) gel(x,lx++) = stoi(k);
3327 0 : setlg(x, lx);
3328 0 : return gerepileuptoint(av, ZV_prod(x));
3329 : }
3330 :
3331 : GEN
3332 105 : mpprimorial(long n)
3333 : {
3334 105 : pari_sp av = avma;
3335 105 : if (n <= 12) switch(n)
3336 : {
3337 14 : case 0: case 1: return gen_1;
3338 7 : case 2: return gen_2;
3339 14 : case 3: case 4: return utoipos(6);
3340 14 : case 5: case 6: return utoipos(30);
3341 28 : case 7: case 8: case 9: case 10: return utoipos(210);
3342 14 : case 11: case 12: return utoipos(2310);
3343 7 : default: pari_err_DOMAIN("primorial", "argument","<",gen_0,stoi(n));
3344 : }
3345 7 : return gerepileuptoint(av, zv_prod_Z(primes_upto_zv(n)));
3346 : }
3347 :
3348 : GEN
3349 503726 : mpfact(long n)
3350 : {
3351 503726 : pari_sp av = avma;
3352 : GEN a, v;
3353 : long k;
3354 503726 : if (n <= 12) switch(n)
3355 : {
3356 423129 : case 0: case 1: return gen_1;
3357 36943 : case 2: return gen_2;
3358 3409 : case 3: return utoipos(6);
3359 4257 : case 4: return utoipos(24);
3360 2873 : case 5: return utoipos(120);
3361 2528 : case 6: return utoipos(720);
3362 2434 : case 7: return utoipos(5040);
3363 2423 : case 8: return utoipos(40320);
3364 2444 : case 9: return utoipos(362880);
3365 2708 : case 10:return utoipos(3628800);
3366 1437 : case 11:return utoipos(39916800);
3367 591 : case 12:return utoipos(479001600);
3368 0 : default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
3369 : }
3370 18550 : v = cgetg(expu(n) + 2, t_VEC);
3371 18560 : for (k = 1;; k++)
3372 86525 : {
3373 105085 : long m = n >> (k-1), l;
3374 105085 : if (m <= 2) break;
3375 86537 : l = (1 + (n >> k)) | 1;
3376 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3377 86537 : a = mulu_interval_step(l, m, 2);
3378 86399 : gel(v,k) = k == 1? a: powiu(a, k);
3379 : }
3380 86481 : a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
3381 18530 : a = shifti(a, factorial_lval(n, 2));
3382 18528 : return gerepileuptoint(av, a);
3383 : }
3384 :
3385 : ulong
3386 56718 : factorial_Fl(long n, ulong p)
3387 : {
3388 : long k;
3389 : ulong v;
3390 56718 : if (p <= (ulong)n) return 0;
3391 56718 : v = Fl_powu(2, factorial_lval(n, 2), p);
3392 56771 : for (k = 1;; k++)
3393 142366 : {
3394 199137 : long m = n >> (k-1), l, i;
3395 199137 : ulong a = 1;
3396 199137 : if (m <= 2) break;
3397 142390 : l = (1 + (n >> k)) | 1;
3398 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3399 779803 : for (i=l; i<=m; i+=2)
3400 637434 : a = Fl_mul(a, i, p);
3401 142369 : v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
3402 : }
3403 56747 : return v;
3404 : }
3405 :
3406 : GEN
3407 158 : factorial_Fp(long n, GEN p)
3408 : {
3409 158 : pari_sp av = avma;
3410 : long k;
3411 158 : GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
3412 158 : for (k = 1;; k++)
3413 344 : {
3414 502 : long m = n >> (k-1), l, i;
3415 502 : GEN a = gen_1;
3416 502 : if (m <= 2) break;
3417 344 : l = (1 + (n >> k)) | 1;
3418 : /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
3419 850 : for (i=l; i<=m; i+=2)
3420 506 : a = Fp_mulu(a, i, p);
3421 344 : v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
3422 344 : v = gerepileuptoint(av, v);
3423 : }
3424 158 : return v;
3425 : }
3426 :
3427 : /*******************************************************************/
3428 : /** LUCAS & FIBONACCI **/
3429 : /*******************************************************************/
3430 : static void
3431 56 : lucas(ulong n, GEN *a, GEN *b)
3432 : {
3433 : GEN z, t, zt;
3434 56 : if (!n) { *a = gen_2; *b = gen_1; return; }
3435 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
3436 49 : switch(n & 3) {
3437 14 : case 0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
3438 14 : case 1: *a = subiu(zt,1); *b = addiu(sqri(t),2); break;
3439 7 : case 2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
3440 14 : case 3: *a = addiu(zt,1); *b = subiu(sqri(t),2);
3441 : }
3442 49 : }
3443 :
3444 : GEN
3445 7 : fibo(long n)
3446 : {
3447 7 : pari_sp av = avma;
3448 : GEN a, b;
3449 7 : if (!n) return gen_0;
3450 7 : lucas((ulong)(labs(n)-1), &a, &b);
3451 7 : a = diviuexact(addii(shifti(a,1),b), 5);
3452 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
3453 7 : return gerepileuptoint(av, a);
3454 : }
3455 :
3456 : /*******************************************************************/
3457 : /* CONTINUED FRACTIONS */
3458 : /*******************************************************************/
3459 : static GEN
3460 2830744 : icopy_lg(GEN x, long l)
3461 : {
3462 2830744 : long lx = lgefint(x);
3463 : GEN y;
3464 :
3465 2830744 : if (lx >= l) return icopy(x);
3466 49 : y = cgeti(l); affii(x, y); return y;
3467 : }
3468 :
3469 : /* continued fraction of a/b. If y != NULL, stop when partial quotients
3470 : * differ from y */
3471 : static GEN
3472 2831094 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
3473 : {
3474 : GEN z, c;
3475 2831094 : ulong i, l, ly = lgefint(b);
3476 :
3477 : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
3478 2831094 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
3479 2831094 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
3480 2831094 : if (l > LGBITS) l = LGBITS;
3481 :
3482 2831094 : z = cgetg(l,t_VEC);
3483 2831094 : l--;
3484 2831094 : if (y) {
3485 350 : pari_sp av = avma;
3486 350 : if (l >= (ulong)lg(y)) l = lg(y)-1;
3487 25209 : for (i = 1; i <= l; i++)
3488 : {
3489 24985 : GEN q = gel(y,i);
3490 24985 : gel(z,i) = q;
3491 24985 : c = b; if (!gequal1(q)) c = mulii(q, b);
3492 24985 : c = subii(a, c);
3493 24985 : if (signe(c) < 0)
3494 : { /* partial quotient too large */
3495 96 : c = addii(c, b);
3496 96 : if (signe(c) >= 0) i++; /* by 1 */
3497 96 : break;
3498 : }
3499 24889 : if (cmpii(c, b) >= 0)
3500 : { /* partial quotient too small */
3501 30 : c = subii(c, b);
3502 30 : if (cmpii(c, b) < 0) {
3503 : /* by 1. If next quotient is 1 in y, add 1 */
3504 12 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
3505 12 : i++;
3506 : }
3507 30 : break;
3508 : }
3509 24859 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
3510 24859 : a = b; b = c;
3511 : }
3512 : } else {
3513 2830744 : a = icopy_lg(a, ly);
3514 2830744 : b = icopy(b);
3515 23447857 : for (i = 1; i <= l; i++)
3516 : {
3517 23447539 : gel(z,i) = truedvmdii(a,b,&c);
3518 23447539 : if (c == gen_0) { i++; break; }
3519 20617113 : affii(c, a); cgiv(c); c = a;
3520 20617113 : a = b; b = c;
3521 : }
3522 : }
3523 2831094 : i--;
3524 2831094 : if (i > 1 && gequal1(gel(z,i)))
3525 : {
3526 101 : cgiv(gel(z,i)); --i;
3527 101 : gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
3528 : }
3529 2831094 : setlg(z,i+1); return z;
3530 : }
3531 :
3532 : static GEN
3533 0 : sersfcont(GEN a, GEN b, long k)
3534 : {
3535 0 : long i, l = typ(a) == t_POL? lg(a): 3;
3536 : GEN y, c;
3537 0 : if (lg(b) > l) l = lg(b);
3538 0 : if (k > 0 && l > k+1) l = k+1;
3539 0 : y = cgetg(l,t_VEC);
3540 0 : for (i=1; i<l; i++)
3541 : {
3542 0 : gel(y,i) = poldivrem(a,b,&c);
3543 0 : if (gequal0(c)) { i++; break; }
3544 0 : a = b; b = c;
3545 : }
3546 0 : setlg(y, i); return y;
3547 : }
3548 :
3549 : GEN
3550 2831759 : gboundcf(GEN x, long k)
3551 : {
3552 : pari_sp av;
3553 2831759 : long tx = typ(x), e;
3554 : GEN y, a, b, c;
3555 :
3556 2831759 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
3557 2831752 : if (is_scalar_t(tx))
3558 : {
3559 2831752 : if (gequal0(x)) return mkvec(gen_0);
3560 2831647 : switch(tx)
3561 : {
3562 896 : case t_INT: return mkveccopy(x);
3563 357 : case t_REAL:
3564 357 : av = avma;
3565 357 : c = mantissa_real(x,&e);
3566 357 : if (e < 0) pari_err_PREC("gboundcf");
3567 350 : y = int2n(e);
3568 350 : a = Qsfcont(c,y, NULL, k);
3569 350 : b = addsi(signe(x), c);
3570 350 : return gerepilecopy(av, Qsfcont(b,y, a, k));
3571 :
3572 2830394 : case t_FRAC:
3573 2830394 : av = avma;
3574 2830394 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
3575 : }
3576 0 : pari_err_TYPE("gboundcf",x);
3577 : }
3578 :
3579 0 : switch(tx)
3580 : {
3581 0 : case t_POL: return mkveccopy(x);
3582 0 : case t_SER:
3583 0 : av = avma;
3584 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
3585 0 : case t_RFRAC:
3586 0 : av = avma;
3587 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
3588 : }
3589 0 : pari_err_TYPE("gboundcf",x);
3590 : return NULL; /* LCOV_EXCL_LINE */
3591 : }
3592 :
3593 : static GEN
3594 14 : sfcont2(GEN b, GEN x, long k)
3595 : {
3596 14 : pari_sp av = avma;
3597 14 : long lb = lg(b), tx = typ(x), i;
3598 : GEN y,p1;
3599 :
3600 14 : if (k)
3601 : {
3602 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
3603 0 : lb = k+1;
3604 : }
3605 7 : y = cgetg(lb,t_VEC);
3606 7 : if (lb==1) return y;
3607 7 : if (is_scalar_t(tx))
3608 : {
3609 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
3610 : }
3611 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
3612 :
3613 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
3614 7 : for (i = 1;;)
3615 : {
3616 35 : if (tx == t_REAL)
3617 : {
3618 35 : long e = expo(x);
3619 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
3620 35 : gel(y,i) = floorr(x);
3621 35 : p1 = subri(x, gel(y,i));
3622 : }
3623 : else
3624 : {
3625 0 : gel(y,i) = gfloor(x);
3626 0 : p1 = gsub(x, gel(y,i));
3627 : }
3628 35 : if (++i >= lb) break;
3629 28 : if (gequal0(p1)) break;
3630 28 : x = gdiv(gel(b,i),p1);
3631 : }
3632 7 : setlg(y,i);
3633 7 : return gerepilecopy(av,y);
3634 : }
3635 :
3636 : GEN
3637 112 : gcf(GEN x) { return gboundcf(x,0); }
3638 : GEN
3639 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
3640 : GEN
3641 49 : contfrac0(GEN x, GEN b, long nmax)
3642 : {
3643 : long tb;
3644 :
3645 49 : if (!b) return gboundcf(x,nmax);
3646 28 : tb = typ(b);
3647 28 : if (tb == t_INT) return gboundcf(x,itos(b));
3648 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
3649 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
3650 14 : return sfcont2(b,x,nmax);
3651 : }
3652 :
3653 : GEN
3654 252 : contfracpnqn(GEN x, long n)
3655 : {
3656 252 : pari_sp av = avma;
3657 252 : long i, lx = lg(x);
3658 : GEN M,A,B, p0,p1, q0,q1;
3659 :
3660 252 : if (lx == 1)
3661 : {
3662 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
3663 21 : if (n >= 0) return cgetg(1,t_MAT);
3664 7 : return matid(2);
3665 : }
3666 224 : switch(typ(x))
3667 : {
3668 182 : case t_VEC: case t_COL: A = x; B = NULL; break;
3669 42 : case t_MAT:
3670 42 : switch(lgcols(x))
3671 : {
3672 0 : case 2: A = row(x,1); B = NULL; break;
3673 35 : case 3: A = row(x,2); B = row(x,1); break;
3674 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
3675 : return NULL; /*LCOV_EXCL_LINE*/
3676 : }
3677 35 : break;
3678 0 : default: pari_err_TYPE("pnqn",x);
3679 : return NULL; /*LCOV_EXCL_LINE*/
3680 : }
3681 217 : p1 = gel(A,1);
3682 217 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
3683 217 : if (n >= 0)
3684 : {
3685 182 : lx = minss(lx, n+2);
3686 182 : if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
3687 : }
3688 35 : else if (lx == 2)
3689 7 : return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
3690 : /* lx >= 3 */
3691 119 : p0 = gen_1;
3692 119 : q0 = gen_0; /* p[-1], q[-1] */
3693 119 : M = cgetg(lx, t_MAT);
3694 119 : gel(M,1) = mkcol2(p1,q1);
3695 399 : for (i=2; i<lx; i++)
3696 : {
3697 280 : GEN a = gel(A,i), p2,q2;
3698 280 : if (B) {
3699 84 : GEN b = gel(B,i);
3700 84 : p0 = gmul(b,p0);
3701 84 : q0 = gmul(b,q0);
3702 : }
3703 280 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
3704 280 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
3705 280 : gel(M,i) = mkcol2(p1,q1);
3706 : }
3707 119 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
3708 119 : return gerepilecopy(av, M);
3709 : }
3710 : GEN
3711 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
3712 : /* x = [a0, ..., an] from gboundcf, n >= 0;
3713 : * return [[p0, ..., pn], [q0,...,qn]] */
3714 : GEN
3715 609308 : ZV_allpnqn(GEN x)
3716 : {
3717 609308 : long i, lx = lg(x);
3718 609308 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
3719 :
3720 609308 : gel(v,1) = P = cgetg(lx, t_VEC);
3721 609308 : gel(v,2) = Q = cgetg(lx, t_VEC);
3722 609308 : p0 = gen_1; q0 = gen_0;
3723 609308 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
3724 2092209 : for (i=2; i<lx; i++)
3725 : {
3726 1482901 : GEN a = gel(x,i);
3727 1482901 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
3728 1482901 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
3729 : }
3730 609308 : return v;
3731 : }
3732 :
3733 : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
3734 : static GEN
3735 42 : mod_to_frac(GEN x, GEN N, GEN B)
3736 : {
3737 : GEN a, b, A;
3738 42 : if (B) A = divii(shifti(N, -1), B);
3739 : else
3740 : {
3741 14 : A = sqrti(shifti(N, -1));
3742 14 : B = A;
3743 : }
3744 42 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
3745 28 : return equali1(b)? a: mkfrac(a,b);
3746 : }
3747 :
3748 : static GEN
3749 70 : mod_to_rfrac(GEN x, GEN N, long B)
3750 : {
3751 : GEN a, b;
3752 70 : long A, d = degpol(N);
3753 70 : if (B >= 0) A = d-1 - B;
3754 : else
3755 : {
3756 42 : B = d >> 1;
3757 42 : A = odd(d)? B : B-1;
3758 : }
3759 70 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
3760 70 : if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
3761 56 : return gdiv(a,b);
3762 : }
3763 :
3764 : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
3765 : * of the continued fraction of x with b <= k maximal */
3766 : static GEN
3767 7 : bestappr_frac(GEN x, GEN k)
3768 : {
3769 : pari_sp av;
3770 : GEN p0, p1, p, q0, q1, q, a, y;
3771 :
3772 7 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
3773 0 : av = avma; y = x;
3774 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
3775 0 : q1 = gen_0; q0 = gen_1;
3776 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
3777 : for(;;)
3778 : {
3779 0 : x = ginv(x); /* > 1 */
3780 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
3781 0 : if (cmpii(a,k) > 0)
3782 : { /* next partial quotient will overflow limits */
3783 : GEN n, d;
3784 0 : a = divii(subii(k, q1), q0);
3785 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3786 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3787 : /* compare |y-p0/q0|, |y-p1/q1| */
3788 0 : n = gel(y,1);
3789 0 : d = gel(y,2);
3790 0 : if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
3791 : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
3792 0 : { p1 = p0; q1 = q0; }
3793 0 : break;
3794 : }
3795 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3796 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3797 :
3798 0 : if (cmpii(q0,k) > 0) break;
3799 0 : x = gsub(x,a); /* 0 <= x < 1 */
3800 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
3801 :
3802 : }
3803 0 : return gerepileupto(av, gdiv(p1,q1));
3804 : }
3805 : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
3806 : * of the continued fraction of x with b <= k maximal */
3807 : static GEN
3808 1228551 : bestappr_real(GEN x, GEN k)
3809 : {
3810 1228551 : pari_sp av = avma;
3811 1228551 : GEN kr, p0, p1, p, q0, q1, q, a, y = x;
3812 :
3813 1228551 : p1 = gen_1; a = p0 = floorr(x);
3814 1228448 : q1 = gen_0; q0 = gen_1;
3815 1228448 : x = subri(x,a); /* 0 <= x < 1 */
3816 1228462 : if (!signe(x)) { cgiv(x); return a; }
3817 1114764 : kr = itor(k, realprec(x));
3818 : for(;;)
3819 1198322 : {
3820 : long d;
3821 2313130 : x = invr(x); /* > 1 */
3822 2313022 : if (cmprr(x,kr) > 0)
3823 : { /* next partial quotient will overflow limits */
3824 1099891 : a = divii(subii(k, q1), q0);
3825 1099839 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3826 1099849 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3827 : /* compare |y-p0/q0|, |y-p1/q1| */
3828 1099797 : if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
3829 : mulir(q0, subri(mulir(q1,y), p1))) < 0)
3830 131678 : { p1 = p0; q1 = q0; }
3831 1099922 : break;
3832 : }
3833 1213194 : d = nbits2prec(expo(x) + 1);
3834 1213193 : if (d > realprec(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
3835 :
3836 1213003 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
3837 1212956 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3838 1212984 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3839 :
3840 1212976 : if (cmpii(q0,k) > 0) break;
3841 1206771 : x = subri(x,a); /* 0 <= x < 1 */
3842 1206773 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
3843 : }
3844 1114780 : if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
3845 1114780 : return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
3846 : }
3847 :
3848 : /* k t_INT or NULL */
3849 : static GEN
3850 2143702 : bestappr_Q(GEN x, GEN k)
3851 : {
3852 2143702 : long lx, tx = typ(x), i;
3853 : GEN a, y;
3854 :
3855 2143702 : switch(tx)
3856 : {
3857 154 : case t_INT: return icopy(x);
3858 7 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
3859 1461322 : case t_REAL:
3860 1461322 : if (!signe(x)) return gen_0;
3861 : /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
3862 1228542 : i = bit_prec(x); if (i <= expo(x)) return NULL;
3863 1228553 : return bestappr_real(x, k? k: int2n(i));
3864 :
3865 28 : case t_INTMOD: {
3866 28 : pari_sp av = avma;
3867 28 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
3868 21 : return gerepilecopy(av, a);
3869 : }
3870 14 : case t_PADIC: {
3871 14 : pari_sp av = avma;
3872 14 : long v = valp(x);
3873 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
3874 7 : if (v) a = gmul(a, powis(gel(x,2), v));
3875 7 : return gerepilecopy(av, a);
3876 : }
3877 :
3878 294 : case t_COMPLEX: {
3879 294 : pari_sp av = avma;
3880 294 : y = cgetg(3, t_COMPLEX);
3881 294 : gel(y,2) = bestappr(gel(x,2), k);
3882 294 : gel(y,1) = bestappr(gel(x,1), k);
3883 294 : if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
3884 91 : return y;
3885 : }
3886 0 : case t_SER:
3887 0 : if (ser_isexactzero(x)) return gcopy(x);
3888 : /* fall through */
3889 : case t_POLMOD: case t_POL: case t_RFRAC:
3890 : case t_VEC: case t_COL: case t_MAT:
3891 681883 : y = cgetg_copy(x, &lx);
3892 681920 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
3893 2758794 : for (; i<lx; i++)
3894 : {
3895 2076883 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
3896 2076874 : gel(y,i) = a;
3897 : }
3898 681911 : if (tx == t_POL) return normalizepol(y);
3899 681897 : if (tx == t_SER) return normalizeser(y);
3900 681897 : return y;
3901 : }
3902 0 : pari_err_TYPE("bestappr_Q",x);
3903 : return NULL; /* LCOV_EXCL_LINE */
3904 : }
3905 :
3906 : static GEN
3907 56 : bestappr_ser(GEN x, long B)
3908 : {
3909 56 : long dN, v = valser(x), lx = lg(x);
3910 : GEN t;
3911 56 : x = normalizepol(ser2pol_i(x, lx));
3912 56 : dN = lx-2;
3913 56 : if (v > 0)
3914 : {
3915 14 : x = RgX_shift_shallow(x, v);
3916 14 : dN += v;
3917 : }
3918 42 : else if (v < 0)
3919 : {
3920 7 : if (B >= 0) B = maxss(B+v, 0);
3921 : }
3922 56 : t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
3923 56 : if (!t) return NULL;
3924 42 : if (v < 0)
3925 : {
3926 : GEN a, b;
3927 : long vx;
3928 7 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
3929 : /* t_RFRAC */
3930 7 : vx = varn(x);
3931 7 : a = gel(t,1);
3932 7 : b = gel(t,2);
3933 7 : v -= RgX_valrem(b, &b);
3934 7 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
3935 7 : if (v < 0) b = RgX_shift(b, -v);
3936 0 : else if (v > 0) {
3937 0 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
3938 0 : a = RgX_shift(a, v);
3939 : }
3940 7 : t = mkrfraccopy(a, b);
3941 : }
3942 42 : return t;
3943 : }
3944 : static GEN bestappr_RgX(GEN x, long B);
3945 : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
3946 : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
3947 : static GEN
3948 77 : bestappr_RgX(GEN x, long B)
3949 : {
3950 77 : long i, lx, tx = typ(x);
3951 : GEN y, t;
3952 77 : switch(tx)
3953 : {
3954 0 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
3955 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
3956 0 : return gcopy(x);
3957 :
3958 14 : case t_RFRAC: {
3959 14 : pari_sp av = avma;
3960 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
3961 7 : x = rfrac_to_ser_i(x, 2*B+1);
3962 7 : t = bestappr_ser(x, B); if (!t) return NULL;
3963 0 : return gerepileupto(av, t);
3964 : }
3965 14 : case t_POLMOD: {
3966 14 : pari_sp av = avma;
3967 14 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
3968 14 : return gerepileupto(av, t);
3969 : }
3970 49 : case t_SER: {
3971 49 : pari_sp av = avma;
3972 49 : t = bestappr_ser(x, B); if (!t) return NULL;
3973 42 : return gerepileupto(av, t);
3974 : }
3975 :
3976 0 : case t_VEC: case t_COL: case t_MAT:
3977 0 : y = cgetg_copy(x, &lx);
3978 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
3979 0 : for (; i<lx; i++)
3980 : {
3981 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
3982 0 : gel(y,i) = t;
3983 : }
3984 0 : return y;
3985 : }
3986 0 : pari_err_TYPE("bestappr_RgX",x);
3987 : return NULL; /* LCOV_EXCL_LINE */
3988 : }
3989 :
3990 : /* allow k = NULL: maximal accuracy */
3991 : GEN
3992 66804 : bestappr(GEN x, GEN k)
3993 : {
3994 66804 : pari_sp av = avma;
3995 66804 : if (k) { /* replace by floor(k) */
3996 66482 : switch(typ(k))
3997 : {
3998 4368 : case t_INT:
3999 4368 : break;
4000 62114 : case t_REAL: case t_FRAC:
4001 62114 : k = floor_safe(k); /* left on stack for efficiency */
4002 62114 : if (!signe(k)) k = gen_1;
4003 62114 : break;
4004 0 : default:
4005 0 : pari_err_TYPE("bestappr [bound type]", k);
4006 0 : break;
4007 : }
4008 322 : }
4009 66804 : x = bestappr_Q(x, k);
4010 66803 : if (!x) { set_avma(av); return cgetg(1,t_VEC); }
4011 66789 : return x;
4012 : }
4013 : GEN
4014 77 : bestapprPade(GEN x, long B)
4015 : {
4016 77 : pari_sp av = avma;
4017 77 : GEN t = bestappr_RgX(x, B);
4018 77 : if (!t) { set_avma(av); return cgetg(1,t_VEC); }
4019 63 : return t;
4020 : }
|
