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. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 :
14 : /*********************************************************************/
15 : /** **/
16 : /** ARITHMETIC FUNCTIONS **/
17 : /** (first part) **/
18 : /** **/
19 : /*********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : /******************************************************************/
24 : /* */
25 : /* GENERATOR of (Z/mZ)* */
26 : /* */
27 : /******************************************************************/
28 : static GEN
29 44 : remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
30 : static ulong
31 22466 : u_remove2(ulong q) { return q >> vals(q); }
32 : GEN
33 44 : odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
34 : static GEN
35 22466 : u_odd_prime_divisors(ulong q) { return gel(factoru(u_remove2(q)), 1); }
36 : /* p odd prime, q=(p-1)/2; L0 list of (some) divisors of q = (p-1)/2 or NULL
37 : * (all prime divisors of q); return the q/l, l in L0 */
38 : static GEN
39 35 : is_gener_expo(GEN p, GEN L0)
40 : {
41 35 : GEN L, q = shifti(p,-1);
42 : long i, l;
43 35 : if (L0) {
44 14 : l = lg(L0);
45 14 : L = cgetg(l, t_VEC);
46 : } else {
47 21 : L0 = L = odd_prime_divisors(q);
48 21 : l = lg(L);
49 : }
50 35 : for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
51 35 : return L;
52 : }
53 : static GEN
54 22195 : u_is_gener_expo(ulong p, GEN L0)
55 : {
56 22195 : const ulong q = p >> 1;
57 : long i, l;
58 : GEN L;
59 22195 : if (L0) {
60 2858 : l = lg(L0);
61 2858 : L = cgetg(l, t_VECSMALL);
62 : } else {
63 19337 : L0 = L = u_odd_prime_divisors(q);
64 19337 : l = lg(L);
65 : }
66 22195 : for (i=1; i<l; i++) L[i] = q / uel(L0,i);
67 22195 : return L;
68 : }
69 :
70 : int
71 53268 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
72 : {
73 : long i;
74 53268 : if (krouu(x, p) >= 0) return 0;
75 67274 : for (i=lg(L)-1; i; i--)
76 : {
77 43527 : ulong t = Fl_powu(x, uel(L,i), p);
78 43527 : if (t == p_1 || t == 1) return 0;
79 : }
80 23747 : return 1;
81 : }
82 : /* assume p prime */
83 : ulong
84 53381 : pgener_Fl_local(ulong p, GEN L0)
85 : {
86 53381 : const pari_sp av = avma;
87 53381 : const ulong p_1 = p-1;
88 : long x;
89 : GEN L;
90 53381 : if (p <= 19) switch(p)
91 : { /* quick trivial cases */
92 42 : case 2: return 1;
93 : case 7:
94 4195 : case 17: return 3;
95 26985 : default: return 2;
96 : }
97 22159 : L = u_is_gener_expo(p,L0);
98 22159 : for (x=2;;x++) { if (is_gener_Fl(x,p,p_1,L)) break; }
99 22159 : avma = av; return x;
100 : }
101 : ulong
102 43692 : pgener_Fl(ulong p) { return pgener_Fl_local(p, NULL); }
103 :
104 : /* L[i] = set of (p-1)/2l, l ODD prime divisor of p-1 (l=2 can be included,
105 : * but wasteful) */
106 : int
107 136 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
108 : {
109 136 : long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
110 136 : if (t >= 0) return 0;
111 368 : for (i = lg(L)-1; i; i--)
112 : {
113 247 : GEN t = Fp_pow(x, gel(L,i), p);
114 247 : if (equalii(t, p_1) || equali1(t)) return 0;
115 : }
116 121 : return 1;
117 : }
118 :
119 : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
120 : GEN
121 7728 : pgener_Fp_local(GEN p, GEN L0)
122 : {
123 7728 : pari_sp av0 = avma;
124 : GEN x, p_1, L;
125 7728 : if (lgefint(p) == 3)
126 : {
127 : ulong z;
128 7699 : if (p[2] == 2) return gen_1;
129 4920 : if (L0) L0 = ZV_to_nv(L0);
130 4920 : z = pgener_Fl_local(uel(p,2), L0);
131 4920 : avma = av0; return utoipos(z);
132 : }
133 29 : p_1 = subis(p,1); L = is_gener_expo(p, L0);
134 29 : x = utoipos(2);
135 32 : for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
136 29 : avma = av0; return utoipos(uel(x,2));
137 : }
138 :
139 : GEN
140 7203 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
141 :
142 : ulong
143 3865 : pgener_Zl(ulong p)
144 : {
145 3865 : if (p == 2) pari_err_DOMAIN("pgener_Zl","p","=",gen_2,gen_2);
146 : /* only p < 2^32 such that znprimroot(p) != znprimroot(p^2) */
147 3865 : if (p == 40487) return 10;
148 : #ifndef LONG_IS_64BIT
149 547 : return pgener_Fl(p);
150 : #else
151 3318 : if (p < (1UL<<32)) return pgener_Fl(p);
152 : else
153 : {
154 36 : const pari_sp av = avma;
155 36 : const ulong p_1 = p-1;
156 : long x ;
157 36 : GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
158 108 : for (x=2;;x++)
159 108 : if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2))) break;
160 72 : avma = av; return x;
161 : }
162 : #endif
163 : }
164 :
165 : /* p prime. Return a primitive root modulo p^e, e > 1 */
166 : GEN
167 3871 : pgener_Zp(GEN p)
168 : {
169 3871 : if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
170 : else
171 : {
172 6 : const pari_sp av = avma;
173 6 : GEN p_1 = subis(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
174 6 : GEN x = utoipos(2);
175 12 : for (;; x[2]++)
176 18 : if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
177 12 : avma = av; return utoipos(uel(x,2));
178 : }
179 : }
180 :
181 : static GEN
182 287 : gener_Zp(GEN q, GEN F)
183 : {
184 287 : GEN p = NULL;
185 287 : long e = 0;
186 287 : if (F)
187 : {
188 14 : GEN P = gel(F,1), E = gel(F,2);
189 14 : long i, l = lg(P);
190 42 : for (i = 1; i < l; i++)
191 : {
192 28 : p = gel(P,i);
193 28 : if (absequaliu(p, 2)) continue;
194 14 : if (i < l-1) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
195 14 : e = itos(gel(E,i));
196 : }
197 14 : if (!p) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
198 : }
199 : else
200 273 : e = Z_isanypower(q, &p);
201 287 : return e > 1? pgener_Zp(p): pgener_Fp(q);
202 : }
203 :
204 : GEN
205 357 : znprimroot(GEN N)
206 : {
207 357 : pari_sp av = avma;
208 : GEN x, n, F;
209 :
210 357 : if ((F = check_arith_non0(N,"znprimroot")))
211 : {
212 14 : F = clean_Z_factor(F);
213 14 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
214 : }
215 350 : if (signe(N) < 0) N = absi(N);
216 350 : if (abscmpiu(N, 4) <= 0) { avma = av; return mkintmodu(N[2]-1,N[2]); }
217 301 : switch(mod4(N))
218 : {
219 : case 0: /* N = 0 mod 4 */
220 14 : pari_err_DOMAIN("znprimroot", "argument","=",N,N);
221 0 : x = NULL; break;
222 : case 2: /* N = 2 mod 4 */
223 21 : n = shifti(N,-1); /* becomes odd */
224 21 : x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
225 21 : break;
226 : default: /* N odd */
227 266 : x = gener_Zp(N,F);
228 266 : break;
229 : }
230 287 : return gerepilecopy(av, mkintmod(x, N));
231 : }
232 :
233 : /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
234 : GEN
235 0 : rootsof1_Fp(GEN n, GEN p)
236 : {
237 0 : pari_sp av = avma;
238 0 : GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
239 0 : GEN z = pgener_Fp_local(p, L);
240 0 : z = Fp_pow(z, diviiexact(subis(p,1), n), p); /* prim. n-th root of 1 */
241 0 : return gerepileuptoint(av, z);
242 : }
243 :
244 : GEN
245 301 : rootsof1u_Fp(ulong n, GEN p)
246 : {
247 301 : pari_sp av = avma;
248 301 : GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
249 301 : z = pgener_Fp_local(p, Flv_to_ZV(L));
250 301 : z = Fp_pow(z, diviuexact(subis(p,1), n), p); /* prim. n-th root of 1 */
251 301 : return gerepileuptoint(av, z);
252 : }
253 :
254 : ulong
255 2828 : rootsof1_Fl(ulong n, ulong p)
256 : {
257 2828 : pari_sp av = avma;
258 2828 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
259 2828 : ulong z = pgener_Fl_local(p, L);
260 2828 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
261 2828 : avma = av; return z;
262 : }
263 :
264 : /*********************************************************************/
265 : /** **/
266 : /** INVERSE TOTIENT FUNCTION **/
267 : /** **/
268 : /*********************************************************************/
269 : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
270 : * primes. Return factor(N) */
271 : GEN
272 350651 : Z_factor_listP(GEN N, GEN L)
273 : {
274 350651 : long i, k, l = lg(L);
275 350651 : GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
276 1346688 : for (i = k = 1; i < l; i++)
277 : {
278 996037 : GEN p = gel(L,i);
279 996037 : long v = Z_pvalrem(N, p, &N);
280 996037 : if (v)
281 : {
282 792176 : gel(P,k) = p;
283 792176 : gel(E,k) = utoipos(v);
284 792176 : k++;
285 : }
286 : }
287 350651 : setlg(P, k);
288 350651 : setlg(E, k); return mkmat2(P,E);
289 : }
290 :
291 : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
292 : * L is a list of primes containing all prime divisors of n. */
293 : static long
294 621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
295 : {
296 621565 : pari_sp av = avma, av2;
297 : GEN k, D;
298 : long i, v;
299 621565 : if (m && mod2(n))
300 : {
301 270914 : if (!equali1(n)) return 0;
302 69986 : if (px) *px = gen_1;
303 69986 : return 1;
304 : }
305 350651 : D = divisors(Z_factor_listP(shifti(n, -1), L));
306 : /* loop through primes p > m, d = p-1 | n */
307 350651 : av2 = avma;
308 350651 : if (!m)
309 : { /* special case p = 2, d = 1 */
310 69986 : k = n;
311 69986 : for (v = 1;; v++) {
312 69986 : if (istotient_i(k, gen_2, L, px)) {
313 69986 : if (px) *px = shifti(*px, v);
314 69986 : return 1;
315 : }
316 0 : if (mod2(k)) break;
317 0 : k = shifti(k,-1);
318 0 : }
319 0 : avma = av2;
320 : }
321 1099462 : for (i = 1; i < lg(D); ++i)
322 : {
323 1001588 : GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
324 1001588 : if (m && cmpii(d, m) < 0) continue;
325 677782 : p = addiu(d, 1);
326 677782 : if (!isprime(p)) continue;
327 442064 : k = diviiexact(n, d);
328 481593 : for (v = 1;; v++) {
329 : GEN r;
330 481593 : if (istotient_i(k, p, L, px)) {
331 182791 : if (px) *px = mulii(*px, powiu(p, v));
332 182791 : return 1;
333 : }
334 298802 : k = dvmdii(k, p, &r);
335 298802 : if (r != gen_0) break;
336 39529 : }
337 259273 : avma = av2;
338 : }
339 97874 : avma = av; return 0;
340 : }
341 :
342 : /* find x such that phi(x) = n */
343 : long
344 70000 : istotient(GEN n, GEN *px)
345 : {
346 70000 : pari_sp av = avma;
347 70000 : if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
348 70000 : if (signe(n) < 1) return 0;
349 70000 : if (mod2(n))
350 : {
351 14 : if (!equali1(n)) return 0;
352 14 : if (px) *px = gen_1;
353 14 : return 1;
354 : }
355 69986 : if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
356 : {
357 69986 : if (!px) avma = av;
358 : else
359 69986 : *px = gerepileuptoint(av, *px);
360 69986 : return 1;
361 : }
362 0 : avma = av; return 0;
363 : }
364 :
365 : /*********************************************************************/
366 : /** **/
367 : /** INTEGRAL LOGARITHM **/
368 : /** **/
369 : /*********************************************************************/
370 :
371 : /* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
372 : * e = floor(log_y B). Set *ptq = y^e if non-NULL */
373 : long
374 443472 : logintall(GEN B, GEN y, GEN *ptq)
375 : {
376 : pari_sp av;
377 443472 : long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
378 : GEN q, pow2;
379 :
380 443472 : if (equaliu(y,2))
381 : {
382 3430 : if (ptq) *ptq = int2n(eB);
383 3430 : return eB;
384 : }
385 440042 : av = avma;
386 440042 : ey = expi(y);
387 : /* eB/(ey+1) - 1 < e <= eB/ey */
388 440042 : emax = eB/ey;
389 440042 : if (emax <= 13) /* e small, be naive */
390 : {
391 419056 : GEN r = y, r2 = gen_1;
392 1035260 : for (e=1;; e++)
393 : { /* here, r = y^e, r2 = y^(e-1) */
394 1035260 : long fl = cmpii(r, B);
395 1035260 : if (fl >= 0)
396 : {
397 419056 : if (fl) { e--; cgiv(r); r = r2; }
398 419056 : if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
399 419056 : return e;
400 : }
401 616204 : r2 = r; r = mulii(r,y);
402 616204 : }
403 : }
404 : /* e >= 13 ey / (ey+1) >= 6.5 */
405 :
406 : /* binary splitting: compute bits of e one by one */
407 : /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
408 20986 : pow2 = new_chunk((long)log2(eB)+2);
409 20986 : gel(pow2,0) = y;
410 20986 : for (i=0, q=y;; )
411 : {
412 109572 : GEN r = gel(pow2,i); /* r = y^2^i */
413 109572 : long fl = cmpii(r,B);
414 109572 : if (!fl)
415 : {
416 0 : e = 1L<<i;
417 0 : if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
418 0 : return e;
419 : }
420 109572 : if (fl > 0) { i--; break; }
421 92347 : q = r;
422 92347 : if (1L<<(i+1) > emax) break;
423 88586 : gel(pow2,++i) = sqri(q);
424 88586 : }
425 :
426 20986 : for (e = 1L<<i;;)
427 : { /* y^e = q < B < r = q * y^(2^i) */
428 92326 : pari_sp av2 = avma;
429 : long fl;
430 : GEN r;
431 92326 : if (--i < 0) break;
432 71347 : r = mulii(q, gel(pow2,i));
433 71347 : fl = cmpii(r, B);
434 71347 : if (fl > 0) avma = av2;
435 : else
436 : {
437 45791 : e += (1L<<i);
438 45791 : q = r;
439 45791 : if (!fl) break; /* B = r */
440 : }
441 71340 : }
442 20986 : if (ptq) *ptq = gerepileuptoint(av, q); else avma = av;
443 20986 : return e;
444 : }
445 :
446 : long
447 42 : logint0(GEN B, GEN y, GEN *ptq)
448 : {
449 42 : if (typ(B) != t_INT) pari_err_TYPE("logint",B);
450 42 : if (signe(B) <= 0) pari_err_DOMAIN("logint", "x" ,"<=", gen_0, B);
451 42 : if (typ(y) != t_INT) pari_err_TYPE("logint",y);
452 42 : if (cmpis(y, 2) < 0) pari_err_DOMAIN("logint", "b" ,"<=", gen_1, y);
453 42 : return logintall(B,y,ptq);
454 : }
455 :
456 : /*********************************************************************/
457 : /** **/
458 : /** INTEGRAL SQUARE ROOT **/
459 : /** **/
460 : /*********************************************************************/
461 : GEN
462 38431 : sqrtint(GEN a)
463 : {
464 38431 : if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
465 38431 : switch (signe(a))
466 : {
467 38417 : case 1: return sqrti(a);
468 7 : case 0: return gen_0;
469 7 : default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
470 : }
471 0 : return NULL; /* not reached */
472 : }
473 :
474 : /*********************************************************************/
475 : /** **/
476 : /** PERFECT SQUARE **/
477 : /** **/
478 : /*********************************************************************/
479 : static int
480 8762681 : carremod(ulong A)
481 : {
482 8762681 : const int carresmod64[]={
483 : 1,1,0,0,1,0,0,0,0,1, 0,0,0,0,0,0,1,1,0,0, 0,0,0,0,0,1,0,0,0,0,
484 : 0,0,0,1,0,0,1,0,0,0, 0,1,0,0,0,0,0,0,0,1, 0,0,0,0,0,0,0,1,0,0, 0,0,0,0};
485 8762681 : const int carresmod63[]={
486 : 1,1,0,0,1,0,0,1,0,1, 0,0,0,0,0,0,1,0,1,0, 0,0,1,0,0,1,0,0,1,0,
487 : 0,0,0,0,0,0,1,1,0,0, 0,0,0,1,0,0,1,0,0,1, 0,0,0,0,0,0,0,0,1,0, 0,0,0};
488 8762681 : const int carresmod65[]={
489 : 1,1,0,0,1,0,0,0,0,1, 1,0,0,0,1,0,1,0,0,0, 0,0,0,0,0,1,1,0,0,1,
490 : 1,0,0,0,0,1,1,0,0,1, 1,0,0,0,0,0,0,0,0,1, 0,1,0,0,0,1,1,0,0,0, 0,1,0,0,1};
491 8762681 : const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
492 17525362 : return (carresmod64[A & 0x3fUL]
493 3227154 : && carresmod63[A % 63UL]
494 1938360 : && carresmod65[A % 65UL]
495 10350123 : && carresmod11[A % 11UL]);
496 : }
497 :
498 : /* emulate Z_issquareall on single-word integers */
499 : long
500 8596364 : uissquareall(ulong A, ulong *sqrtA)
501 : {
502 8596364 : if (!A) { *sqrtA = 0; return 1; }
503 8596364 : if (carremod(A))
504 : {
505 1496940 : ulong a = usqrt(A);
506 1496936 : if (a * a == A) { *sqrtA = a; return 1; }
507 : }
508 7189223 : return 0;
509 : }
510 : long
511 115576 : uissquare(ulong A)
512 : {
513 115576 : if (!A) return 1;
514 115576 : if (carremod(A))
515 : {
516 2335 : ulong a = usqrt(A);
517 2335 : if (a * a == A) return 1;
518 : }
519 113265 : return 0;
520 : }
521 :
522 : long
523 4606346 : Z_issquareall(GEN x, GEN *pt)
524 : {
525 : pari_sp av;
526 : GEN y, r;
527 :
528 4606346 : switch(signe(x))
529 : {
530 1856461 : case -1: return 0;
531 686 : case 0: if (pt) *pt=gen_0; return 1;
532 : }
533 2749199 : if (lgefint(x) == 3)
534 : {
535 2698599 : ulong u = uel(x,2), a;
536 2698599 : if (!pt) return uissquare(u);
537 2583023 : if (!uissquareall(u, &a)) return 0;
538 1376101 : *pt = utoipos(a); return 1;
539 : }
540 50600 : if (!carremod(umodiu(x, 64*63*65*11))) return 0;
541 6327 : av = avma; y = sqrtremi(x, &r);
542 6327 : if (r != gen_0) { avma = av; return 0; }
543 5740 : if (pt) { *pt = y; avma = (pari_sp)y; } else avma = av;
544 5740 : return 1;
545 : }
546 :
547 : /* a t_INT, p prime */
548 : long
549 0 : Zp_issquare(GEN a, GEN p)
550 : {
551 : long v;
552 : GEN ap;
553 :
554 0 : if (!signe(a) || gequal1(a)) return 1;
555 0 : v = Z_pvalrem(a, p, &ap);
556 0 : if (v&1) return 0;
557 0 : return absequaliu(p, 2)? umodiu(ap, 8) == 1
558 0 : : kronecker(ap,p) == 1;
559 : }
560 :
561 : static long
562 1540 : polissquareall(GEN x, GEN *pt)
563 : {
564 : pari_sp av;
565 : long v;
566 : GEN y, a, b, p;
567 :
568 1540 : if (!signe(x))
569 : {
570 7 : if (pt) *pt = gcopy(x);
571 7 : return 1;
572 : }
573 1533 : if (odd(degpol(x))) return 0; /* odd degree */
574 987 : av = avma;
575 987 : v = RgX_valrem(x, &x);
576 987 : if (v & 1) { avma = av; return 0; }
577 980 : a = gel(x,2); /* test constant coeff */
578 980 : if (!pt)
579 63 : { if (!issquare(a)) { avma = av; return 0; } }
580 : else
581 917 : { if (!issquareall(a,&b)) { avma = av; return 0; } }
582 980 : if (!degpol(x)) { /* constant polynomial */
583 77 : if (!pt) { avma = av; return 1; }
584 35 : y = scalarpol(b, varn(x)); goto END;
585 : }
586 903 : p = characteristic(x);
587 903 : if (signe(p) && !mod2(p))
588 : {
589 : long i, lx;
590 35 : if (!absequaliu(p,2)) pari_err_IMPL("issquare for even characteristic != 2");
591 28 : x = gmul(x, mkintmod(gen_1, gen_2));
592 28 : lx = lg(x);
593 28 : if ((lx-3) & 1) { avma = av; return 0; }
594 49 : for (i = 3; i < lx; i+=2)
595 28 : if (!gequal0(gel(x,i))) { avma = av; return 0; }
596 21 : if (pt) {
597 14 : y = cgetg((lx+3) / 2, t_POL);
598 49 : for (i = 2; i < lx; i+=2)
599 35 : if (!issquareall(gel(x,i), &gel(y, (i+2)>>1))) { avma = av; return 0; }
600 14 : y[1] = evalsigne(1) | evalvarn(varn(x));
601 14 : goto END;
602 : } else {
603 21 : for (i = 2; i < lx; i+=2)
604 14 : if (!issquare(gel(x,i))) { avma = av; return 0; }
605 7 : avma = av; return 1;
606 : }
607 : }
608 : else
609 : {
610 868 : long m = 1;
611 868 : x = RgX_Rg_div(x,a);
612 : /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
613 868 : if (!signe(p)) x = RgX_deflate_max(x,&m);
614 868 : y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
615 1127 : if (!RgX_equal(RgX_sqr(y), x)) { avma = av; return 0; }
616 609 : if (!pt) { avma = av; return 1; }
617 609 : if (!gequal1(a)) y = gmul(b, y);
618 609 : if (m != 1) y = RgX_inflate(y,m);
619 : }
620 : END:
621 658 : if (v) y = RgX_shift_shallow(y, v>>1);
622 658 : *pt = gerepilecopy(av, y); return 1;
623 : }
624 :
625 : /* b unit mod p */
626 : static int
627 294 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
628 : {
629 294 : if (d == 1)
630 : { /* mod p: faster */
631 210 : if (!Fp_ispower(b, K, p)) return 0;
632 210 : if (pt) *pt = Fp_sqrtn(b, K, p, NULL);
633 : }
634 : else
635 : { /* mod p^{2 +} */
636 84 : if (!ispower(cvtop(b, p, d), K, pt)) return 0;
637 63 : if (pt) *pt = gtrunc(*pt);
638 : }
639 273 : return 1;
640 : }
641 :
642 : /* We're studying whether a mod (q*p^e) is a K-th power, (q,p) = 1.
643 : * Decide mod p^e, then reduce a mod q unless q = NULL. */
644 : static int
645 434 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
646 : {
647 : GEN t, A;
648 434 : long v = Z_pvalrem(*pa, p, &A), d = e - v;
649 434 : if (d <= 0) t = gen_0;
650 : else
651 : {
652 : ulong r;
653 378 : v = udivui_rem(v, K, &r);
654 378 : if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
655 273 : if (L && v) t = mulii(t, powiu(p, v));
656 : }
657 329 : if (q) *pa = modii(*pa, q);
658 329 : if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
659 329 : return 1;
660 : }
661 : long
662 322 : Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
663 : {
664 : GEN L, N;
665 : pari_sp av;
666 : long e, i, l;
667 : ulong pp;
668 : forprime_t S;
669 :
670 322 : if (!signe(a))
671 : {
672 7 : if (pt) {
673 7 : GEN t = cgetg(3, t_INTMOD);
674 7 : gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
675 : }
676 7 : return 1;
677 : }
678 : /* a != 0 */
679 315 : av = avma;
680 :
681 315 : if (typ(q) != t_INT) /* integer factorization */
682 : {
683 0 : GEN P = gel(q,1), E = gel(q,2);
684 0 : l = lg(P);
685 0 : L = pt? vectrunc_init(l): NULL;
686 0 : for (i = 1; i < l; i++)
687 : {
688 0 : GEN p = gel(P,i);
689 0 : long e = itos(gel(E,i));
690 0 : if (!handle_pe(&a, NULL, L, K, p, e)) { avma = av; return 0; }
691 : }
692 0 : goto END;
693 : }
694 315 : if (!mod2(K)
695 189 : && kronecker(a, shifti(q,-vali(q))) == -1) { avma = av; return 0; }
696 308 : L = pt? vectrunc_init(expi(q)+1): NULL;
697 308 : u_forprime_init(&S, 2, tridiv_bound(q));
698 308 : while ((pp = u_forprime_next(&S)))
699 : {
700 : int stop;
701 883204 : e = Z_lvalrem_stop(&q, pp, &stop);
702 883204 : if (!e) continue;
703 210 : if (!handle_pe(&a, q, L, K, utoipos(pp), e)) { avma = av; return 0; }
704 168 : if (stop)
705 : {
706 133 : if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
707 133 : goto END;
708 : }
709 : }
710 154 : l = lg(primetab);
711 154 : for (i = 1; i < l; i++)
712 : {
713 0 : GEN p = gel(primetab,i);
714 0 : e = Z_pvalrem(q, p, &q);
715 0 : if (!e) continue;
716 0 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
717 0 : if (is_pm1(q)) goto END;
718 : }
719 154 : N = gcdii(a,q);
720 154 : if (!is_pm1(N))
721 : {
722 112 : if (ifac_isprime(N))
723 : {
724 70 : e = Z_pvalrem(q, N, &q);
725 70 : if (!handle_pe(&a, q, L, K, N, e)) { avma = av; return 0; }
726 : }
727 : else
728 : {
729 42 : GEN part = ifac_start(N, 0);
730 : for(;;)
731 : {
732 : long e;
733 : GEN p;
734 84 : if (!ifac_next(&part, &p, &e)) break;
735 42 : e = Z_pvalrem(q, p, &q);
736 42 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
737 42 : }
738 : }
739 : }
740 84 : if (!is_pm1(q))
741 : {
742 84 : if (ifac_isprime(q))
743 : {
744 28 : if (!handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
745 : }
746 : else
747 : {
748 56 : GEN part = ifac_start(q, 0);
749 : for(;;)
750 : {
751 : long e;
752 : GEN p;
753 140 : if (!ifac_next(&part, &p, &e)) break;
754 98 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
755 84 : }
756 : }
757 : }
758 : END:
759 203 : if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
760 203 : return 1;
761 : }
762 :
763 : static long
764 126 : polmodispower(GEN x, GEN K, GEN *pt)
765 : {
766 126 : pari_sp av = avma;
767 126 : GEN p = NULL, T = NULL;
768 126 : if (Rg_is_FpXQ(x, &T,&p) && p)
769 : {
770 112 : x = liftall_shallow(x);
771 112 : if (!Fq_ispower(x, K, T, p)) { avma = av; return 0; }
772 105 : if (!pt) { avma = av; return 1; }
773 98 : x = Fq_sqrtn(x, K, T,p, NULL);
774 98 : if (typ(x) == t_INT)
775 7 : x = Fp_to_mod(x,p);
776 : else
777 91 : x = mkpolmod(FpX_to_mod(x,p), FpX_to_mod(T,p));
778 98 : *pt = gerepilecopy(av, x); return 1;
779 : }
780 14 : pari_err_IMPL("ispower for general t_POLMOD");
781 0 : return 0;
782 : }
783 :
784 : long
785 146111 : issquareall(GEN x, GEN *pt)
786 : {
787 146111 : long tx = typ(x);
788 : GEN F;
789 : pari_sp av;
790 :
791 146111 : if (!pt) return issquare(x);
792 3276 : switch(tx)
793 : {
794 1386 : case t_INT: return Z_issquareall(x, pt);
795 238 : case t_FRAC: av = avma;
796 238 : F = cgetg(3, t_FRAC);
797 238 : if ( !Z_issquareall(gel(x,1), &gel(F,1))
798 105 : || !Z_issquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
799 105 : *pt = F; return 1;
800 :
801 : case t_POLMOD:
802 21 : return polmodispower(x, gen_2, pt);
803 1463 : case t_POL: return polissquareall(x,pt);
804 7 : case t_RFRAC: av = avma;
805 7 : F = cgetg(3, t_RFRAC);
806 7 : if ( !issquareall(gel(x,1), &gel(F,1))
807 7 : || !polissquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
808 7 : *pt = F; return 1;
809 :
810 : case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
811 63 : if (!issquare(x)) return 0;
812 63 : *pt = gsqrt(x, DEFAULTPREC); return 1;
813 :
814 : case t_INTMOD:
815 63 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
816 :
817 35 : case t_FFELT: return FF_issquareall(x, pt);
818 :
819 : }
820 0 : pari_err_TYPE("issquareall",x);
821 0 : return 0; /* not reached */
822 : }
823 :
824 : long
825 143108 : issquare(GEN x)
826 : {
827 : pari_sp av;
828 : GEN a, p;
829 : long i, v;
830 :
831 143108 : switch(typ(x))
832 : {
833 : case t_INT:
834 142716 : return Z_issquare(x);
835 :
836 : case t_REAL:
837 14 : return (signe(x)>=0);
838 :
839 : case t_INTMOD:
840 70 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, NULL);
841 :
842 : case t_FRAC:
843 21 : return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
844 :
845 7 : case t_FFELT: return FF_issquareall(x, NULL);
846 :
847 : case t_COMPLEX:
848 63 : return 1;
849 :
850 : case t_PADIC:
851 105 : a = gel(x,4); if (!signe(a)) return 1;
852 105 : if (valp(x)&1) return 0;
853 91 : p = gel(x,2);
854 91 : if (!absequaliu(p, 2)) return (kronecker(a,p) != -1);
855 :
856 42 : v = precp(x); /* here p=2, a is odd */
857 42 : if ((v>=3 && mod8(a) != 1 ) ||
858 0 : (v==2 && mod4(a) != 1)) return 0;
859 21 : return 1;
860 :
861 : case t_POLMOD:
862 14 : return polmodispower(x, gen_2, NULL);
863 :
864 : case t_POL:
865 70 : return polissquareall(x,NULL);
866 :
867 : case t_SER:
868 21 : if (!signe(x)) return 1;
869 14 : if (valp(x)&1) return 0;
870 7 : return issquare(gel(x,2));
871 :
872 : case t_RFRAC:
873 7 : av = avma; i = issquare(gmul(gel(x,1),gel(x,2)));
874 7 : avma = av; return i;
875 : }
876 0 : pari_err_TYPE("issquare",x);
877 0 : return 0; /* not reached */
878 : }
879 0 : GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
880 0 : GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
881 :
882 : long
883 1386 : ispolygonal(GEN x, GEN S, GEN *N)
884 : {
885 1386 : pari_sp av = avma;
886 : GEN D, d, n;
887 1386 : if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
888 1386 : if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
889 1386 : if (abscmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
890 1386 : if (signe(x) < 0) return 0;
891 1386 : if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
892 1260 : if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
893 : /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
894 1134 : if (abscmpiu(S, 1<<16) < 0) /* common case ! */
895 : {
896 441 : ulong s = S[2], r;
897 504 : if (s == 4) return Z_issquareall(x, N);
898 378 : if (s == 3)
899 0 : D = addiu(shifti(x, 3), 1);
900 : else
901 378 : D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
902 378 : if (!Z_issquareall(D, &d)) { avma = av; return 0; }
903 378 : if (s == 3)
904 0 : d = subiu(d, 1);
905 : else
906 378 : d = addiu(d, s - 4);
907 378 : n = diviu_rem(d, 2*s - 4, &r);
908 378 : if (r) { avma = av; return 0; }
909 : }
910 : else
911 : {
912 693 : GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
913 693 : D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
914 693 : if (!Z_issquareall(D, &d)) { avma = av; return 0; }
915 693 : d = addii(d, S_4);
916 693 : n = dvmdii(shifti(d,-1), S_2, &r);
917 693 : if (r != gen_0) { avma = av; return 0; }
918 : }
919 1071 : if (N) *N = gerepileuptoint(av, n); else avma = av;
920 1071 : return 1;
921 : }
922 :
923 : /*********************************************************************/
924 : /** **/
925 : /** PERFECT POWER **/
926 : /** **/
927 : /*********************************************************************/
928 : static long
929 497 : polispower(GEN x, GEN K, GEN *pt)
930 : {
931 : pari_sp av;
932 497 : long v, d, k = itos(K);
933 : GEN y, a, b;
934 :
935 497 : if (!signe(x))
936 : {
937 7 : if (pt) *pt = gcopy(x);
938 7 : return 1;
939 : }
940 490 : if (degpol(x) % k) return 0; /* degree not multiple of k */
941 483 : av = avma;
942 483 : y = NULL; /*-Wall*/
943 483 : v = RgX_valrem(x, &x);
944 483 : if (v % k) return 0;
945 476 : v /= k;
946 476 : a = gel(x,2); b = NULL;
947 476 : if (!ispower(a, K, &b)) { avma = av; return 0; }
948 462 : d = degpol(x);
949 462 : if (d)
950 : {
951 322 : GEN p = characteristic(x);
952 322 : a = leading_coeff(x);
953 350 : if (!ispower(a, K, &b)) { avma = av; return 0; }
954 322 : x = RgX_normalize(x);
955 322 : if (signe(p))
956 : {
957 168 : GEN T0, T = NULL;
958 168 : if (!BPSW_isprime(p))
959 0 : pari_err_IMPL("ispower in non-prime characteristic");
960 168 : if (RgX_is_FpXQX(x,&T,&p))
961 : { /* over Fq */
962 168 : T0 = T;
963 168 : if (T && typ(T) == t_FFELT) T = FF_mod(T);
964 168 : x = RgX_to_FqX(x,T,p);
965 168 : if (!FqX_ispower(x, k, T,p, pt)) { avma = av; return 0; }
966 140 : if (pt)
967 : {
968 140 : y = *pt;
969 140 : if (!T) y = FpX_to_mod(y, p);
970 105 : else if (typ(T0) == t_FFELT)
971 70 : y = FqX_to_FFX(y, T0);
972 : else
973 : {
974 35 : T = FpX_to_mod(T, p);
975 35 : y = gmul(y, gmodulsg(1,T));
976 : }
977 : }
978 140 : goto END;
979 : }
980 0 : if (cmpii(p,K) <= 0)
981 0 : pari_err_IMPL("ispower(general t_POL) in small characteristic");
982 : }
983 154 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
984 154 : if (!RgX_equal(powgi(y, K), x)) { avma = av; return 0; }
985 : }
986 : else
987 140 : y = pol_1(varn(x));
988 : END:
989 434 : if (pt)
990 : {
991 434 : if (!gequal1(a))
992 : {
993 14 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
994 14 : y = gmul(b,y);
995 : }
996 434 : if (v) y = RgX_shift_shallow(y, v);
997 434 : *pt = gerepilecopy(av, y);
998 : }
999 0 : else avma = av;
1000 434 : return 1;
1001 : }
1002 :
1003 : long
1004 1498 : Z_ispowerall(GEN x, ulong k, GEN *pt)
1005 : {
1006 1498 : long s = signe(x);
1007 : ulong mask;
1008 1498 : if (!s) { if (pt) *pt = gen_0; return 1; }
1009 1498 : if (s > 0) {
1010 1344 : if (k == 2) return Z_issquareall(x, pt);
1011 973 : if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
1012 203 : if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
1013 189 : if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
1014 182 : return is_kth_power(x, k, pt);
1015 : }
1016 154 : if (!odd(k)) return 0;
1017 140 : if (Z_ispowerall(absi(x), k, pt))
1018 : {
1019 140 : if (pt) *pt = negi(*pt);
1020 140 : return 1;
1021 : };
1022 0 : return 0;
1023 : }
1024 :
1025 : /* is x a K-th power mod p ? Assume p prime. */
1026 : int
1027 231 : Fp_ispower(GEN x, GEN K, GEN p)
1028 : {
1029 231 : pari_sp av = avma;
1030 : GEN p_1;
1031 : long r;
1032 231 : x = modii(x, p);
1033 231 : if (!signe(x) || equali1(x)) { avma = av; return 1; }
1034 : /* implies p > 2 */
1035 105 : p_1 = subiu(p,1);
1036 105 : K = gcdii(K, p_1);
1037 105 : if (absequaliu(K, 2)) { r = kronecker(x,p); avma = av; return (r > 0); }
1038 35 : x = Fp_pow(x, diviiexact(p_1,K), p);
1039 35 : avma = av; return equali1(x);
1040 : }
1041 :
1042 : /* x unit defined modulo 2^e, e > 0, p prime */
1043 : static int
1044 2373 : U2_issquare(GEN x, long e)
1045 : {
1046 2373 : long r = signe(x)>=0?mod8(x):8-mod8(x);
1047 2373 : if (e==1) return 1;
1048 2373 : if (e==2) return (r&3L) == 1;
1049 2009 : return r == 1;
1050 : }
1051 : /* x unit defined modulo p^e, e > 0, p prime */
1052 : static int
1053 4690 : Up_issquare(GEN x, GEN p, long e)
1054 4690 : { return (absequaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
1055 :
1056 : long
1057 2548 : Zn_issquare(GEN d, GEN fn)
1058 : {
1059 : long j, np;
1060 2548 : if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
1061 2548 : if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
1062 : /* integer factorization */
1063 2548 : np = nbrows(fn);
1064 5320 : for (j = 1; j <= np; ++j)
1065 : {
1066 4970 : GEN r, p = gcoeff(fn, j, 1);
1067 4970 : long e = itos(gcoeff(fn, j, 2));
1068 4970 : long v = Z_pvalrem(d,p,&r);
1069 4970 : if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
1070 : }
1071 350 : return 1;
1072 : }
1073 :
1074 : static long
1075 1113 : Qp_ispower(GEN x, GEN K, GEN *pt)
1076 : {
1077 1113 : pari_sp av = avma;
1078 1113 : GEN z = Qp_sqrtn(x, K, NULL);
1079 1113 : if (!z) { avma = av; return 0; }
1080 819 : if (pt) *pt = z;
1081 819 : return 1;
1082 : }
1083 :
1084 : long
1085 7002947 : ispower(GEN x, GEN K, GEN *pt)
1086 : {
1087 : GEN z;
1088 :
1089 7002947 : if (!K) return gisanypower(x, pt);
1090 2786 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
1091 2786 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
1092 2786 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
1093 2758 : switch(typ(x)) {
1094 : case t_INT:
1095 714 : return Z_ispowerall(x, itou(K), pt);
1096 : case t_FRAC:
1097 : {
1098 21 : GEN a = gel(x,1), b = gel(x,2);
1099 21 : ulong k = itou(K);
1100 21 : if (pt) {
1101 14 : z = cgetg(3, t_FRAC);
1102 14 : if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
1103 14 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1104 : }
1105 0 : avma = (pari_sp)(z + 3); return 0;
1106 : }
1107 7 : return Z_ispower(a, k) && Z_ispower(b, k);
1108 : }
1109 : case t_INTMOD:
1110 189 : return Zn_ispower(gel(x,2), gel(x,1), K, pt);
1111 : case t_FFELT:
1112 112 : return FF_ispower(x, K, pt);
1113 :
1114 : case t_PADIC:
1115 1113 : return Qp_ispower(x, K, pt);
1116 : case t_POLMOD:
1117 91 : return polmodispower(x, K, pt);
1118 : case t_POL:
1119 490 : return polispower(x, K, pt);
1120 : case t_RFRAC: {
1121 7 : GEN a = gel(x,1), b = gel(x,2);
1122 7 : if (pt) {
1123 7 : z = cgetg(3, t_RFRAC);
1124 7 : if (ispower(a, K, &a) && polispower(b, K, &b)) {
1125 7 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1126 : }
1127 0 : avma = (pari_sp)(z + 3); return 0;
1128 : }
1129 0 : return (ispower(a, K, NULL) && polispower(b, K, NULL));
1130 : }
1131 : case t_REAL:
1132 7 : if (signe(x) < 0 && !mpodd(K)) return 0;
1133 : case t_COMPLEX:
1134 14 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1135 14 : return 1;
1136 :
1137 : case t_SER:
1138 7 : if (signe(x) && (!dvdsi(valp(x), K) || !ispower(gel(x,2), K, NULL)))
1139 0 : return 0;
1140 7 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1141 7 : return 1;
1142 : }
1143 0 : pari_err_TYPE("ispower",x);
1144 0 : return 0; /* not reached */
1145 : }
1146 :
1147 : long
1148 7000161 : gisanypower(GEN x, GEN *pty)
1149 : {
1150 7000161 : long tx = typ(x);
1151 : ulong k, h;
1152 7000161 : if (tx == t_INT) return Z_isanypower(x, pty);
1153 14 : if (tx == t_FRAC)
1154 : {
1155 14 : pari_sp av = avma;
1156 14 : GEN fa, P, E, a = gel(x,1), b = gel(x,2);
1157 : long i, j, p, e;
1158 14 : int sw = (abscmpii(a, b) > 0);
1159 :
1160 14 : if (sw) swap(a, b);
1161 14 : k = Z_isanypower(a, pty? &a: NULL);
1162 14 : if (!k)
1163 : { /* a = -1,1 or not a pure power */
1164 7 : if (!is_pm1(a)) { avma = av; return 0; }
1165 7 : if (signe(a) < 0) b = negi(b);
1166 7 : k = Z_isanypower(b, pty? &b: NULL);
1167 7 : if (!k || !pty) { avma = av; return k; }
1168 7 : *pty = gerepileupto(av, ginv(b));
1169 7 : return k;
1170 : }
1171 7 : fa = factoru(k);
1172 7 : P = gel(fa,1);
1173 7 : E = gel(fa,2); h = k;
1174 14 : for (i = lg(P) - 1; i > 0; i--)
1175 : {
1176 7 : p = P[i];
1177 7 : e = E[i];
1178 21 : for (j = 0; j < e; j++)
1179 14 : if (!is_kth_power(b, p, &b)) break;
1180 7 : if (j < e) k /= upowuu(p, e - j);
1181 : }
1182 7 : if (k == 1) { avma = av; return 0; }
1183 7 : if (!pty) { avma = av; return k; }
1184 0 : if (k != h) a = powiu(a, h/k);
1185 0 : *pty = gerepilecopy(av, mkfrac(a, b));
1186 0 : return k;
1187 : }
1188 0 : pari_err_TYPE("gisanypower", x);
1189 0 : return 0; /* not reached */
1190 : }
1191 :
1192 : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
1193 : * No need to optimize for 2,3,5,7 powers (done before) */
1194 : static long
1195 505729 : split_exponent(ulong e, GEN *x)
1196 : {
1197 : GEN fa, P, E;
1198 505729 : long i, j, l, k = 1;
1199 505729 : if (e == 1) return 1;
1200 14 : fa = factoru(e);
1201 14 : P = gel(fa,1);
1202 14 : E = gel(fa,2); l = lg(P);
1203 28 : for (i = 1; i < l; i++)
1204 : {
1205 14 : ulong p = P[i];
1206 28 : for (j = 0; j < E[i]; j++)
1207 : {
1208 : GEN y;
1209 14 : if (!is_kth_power(*x, p, &y)) break;
1210 14 : k *= p; *x = y;
1211 : }
1212 : }
1213 14 : return k;
1214 : }
1215 :
1216 : static long
1217 864316 : Z_isanypower_nosmalldiv(GEN *px)
1218 : { /* any prime divisor of x is > 102 */
1219 864316 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
1220 864316 : const double LOG103 = 4.6347; /* lower bound for log(103) */
1221 : forprime_t T;
1222 864316 : ulong mask = 7, e2;
1223 : long k, ex;
1224 864316 : GEN y, x = *px;
1225 :
1226 864316 : k = 1;
1227 864316 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
1228 864316 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
1229 864316 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
1230 864316 : if (u_forprime_init(&T, 11, e2))
1231 : {
1232 16814 : GEN logx = NULL;
1233 16814 : const ulong Q = 30011; /* prime */
1234 : ulong p, xmodQ;
1235 16814 : double dlogx = 0;
1236 : /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
1237 : * for large p the modular checks are no longer competitively fast */
1238 33649 : while ( (ex = is_pth_power(x, &y, &T, 30)) )
1239 : {
1240 21 : k *= ex; x = y;
1241 21 : e2 = (ulong)((expi(x) + 1) / LOG2_103);
1242 21 : u_forprime_restrict(&T, e2);
1243 : }
1244 16814 : if (DEBUGLEVEL>4) err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x));
1245 16814 : xmodQ = umodiu(x, Q);
1246 : /* test Q | x, just in case */
1247 16814 : if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
1248 : /* x^(1/p) < 2^31 */
1249 16800 : p = T.p;
1250 16800 : if (p <= e2)
1251 : {
1252 16786 : logx = logr_abs( itor(x, DEFAULTPREC) );
1253 16786 : dlogx = rtodbl(logx);
1254 16786 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1255 : }
1256 151480 : while (p && p <= e2)
1257 : { /* is x a p-th power ? By computing y = round(x^(1/p)).
1258 : * Check whether y^p = x, first mod Q, then exactly. */
1259 117880 : pari_sp av = avma;
1260 : long e;
1261 117880 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
1262 117880 : ulong ymodQ = umodiu(y,Q);
1263 117880 : if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
1264 21 : || !equalii(powiu(y, p), x)) avma = av;
1265 : else
1266 : {
1267 21 : k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
1268 21 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1269 21 : u_forprime_restrict(&T, e2);
1270 21 : continue; /* if success, retry same p */
1271 : }
1272 117859 : p = u_forprime_next(&T);
1273 : }
1274 : }
1275 864302 : *px = x; return k;
1276 : }
1277 :
1278 : static ulong tinyprimes[] = {
1279 : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
1280 : 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
1281 : 157, 163, 167, 173, 179, 181, 191, 193, 197, 199
1282 : };
1283 :
1284 : /* disregard the sign of x, caller will take care of x < 0 */
1285 : static long
1286 7000630 : Z_isanypower_aux(GEN x, GEN *pty)
1287 : {
1288 : long ex, v, i, l, k;
1289 : GEN y, P, E;
1290 7000630 : ulong mask, e = 0;
1291 :
1292 7000630 : if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
1293 :
1294 7000616 : if (signe(x) < 0) x = negi(x);
1295 7000616 : k = l = 1;
1296 7000616 : P = cgetg(26 + 1, t_VECSMALL);
1297 7000616 : E = cgetg(26 + 1, t_VECSMALL);
1298 : /* trial division */
1299 122926272 : for(i = 0; i < 26; i++)
1300 : {
1301 60128943 : ulong p = tinyprimes[i];
1302 : int stop;
1303 60128943 : v = Z_lvalrem_stop(&x, p, &stop);
1304 60128943 : if (v)
1305 : {
1306 7922439 : P[l] = p;
1307 7922439 : E[l] = v; l++;
1308 13588862 : e = ugcd(e, v); if (e == 1) goto END;
1309 : }
1310 54710628 : if (stop) {
1311 248108 : if (is_pm1(x)) k = e;
1312 248108 : goto END;
1313 : }
1314 : }
1315 :
1316 1334193 : if (e)
1317 : { /* Bingo. Result divides e */
1318 : long v3, v5, v7;
1319 505715 : ulong e2 = e;
1320 505715 : v = u_lvalrem(e2, 2, &e2);
1321 505715 : if (v)
1322 : {
1323 375277 : for (i = 0; i < v; i++)
1324 : {
1325 374185 : if (!Z_issquareall(x, &y)) break;
1326 1302 : k <<= 1; x = y;
1327 : }
1328 : }
1329 505715 : mask = 0;
1330 505715 : v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
1331 505715 : v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
1332 505715 : v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
1333 1011507 : while ( (ex = is_357_power(x, &y, &mask)) ) {
1334 77 : x = y;
1335 77 : switch(ex)
1336 : {
1337 28 : case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
1338 28 : case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
1339 21 : case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
1340 : }
1341 : }
1342 505715 : k *= split_exponent(e2, &x);
1343 : }
1344 : else
1345 828478 : k = Z_isanypower_nosmalldiv(&x);
1346 : END:
1347 7000616 : if (pty && k != 1)
1348 : {
1349 7959 : if (e)
1350 : { /* add missing small factors */
1351 6881 : y = powuu(P[1], E[1] / k);
1352 6881 : for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
1353 6881 : x = equali1(x)? y: mulii(x,y);
1354 : }
1355 7959 : *pty = x;
1356 : }
1357 7000616 : return k == 1? 0: k;
1358 : }
1359 :
1360 : long
1361 7000630 : Z_isanypower(GEN x, GEN *pty)
1362 : {
1363 7000630 : pari_sp av = avma;
1364 7000630 : long k = Z_isanypower_aux(x, pty);
1365 7000630 : if (!k) { avma = av; return 0; }
1366 8001 : if (signe(x) < 0)
1367 : {
1368 42 : long v = vals(k);
1369 42 : if (v)
1370 : {
1371 : GEN y;
1372 28 : k >>= v;
1373 28 : if (k == 1) { avma = av; return 0; }
1374 21 : if (!pty) { avma = av; return k; }
1375 14 : y = *pty;
1376 14 : y = powiu(y, 1<<v);
1377 14 : togglesign(y);
1378 14 : *pty = gerepileuptoint(av, y);
1379 14 : return k;
1380 : }
1381 14 : if (pty) togglesign_safe(pty);
1382 : }
1383 7973 : if (pty) *pty = gerepilecopy(av, *pty); else avma = av;
1384 7973 : return k;
1385 : }
1386 :
1387 : /* Faster than */
1388 : /* !cmpii(n, int2n(vali(n))) */
1389 : /* !cmpis(shifti(n, -vali(n)), 1) */
1390 : /* expi(n) == vali(n) */
1391 : /* hamming(n) == 1 */
1392 : /* even for single-word values, and much faster for multiword values. */
1393 : /* If all you have is a word, you can just use n & !(n & (n-1)). */
1394 : long
1395 29745 : Z_ispow2(GEN n)
1396 : {
1397 : GEN xp;
1398 : long i, lx;
1399 : ulong u;
1400 29745 : if (signe(n) != 1) return 0;
1401 29738 : xp = int_LSW(n);
1402 29738 : lx = lgefint(n);
1403 29738 : u = *xp;
1404 29770 : for (i = 3; i < lx; ++i)
1405 : {
1406 27308 : if (u) return 0;
1407 32 : xp = int_nextW(xp);
1408 32 : u = *xp;
1409 : }
1410 2462 : return !(u & (u-1)); /* faster than hamming_word(u) == 1 */
1411 : }
1412 :
1413 : static long
1414 841628 : isprimepower_i(GEN n, GEN *pt, long flag)
1415 : {
1416 841628 : pari_sp av = avma;
1417 : long i, v;
1418 :
1419 841628 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1420 841628 : if (signe(n) <= 0) return 0;
1421 :
1422 841628 : if (lgefint(n) == 3)
1423 : {
1424 : ulong p;
1425 541061 : v = uisprimepower(n[2], &p);
1426 541061 : if (v)
1427 : {
1428 54849 : if (pt) *pt = utoipos(p);
1429 54849 : return v;
1430 : }
1431 486212 : return 0;
1432 : }
1433 1662027 : for (i = 0; i < 26; i++)
1434 : {
1435 1626189 : ulong p = tinyprimes[i];
1436 1626189 : v = Z_lvalrem(n, p, &n);
1437 1626190 : if (v)
1438 : {
1439 264730 : avma = av;
1440 264730 : if (!is_pm1(n)) return 0;
1441 333 : if (pt) *pt = utoipos(p);
1442 334 : return v;
1443 : }
1444 : }
1445 : /* p | n => p >= 103 */
1446 35838 : v = Z_isanypower_nosmalldiv(&n); /* expensive */
1447 35838 : if (!(flag? isprime(n): BPSW_psp(n))) { avma = av; return 0; }
1448 5526 : if (pt) *pt = gerepilecopy(av, n); else avma = av;
1449 5526 : return v;
1450 : }
1451 : long
1452 840098 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
1453 : long
1454 1530 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1455 :
1456 : long
1457 541698 : uisprimepower(ulong n, ulong *pp)
1458 : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1459 : * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1460 541698 : const ulong CUTOFF = 200UL;
1461 541698 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1462 541698 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1463 : #ifdef LONG_IS_64BIT
1464 : /* primes preceeding the appropriate root of ULONG_MAX. */
1465 481464 : const ulong ROOT9 = 137;
1466 481464 : const ulong ROOT8 = 251;
1467 481464 : const ulong ROOT7 = 563;
1468 481464 : const ulong ROOT5 = 7129;
1469 481464 : const ulong ROOT4 = 65521;
1470 : #else
1471 60234 : const ulong ROOT9 = 11;
1472 60234 : const ulong ROOT8 = 13;
1473 60234 : const ulong ROOT7 = 23;
1474 60234 : const ulong ROOT5 = 83;
1475 60234 : const ulong ROOT4 = 251;
1476 : #endif
1477 : ulong mask;
1478 : long v, i;
1479 : int e;
1480 541698 : if (n < 2) return 0;
1481 541684 : if (!odd(n)) {
1482 270690 : if (n & (n-1)) return 0;
1483 886 : *pp = 2; return vals(n);
1484 : }
1485 270994 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1486 3653596 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1487 : {
1488 3594527 : ulong p = tinyprimes[i];
1489 3594527 : if (n % p == 0)
1490 : {
1491 211519 : v = u_lvalrem(n, p, &n);
1492 211519 : if (n == 1) { *pp = p; return v; }
1493 209438 : return 0;
1494 : }
1495 : }
1496 : /* p | n => p >= CUTOFF */
1497 :
1498 59069 : if (n < CUTOFF3)
1499 : {
1500 46354 : if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1501 0 : if (uissquareall(n, &n)) { *pp = n; return 2; }
1502 0 : return 0;
1503 : }
1504 :
1505 : /* Check for squares, fourth powers, and eighth powers as appropriate. */
1506 12715 : v = 1;
1507 12715 : if (uissquareall(n, &n)) {
1508 0 : v <<= 1;
1509 0 : if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1510 0 : v <<= 1;
1511 0 : if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1512 : }
1513 : }
1514 :
1515 12715 : if (CUTOFF > ROOT5) mask = 1;
1516 : else
1517 : {
1518 12714 : const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1519 12714 : if (n < CUTOFF5) mask = 1; else mask = 3;
1520 12714 : if (CUTOFF <= ROOT7)
1521 : {
1522 12714 : const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1523 12714 : if (n >= CUTOFF7) mask = 7;
1524 : }
1525 : }
1526 :
1527 12715 : if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1528 12715 : if ((e = uis_357_power(n, &n, &mask))) v *= e;
1529 :
1530 12715 : if (uisprime_101(n)) { *pp = n; return v; }
1531 6984 : return 0;
1532 : }
1533 :
1534 : /*********************************************************************/
1535 : /** **/
1536 : /** KRONECKER SYMBOL **/
1537 : /** **/
1538 : /*********************************************************************/
1539 : /* t = 3,5 mod 8 ? (= 2 not a square mod t) */
1540 : static int
1541 533403418 : ome(long t)
1542 : {
1543 533403418 : switch(t & 7)
1544 : {
1545 : case 3:
1546 304064181 : case 5: return 1;
1547 229339237 : default: return 0;
1548 : }
1549 : }
1550 : /* t a t_INT, is t = 3,5 mod 8 ? */
1551 : static int
1552 4412292 : gome(GEN t)
1553 4412292 : { return signe(t)? ome( mod2BIL(t) ): 0; }
1554 :
1555 : /* assume y odd, return kronecker(x,y) * s */
1556 : static long
1557 401424969 : krouu_s(ulong x, ulong y, long s)
1558 : {
1559 401424969 : ulong x1 = x, y1 = y, z;
1560 2178370899 : while (x1)
1561 : {
1562 1375590935 : long r = vals(x1);
1563 1375677878 : if (r)
1564 : {
1565 749012458 : if (odd(r) && ome(y1)) s = -s;
1566 748855541 : x1 >>= r;
1567 : }
1568 1375520961 : if (x1 & y1 & 2) s = -s;
1569 1375520961 : z = y1 % x1; y1 = x1; x1 = z;
1570 : }
1571 401354995 : return (y1 == 1)? s: 0;
1572 : }
1573 :
1574 : long
1575 5057786 : kronecker(GEN x, GEN y)
1576 : {
1577 5057786 : pari_sp av = avma;
1578 5057786 : long s = 1, r;
1579 : ulong xu, yu;
1580 :
1581 5057786 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
1582 5057786 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
1583 5057786 : switch (signe(y))
1584 : {
1585 0 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
1586 0 : case 0: return is_pm1(x);
1587 : }
1588 5057786 : r = vali(y);
1589 5057786 : if (r)
1590 : {
1591 11943 : if (!mpodd(x)) { avma = av; return 0; }
1592 10228 : if (odd(r) && gome(x)) s = -s;
1593 10228 : y = shifti(y,-r);
1594 : }
1595 5056071 : x = modii(x,y);
1596 10186886 : while (lgefint(x) > 3) /* x < y */
1597 : {
1598 : GEN z;
1599 74744 : r = vali(x);
1600 74744 : if (r)
1601 : {
1602 40926 : if (odd(r) && gome(y)) s = -s;
1603 40926 : x = shifti(x,-r);
1604 : }
1605 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1606 74744 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
1607 74744 : z = remii(y,x); y = x; x = z;
1608 74744 : if (gc_needed(av,2))
1609 : {
1610 0 : if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
1611 0 : gerepileall(av, 2, &x, &y);
1612 : }
1613 : }
1614 5056071 : xu = itou(x);
1615 5056071 : if (!xu) return is_pm1(y)? s: 0;
1616 5035446 : r = vals(xu);
1617 5035446 : if (r)
1618 : {
1619 3493837 : if (odd(r) && gome(y)) s = -s;
1620 3493837 : xu >>= r;
1621 : }
1622 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1623 5035446 : if (xu & mod2BIL(y) & 2) s = -s;
1624 5035446 : yu = umodiu(y, xu);
1625 5035446 : avma = av; return krouu_s(yu, xu, s);
1626 : }
1627 :
1628 : long
1629 27699 : krois(GEN x, long y)
1630 : {
1631 : ulong yu;
1632 27699 : long s = 1;
1633 :
1634 27699 : if (y <= 0)
1635 : {
1636 0 : if (y == 0) return is_pm1(x);
1637 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
1638 : }
1639 : else
1640 27699 : yu = (ulong)y;
1641 27699 : if (!odd(yu))
1642 : {
1643 : long r;
1644 13335 : if (!mpodd(x)) return 0;
1645 9961 : r = vals(yu); yu >>= r;
1646 9961 : if (odd(r) && gome(x)) s = -s;
1647 : }
1648 24325 : return krouu_s(umodiu(x, yu), yu, s);
1649 : }
1650 : /* assume y != 0 */
1651 : long
1652 283954426 : kroiu(GEN x, ulong y)
1653 : {
1654 : long r;
1655 283954426 : if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
1656 1800022 : if (!mpodd(x)) return 0;
1657 1764042 : r = vals(y); y >>= r;
1658 1764042 : return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
1659 : }
1660 :
1661 : /* assume y > 0, odd, return s * kronecker(x,y) */
1662 : static long
1663 31309 : krouodd(ulong x, GEN y, long s)
1664 : {
1665 : long r;
1666 31309 : if (lgefint(y) == 3) return krouu_s(x, y[2], s);
1667 13046 : if (!x) return 0; /* y != 1 */
1668 13046 : r = vals(x);
1669 13046 : if (r)
1670 : {
1671 8426 : if (odd(r) && gome(y)) s = -s;
1672 8426 : x >>= r;
1673 : }
1674 : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1675 13046 : if (x & mod2BIL(y) & 2) s = -s;
1676 13046 : return krouu_s(umodiu(y,x), x, s);
1677 : }
1678 :
1679 : long
1680 31201 : krosi(long x, GEN y)
1681 : {
1682 31201 : const pari_sp av = avma;
1683 31201 : long s = 1, r;
1684 31201 : switch (signe(y))
1685 : {
1686 0 : case -1: y = negi(y); if (x < 0) s = -1; break;
1687 0 : case 0: return (x==1 || x==-1);
1688 : }
1689 31201 : r = vali(y);
1690 31201 : if (r)
1691 : {
1692 0 : if (!odd(x)) { avma = av; return 0; }
1693 0 : if (odd(r) && ome(x)) s = -s;
1694 0 : y = shifti(y,-r);
1695 : }
1696 31201 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
1697 31201 : s = krouodd((ulong)x, y, s);
1698 31201 : avma = av; return s;
1699 : }
1700 :
1701 : long
1702 108 : kroui(ulong x, GEN y)
1703 : {
1704 108 : const pari_sp av = avma;
1705 108 : long s = 1, r;
1706 108 : switch (signe(y))
1707 : {
1708 0 : case -1: y = negi(y); break;
1709 0 : case 0: return x==1UL;
1710 : }
1711 108 : r = vali(y);
1712 108 : if (r)
1713 : {
1714 0 : if (!odd(x)) { avma = av; return 0; }
1715 0 : if (odd(r) && ome(x)) s = -s;
1716 0 : y = shifti(y,-r);
1717 : }
1718 108 : s = krouodd(x, y, s);
1719 108 : avma = av; return s;
1720 : }
1721 :
1722 : long
1723 46125998 : kross(long x, long y)
1724 : {
1725 : ulong yu;
1726 46125998 : long s = 1;
1727 :
1728 46125998 : if (y <= 0)
1729 : {
1730 385 : if (y == 0) return (labs(x)==1);
1731 385 : yu = (ulong)-y; if (x < 0) s = -1;
1732 : }
1733 : else
1734 46125613 : yu = (ulong)y;
1735 46125998 : if (!odd(yu))
1736 : {
1737 : long r;
1738 4745297 : if (!odd(x)) return 0;
1739 4744618 : r = vals(yu); yu >>= r;
1740 4744618 : if (odd(r) && ome(x)) s = -s;
1741 : }
1742 46125319 : x %= (long)yu; if (x < 0) x += yu;
1743 46125319 : return krouu_s((ulong)x, yu, s);
1744 : }
1745 :
1746 : long
1747 66186207 : krouu(ulong x, ulong y)
1748 : {
1749 : long r;
1750 66186207 : if (odd(y)) return krouu_s(x, y, 1);
1751 1926 : if (!odd(x)) return 0;
1752 1926 : r = vals(y); y >>= r;
1753 1926 : return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
1754 : }
1755 :
1756 : /*********************************************************************/
1757 : /** **/
1758 : /** HILBERT SYMBOL **/
1759 : /** **/
1760 : /*********************************************************************/
1761 : /* x,y are t_INT or t_REAL */
1762 : static long
1763 9982 : mphilbertoo(GEN x, GEN y)
1764 : {
1765 9982 : long sx = signe(x), sy = signe(y);
1766 9982 : if (!sx || !sy) return 0;
1767 9982 : return (sx < 0 && sy < 0)? -1: 1;
1768 : }
1769 :
1770 : long
1771 53144 : hilbertii(GEN x, GEN y, GEN p)
1772 : {
1773 : pari_sp av;
1774 : long oddvx, oddvy, z;
1775 :
1776 53144 : if (!p) return mphilbertoo(x,y);
1777 43183 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
1778 43183 : if (!signe(x) || !signe(y)) return 0;
1779 43162 : av = avma;
1780 43162 : oddvx = odd(Z_pvalrem(x,p,&x));
1781 43162 : oddvy = odd(Z_pvalrem(y,p,&y));
1782 : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
1783 43162 : if (absequaliu(p, 2))
1784 : {
1785 10689 : z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
1786 10689 : if (oddvx && gome(y)) z = -z;
1787 10689 : if (oddvy && gome(x)) z = -z;
1788 : }
1789 : else
1790 : {
1791 32473 : z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
1792 32473 : if (oddvx && kronecker(y,p) < 0) z = -z;
1793 32473 : if (oddvy && kronecker(x,p) < 0) z = -z;
1794 : }
1795 43162 : avma = av; return z;
1796 : }
1797 :
1798 : static void
1799 196 : err_prec(void) { pari_err_PREC("hilbert"); }
1800 : static void
1801 161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
1802 : static void
1803 56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
1804 :
1805 : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
1806 : * Return lift(x) provided it's p-adic accuracy is large enough to decide
1807 : * hilbert()'s value [ problem at p = 2 ] */
1808 : static GEN
1809 420 : lift_intmod(GEN x, GEN *pp)
1810 : {
1811 420 : GEN p = *pp, N = gel(x,1);
1812 420 : x = gel(x,2);
1813 420 : if (!p)
1814 : {
1815 266 : *pp = p = N;
1816 266 : switch(itos_or_0(p))
1817 : {
1818 : case 2:
1819 126 : case 4: err_prec();
1820 : }
1821 140 : return x;
1822 : }
1823 154 : if (!signe(p)) err_oo(N);
1824 112 : if (absequaliu(p,2))
1825 42 : { if (vali(N) <= 2) err_prec(); }
1826 : else
1827 70 : { if (!dvdii(N,p)) err_p(N,p); }
1828 28 : if (!signe(x)) err_prec();
1829 21 : return x;
1830 : }
1831 : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
1832 : * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
1833 : * to decide hilbert()'s value [ problem at p = 2 ]*/
1834 : static GEN
1835 210 : lift_padic(GEN x, GEN *pp)
1836 : {
1837 210 : GEN p = *pp, q = gel(x,2), y = gel(x,4);
1838 210 : if (!p) *pp = p = q;
1839 147 : else if (!equalii(p,q)) err_p(p, q);
1840 105 : if (absequaliu(p,2) && precp(x) <= 2) err_prec();
1841 70 : if (!signe(y)) err_prec();
1842 70 : return odd(valp(x))? mulii(p,y): y;
1843 : }
1844 :
1845 : long
1846 658 : hilbert(GEN x, GEN y, GEN p)
1847 : {
1848 658 : pari_sp av = avma;
1849 658 : long tx = typ(x), ty = typ(y), z;
1850 :
1851 658 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
1852 658 : if (tx == t_REAL)
1853 : {
1854 77 : if (p && signe(p)) err_oo(p);
1855 63 : switch (ty)
1856 : {
1857 : case t_INT:
1858 7 : case t_REAL: return mphilbertoo(x,y);
1859 0 : case t_FRAC: return mphilbertoo(x,gel(y,1));
1860 56 : default: pari_err_TYPE2("hilbert",x,y);
1861 : }
1862 : }
1863 581 : if (ty == t_REAL)
1864 : {
1865 14 : if (p && signe(p)) err_oo(p);
1866 14 : switch (tx)
1867 : {
1868 : case t_INT:
1869 14 : case t_REAL: return mphilbertoo(x,y);
1870 0 : case t_FRAC: return mphilbertoo(gel(x,1),y);
1871 0 : default: pari_err_TYPE2("hilbert",x,y);
1872 : }
1873 : }
1874 567 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
1875 364 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
1876 :
1877 308 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
1878 245 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
1879 :
1880 168 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
1881 168 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
1882 :
1883 168 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
1884 168 : if (p && !signe(p)) p = NULL;
1885 168 : z = hilbertii(x,y,p); avma = av; return z;
1886 : }
1887 :
1888 : /*******************************************************************/
1889 : /* */
1890 : /* SQUARE ROOT MODULO p */
1891 : /* */
1892 : /*******************************************************************/
1893 :
1894 : static ulong
1895 22208927 : Fl_2gener_pre_all(long e, ulong p, ulong pi, ulong *pt_m)
1896 : {
1897 : ulong y, m;
1898 : long k, i;
1899 22208927 : ulong q = (p-1) >> e; /* q = (p-1)/2^oo is odd */
1900 39582997 : for (k=2; ; k++)
1901 : { /* loop terminates for k < p (even if p composite) */
1902 39582997 : i = krouu(k, p);
1903 39582997 : if (i >= 0)
1904 : {
1905 17374077 : if (i) continue;
1906 7 : pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1907 : }
1908 22208920 : y = m = Fl_powu_pre(k, q, p, pi);
1909 60670455 : for (i=1; i<e; i++)
1910 38461535 : if ((m = Fl_sqr_pre(m, p, pi)) == 1) break;
1911 22208920 : if (i == e) break; /* success */
1912 17374070 : }
1913 22208920 : *pt_m = m;
1914 22208920 : return y;
1915 : }
1916 :
1917 : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
1918 : static ulong
1919 53892830 : Fl_sqrt_i(ulong a, ulong p, ulong pi, ulong y, ulong m)
1920 : {
1921 : long i, e, k;
1922 : ulong p1, q, v, w;
1923 :
1924 53892830 : if (!a) return 0;
1925 53046221 : p1 = p - 1; e = vals(p1);
1926 53057593 : if (e == 0) /* p = 2 */
1927 : {
1928 417354 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1929 417347 : return ((a & 1) == 0)? 0: 1;
1930 : }
1931 52640239 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
1932 52640239 : if (e == 1) y = p1;
1933 22208927 : else if (y==0) y = Fl_2gener_pre_all(e, p, pi, &m);
1934 52640232 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
1935 52614980 : if (!p1) return 0;
1936 52614980 : v = Fl_mul_pre(a, p1, p, pi);
1937 52619383 : w = Fl_mul_pre(v, p1, p, pi);
1938 123826761 : while (w != 1)
1939 : { /* a*w = v^2, y primitive 2^e-th root of 1
1940 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1941 18686035 : p1 = Fl_sqr_pre(w,p,pi);
1942 18686035 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
1943 18686035 : if (k == e) return ~0UL;
1944 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1945 18590563 : p1 = y;
1946 18590563 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
1947 18590563 : y = Fl_sqr_pre(p1, p, pi); e = k;
1948 18590563 : w = Fl_mul_pre(y, w, p, pi);
1949 18590563 : v = Fl_mul_pre(v, p1, p, pi);
1950 : }
1951 52522366 : p1 = p - v; if (v > p1) v = p1;
1952 52522366 : return v;
1953 : }
1954 :
1955 : ulong
1956 50431534 : Fl_sqrt(ulong a, ulong p)
1957 : {
1958 50431534 : ulong pi = get_Fl_red(p);
1959 50421962 : return Fl_sqrt_i(a, p, pi, 0, 0);
1960 : }
1961 :
1962 : ulong
1963 3471489 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
1964 : {
1965 3471489 : return Fl_sqrt_i(a, p, pi, 0, 0);
1966 : }
1967 :
1968 : static ulong
1969 63265 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
1970 : {
1971 : ulong x, y, m;
1972 63265 : ulong le1 = upowuu(l, e-1);
1973 100468 : for (x = 2; ; x++)
1974 : {
1975 100468 : y = Fl_powu_pre(x, r, p, pi);
1976 100468 : if (y==1) continue;
1977 78231 : m = Fl_powu_pre(y, le1, p, pi);
1978 78231 : if (m != 1) break;
1979 37203 : }
1980 63265 : *pt_m = m;
1981 63265 : return y;
1982 : }
1983 :
1984 : /* solve x^l = a , l prime in G of order q.
1985 : *
1986 : * q = (l^e)*r, e >= 1, (r,l) = 1
1987 : * y generates the l-Sylow of G
1988 : * m = y^(l^(e-1)) != 1 */
1989 : static ulong
1990 115539 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
1991 : {
1992 : ulong p1, v, w, z, dl, zm;
1993 : ulong r, e, u2;
1994 115539 : if (a==0) return a;
1995 115535 : e = u_lvalrem(p-1, l, &r);
1996 115534 : u2 = Fl_inv(l%r, r);
1997 115537 : v = Fl_powu_pre(a, u2, p,pi);
1998 115537 : w = Fl_powu_pre(v, l, p,pi);
1999 115535 : w = Fl_mul_pre(w, Fl_inv(a, p),p,pi);
2000 115537 : if (w==1) return v;
2001 63265 : if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
2002 145502 : while (w!=1)
2003 : {
2004 67688 : ulong k = 0;
2005 67688 : p1 = w;
2006 : do
2007 : {
2008 101192 : z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
2009 101192 : k++;
2010 101192 : } while (p1!=1);
2011 67688 : if (k==e) return ~0UL;
2012 18972 : dl = 0; zm = 1;
2013 66838 : while (z!=zm)
2014 : {
2015 28894 : zm = Fl_mul_pre(zm, m, p, pi); dl++;
2016 : }
2017 18972 : dl = Fl_neg(dl, l);
2018 18972 : p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
2019 18972 : m = Fl_powu_pre(m, dl, p, pi);
2020 18972 : e = k;
2021 18972 : v = Fl_mul_pre(p1,v,p,pi);
2022 18972 : y = Fl_powu_pre(p1,l,p,pi);
2023 18972 : w = Fl_mul_pre(y,w,p,pi);
2024 : }
2025 14549 : return v;
2026 : }
2027 :
2028 : ulong
2029 115541 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
2030 : {
2031 115541 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2032 : }
2033 :
2034 : ulong
2035 0 : Fl_sqrtl(ulong a, ulong l, ulong p)
2036 : {
2037 0 : ulong pi = get_Fl_red(p);
2038 0 : return Fl_sqrtl_i(a, l, p, pi, 0, 0);
2039 : }
2040 :
2041 : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
2042 : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
2043 : * and in average is about 2 or 2.5 times worse. But try both algorithms for
2044 : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
2045 : *
2046 : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
2047 : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
2048 : * If (a|p)=1, then sqrt(a) is in F_p.
2049 : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
2050 :
2051 : /* compute y^2, y = y[1] + y[2] X */
2052 : static GEN
2053 449 : sqrt_Cipolla_sqr(void *data, GEN y)
2054 : {
2055 449 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
2056 449 : GEN u2 = sqri(u), v2 = sqri(v);
2057 449 : v = subii(sqri(addii(v,u)), addii(u2,v2));
2058 449 : u = addii(u2, mulii(v2,n));
2059 : /* NOT mkvec2: must be gerepileupto-able */
2060 449 : retmkvec2(modii(u,p), modii(v,p));
2061 : }
2062 : /* compute (t+X) y^2 */
2063 : static GEN
2064 23 : sqrt_Cipolla_msqr(void *data, GEN y)
2065 : {
2066 23 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
2067 23 : ulong t = gt[2];
2068 23 : GEN d = addii(u, mului(t,v)), d2= sqri(d);
2069 23 : GEN b = remii(mulii(a,v), p);
2070 23 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
2071 23 : v = subii(d2, mulii(b,v));
2072 : /* NOT mkvec2: must be gerepileupto-able */
2073 23 : retmkvec2(modii(u,p), modii(v,p));
2074 : }
2075 : /* assume a reduced mod p [ otherwise correct but inefficient ] */
2076 : static GEN
2077 8 : sqrt_Cipolla(GEN a, GEN p)
2078 : {
2079 : pari_sp av1;
2080 : GEN u, v, n, y, pov2;
2081 : ulong t;
2082 :
2083 8 : if (kronecker(a, p) < 0) return NULL;
2084 8 : pov2 = shifti(p,-1);
2085 8 : if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
2086 :
2087 8 : av1 = avma;
2088 41 : for(t=1; ; t++)
2089 : {
2090 41 : n = subsi((long)(t*t), a);
2091 41 : if (kronecker(n, p) < 0) break;
2092 33 : avma = av1;
2093 33 : }
2094 :
2095 : /* compute (t+X)^((p-1)/2) =: u+vX */
2096 8 : u = utoipos(t);
2097 8 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
2098 : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
2099 : /* Now u+vX = (t+X)^((p-1)/2); thus
2100 : * (u+vX)(t+X) = sqrt(a) + 0 X
2101 : * Whence,
2102 : * sqrt(a) = (u+vt)t - v*a
2103 : * 0 = (u+vt)
2104 : * Thus a square root is v*a */
2105 :
2106 8 : v = Fp_mul(gel(y, 2), a, p);
2107 8 : if (cmpii(v,pov2) > 0) v = subii(p,v);
2108 8 : return v;
2109 : }
2110 :
2111 : #define sqrmod(x,p) (remii(sqri(x),p))
2112 :
2113 : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
2114 : GEN
2115 2338164 : Fp_sqrt(GEN a, GEN p)
2116 : {
2117 2338164 : pari_sp av = avma, av1;
2118 : long i, k, e;
2119 : GEN p1, q, v, y, w, m;
2120 :
2121 2338164 : if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
2122 2338164 : if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
2123 2338164 : if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
2124 2338164 : if (lgefint(p) == 3)
2125 : {
2126 2328056 : ulong pp = uel(p,2), u = Fl_sqrt(umodiu(a, pp), pp);
2127 2328042 : if (u == ~0UL) return NULL;
2128 2327993 : return utoi(u);
2129 : }
2130 :
2131 10108 : p1 = addsi(-1,p); e = vali(p1);
2132 10108 : a = modii(a, p);
2133 :
2134 : /* On average, the algorithm of Cipolla is better than the algorithm of
2135 : * Tonelli and Shanks if and only if e(e-1)>8*log2(n)+20
2136 : * see LNCS 2286 pp 430 [GTL] */
2137 10108 : if (e*(e-1) > 20 + 8 * expi(p))
2138 : {
2139 8 : v = sqrt_Cipolla(a,p);
2140 8 : if (!v) { avma = av; return NULL; }
2141 8 : return gerepileuptoint(av,v);
2142 : }
2143 :
2144 10100 : if (e == 0) /* p = 2 */
2145 : {
2146 0 : avma = av;
2147 0 : if (!absequaliu(p,2)) pari_err_PRIME("Fp_sqrt [modulus]",p);
2148 0 : if (!signe(a) || !mod2(a)) return gen_0;
2149 0 : return gen_1;
2150 : }
2151 10100 : q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
2152 10100 : if (e == 1) y = p1;
2153 : else /* look for an odd power of a primitive root */
2154 8330 : for (k=2; ; k++)
2155 : { /* loop terminates for k < p (even if p composite) */
2156 :
2157 8330 : i = krosi(k,p);
2158 8330 : if (i >= 0)
2159 : {
2160 3137 : if (i) continue;
2161 0 : pari_err_PRIME("Fp_sqrt [modulus]",p);
2162 : }
2163 5193 : av1 = avma;
2164 5193 : y = m = Fp_pow(utoipos((ulong)k),q,p);
2165 12825 : for (i=1; i<e; i++)
2166 7632 : if (gequal1(m = sqrmod(m,p))) break;
2167 5193 : if (i == e) break; /* success */
2168 0 : avma = av1;
2169 3137 : }
2170 :
2171 10100 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ [(q-1)/2] */
2172 10100 : if (!signe(p1)) { avma=av; return gen_0; }
2173 10093 : v = Fp_mul(a, p1, p);
2174 10093 : w = Fp_mul(v, p1, p);
2175 23995 : while (!equali1(w))
2176 : { /* a*w = v^2, y primitive 2^e-th root of 1
2177 : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
2178 3819 : p1 = sqrmod(w,p);
2179 3819 : for (k=1; !equali1(p1) && k < e; k++) p1 = sqrmod(p1,p);
2180 3819 : if (k == e) { avma=av; return NULL; } /* p composite or (a/p) != 1 */
2181 : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
2182 3809 : p1 = y;
2183 3809 : for (i=1; i < e-k; i++) p1 = sqrmod(p1,p);
2184 3809 : y = sqrmod(p1, p); e = k;
2185 3809 : w = Fp_mul(y, w, p);
2186 3809 : v = Fp_mul(v, p1, p);
2187 3809 : if (gc_needed(av,1))
2188 : {
2189 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
2190 0 : gerepileall(av,3, &y,&w,&v);
2191 : }
2192 : }
2193 10083 : av1 = avma;
2194 10083 : p1 = subii(p,v); if (cmpii(v,p1) > 0) v = p1; else avma = av1;
2195 10083 : return gerepileuptoint(av, v);
2196 : }
2197 :
2198 : /*********************************************************************/
2199 : /** **/
2200 : /** GCD & BEZOUT **/
2201 : /** **/
2202 : /*********************************************************************/
2203 :
2204 : GEN
2205 4077004 : lcmii(GEN x, GEN y)
2206 : {
2207 : pari_sp av;
2208 : GEN a, b;
2209 4077004 : if (!signe(x) || !signe(y)) return gen_0;
2210 4077004 : av = avma;
2211 4077004 : a = gcdii(x,y); if (!equali1(a)) y = diviiexact(y,a);
2212 4077004 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
2213 : }
2214 :
2215 : /*********************************************************************/
2216 : /** **/
2217 : /** CHINESE REMAINDERS **/
2218 : /** **/
2219 : /*********************************************************************/
2220 :
2221 : /* Chinese Remainder Theorem. x and y must have the same type (integermod,
2222 : * polymod, or polynomial/vector/matrix recursively constructed with these
2223 : * as coefficients). Creates (with the same type) a z in the same residue
2224 : * class as x and the same residue class as y, if it is possible.
2225 : *
2226 : * We also allow (during recursion) two identical objects even if they are
2227 : * not integermod or polymod. For example:
2228 : *
2229 : * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
2230 : * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
2231 : * ? chinese(x, y)
2232 : * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
2233 :
2234 : static GEN
2235 101817 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
2236 : {
2237 101817 : GEN z = gassoc_proto(f,x,NULL);
2238 101810 : if (z == gen_1) retmkintmod(gen_0,gen_1);
2239 101796 : return z;
2240 : }
2241 :
2242 : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
2243 : * call chinese: makes Mod(0,1) a better "neutral" element */
2244 : static GEN
2245 21 : chinese_intpol(GEN x,GEN y)
2246 : {
2247 21 : pari_sp av = avma;
2248 21 : GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
2249 21 : return gerepileupto(av, chinese(z, y));
2250 : }
2251 :
2252 : GEN
2253 49 : chinese1(GEN x) { return gen_chinese(x,chinese); }
2254 :
2255 : GEN
2256 16520 : chinese(GEN x, GEN y)
2257 : {
2258 : pari_sp av;
2259 16520 : long tx = typ(x), ty;
2260 : GEN z,p1,p2,d,u,v;
2261 :
2262 16520 : if (!y) return chinese1(x);
2263 16471 : if (gidentical(x,y)) return gcopy(x);
2264 16464 : ty = typ(y);
2265 16464 : if (tx == ty) switch(tx)
2266 : {
2267 : case t_POLMOD:
2268 : {
2269 28 : GEN A = gel(x,1), B = gel(y,1);
2270 28 : GEN a = gel(x,2), b = gel(y,2);
2271 28 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
2272 28 : if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
2273 28 : av = avma;
2274 28 : d = RgX_extgcd(A,B,&u,&v);
2275 28 : p2 = gsub(b, a);
2276 28 : if (!gequal0(gmod(p2, d))) break;
2277 28 : p1 = gdiv(A,d);
2278 28 : p2 = gadd(a, gmul(gmul(u,p1), p2));
2279 :
2280 28 : z = cgetg(3, t_POLMOD);
2281 28 : gel(z,1) = gmul(p1,B);
2282 28 : gel(z,2) = gmod(p2,gel(z,1));
2283 28 : return gerepileupto(av, z);
2284 : }
2285 : case t_INTMOD:
2286 : {
2287 16401 : GEN A = gel(x,1), B = gel(y,1);
2288 16401 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
2289 16401 : z = cgetg(3,t_INTMOD);
2290 16401 : Z_chinese_pre(A, B, &C, &U, &d);
2291 16401 : c = Z_chinese_post(a, b, C, U, d);
2292 16401 : if (!c) pari_err_OP("chinese", x,y);
2293 16401 : gel(z,1) = icopy_avma(C, (pari_sp)z);
2294 16401 : gel(z,2) = icopy_avma(c, (pari_sp)gel(z,1));
2295 16401 : avma = (pari_sp)gel(z,2); return z;
2296 : }
2297 : case t_POL:
2298 : {
2299 7 : long i, lx = lg(x), ly = lg(y);
2300 7 : if (varn(x) != varn(y)) break;
2301 7 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
2302 7 : z = cgetg(lx, t_POL); z[1] = x[1];
2303 7 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2304 7 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
2305 7 : return z;
2306 : }
2307 :
2308 : case t_VEC: case t_COL: case t_MAT:
2309 : {
2310 : long i, lx;
2311 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
2312 7 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2313 7 : return z;
2314 : }
2315 : }
2316 21 : if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
2317 7 : if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
2318 0 : pari_err_OP("chinese",x,y);
2319 0 : return NULL; /* not reached */
2320 : }
2321 :
2322 : /* init chinese(Mod(.,A), Mod(.,B)) */
2323 : void
2324 436786 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
2325 : {
2326 436786 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
2327 436786 : GEN t = diviiexact(A,d);
2328 436786 : *pU = mulii(u, t);
2329 436786 : *pC = mulii(t, B);
2330 436786 : if (pd) *pd = d;
2331 436786 : }
2332 : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
2333 : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
2334 : * If d not NULL, check wether a = b mod d. */
2335 : GEN
2336 589433 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
2337 : {
2338 : GEN b_a;
2339 589433 : if (!signe(a))
2340 : {
2341 446000 : if (d && remii(b, d) != gen_0) return NULL;
2342 446000 : return Fp_mul(b, U, C);
2343 : }
2344 143433 : b_a = subii(b,a);
2345 143433 : if (d && remii(b_a, d) != gen_0) return NULL;
2346 143433 : return modii(addii(a, mulii(U, b_a)), C);
2347 : }
2348 : GEN
2349 2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
2350 : {
2351 2142 : pari_sp av = avma;
2352 2142 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
2353 2142 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
2354 : }
2355 : GEN
2356 418187 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
2357 : {
2358 418187 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
2359 418187 : return Z_chinese_post(a,b, *pC, U, NULL);
2360 : }
2361 :
2362 : /* return lift(chinese(a mod A, b mod B))
2363 : * assume(A,B)=1, a,b,A,B integers and C = A*B */
2364 : GEN
2365 83153 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
2366 : {
2367 83153 : pari_sp av = avma;
2368 83153 : GEN U = mulii(Fp_inv(A,B), A);
2369 83153 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2370 : }
2371 :
2372 : /* chinese1 for coprime moduli in Z */
2373 : static GEN
2374 69228 : chinese1_coprime_Z_aux(GEN x, GEN y)
2375 : {
2376 69228 : GEN z = cgetg(3, t_INTMOD);
2377 69228 : GEN A = gel(x,1), a = gel(x, 2);
2378 69228 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
2379 69228 : pari_sp av = avma;
2380 69228 : GEN U = mulii(Fp_inv(A,B), A);
2381 69228 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2382 69228 : gel(z,1) = C; return z;
2383 : }
2384 : GEN
2385 101768 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
2386 :
2387 : /*********************************************************************/
2388 : /** **/
2389 : /** MODULAR EXPONENTIATION **/
2390 : /** **/
2391 : /*********************************************************************/
2392 :
2393 : /* xa, ya = t_VECSMALL */
2394 : GEN
2395 161905 : ZV_producttree(GEN xa)
2396 : {
2397 161905 : long n = lg(xa)-1;
2398 161905 : long m = n==1 ? 1: expu(n-1)+1;
2399 161905 : GEN T = cgetg(m+1, t_VEC), t;
2400 : long i, j, k;
2401 161905 : t = cgetg(((n+1)>>1)+1, t_VEC);
2402 161900 : if (typ(xa)==t_VECSMALL)
2403 : {
2404 189226 : for (j=1, k=1; k<n; j++, k+=2)
2405 136788 : gel(t, j) = muluu(xa[k], xa[k+1]);
2406 52438 : if (k==n) gel(t, j) = utoi(xa[k]);
2407 : } else {
2408 292083 : for (j=1, k=1; k<n; j++, k+=2)
2409 182619 : gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
2410 109464 : if (k==n) gel(t, j) = icopy(gel(xa,k));
2411 : }
2412 161904 : gel(T,1) = t;
2413 289212 : for (i=2; i<=m; i++)
2414 : {
2415 127308 : GEN u = gel(T, i-1);
2416 127308 : long n = lg(u)-1;
2417 127308 : t = cgetg(((n+1)>>1)+1, t_VEC);
2418 336927 : for (j=1, k=1; k<n; j++, k+=2)
2419 209619 : gel(t, j) = mulii(gel(u, k), gel(u, k+1));
2420 127308 : if (k==n) gel(t, j) = gel(u, k);
2421 127308 : gel(T, i) = t;
2422 : }
2423 161904 : return T;
2424 : }
2425 :
2426 : static GEN
2427 870677 : Z_ZV_mod_tree(GEN P, GEN xa, GEN T)
2428 : {
2429 : long i,j,k;
2430 870677 : long m = lg(T)-1, n = lg(xa)-1;
2431 : GEN t;
2432 870677 : GEN Tp = cgetg(m+1, t_VEC);
2433 870666 : gel(Tp, m) = mkvec(P);
2434 1542835 : for (i=m-1; i>=1; i--)
2435 : {
2436 672184 : GEN u = gel(T, i);
2437 672184 : GEN v = gel(Tp, i+1);
2438 672184 : long n = lg(u)-1;
2439 672184 : t = cgetg(n+1, t_VEC);
2440 1650575 : for (j=1, k=1; k<n; j++, k+=2)
2441 : {
2442 978385 : gel(t, k) = modii(gel(v, j), gel(u, k));
2443 978388 : gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
2444 : }
2445 672190 : if (k==n) gel(t, k) = gel(v, j);
2446 672190 : gel(Tp, i) = t;
2447 : }
2448 : {
2449 870651 : GEN u = gel(T, i+1);
2450 870651 : GEN v = gel(Tp, i+1);
2451 870651 : long l = lg(u)-1;
2452 870651 : if (typ(xa)==t_VECSMALL)
2453 : {
2454 708744 : GEN R = cgetg(n+1, t_VECSMALL);
2455 2186270 : for (j=1, k=1; j<=l; j++, k+=2)
2456 : {
2457 1477500 : uel(R,k) = umodiu(gel(v, j), xa[k]);
2458 1477509 : if (k < n)
2459 1226693 : uel(R,k+1) = umodiu(gel(v, j), xa[k+1]);
2460 : }
2461 708770 : return R;
2462 : }
2463 : else
2464 : {
2465 161907 : GEN R = cgetg(n+1, t_VEC);
2466 533438 : for (j=1, k=1; j<=l; j++, k+=2)
2467 : {
2468 371526 : gel(R,k) = modii(gel(v, j), gel(xa,k));
2469 371530 : if (k < n)
2470 319413 : gel(R,k+1) = modii(gel(v, j), gel(xa,k+1));
2471 : }
2472 161912 : return R;
2473 : }
2474 : }
2475 : }
2476 :
2477 : static GEN
2478 2272203 : ZV_polint_tree(GEN T, GEN R, GEN xa, GEN ya)
2479 : {
2480 2272203 : long m = lg(T)-1, n = lg(ya)-1;
2481 : long i,j,k;
2482 2272203 : GEN Tp = cgetg(m+1, t_VEC);
2483 2263469 : GEN M = gel(T, 1);
2484 2263469 : GEN t = cgetg(lg(M), t_VEC);
2485 2315631 : if (typ(xa)==t_VECSMALL)
2486 : {
2487 18009899 : for (j=1, k=1; k<n; j++, k+=2)
2488 : {
2489 15855872 : pari_sp av = avma;
2490 15855872 : GEN a = mului(ya[k], gel(R,k)), b = mului(ya[k+1], gel(R,k+1));
2491 15667788 : GEN tj = modii(addii(mului(xa[k],b), mului(xa[k+1],a)), gel(M,j));
2492 15667609 : gel(t, j) = gerepileuptoint(av, tj);
2493 : }
2494 2154027 : if (k==n) gel(t, j) = modii(mului(ya[k], gel(R,k)), gel(M, j));
2495 : } else
2496 : {
2497 292083 : for (j=1, k=1; k<n; j++, k+=2)
2498 : {
2499 182619 : pari_sp av = avma;
2500 182619 : GEN a = mulii(gel(ya,k), gel(R,k)), b = mulii(gel(ya,k+1), gel(R,k+1));
2501 182619 : GEN tj = modii(addii(mulii(gel(xa,k),b), mulii(gel(xa,k+1),a)), gel(M,j));
2502 182619 : gel(t, j) = gerepileuptoint(av, tj);
2503 : }
2504 109464 : if (k==n) gel(t, j) = modii(mulii(gel(ya,k), gel(R,k)), gel(M, j));
2505 : }
2506 2262828 : gel(Tp, 1) = t;
2507 8455867 : for (i=2; i<=m; i++)
2508 : {
2509 6199585 : GEN u = gel(T, i-1), M = gel(T, i);
2510 6199585 : GEN t = cgetg(lg(M), t_VEC);
2511 6275312 : GEN v = gel(Tp, i-1);
2512 6275312 : long n = lg(v)-1;
2513 20620535 : for (j=1, k=1; k<n; j++, k+=2)
2514 : {
2515 14427496 : pari_sp av = avma;
2516 57709984 : gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
2517 43282488 : mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
2518 : }
2519 6193039 : if (k==n) gel(t, j) = gel(v, k);
2520 6193039 : gel(Tp, i) = t;
2521 : }
2522 2256282 : return gmael(Tp,m,1);
2523 : }
2524 :
2525 : static GEN
2526 2130399 : ZV_polint_center_tree(GEN T, GEN R, GEN xa, GEN ya, GEN m2)
2527 : {
2528 2130399 : GEN mod = gmael(T, lg(T)-1, 1);
2529 2130399 : GEN a = ZV_polint_tree(T, R, xa, ya);
2530 2099299 : return Fp_center(a, mod, m2);
2531 : }
2532 :
2533 : static GEN
2534 40803 : ncV_polint_center_tree(GEN T, GEN R, GEN xa, GEN Va, GEN m2)
2535 : {
2536 40803 : long i, j, l = lg(gel(Va,1)), n = lg(xa);
2537 40803 : GEN V = cgetg(l, t_COL);
2538 2165667 : for(i=1; i < l; i++)
2539 : {
2540 2124665 : pari_sp av = avma;
2541 2124665 : GEN ya = cgetg(n, t_VECSMALL);
2542 34929493 : for(j=1; j < n; j++)
2543 32795733 : ya[j] = mael(Va,j,i);
2544 2133760 : gel(V,i) = gerepilecopy(av, ZV_polint_center_tree(T, R, xa, ya, m2));
2545 : }
2546 41002 : return V;
2547 : }
2548 :
2549 : GEN
2550 38014 : nmV_polint_center_tree_worker(GEN Va, GEN T, GEN R, GEN xa, GEN m2)
2551 : {
2552 38014 : return ncV_polint_center_tree(T, R, xa, Va, m2);
2553 : }
2554 :
2555 : static GEN
2556 2090 : nmV_polint_center_tree(GEN T, GEN R, GEN xa, GEN Ma, GEN m2)
2557 : {
2558 2090 : long i, j, l = lg(gel(Ma,1)), n = lg(xa);
2559 2090 : long pending = 0, workid, cnt = 0;
2560 : struct pari_mt pt;
2561 : GEN worker, done, va, M;
2562 2090 : GEN ya = cgetg(n, t_VEC);
2563 2090 : worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, xa, m2));
2564 2090 : va = mkvec(gen_0);
2565 2090 : M = cgetg(l, t_MAT);
2566 2090 : if (DEBUGLEVEL>2) err_printf("Start parallel Chinese remainder: ");
2567 2090 : mt_queue_start(&pt, worker);
2568 43765 : for (i=1; i<l || pending; i++)
2569 : {
2570 611825 : for(j=1; j < n; j++)
2571 570150 : gel(ya,j) = gmael(Ma,j,i);
2572 41675 : gel(va, 1) = ya;
2573 41675 : mt_queue_submit(&pt, i, i<l? va: NULL);
2574 41675 : done = mt_queue_get(&pt, &workid, &pending);
2575 41675 : if (done)
2576 : {
2577 38256 : gel(M,workid) = done;
2578 38256 : if (DEBUGLEVEL>2) err_printf("%ld%% ",(++cnt)*100/(l-1));
2579 : }
2580 : }
2581 2090 : if (DEBUGLEVEL>2) err_printf("\n");
2582 2090 : mt_queue_end(&pt);
2583 2090 : return M;
2584 : }
2585 :
2586 : GEN
2587 0 : Z_ZV_mod(GEN P, GEN xa)
2588 : {
2589 0 : pari_sp av = avma;
2590 0 : GEN T = ZV_producttree(xa);
2591 0 : return gerepileuptoleaf(av, Z_ZV_mod_tree(P, xa, T));
2592 : }
2593 :
2594 : GEN
2595 0 : Z_nv_mod(GEN P, GEN xa)
2596 : {
2597 0 : return Z_ZV_mod(P, xa);
2598 : }
2599 :
2600 : GEN
2601 60827 : ZX_nv_mod_tree(GEN P, GEN xa, GEN T)
2602 : {
2603 60827 : long i, j, l = lg(P), n = lg(xa)-1;
2604 60827 : GEN V = cgetg(n+1, t_VEC);
2605 263192 : for (j=1; j <= n; j++)
2606 : {
2607 202350 : gel(V, j) = cgetg(l, t_VECSMALL);
2608 202349 : mael(V, j, 1) = P[1]&VARNBITS;
2609 : }
2610 769607 : for (i=2; i < l; i++)
2611 : {
2612 708774 : GEN v = Z_ZV_mod_tree(gel(P, i), xa, T);
2613 3412946 : for (j=1; j <= n; j++)
2614 2704181 : mael(V, j, i) = v[j];
2615 : }
2616 263197 : for (j=1; j <= n; j++)
2617 202366 : (void) Flx_renormalize(gel(V, j), l);
2618 60831 : return V;
2619 : }
2620 :
2621 : static GEN
2622 161910 : ZV_sqr(GEN z)
2623 : {
2624 161910 : long i,l = lg(z);
2625 161910 : GEN x = cgetg(l, t_VEC);
2626 161911 : if (typ(z)==t_VECSMALL)
2627 52447 : for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
2628 : else
2629 109464 : for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
2630 161909 : return x;
2631 : }
2632 :
2633 : static GEN
2634 1062370 : ZT_sqr(GEN z)
2635 : {
2636 1062370 : if (typ(z) == t_INT)
2637 611252 : return sqri(z);
2638 : else
2639 : {
2640 451118 : long i,l = lg(z);
2641 451118 : GEN x = cgetg(l, t_VEC);
2642 451120 : for (i=1; i<l; i++) gel(x,i) = ZT_sqr(gel(z,i));
2643 451124 : return x;
2644 : }
2645 : }
2646 :
2647 : static GEN
2648 161910 : ZV_invdivexact(GEN y, GEN x)
2649 : {
2650 161910 : long i, l = lg(y);
2651 161910 : GEN z = cgetg(l,t_VEC);
2652 161912 : if (typ(x)==t_VECSMALL)
2653 342974 : for (i=1; i<l; i++)
2654 : {
2655 290526 : pari_sp av = avma;
2656 290526 : ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
2657 290524 : avma = av;
2658 290524 : gel(z,i) = utoi(a);
2659 : }
2660 : else
2661 509880 : for (i=1; i<l; i++)
2662 400416 : gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
2663 161912 : return z;
2664 : }
2665 :
2666 : static GEN
2667 161908 : ZV_chinesetree(GEN T, GEN xa)
2668 : {
2669 161908 : GEN T2 = ZT_sqr(T), xa2 = ZV_sqr(xa);
2670 161910 : GEN mod = gmael(T,lg(T)-1,1);
2671 161910 : return ZV_invdivexact(Z_ZV_mod_tree(mod, xa2, T2), xa);
2672 : }
2673 :
2674 : static GEN
2675 161901 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
2676 : {
2677 161901 : if (!pt_mod)
2678 2090 : return gerepileupto(av, a);
2679 : else
2680 : {
2681 159811 : GEN mod = gmael(T, lg(T)-1, 1);
2682 159811 : gerepileall(av, 2, &a, &mod);
2683 159817 : *pt_mod = mod;
2684 159817 : return a;
2685 : }
2686 : }
2687 :
2688 : GEN
2689 47559 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN *pt_mod)
2690 : {
2691 47559 : pari_sp av = avma;
2692 47559 : GEN R = ZV_chinesetree(T, P);
2693 47562 : GEN a = ZV_polint_tree(T, R, P, A);
2694 47555 : return gc_chinese(av, T, a, pt_mod);
2695 : }
2696 :
2697 : GEN
2698 109464 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
2699 : {
2700 109464 : pari_sp av = avma;
2701 109464 : GEN T = ZV_producttree(P);
2702 109464 : GEN R = ZV_chinesetree(T, P);
2703 109464 : GEN a = ZV_polint_tree(T, R, P, A);
2704 109464 : return gc_chinese(av, T, a, pt_mod);
2705 : }
2706 :
2707 : GEN
2708 2794 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2709 : {
2710 2794 : pari_sp av = avma;
2711 2794 : GEN T = ZV_producttree(P);
2712 2794 : GEN R = ZV_chinesetree(T, P);
2713 2794 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2714 2794 : GEN a = ncV_polint_center_tree(T, R, P, A, m2);
2715 2794 : return gc_chinese(av, T, a, pt_mod);
2716 : }
2717 :
2718 : GEN
2719 2090 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
2720 : {
2721 2090 : pari_sp av = avma;
2722 2090 : GEN T = ZV_producttree(P);
2723 2090 : GEN R = ZV_chinesetree(T, P);
2724 2090 : GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
2725 2090 : GEN a = nmV_polint_center_tree(T, R, P, A, m2);
2726 2090 : return gc_chinese(av, T, a, pt_mod);
2727 : }
2728 :
2729 : /**********************************************************************
2730 : ** **
2731 : ** Powering over (Z/NZ)^*, small N **
2732 : ** **
2733 : **********************************************************************/
2734 :
2735 : ulong
2736 80311218 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
2737 : {
2738 : ulong y, z, n;
2739 80311218 : if (n0 <= 1)
2740 : { /* frequent special cases */
2741 9176488 : if (n0 == 1) return x;
2742 3485951 : if (n0 == 0) return 1;
2743 : }
2744 71134730 : if (x <= 1) return x; /* 0 or 1 */
2745 70117371 : y = 1; z = x; n = n0;
2746 : for(;;)
2747 : {
2748 819507426 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
2749 819412634 : n>>=1; if (!n) return y;
2750 749309309 : z = Fl_sqr_pre(z,p,pi);
2751 749390055 : }
2752 : }
2753 :
2754 : ulong
2755 17002709 : Fl_powu(ulong x, ulong n0, ulong p)
2756 : {
2757 : ulong y, z, n;
2758 17002709 : if (n0 <= 2)
2759 : { /* frequent special cases */
2760 8460050 : if (n0 == 2) return Fl_sqr(x,p);
2761 1951646 : if (n0 == 1) return x;
2762 2638 : if (n0 == 0) return 1;
2763 : }
2764 8542659 : if (x <= 1) return x; /* 0 or 1 */
2765 8506556 : if (!SMALL_ULONG(p))
2766 4433006 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
2767 4073550 : y = 1; z = x; n = n0;
2768 : for(;;)
2769 : {
2770 52542855 : if (n&1) y = Fl_mul(y,z,p);
2771 52564515 : n>>=1; if (!n) return y;
2772 48483909 : z = Fl_sqr(z,p);
2773 48469305 : }
2774 : }
2775 :
2776 : /* Reduce data dependency to maximize internal parallelism */
2777 : GEN
2778 8285426 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
2779 : {
2780 : long i, k;
2781 8285426 : GEN powers = cgetg(n + 2, t_VECSMALL);
2782 8295960 : powers[1] = 1; if (n == 0) return powers;
2783 8295960 : powers[2] = x;
2784 43857389 : for (i = 3, k=2; i <= n; i+=2, k++)
2785 : {
2786 35576024 : powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
2787 35563630 : powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
2788 : }
2789 8281365 : if (i==n+1)
2790 7344316 : powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
2791 8281520 : return powers;
2792 : }
2793 :
2794 : GEN
2795 1520 : Fl_powers(ulong x, long n, ulong p)
2796 : {
2797 1520 : return Fl_powers_pre(x, n, p, get_Fl_red(p));
2798 : }
2799 :
2800 : /**********************************************************************
2801 : ** **
2802 : ** Powering over (Z/NZ)^*, large N **
2803 : ** **
2804 : **********************************************************************/
2805 :
2806 : static GEN
2807 388111 : Fp_dblsqr(GEN x, GEN N)
2808 : {
2809 388111 : GEN z = shifti(Fp_sqr(x, N), 1);
2810 388111 : return cmpii(z, N) >= 0? subii(z, N): z;
2811 : }
2812 :
2813 : typedef struct muldata {
2814 : GEN (*sqr)(void * E, GEN x);
2815 : GEN (*mul)(void * E, GEN x, GEN y);
2816 : GEN (*mul2)(void * E, GEN x);
2817 : } muldata;
2818 :
2819 : /* modified Barrett reduction with one fold */
2820 : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
2821 :
2822 : static GEN
2823 5248 : Fp_invmBarrett(GEN p, long s)
2824 : {
2825 5248 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
2826 5248 : return mkvec2(Q,R);
2827 : }
2828 :
2829 : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
2830 : * a = r (mod N) */
2831 : static GEN
2832 317929 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
2833 : {
2834 317929 : pari_sp av = avma;
2835 317929 : GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
2836 317929 : long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
2837 317929 : GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
2838 317929 : GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
2839 317929 : GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
2840 317929 : GEN r = subii(A, mulii(q, N));
2841 317929 : GEN sr= subii(r,N); /* 0 <= r < 4*N */
2842 317929 : if (signe(sr)<0) return gerepileuptoint(av, r);
2843 305303 : r=sr; sr = subii(r,N); /* 0 <= r < 3*N */
2844 305303 : if (signe(sr)<0) return gerepileuptoint(av, r);
2845 17459 : r=sr; sr = subii(r,N); /* 0 <= r < 2*N */
2846 17459 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
2847 : }
2848 :
2849 : /* Montgomery reduction */
2850 :
2851 : INLINE ulong
2852 253908 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
2853 :
2854 : struct montred
2855 : {
2856 : GEN N;
2857 : ulong inv;
2858 : };
2859 :
2860 : /* Montgomery reduction */
2861 : static GEN
2862 11895234 : _sqr_montred(void * E, GEN x)
2863 : {
2864 11895234 : struct montred * D = (struct montred *) E;
2865 11895234 : return red_montgomery(sqri(x), D->N, D->inv);
2866 : }
2867 :
2868 : /* Montgomery reduction */
2869 : static GEN
2870 1251024 : _mul_montred(void * E, GEN x, GEN y)
2871 : {
2872 1251024 : struct montred * D = (struct montred *) E;
2873 1251024 : return red_montgomery(mulii(x, y), D->N, D->inv);
2874 : }
2875 :
2876 : static GEN
2877 816959 : _mul2_montred(void * E, GEN x)
2878 : {
2879 816959 : struct montred * D = (struct montred *) E;
2880 816959 : GEN z = shifti(_sqr_montred(E, x), 1);
2881 816959 : long l = lgefint(D->N);
2882 816959 : while (lgefint(z) > l) z = subii(z, D->N);
2883 816959 : return z;
2884 : }
2885 :
2886 : static GEN
2887 3638726 : _sqr_remii(void* N, GEN x)
2888 3638726 : { return remii(sqri(x), (GEN) N); }
2889 :
2890 : static GEN
2891 243947 : _mul_remii(void* N, GEN x, GEN y)
2892 243947 : { return remii(mulii(x, y), (GEN) N); }
2893 :
2894 : static GEN
2895 383748 : _mul2_remii(void* N, GEN x)
2896 383748 : { return Fp_dblsqr(x, (GEN) N); }
2897 :
2898 : struct redbarrett
2899 : {
2900 : GEN iM, N;
2901 : long s;
2902 : };
2903 :
2904 : static GEN
2905 309231 : _sqr_remiibar(void *E, GEN x)
2906 : {
2907 309231 : struct redbarrett * D = (struct redbarrett *) E;
2908 309231 : return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
2909 : }
2910 :
2911 : static GEN
2912 8698 : _mul_remiibar(void *E, GEN x, GEN y)
2913 : {
2914 8698 : struct redbarrett * D = (struct redbarrett *) E;
2915 8698 : return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
2916 : }
2917 :
2918 : static GEN
2919 4363 : _mul2_remiibar(void *E, GEN x)
2920 : {
2921 4363 : struct redbarrett * D = (struct redbarrett *) E;
2922 4363 : return Fp_dblsqr(x, D->N);
2923 : }
2924 :
2925 : static long
2926 329794 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
2927 : {
2928 329794 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
2929 : {
2930 5248 : struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
2931 5248 : D->sqr = &_sqr_remiibar;
2932 5248 : D->mul = &_mul_remiibar;
2933 5248 : D->mul2 = &_mul2_remiibar;
2934 5248 : E->N = N;
2935 5248 : E->s = 1+(expi(N)>>1);
2936 5248 : E->iM = Fp_invmBarrett(N, E->s);
2937 5248 : *pt_E = (void*) E;
2938 5248 : return 0;
2939 : }
2940 324546 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
2941 : {
2942 253909 : struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
2943 253909 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
2944 253909 : D->sqr = &_sqr_montred;
2945 253909 : D->mul = &_mul_montred;
2946 253909 : D->mul2 = &_mul2_montred;
2947 253909 : E->N = N;
2948 253909 : E->inv = init_montdata(N);
2949 253908 : *pt_E = (void*) E;
2950 253908 : return 1;
2951 : }
2952 : else
2953 : {
2954 70639 : D->sqr = &_sqr_remii;
2955 70639 : D->mul = &_mul_remii;
2956 70639 : D->mul2 = &_mul2_remii;
2957 70639 : *pt_E = (void*) N;
2958 70639 : return 0;
2959 : }
2960 : }
2961 :
2962 : GEN
2963 877607 : Fp_powu(GEN A, ulong k, GEN N)
2964 : {
2965 877607 : long lN = lgefint(N), sA;
2966 : int base_is_2, use_montgomery;
2967 : muldata D;
2968 : void *E;
2969 : pari_sp av;
2970 :
2971 877607 : if (lN == 3) {
2972 117000 : ulong n = uel(N,2);
2973 117000 : return utoi( Fl_powu(umodiu(A, n), k, n) );
2974 : }
2975 760607 : if (k <= 2)
2976 : { /* frequent special cases */
2977 523980 : if (k == 2) return Fp_sqr(A,N);
2978 119477 : if (k == 1) return A;
2979 0 : if (k == 0) return gen_1;
2980 : }
2981 236627 : sA = signe(A)==-1 && odd(k);
2982 236627 : base_is_2 = 0;
2983 236627 : if (lgefint(A) == 3) switch(A[2])
2984 : {
2985 492 : case 1: return sA ? gen_m1 : gen_1;
2986 36828 : case 2: base_is_2 = 1; break;
2987 : }
2988 :
2989 : /* TODO: Move this out of here and use for general modular computations */
2990 236135 : av = avma;
2991 236135 : use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
2992 236135 : if (base_is_2)
2993 36828 : A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
2994 : else
2995 199307 : A = gen_powu_i(A, k, E, D.sqr, D.mul);
2996 236135 : if (use_montgomery)
2997 : {
2998 209799 : A = red_montgomery(A, N, ((struct montred *) E)->inv);
2999 209799 : if (cmpii(A, N) >= 0) A = subii(A,N);
3000 209799 : if (sA) A = subii(N, A);
3001 : }
3002 236135 : return gerepileuptoint(av, A);
3003 : }
3004 :
3005 : GEN
3006 21000 : Fp_pows(GEN A, long k, GEN N)
3007 : {
3008 21000 : if (lgefint(N) == 3) {
3009 6973 : ulong n = N[2];
3010 6973 : ulong a = umodiu(A, n);
3011 6973 : if (k < 0) {
3012 14 : a = Fl_inv(a, n);
3013 14 : k = -k;
3014 : }
3015 6973 : return utoi( Fl_powu(a, (ulong)k, n) );
3016 : }
3017 14027 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
3018 14027 : return Fp_powu(A, (ulong)k, N);
3019 : }
3020 :
3021 : /* A^K mod N */
3022 : GEN
3023 2476626 : Fp_pow(GEN A, GEN K, GEN N)
3024 : {
3025 2476626 : pari_sp av = avma;
3026 2476626 : long t,s, lN = lgefint(N), sA;
3027 : int base_is_2, use_montgomery;
3028 : GEN y;
3029 : muldata D;
3030 : void *E;
3031 :
3032 2476626 : s = signe(K);
3033 2476626 : if (!s)
3034 : {
3035 14603 : t = signe(remii(A,N)); avma = av;
3036 14603 : return t? gen_1: gen_0;
3037 : }
3038 2462023 : if (lN == 3)
3039 : {
3040 2218107 : ulong k, n = N[2], a = umodiu(A, n);
3041 2218107 : if (s < 0) a = Fl_inv(a, n);
3042 2218093 : if (a <= 1) return utoi(a); /* 0 or 1 */
3043 2112300 : if (lgefint(K) > 3)
3044 : { /* silly case : huge exponent, small modulus */
3045 21 : pari_warn(warner, "Mod(a,b)^n with n >> b : wasteful");
3046 21 : if (s > 0)
3047 : {
3048 14 : ulong d = ugcd(a, n);
3049 14 : if (d != 1)
3050 : { /* write n = n1 n2, with n2 maximal such that (n1,a) = 1 */
3051 7 : ulong n1 = ucoprime_part(n, d), n2 = n/n1;
3052 :
3053 7 : k = umodiu(K, eulerphiu(n1));
3054 : /* CRT: = a^K (mod n1), = 0 (mod n2)*/
3055 7 : return utoi( Fl_mul(Fl_powu(a, k, n1), n2 * Fl_inv(n2,n1), n) );
3056 : }
3057 : /* gcd(a,n) = 1 */
3058 7 : k = umodiu(K, eulerphiu(n));
3059 : }
3060 : else
3061 7 : k = umodiu(negi(K), eulerphiu(n));
3062 : }
3063 : else
3064 2112279 : k = uel(K,2);
3065 2112293 : return utoi(Fl_powu(a, k, n));
3066 : }
3067 :
3068 243916 : if (s < 0) y = Fp_inv(A,N);
3069 : else
3070 : {
3071 243537 : y = modii(A,N);
3072 243536 : if (!signe(y)) { avma = av; return gen_0; }
3073 : }
3074 243908 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
3075 :
3076 93720 : base_is_2 = 0;
3077 93720 : sA = signe(y)==-1 && mod2(K);
3078 93720 : if (lgefint(y) == 3) switch(y[2])
3079 : {
3080 61 : case 1: return sA ? gen_m1 : gen_1;
3081 64810 : case 2: base_is_2 = 1; break;
3082 : }
3083 :
3084 : /* TODO: Move this out of here and use for general modular computations */
3085 93659 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
3086 93660 : if (base_is_2)
3087 64811 : y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
3088 : else
3089 28849 : y = gen_pow_i(y, K, E, D.sqr, D.mul);
3090 93661 : if (use_montgomery)
3091 : {
3092 44110 : y = red_montgomery(y, N, ((struct montred *) E)->inv);
3093 44110 : if (cmpii(y,N) >= 0) y = subii(y,N);
3094 44110 : if (sA) y = subii(N, y);
3095 : }
3096 93661 : return gerepileuptoint(av,y);
3097 : }
3098 :
3099 : static GEN
3100 168089 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
3101 :
3102 : static GEN
3103 948 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
3104 :
3105 : static GEN
3106 104 : _Fp_one(void *E) { (void) E; return gen_1; }
3107 :
3108 : GEN
3109 1624 : Fp_powers(GEN x, long n, GEN p)
3110 : {
3111 1624 : if (lgefint(p) == 3)
3112 1520 : return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
3113 104 : return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
3114 : }
3115 :
3116 : static GEN
3117 595811 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
3118 :
3119 : static GEN
3120 960 : _Fp_rand(void *E) { return addis(randomi(subis((GEN)E,1)),1); }
3121 :
3122 : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
3123 :
3124 : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
3125 : equalii,equali1,Fp_easylog};
3126 :
3127 : static GEN
3128 935783 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
3129 :
3130 : static GEN
3131 944352 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
3132 :
3133 : static GEN
3134 103113 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
3135 :
3136 : static GEN
3137 1181019 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
3138 :
3139 : static GEN
3140 18323 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
3141 :
3142 : static int
3143 526253 : _Fp_equal0(GEN x) { return signe(x)==0; }
3144 :
3145 : static GEN
3146 82540 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
3147 :
3148 : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
3149 : _Fp_inv,_Fp_equal0,_Fp_s};
3150 :
3151 4557 : const struct bb_field *get_Fp_field(void **E, GEN p)
3152 : {
3153 4557 : *E = (void*)p; return &Fp_field;
3154 : }
3155 :
3156 : /*********************************************************************/
3157 : /** **/
3158 : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
3159 : /** **/
3160 : /*********************************************************************/
3161 : ulong
3162 5229 : Fl_order(ulong a, ulong o, ulong p)
3163 : {
3164 5229 : pari_sp av = avma;
3165 : GEN m, P, E;
3166 : long i;
3167 5229 : if (!o) o = p-1;
3168 5229 : m = factoru(o);
3169 5229 : P = gel(m,1);
3170 5229 : E = gel(m,2);
3171 13076 : for (i = lg(P)-1; i; i--)
3172 : {
3173 7847 : ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
3174 7847 : if (y == 1) o = t;
3175 7301 : else for (j = 1; j < e; j++)
3176 : {
3177 2261 : y = Fl_powu(y, l, p);
3178 2261 : if (y == 1) { o = t * upowuu(l, j); break; }
3179 : }
3180 : }
3181 5229 : avma = av; return o;
3182 : }
3183 :
3184 : /*Find the exact order of a assuming a^o==1*/
3185 : GEN
3186 10443 : Fp_order(GEN a, GEN o, GEN p) {
3187 10443 : if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
3188 : {
3189 21 : ulong pp = p[2], oo = (o && lgefint(o)==3)? o[2]: pp-1;
3190 21 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
3191 : }
3192 10422 : return gen_order(a, o, (void*)p, &Fp_star);
3193 : }
3194 : GEN
3195 56 : Fp_factored_order(GEN a, GEN o, GEN p)
3196 56 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
3197 :
3198 : /* return order of a mod p^e, e > 0, pe = p^e */
3199 : static GEN
3200 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
3201 : {
3202 : GEN ap, op;
3203 70 : if (absequaliu(p, 2))
3204 : {
3205 56 : if (e == 1) return gen_1;
3206 56 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
3207 49 : if (mod4(a) == 1)
3208 14 : op = gen_1;
3209 : else {
3210 35 : op = gen_2;
3211 35 : a = Fp_sqr(a, pe);
3212 : }
3213 : } else {
3214 14 : ap = (e == 1)? a: remii(a,p);
3215 14 : op = Fp_order(ap, subis(p,1), p);
3216 14 : if (e == 1) return op;
3217 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
3218 : }
3219 49 : if (equali1(a)) return op;
3220 7 : return mulii(op, powiu(p, e - Z_pval(subis(a,1), p)));
3221 : }
3222 :
3223 : GEN
3224 63 : znorder(GEN x, GEN o)
3225 : {
3226 63 : pari_sp av = avma;
3227 : GEN b, a;
3228 :
3229 63 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
3230 56 : b = gel(x,1); a = gel(x,2);
3231 56 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
3232 49 : if (!o)
3233 : {
3234 35 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
3235 35 : long i, l = lg(P);
3236 35 : o = gen_1;
3237 70 : for (i = 1; i < l; i++)
3238 : {
3239 35 : GEN p = gel(P,i);
3240 35 : long e = itos(gel(E,i));
3241 :
3242 35 : if (l == 2)
3243 35 : o = Zp_order(a, p, e, b);
3244 : else {
3245 0 : GEN pe = powiu(p,e);
3246 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
3247 : }
3248 : }
3249 35 : return gerepileuptoint(av, o);
3250 : }
3251 14 : return Fp_order(a, o, b);
3252 : }
3253 : GEN
3254 0 : order(GEN x) { return znorder(x, NULL); }
3255 :
3256 : /*********************************************************************/
3257 : /** **/
3258 : /** DISCRETE LOGARITHM in (Z/nZ)* **/
3259 : /** **/
3260 : /*********************************************************************/
3261 : static GEN
3262 56805 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
3263 : {
3264 56805 : pari_sp av = avma;
3265 : GEN h1, h2, F, G;
3266 56805 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return NULL;
3267 34174 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
3268 : {
3269 231 : GEN M = cgetg(3, t_MAT);
3270 231 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
3271 231 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
3272 231 : return gerepileupto(av, M);
3273 : }
3274 33943 : avma = av; return NULL;
3275 : }
3276 :
3277 : static GEN
3278 56805 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
3279 : {
3280 : GEN rel;
3281 : do
3282 : {
3283 56805 : (*e)++; *g = Fp_mul(*g, b, p);
3284 56805 : rel = Fp_log_halfgcd(bnd, C, *g, p);
3285 56805 : } while (!rel);
3286 231 : return rel;
3287 : }
3288 :
3289 : struct Fp_log_rel
3290 : {
3291 : GEN rel;
3292 : long *sieve;
3293 : ulong prmax;
3294 : long nbrel, nbmax;
3295 : };
3296 :
3297 : /* add u^e */
3298 : static long
3299 37401 : addifsmooth1(struct Fp_log_rel *r, GEN h, long u, long e)
3300 : {
3301 37401 : pari_sp av = avma;
3302 : GEN z;
3303 37401 : if ((z = Z_issmooth_fact(h, r->prmax)))
3304 : {
3305 2485 : long off = r->prmax+1;
3306 2485 : GEN F = cgetg(3, t_MAT);
3307 2485 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
3308 2485 : gel(F,2) = vecsmall_append(gel(z,2), e);
3309 2485 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3310 : }
3311 37401 : return r->nbrel==r->nbmax;
3312 : }
3313 :
3314 : /* add u^-1 v^-1 */
3315 : static long
3316 253358 : addifsmooth2(struct Fp_log_rel *r, GEN h, long u, long v)
3317 : {
3318 253358 : pari_sp av = avma;
3319 : GEN z;
3320 253358 : if ((z = Z_issmooth_fact(h, r->prmax)))
3321 : {
3322 95011 : long off = r->prmax+1;
3323 95011 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
3324 95011 : GEN F = cgetg(3, t_MAT);
3325 95011 : gel(F,1) = vecsmall_concat(gel(z,1), P);
3326 95011 : gel(F,2) = vecsmall_concat(gel(z,2), E);
3327 95011 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
3328 : }
3329 253358 : return r->nbrel==r->nbmax;
3330 : }
3331 :
3332 : /*
3333 : Let p=C^2+c
3334 : Solve h = (C+x)*(C+a)-p = 0 [mod l]
3335 : h= -c+x*(C+a)+C*a = 0 [mod l]
3336 : x = (c-C*a)/(C+a) [mod l]
3337 : h = -c+C*(x+a)+a*x
3338 : */
3339 :
3340 : static void
3341 36862 : Fp_log_sieve_h(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, long a, GEN pr, GEN sz)
3342 : {
3343 36862 : long th = expi(C), n = lg(pr)-1;
3344 : long i,j;
3345 36862 : if (addifsmooth1(r, addis(C,a), a, -1)) return;
3346 25526508 : for(j=0; j<=a; j++)
3347 25489646 : r->sieve[j]=0;
3348 16626820 : for(i=1; i<=n; i++)
3349 : {
3350 16589958 : ulong li = pr[i], s = sz[i], al = a % li;
3351 16589958 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
3352 16589958 : if (!iv) continue;
3353 16504691 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
3354 65360344 : for(j = u; j<=a; j+=li)
3355 48855653 : r->sieve[j] += s;
3356 : }
3357 36862 : th = th - expu(th)-1;
3358 25465797 : for(j=0; j<a; j++)
3359 25428984 : if (r->sieve[j]>=th)
3360 : {
3361 253358 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
3362 253358 : if (addifsmooth2(r, h, a, j)) return;
3363 : }
3364 : /* j = a */
3365 36813 : if (r->sieve[a]>=th)
3366 : {
3367 539 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
3368 539 : if (addifsmooth1(r, h, a, -2)) return;
3369 : }
3370 : }
3371 :
3372 : static GEN
3373 1367 : _psi(void*E, GEN y)
3374 : {
3375 1367 : GEN lx = (GEN) E;
3376 1367 : long prec = realprec(lx);
3377 1367 : GEN ly = glog(y, prec);
3378 1367 : GEN u = gdiv(lx, ly);
3379 1367 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
3380 : }
3381 :
3382 : static GEN
3383 49 : opt_param(GEN x, long prec)
3384 : {
3385 49 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
3386 : }
3387 :
3388 : static GEN
3389 49 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
3390 : {
3391 49 : pari_sp av = avma;
3392 49 : long lM = lg(M)-1, nbcol = lM;
3393 : for (;;)
3394 : {
3395 84 : GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
3396 84 : long i, f=0;
3397 84 : long l = lg(K), lm = lgefint(m);
3398 84 : GEN idx = diviiexact(subis(p,1),m), g;
3399 : pari_timer ti;
3400 84 : if (DEBUGLEVEL) timer_start(&ti);
3401 133 : for(i=1; i<l; i++)
3402 133 : if (signe(gel(K,i)))
3403 84 : break;
3404 84 : g = Fp_pow(utoi(i), idx, p);
3405 84 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
3406 129185 : for(i=1; i<l; i++)
3407 : {
3408 129101 : GEN k = gel(K,i);
3409 129101 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
3410 129101 : if (signe(k)==0 || !equalii(Fp_pow(g, k, p),
3411 : Fp_pow(j, idx, p)))
3412 80234 : gel(K,i) = cgetineg(lm);
3413 : else
3414 48867 : f++;
3415 : }
3416 84 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld logs", f);
3417 133 : if(f > (nbg>>1)) return gerepileupto(av, K);
3418 9198 : for(i=1; i<=nbcol; i++)
3419 : {
3420 9163 : long a = 1+random_Fl(lM);
3421 9163 : swap(gel(M,a),gel(M,i));
3422 : }
3423 35 : if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
3424 35 : }
3425 : }
3426 :
3427 : static GEN
3428 98 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
3429 : {
3430 98 : pari_sp av=avma;
3431 98 : GEN aa = gen_1;
3432 98 : long AV = 0;
3433 : for(;;)
3434 : {
3435 231 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
3436 231 : GEN F = gel(A,1), E = gel(A,2);
3437 231 : GEN Ao = gen_0;
3438 231 : long i, l = lg(F);
3439 1162 : for(i=1; i<l; i++)
3440 : {
3441 1064 : GEN Ki = gel(K,F[i]);
3442 1064 : if (signe(Ki)<0) break;
3443 931 : Ao = addii(Ao, mulis(Ki, E[i]));
3444 : }
3445 329 : if (i==l) return Fp_div(Ao, utoi(AV), m);
3446 133 : aa = gerepileuptoint(av, aa);
3447 133 : }
3448 : }
3449 :
3450 : static GEN
3451 49 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
3452 : {
3453 49 : pari_sp av = avma, av2;
3454 : long i, nbi, nbrow;
3455 : GEN C, c, Ci, ci, pr, sz, l, Ao, Bo, K, d, p_1;
3456 : pari_timer ti;
3457 : struct Fp_log_rel r;
3458 49 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
3459 49 : ulong bnd = 4*bnds;
3460 49 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
3461 :
3462 49 : p_1 = subiu(p,1);
3463 49 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
3464 0 : m = diviiexact(p_1, coprime_part(p_1, m));
3465 49 : pr = primes_upto_zv(bnd);
3466 49 : nbi = lg(pr)-1;
3467 49 : if (DEBUGLEVEL)
3468 : {
3469 0 : err_printf("bnd=%lu Size FB=%ld\n", bnd, nbi);
3470 0 : timer_start(&ti);
3471 : }
3472 49 : C = sqrtremi(p, &c);
3473 49 : av2 = avma;
3474 49 : Ci = cgetg(nbi+1,t_VECSMALL);
3475 49 : ci = cgetg(nbi+1,t_VECSMALL);
3476 49 : sz = cgetg(nbi+1,t_VECSMALL);
3477 12236 : for (i = 1; i <= nbi; ++i)
3478 : {
3479 12187 : ulong lp = pr[i];
3480 12187 : Ci[i] = umodiu(C, lp);
3481 12187 : ci[i] = umodiu(c, lp);
3482 12187 : sz[i] = expu(lp);
3483 : }
3484 49 : r.nbrel = 0;
3485 49 : r.nbmax = 8*nbi;
3486 49 : r.rel = cgetg(r.nbmax+1,t_VEC);
3487 49 : r.sieve = cgetg(r.nbmax+2,t_VECSMALL)+1;
3488 49 : r.prmax = pr[nbi];
3489 36911 : for(i=0; r.nbrel < r.nbmax; i++)
3490 : {
3491 36862 : Fp_log_sieve_h(&r, C, c, Ci, ci, i, pr, sz);
3492 36862 : if (DEBUGLEVEL && (i&127)==0)
3493 0 : err_printf("%ld%% ",100*r.nbrel/(r.nbmax));
3494 : }
3495 49 : nbrow = r.prmax+i;
3496 49 : if (DEBUGLEVEL)
3497 : {
3498 0 : err_printf("\n");
3499 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+i);
3500 : }
3501 49 : setlg(r.rel,r.nbrel+1);
3502 49 : r.rel = gerepileupto(av2, r.rel);
3503 49 : K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
3504 49 : if (DEBUGLEVEL) timer_start(&ti);
3505 49 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
3506 49 : if (DEBUGLEVEL) timer_printf(&ti," log element");
3507 49 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
3508 49 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
3509 49 : d = gcdii(Ao,Bo);
3510 49 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
3511 49 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
3512 49 : return gerepileuptoint(av, l);
3513 : }
3514 :
3515 : static int
3516 194188 : Fp_log_use_index(long e, long p)
3517 : {
3518 194188 : return (e >= 27 && 20*(p+6)<=e*e);
3519 : }
3520 :
3521 : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
3522 : static GEN
3523 213594 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
3524 : {
3525 213594 : pari_sp av = avma;
3526 213594 : GEN p = (GEN)E;
3527 : /* assume a reduced mod p, p not necessarily prime */
3528 213594 : if (equali1(a)) return gen_0;
3529 : /* p > 2 */
3530 145150 : if (equalii(subis(p,1), a)) /* -1 */
3531 : {
3532 : pari_sp av2;
3533 : GEN t;
3534 60371 : ord = get_arith_Z(ord);
3535 60371 : if (mpodd(ord)) { avma = av; return cgetg(1, t_VEC); } /* no solution */
3536 60357 : t = shifti(ord,-1); /* only possible solution */
3537 60357 : av2 = avma;
3538 60357 : if (!equalii(Fp_pow(g, t, p), a)) { avma = av; return cgetg(1, t_VEC); }
3539 60175 : avma = av2; return gerepileuptoint(av, t);
3540 : }
3541 84779 : if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
3542 49 : return Fp_log_index(a, g, ord, p);
3543 84730 : avma = av; return NULL; /* not easy */
3544 : }
3545 :
3546 : GEN
3547 209882 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
3548 : {
3549 209882 : GEN v = get_arith_ZZM(ord);
3550 209854 : GEN F = gmael(v,2,1);
3551 209854 : long lF = lg(F)-1, lmax;
3552 209854 : if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
3553 158614 : lmax = expi(gel(F,lF));
3554 158614 : if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
3555 49 : v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
3556 158614 : return gen_PH_log(a,g,v,(void*)p,&Fp_star);
3557 : }
3558 :
3559 : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
3560 : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
3561 : static GEN
3562 112 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
3563 : {
3564 112 : long l = lg(P) - 1, e = E[l];
3565 112 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
3566 : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
3567 :
3568 112 : if (l == 1) {
3569 84 : hpe = h;
3570 84 : gpe = g;
3571 : } else {
3572 28 : hpe = modii(h, pe);
3573 28 : gpe = modii(g, pe);
3574 : }
3575 112 : if (e == 1) {
3576 28 : hp = hpe;
3577 28 : gp = gpe;
3578 : } else {
3579 84 : hp = remii(hpe, p);
3580 84 : gp = remii(gpe, p);
3581 : }
3582 112 : if (hp == gen_0 || gp == gen_0) return NULL;
3583 91 : if (absequaliu(p, 2))
3584 : {
3585 35 : GEN n = int2n(e);
3586 35 : ogpe = Zp_order(gpe, gen_2, e, n);
3587 35 : a = Fp_log(hpe, gpe, ogpe, n);
3588 35 : if (typ(a) != t_INT) return NULL;
3589 : }
3590 : else
3591 : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
3592 : is trivial */
3593 : /* [order(gp), factor(order(gp))] */
3594 56 : GEN v = Fp_factored_order(gp, subis(p,1), p);
3595 56 : GEN ogp = gel(v,1);
3596 56 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
3597 56 : a = Fp_log(hp, gp, v, p);
3598 56 : if (typ(a) != t_INT) return NULL;
3599 56 : if (e == 1) ogpe = ogp;
3600 : else
3601 : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
3602 : /* use p-adic log: O(log p + e) mul*/
3603 : long vpogpe, vpohpe;
3604 :
3605 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
3606 28 : gpe = Fp_pow(gpe, ogp, pe);
3607 : /* g,h = 1 mod p; compute b s.t. h = g^b */
3608 :
3609 : /* v_p(order g mod pe) */
3610 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subis(gpe,1), p);
3611 : /* v_p(order h mod pe) */
3612 28 : vpohpe = equali1(hpe)? 0: e - Z_pval(subis(hpe,1), p);
3613 28 : if (vpohpe > vpogpe) return NULL;
3614 :
3615 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
3616 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
3617 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
3618 28 : a = addii(a, mulii(ogp, padic_to_Q(b)));
3619 : }
3620 : }
3621 : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
3622 77 : if (l == 1) return a;
3623 :
3624 28 : N = diviiexact(N, pe); /* make N coprime to p */
3625 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
3626 28 : g = Fp_pow(g, modii(ogpe, phi), N);
3627 28 : setlg(P, l); /* remove last element */
3628 28 : setlg(E, l);
3629 28 : b = znlog_rec(h, g, N, P, E, PHI);
3630 28 : if (!b) return NULL;
3631 28 : return addmulii(a, b, ogpe);
3632 : }
3633 :
3634 : static GEN
3635 84 : get_PHI(GEN P, GEN E)
3636 : {
3637 84 : long i, l = lg(P);
3638 84 : GEN PHI = cgetg(l, t_VEC);
3639 84 : gel(PHI,1) = gen_1;
3640 112 : for (i=1; i<l-1; i++)
3641 : {
3642 28 : GEN t, p = gel(P,i);
3643 28 : long e = E[i];
3644 28 : t = mulii(powiu(p, e-1), subis(p,1));
3645 28 : if (i > 1) t = mulii(t, gel(PHI,i));
3646 28 : gel(PHI,i+1) = t;
3647 : }
3648 84 : return PHI;
3649 : }
3650 :
3651 : GEN
3652 224 : znlog(GEN h, GEN g, GEN o)
3653 : {
3654 224 : pari_sp av = avma;
3655 : GEN N, fa, P, E, x;
3656 224 : switch (typ(g))
3657 : {
3658 : case t_PADIC:
3659 : {
3660 28 : GEN p = gel(g,2);
3661 28 : long v = valp(g);
3662 28 : if (v < 0) pari_err_DIM("znlog");
3663 28 : if (v > 0) {
3664 0 : long k = gvaluation(h, p);
3665 0 : if (k % v) return cgetg(1,t_VEC);
3666 0 : k /= v;
3667 0 : if (!gequal(h, gpowgs(g,k))) { avma = av; return cgetg(1,t_VEC); }
3668 0 : avma = av; return stoi(k);
3669 : }
3670 28 : N = gel(g,3);
3671 28 : g = Rg_to_Fp(g, N);
3672 28 : break;
3673 : }
3674 : case t_INTMOD:
3675 189 : N = gel(g,1);
3676 189 : g = gel(g,2); break;
3677 7 : default: pari_err_TYPE("znlog", g);
3678 0 : return NULL; /* not reached */
3679 : }
3680 217 : if (equali1(N)) { avma = av; return gen_0; }
3681 217 : h = Rg_to_Fp(h, N);
3682 210 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
3683 84 : fa = Z_factor(N);
3684 84 : P = gel(fa,1);
3685 84 : E = vec_to_vecsmall(gel(fa,2));
3686 84 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
3687 84 : if (!x) { avma = av; return cgetg(1,t_VEC); }
3688 49 : return gerepileuptoint(av, x);
3689 : }
3690 :
3691 : GEN
3692 63806 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
3693 : {
3694 63806 : a = modii(a,p);
3695 63806 : if (!signe(a))
3696 : {
3697 48034 : if (zeta) *zeta = gen_1;
3698 48034 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
3699 48027 : return gen_0;
3700 : }
3701 15772 : if (absequaliu(n,2))
3702 : {
3703 2954 : if (zeta) *zeta = addis(p,-1);
3704 2954 : return Fp_sqrt(a,p);
3705 : }
3706 12818 : return gen_Shanks_sqrtn(a,n,addis(p,-1),zeta,(void*)p,&Fp_star);
3707 : }
3708 :
3709 : /*********************************************************************/
3710 : /** **/
3711 : /** FUNDAMENTAL DISCRIMINANTS **/
3712 : /** **/
3713 : /*********************************************************************/
3714 : long
3715 14371 : isfundamental(GEN x) {
3716 14371 : if (typ(x) != t_INT) pari_err_TYPE("isfundamental",x);
3717 14371 : return Z_isfundamental(x);
3718 : }
3719 :
3720 : /* x fundamental ? */
3721 : long
3722 7209 : uposisfundamental(ulong x)
3723 : {
3724 7209 : ulong r = x & 15; /* x mod 16 */
3725 7209 : if (!r) return 0;
3726 6817 : switch(r & 3)
3727 : { /* x mod 4 */
3728 1303 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
3729 2234 : case 1: return uissquarefree(x);
3730 3280 : default: return 0;
3731 : }
3732 : }
3733 : /* -x fundamental ? */
3734 : long
3735 12545 : unegisfundamental(ulong x)
3736 : {
3737 12545 : ulong r = x & 15; /* x mod 16 */
3738 12545 : if (!r) return 0;
3739 12020 : switch(r & 3)
3740 : { /* x mod 4 */
3741 1982 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
3742 6912 : case 3: return uissquarefree(x);
3743 3126 : default: return 0;
3744 : }
3745 : }
3746 : long
3747 1078 : sisfundamental(long x)
3748 1078 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
3749 :
3750 : long
3751 14931 : Z_isfundamental(GEN x)
3752 : {
3753 : long r;
3754 14931 : switch(lgefint(x))
3755 : {
3756 0 : case 2: return 0;
3757 32415 : case 3: return signe(x) < 0? unegisfundamental(x[2])
3758 19494 : : uposisfundamental(x[2]);
3759 : }
3760 2010 : r = mod16(x);
3761 2010 : if (!r) return 0;
3762 1884 : if ((r & 3) == 0)
3763 : {
3764 : pari_sp av;
3765 376 : r >>= 2; /* |x|/4 mod 4 */
3766 376 : if (signe(x) < 0) r = 4-r;
3767 376 : if (r == 1) return 0;
3768 250 : av = avma;
3769 250 : r = Z_issquarefree( shifti(x,-2) );
3770 250 : avma = av; return r;
3771 : }
3772 1508 : r &= 3; /* |x| mod 4 */
3773 1508 : if (signe(x) < 0) r = 4-r;
3774 1508 : return (r==1) ? Z_issquarefree(x) : 0;
3775 : }
3776 :
3777 : GEN
3778 7 : quaddisc(GEN x)
3779 : {
3780 7 : const pari_sp av = avma;
3781 7 : long i,r,tx=typ(x);
3782 : GEN P,E,f,s;
3783 :
3784 7 : if (!is_rational_t(tx)) pari_err_TYPE("quaddisc",x);
3785 7 : f = factor(x);
3786 7 : P = gel(f,1);
3787 7 : E = gel(f,2); s = gen_1;
3788 35 : for (i=1; i<lg(P); i++)
3789 28 : if (odd(mael(E,i,2))) s = mulii(s,gel(P,i));
3790 7 : r = mod4(s); if (gsigne(x) < 0) r = 4-r;
3791 7 : if (r>1) s = shifti(s,2);
3792 7 : return gerepileuptoint(av, s);
3793 : }
3794 :
3795 : /*********************************************************************/
3796 : /** **/
3797 : /** FACTORIAL **/
3798 : /** **/
3799 : /*********************************************************************/
3800 : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
3801 : * first is slower ... ] */
3802 : GEN
3803 993340 : mulu_interval(ulong a, ulong b)
3804 : {
3805 993340 : pari_sp av = avma;
3806 : ulong k, l, N, n;
3807 : long lx;
3808 : GEN x;
3809 :
3810 993340 : if (!a) return gen_0;
3811 993340 : n = b - a + 1;
3812 993340 : if (n < 61)
3813 : {
3814 988313 : x = utoi(a);
3815 988314 : for (k=a+1; k<=b; k++) x = mului(k,x);
3816 988313 : return gerepileuptoint(av, x);
3817 : }
3818 5027 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3819 5027 : N = b + a;
3820 722381 : for (k = a;; k++)
3821 : {
3822 722381 : l = N - k; if (l <= k) break;
3823 717354 : gel(x,lx++) = muluu(k,l);
3824 717354 : }
3825 5027 : if (l == k) gel(x,lx++) = utoipos(k);
3826 5027 : setlg(x, lx);
3827 5027 : return gerepileuptoint(av, ZV_prod(x));
3828 : }
3829 : GEN
3830 245 : muls_interval(long a, long b)
3831 : {
3832 245 : pari_sp av = avma;
3833 245 : long lx, k, l, N, n = b - a + 1;
3834 : GEN x;
3835 :
3836 245 : if (a <= 0 && b >= 0) return gen_0;
3837 98 : if (n < 61)
3838 : {
3839 98 : x = stoi(a);
3840 98 : for (k=a+1; k<=b; k++) x = mulsi(k,x);
3841 98 : return gerepileuptoint(av, x);
3842 : }
3843 0 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3844 0 : N = b + a;
3845 0 : for (k = a;; k++)
3846 : {
3847 0 : l = N - k; if (l <= k) break;
3848 0 : gel(x,lx++) = mulss(k,l);
3849 0 : }
3850 0 : if (l == k) gel(x,lx++) = stoi(k);
3851 0 : setlg(x, lx);
3852 0 : return gerepileuptoint(av, ZV_prod(x));
3853 : }
3854 :
3855 : GEN
3856 2324454 : mpfact(long n)
3857 : {
3858 2324454 : if (n < 2)
3859 : {
3860 1451696 : if (n < 0) pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
3861 1451696 : return gen_1;
3862 : }
3863 872758 : return mulu_interval(2UL, (ulong)n);
3864 : }
3865 :
3866 : /*******************************************************************/
3867 : /** **/
3868 : /** LUCAS & FIBONACCI **/
3869 : /** **/
3870 : /*******************************************************************/
3871 : static void
3872 56 : lucas(ulong n, GEN *a, GEN *b)
3873 : {
3874 : GEN z, t, zt;
3875 112 : if (!n) { *a = gen_2; *b = gen_1; return; }
3876 49 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
3877 49 : switch(n & 3) {
3878 14 : case 0: *a = addsi(-2,sqri(z)); *b = addsi(-1,zt); break;
3879 14 : case 1: *a = addsi(-1,zt); *b = addsi(2,sqri(t)); break;
3880 7 : case 2: *a = addsi(2,sqri(z)); *b = addsi(1,zt); break;
3881 14 : case 3: *a = addsi(1,zt); *b = addsi(-2,sqri(t));
3882 : }
3883 : }
3884 :
3885 : GEN
3886 7 : fibo(long n)
3887 : {
3888 7 : pari_sp av = avma;
3889 : GEN a, b;
3890 7 : if (!n) return gen_0;
3891 7 : lucas((ulong)(labs(n)-1), &a, &b);
3892 7 : a = diviuexact(addii(shifti(a,1),b), 5);
3893 7 : if (n < 0 && !odd(n)) setsigne(a, -1);
3894 7 : return gerepileuptoint(av, a);
3895 : }
3896 :
3897 : /*******************************************************************/
3898 : /* */
3899 : /* CONTINUED FRACTIONS */
3900 : /* */
3901 : /*******************************************************************/
3902 : static GEN
3903 359179 : icopy_lg(GEN x, long l)
3904 : {
3905 359179 : long lx = lgefint(x);
3906 : GEN y;
3907 :
3908 359179 : if (lx >= l) return icopy(x);
3909 35 : y = cgeti(l); affii(x, y); return y;
3910 : }
3911 :
3912 : /* continued fraction of a/b. If y != NULL, stop when partial quotients
3913 : * differ from y */
3914 : static GEN
3915 359468 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
3916 : {
3917 : GEN z, c;
3918 359468 : ulong i, l, ly = lgefint(b);
3919 :
3920 : /* times 1 / log2( (1+sqrt(5)) / 2 ) */
3921 359468 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
3922 359468 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
3923 359468 : if (l > LGBITS) l = LGBITS;
3924 :
3925 359468 : z = cgetg(l,t_VEC);
3926 359468 : l--;
3927 359468 : if (y) {
3928 289 : pari_sp av = avma;
3929 289 : if (l >= (ulong)lg(y)) l = lg(y)-1;
3930 19467 : for (i = 1; i <= l; i++)
3931 : {
3932 19289 : GEN q = gel(y,i);
3933 19289 : gel(z,i) = q;
3934 19289 : c = b; if (!gequal1(q)) c = mulii(q, b);
3935 19289 : c = subii(a, c);
3936 19289 : if (signe(c) < 0)
3937 : { /* partial quotient too large */
3938 110 : c = addii(c, b);
3939 110 : if (signe(c) >= 0) i++; /* by 1 */
3940 110 : break;
3941 : }
3942 19179 : if (cmpii(c, b) >= 0)
3943 : { /* partial quotient too small */
3944 1 : c = subii(c, b);
3945 1 : if (cmpii(c, b) < 0) {
3946 : /* by 1. If next quotient is 1 in y, add 1 */
3947 0 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addis(q,1);
3948 0 : i++;
3949 : }
3950 1 : break;
3951 : }
3952 19178 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
3953 19178 : a = b; b = c;
3954 : }
3955 : } else {
3956 359179 : a = icopy_lg(a, ly);
3957 359179 : b = icopy(b);
3958 1254574 : for (i = 1; i <= l; i++)
3959 : {
3960 1254310 : gel(z,i) = truedvmdii(a,b,&c);
3961 1254310 : if (c == gen_0) { i++; break; }
3962 895395 : affii(c, a); cgiv(c); c = a;
3963 895395 : a = b; b = c;
3964 : }
3965 : }
3966 359468 : i--;
3967 359468 : if (i > 1 && gequal1(gel(z,i)))
3968 : {
3969 85 : cgiv(gel(z,i)); --i;
3970 85 : gel(z,i) = addsi(1, gel(z,i)); /* unclean: leave old z[i] on stack */
3971 : }
3972 359468 : setlg(z,i+1); return z;
3973 : }
3974 :
3975 : static GEN
3976 0 : sersfcont(GEN a, GEN b, long k)
3977 : {
3978 0 : long i, l = typ(a) == t_POL? lg(a): 3;
3979 : GEN y, c;
3980 0 : if (lg(b) > l) l = lg(b);
3981 0 : if (k > 0 && l > k+1) l = k+1;
3982 0 : y = cgetg(l,t_VEC);
3983 0 : for (i=1; i<l; i++)
3984 : {
3985 0 : gel(y,i) = poldivrem(a,b,&c);
3986 0 : if (gequal0(c)) { i++; break; }
3987 0 : a = b; b = c;
3988 : }
3989 0 : setlg(y, i); return y;
3990 : }
3991 :
3992 : GEN
3993 359473 : gboundcf(GEN x, long k)
3994 : {
3995 : pari_sp av;
3996 359473 : long tx = typ(x), e;
3997 : GEN y, a, b, c;
3998 :
3999 359473 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
4000 359466 : if (is_scalar_t(tx))
4001 : {
4002 359466 : if (gequal0(x)) return mkvec(gen_0);
4003 359186 : switch(tx)
4004 : {
4005 0 : case t_INT: return mkveccopy(x);
4006 : case t_REAL:
4007 296 : av = avma;
4008 296 : c = mantissa_real(x,&e);
4009 296 : if (e < 0) pari_err_PREC("gboundcf");
4010 289 : y = int2n(e);
4011 289 : a = Qsfcont(c,y, NULL, k);
4012 289 : b = addsi(signe(x), c);
4013 289 : return gerepilecopy(av, Qsfcont(b,y, a, k));
4014 :
4015 : case t_FRAC:
4016 358890 : av = avma;
4017 358890 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
4018 : }
4019 0 : pari_err_TYPE("gboundcf",x);
4020 : }
4021 :
4022 0 : switch(tx)
4023 : {
4024 0 : case t_POL: return mkveccopy(x);
4025 : case t_SER:
4026 0 : av = avma;
4027 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
4028 : case t_RFRAC:
4029 0 : av = avma;
4030 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
4031 : }
4032 0 : pari_err_TYPE("gboundcf",x);
4033 0 : return NULL; /* not reached */
4034 : }
4035 :
4036 : static GEN
4037 14 : sfcont2(GEN b, GEN x, long k)
4038 : {
4039 14 : pari_sp av = avma;
4040 14 : long lb = lg(b), tx = typ(x), i;
4041 : GEN y,p1;
4042 :
4043 14 : if (k)
4044 : {
4045 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
4046 0 : lb = k+1;
4047 : }
4048 7 : y = cgetg(lb,t_VEC);
4049 7 : if (lb==1) return y;
4050 7 : if (is_scalar_t(tx))
4051 : {
4052 7 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
4053 : }
4054 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
4055 :
4056 7 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
4057 7 : for (i = 1;;)
4058 : {
4059 35 : if (tx == t_REAL)
4060 : {
4061 35 : long e = expo(x);
4062 35 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
4063 35 : gel(y,i) = floorr(x);
4064 35 : p1 = subri(x, gel(y,i));
4065 : }
4066 : else
4067 : {
4068 0 : gel(y,i) = gfloor(x);
4069 0 : p1 = gsub(x, gel(y,i));
4070 : }
4071 35 : if (++i >= lb) break;
4072 28 : if (gequal0(p1)) break;
4073 28 : x = gdiv(gel(b,i),p1);
4074 28 : }
4075 7 : setlg(y,i);
4076 7 : return gerepilecopy(av,y);
4077 : }
4078 :
4079 :
4080 : GEN
4081 0 : gcf(GEN x) { return gboundcf(x,0); }
4082 : GEN
4083 0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
4084 : GEN
4085 49 : contfrac0(GEN x, GEN b, long nmax)
4086 : {
4087 : long tb;
4088 :
4089 49 : if (!b) return gboundcf(x,nmax);
4090 28 : tb = typ(b);
4091 28 : if (tb == t_INT) return gboundcf(x,itos(b));
4092 21 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
4093 21 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
4094 14 : return sfcont2(b,x,nmax);
4095 : }
4096 :
4097 : GEN
4098 126 : contfracpnqn(GEN x, long n)
4099 : {
4100 126 : pari_sp av = avma;
4101 126 : long i, lx = lg(x);
4102 : GEN M,A,B, p0,p1, q0,q1;
4103 :
4104 126 : if (lx == 1)
4105 : {
4106 28 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
4107 21 : if (n >= 0) return cgetg(1,t_MAT);
4108 7 : return matid(2);
4109 : }
4110 98 : switch(typ(x))
4111 : {
4112 56 : case t_VEC: case t_COL: A = x; B = NULL; break;
4113 : case t_MAT:
4114 42 : switch(lgcols(x))
4115 : {
4116 0 : case 2: A = row(x,1); B = NULL; break;
4117 35 : case 3: A = row(x,2); B = row(x,1); break;
4118 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
4119 0 : return NULL; /*not reached*/
4120 : }
4121 35 : break;
4122 0 : default: pari_err_TYPE("pnqn",x);
4123 0 : return NULL; /*not reached*/
4124 : }
4125 91 : p1 = gel(A,1);
4126 91 : q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
4127 91 : if (n >= 0)
4128 : {
4129 56 : lx = minss(lx, n+2);
4130 56 : if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
4131 : }
4132 35 : else if (lx == 2)
4133 7 : return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
4134 : /* lx >= 3 */
4135 56 : p0 = gen_1;
4136 56 : q0 = gen_0; /* p[-1], q[-1] */
4137 56 : M = cgetg(lx, t_MAT);
4138 56 : gel(M,1) = mkcol2(p1,q1);
4139 217 : for (i=2; i<lx; i++)
4140 : {
4141 161 : GEN a = gel(A,i), p2,q2;
4142 161 : if (B) {
4143 84 : GEN b = gel(B,i);
4144 84 : p0 = gmul(b,p0);
4145 84 : q0 = gmul(b,q0);
4146 : }
4147 161 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
4148 161 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
4149 161 : gel(M,i) = mkcol2(p1,q1);
4150 : }
4151 56 : if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
4152 56 : return gerepilecopy(av, M);
4153 : }
4154 : GEN
4155 0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
4156 : /* x = [a0, ..., an] from gboundcf, n >= 0;
4157 : * return [[p0, ..., pn], [q0,...,qn]] */
4158 : GEN
4159 352730 : ZV_allpnqn(GEN x)
4160 : {
4161 352730 : long i, lx = lg(x);
4162 352730 : GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
4163 :
4164 352730 : gel(v,1) = P = cgetg(lx, t_VEC);
4165 352730 : gel(v,2) = Q = cgetg(lx, t_VEC);
4166 352730 : p0 = gen_1; q0 = gen_0;
4167 352730 : gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
4168 1204644 : for (i=2; i<lx; i++)
4169 : {
4170 851914 : GEN a = gel(x,i);
4171 851914 : gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
4172 851914 : gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
4173 : }
4174 352730 : return v;
4175 : }
4176 :
4177 : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
4178 : static GEN
4179 35 : mod_to_frac(GEN x, GEN N, GEN B)
4180 : {
4181 : GEN a, b, A;
4182 35 : if (B)
4183 : {
4184 21 : A = divii(shifti(N, -1), B);
4185 : /* denominator bound useless, don't use it */
4186 21 : if (cmpii(A, B) < 0) B = NULL;
4187 : }
4188 35 : if (!B)
4189 : {
4190 14 : A = sqrti(shifti(N, -1));
4191 14 : B = A;
4192 : }
4193 35 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
4194 21 : return equali1(b)? a: mkfrac(a,b);
4195 : }
4196 :
4197 : static GEN
4198 56 : mod_to_rfrac(GEN x, GEN N, long B)
4199 : {
4200 : GEN a, b;
4201 56 : long A, d = degpol(N);
4202 56 : if (B >= 0)
4203 : {
4204 21 : A = d-1 - B;
4205 : /* denominator bound useless, don't use it */
4206 21 : if (A < B) B = -1;
4207 : }
4208 56 : if (B < 0)
4209 : {
4210 42 : B = d >> 1;
4211 42 : A = odd(d)? B : B-1;
4212 : }
4213 56 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
4214 56 : if (! RgXQ_ratlift(x, N, A, B, &a,&b)) return NULL;
4215 56 : if (degpol(RgX_gcd(a,b)) > 0) return NULL;
4216 49 : return gdiv(a,b);
4217 : }
4218 :
4219 : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
4220 : * of the continued fraction of x with b <= k maximal */
4221 : static GEN
4222 0 : bestappr_frac(GEN x, GEN k)
4223 : {
4224 : pari_sp av;
4225 : GEN p0, p1, p, q0, q1, q, a, y;
4226 :
4227 0 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
4228 0 : av = avma; y = x;
4229 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
4230 0 : q1 = gen_0; q0 = gen_1;
4231 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
4232 : for(;;)
4233 : {
4234 0 : x = ginv(x); /* > 1 */
4235 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
4236 0 : if (cmpii(a,k) > 0)
4237 : { /* next partial quotient will overflow limits */
4238 : GEN n, d;
4239 0 : a = divii(subii(k, q1), q0);
4240 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4241 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4242 : /* compare |y-p0/q0|, |y-p1/q1| */
4243 0 : n = gel(y,1);
4244 0 : d = gel(y,2);
4245 0 : if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
4246 : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
4247 0 : { p1 = p0; q1 = q0; }
4248 0 : break;
4249 : }
4250 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4251 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4252 :
4253 0 : if (cmpii(q0,k) > 0) break;
4254 0 : x = gsub(x,a); /* 0 <= x < 1 */
4255 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
4256 :
4257 0 : }
4258 0 : return gerepileupto(av, gdiv(p1,q1));
4259 : }
4260 : /* bestappr(t_REAL != 0), to maximal accuracy */
4261 : static GEN
4262 21 : bestappr_real_max(GEN x)
4263 : {
4264 21 : pari_sp av = avma;
4265 : GEN p0, p1, p, q0, q1, q, a;
4266 : long e;
4267 21 : p1 = gen_1; a = p0 = floorr(x); q1 = gen_0; q0 = gen_1;
4268 21 : x = subri(x,a); /* 0 <= x < 1 */
4269 21 : e = bit_prec(x) - expo(x);
4270 : for(;;)
4271 : {
4272 : long d;
4273 42 : if (!signe(x) || expi(q0) > e) { p1 = p0; q1 = q0; break; }
4274 22 : x = invr(x); /* > 1 */
4275 22 : d = nbits2prec(expo(x) + 1);
4276 22 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4277 :
4278 21 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4279 21 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4280 21 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4281 21 : x = subri(x,a); /* 0 <= x < 1 */
4282 21 : }
4283 21 : return gerepileupto(av, gdiv(p1,q1));
4284 : }
4285 : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
4286 : * of the continued fraction of x with b <= k maximal */
4287 : static GEN
4288 272599 : bestappr_real(GEN x, GEN k)
4289 : {
4290 272599 : pari_sp av = avma;
4291 : GEN kr, p0, p1, p, q0, q1, q, a, y;
4292 :
4293 272599 : y = x;
4294 272599 : p1 = gen_1; a = p0 = floorr(x);
4295 272599 : q1 = gen_0; q0 = gen_1;
4296 272599 : x = subri(x,a); /* 0 <= x < 1 */
4297 272599 : if (!signe(x)) { cgiv(x); return a; }
4298 268274 : kr = itor(k, realprec(x));
4299 : for(;;)
4300 : {
4301 : long d;
4302 740230 : x = invr(x); /* > 1 */
4303 740230 : if (cmprr(x,kr) > 0)
4304 : { /* next partial quotient will overflow limits */
4305 263704 : a = divii(subii(k, q1), q0);
4306 263704 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4307 263704 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4308 : /* compare |y-p0/q0|, |y-p1/q1| */
4309 263704 : if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
4310 : mulir(q0, subri(mulir(q1,y), p1))) < 0)
4311 3627 : { p1 = p0; q1 = q0; }
4312 263704 : break;
4313 : }
4314 476526 : d = nbits2prec(expo(x) + 1);
4315 476526 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
4316 :
4317 476523 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
4318 476523 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
4319 476523 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
4320 :
4321 476523 : if (cmpii(q0,k) > 0) break;
4322 472445 : x = subri(x,a); /* 0 <= x < 1 */
4323 472445 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
4324 471956 : }
4325 268274 : return gerepileupto(av, gdiv(p1,q1));
4326 : }
4327 :
4328 : /* k t_INT or NULL */
4329 : static GEN
4330 401953 : bestappr_Q(GEN x, GEN k)
4331 : {
4332 401953 : long lx, tx = typ(x), i;
4333 : GEN a, y;
4334 :
4335 401953 : switch(tx)
4336 : {
4337 0 : case t_INT: return icopy(x);
4338 0 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
4339 : case t_REAL:
4340 306577 : if (!signe(x)) return gen_0;
4341 272620 : return k? bestappr_real(x, k): bestappr_real_max(x);
4342 :
4343 : case t_INTMOD: {
4344 21 : pari_sp av = avma;
4345 21 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
4346 14 : return gerepilecopy(av, a);
4347 : }
4348 : case t_PADIC: {
4349 14 : pari_sp av = avma;
4350 14 : long v = valp(x);
4351 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
4352 7 : if (v) a = gmul(a, powis(gel(x,2), v));
4353 7 : return gerepilecopy(av, a);
4354 : }
4355 :
4356 : case t_SER:
4357 0 : if (ser_isexactzero(x)) return gcopy(x);
4358 : /* fall through */
4359 : case t_COMPLEX: case t_POLMOD: case t_POL: case t_RFRAC:
4360 : case t_VEC: case t_COL: case t_MAT:
4361 95341 : y = cgetg_copy(x, &lx);
4362 95341 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
4363 488734 : for (; i<lx; i++)
4364 : {
4365 393393 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
4366 393393 : gel(y,i) = a;
4367 : }
4368 95341 : if (tx == t_POL) return normalizepol(y);
4369 95341 : if (tx == t_SER) return normalize(y);
4370 95341 : return y;
4371 : }
4372 0 : pari_err_TYPE("bestappr_Q",x);
4373 0 : return NULL; /* not reached */
4374 : }
4375 :
4376 : static GEN
4377 49 : bestappr_ser(GEN x, long B)
4378 : {
4379 49 : long v = valp(x), lx = lg(x);
4380 : GEN N, t;
4381 49 : x = normalizepol(ser2pol_i(x, lx));
4382 49 : N = monomial(gen_1, lx-2, varn(x));
4383 49 : t = mod_to_rfrac(x, N, B); if (!t) return NULL;
4384 42 : if (v)
4385 : {
4386 : GEN a, b;
4387 : long vx;
4388 14 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
4389 : /* t_RFRAC */
4390 14 : vx = varn(x);
4391 14 : a = gel(t,1);
4392 14 : b = gel(t,2);
4393 14 : v -= RgX_valrem(b, &b);
4394 14 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
4395 14 : if (v < 0) b = RgX_shift(b, -v);
4396 7 : else if (v > 0) {
4397 7 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
4398 7 : a = RgX_shift(a, v);
4399 : }
4400 14 : t = mkrfraccopy(a, b);
4401 : }
4402 42 : return t;
4403 : }
4404 : static GEN bestappr_RgX(GEN x, long B);
4405 : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
4406 : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
4407 : static GEN
4408 63 : bestappr_RgX(GEN x, long B)
4409 : {
4410 63 : long i, lx, tx = typ(x);
4411 : GEN y, t;
4412 63 : switch(tx)
4413 : {
4414 : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
4415 : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
4416 0 : return gcopy(x);
4417 :
4418 : case t_RFRAC: {
4419 14 : pari_sp av = avma;
4420 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
4421 7 : x = rfractoser(x, varn(gel(x,2)), 2*B+1);
4422 7 : t = bestappr_ser(x, B); if (!t) return NULL;
4423 7 : return gerepileupto(av, t);
4424 : }
4425 : case t_POLMOD: {
4426 7 : pari_sp av = avma;
4427 7 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
4428 7 : return gerepileupto(av, t);
4429 : }
4430 : case t_SER: {
4431 42 : pari_sp av = avma;
4432 42 : t = bestappr_ser(x, B); if (!t) return NULL;
4433 35 : return gerepileupto(av, t);
4434 : }
4435 :
4436 : case t_VEC: case t_COL: case t_MAT:
4437 0 : y = cgetg_copy(x, &lx);
4438 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
4439 0 : for (; i<lx; i++)
4440 : {
4441 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
4442 0 : gel(y,i) = t;
4443 : }
4444 0 : return y;
4445 : }
4446 0 : pari_err_TYPE("bestappr_RgX",x);
4447 0 : return NULL; /* not reached */
4448 : }
4449 :
4450 : /* allow k = NULL: maximal accuracy */
4451 : GEN
4452 8560 : bestappr(GEN x, GEN k)
4453 : {
4454 8560 : pari_sp av = avma;
4455 8560 : if (k) { /* replace by floor(k) */
4456 8525 : switch(typ(k))
4457 : {
4458 : case t_INT:
4459 168 : break;
4460 : case t_REAL: case t_FRAC:
4461 8357 : k = floor_safe(k); /* left on stack for efficiency */
4462 8357 : if (!signe(k)) k = gen_1;
4463 8357 : break;
4464 : default:
4465 0 : pari_err_TYPE("bestappr [bound type]", k);
4466 0 : break;
4467 : }
4468 : }
4469 8560 : x = bestappr_Q(x, k);
4470 8560 : if (!x) { avma = av; return cgetg(1,t_VEC); }
4471 8546 : return x;
4472 : }
4473 : GEN
4474 63 : bestapprPade(GEN x, long B)
4475 : {
4476 63 : pari_sp av = avma;
4477 63 : GEN t = bestappr_RgX(x, B);
4478 63 : if (!t) { avma = av; return cgetg(1,t_VEC); }
4479 56 : return t;
4480 : }
4481 :
4482 : /***********************************************************************/
4483 : /** **/
4484 : /** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
4485 : /** **/
4486 : /***********************************************************************/
4487 :
4488 : static GEN
4489 14 : get_quad(GEN f, GEN pol, long r)
4490 : {
4491 14 : GEN p1 = gcoeff(f,1,2), q1 = gcoeff(f,2,2);
4492 14 : return mkquad(pol, r? subii(p1,q1): p1, q1);
4493 : }
4494 :
4495 : /* replace f by f * [a,1; 1,0] */
4496 : static void
4497 14 : update_f(GEN f, GEN a)
4498 : {
4499 : GEN p1;
4500 14 : p1 = gcoeff(f,1,1);
4501 14 : gcoeff(f,1,1) = addii(mulii(a,p1), gcoeff(f,1,2));
4502 14 : gcoeff(f,1,2) = p1;
4503 :
4504 14 : p1 = gcoeff(f,2,1);
4505 14 : gcoeff(f,2,1) = addii(mulii(a,p1), gcoeff(f,2,2));
4506 14 : gcoeff(f,2,2) = p1;
4507 14 : }
4508 :
4509 : GEN
4510 7 : quadunit(GEN x)
4511 : {
4512 7 : pari_sp av = avma, av2;
4513 : GEN pol, y, a, u, v, sqd, f;
4514 : long r;
4515 :
4516 7 : check_quaddisc_real(x, &r, "quadunit");
4517 7 : pol = quadpoly(x);
4518 7 : sqd = sqrti(x); av2 = avma;
4519 7 : a = shifti(addsi(r,sqd),-1);
4520 7 : f = mkmat2(mkcol2(a, gen_1), mkcol2(gen_1, gen_0)); /* [a,0; 1,0] */
4521 7 : u = stoi(r); v = gen_2;
4522 : for(;;)
4523 : {
4524 : GEN u1, v1;
4525 14 : u1 = subii(mulii(a,v),u);
4526 14 : v1 = divii(subii(x,sqri(u1)),v);
4527 14 : if ( equalii(v,v1) ) {
4528 7 : y = get_quad(f,pol,r);
4529 7 : update_f(f,a);
4530 7 : y = gdiv(get_quad(f,pol,r), gconj(y));
4531 7 : break;
4532 : }
4533 7 : a = divii(addii(sqd,u1), v1);
4534 7 : if ( equalii(u,u1) ) {
4535 0 : y = get_quad(f,pol,r);
4536 0 : y = gdiv(y, gconj(y));
4537 0 : break;
4538 : }
4539 7 : update_f(f,a);
4540 7 : u = u1; v = v1;
4541 7 : if (gc_needed(av2,2))
4542 : {
4543 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit");
4544 0 : gerepileall(av2,4, &a,&f,&u,&v);
4545 : }
4546 7 : }
4547 7 : if (signe(gel(y,3)) < 0) y = gneg(y);
4548 7 : return gerepileupto(av, y);
4549 : }
4550 :
4551 : GEN
4552 21 : quadregulator(GEN x, long prec)
4553 : {
4554 21 : pari_sp av = avma, av2;
4555 : GEN R, rsqd, u, v, sqd;
4556 : long r, Rexpo;
4557 :
4558 21 : check_quaddisc_real(x, &r, "quadregulator");
4559 21 : sqd = sqrti(x);
4560 21 : rsqd = gsqrt(x,prec);
4561 21 : Rexpo = 0; R = real2n(1, prec); /* = 2 */
4562 21 : av2 = avma;
4563 21 : u = stoi(r); v = gen_2;
4564 : for(;;)
4565 : {
4566 70 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
4567 70 : GEN v1 = divii(subii(x,sqri(u1)),v);
4568 70 : if (equalii(v,v1))
4569 : {
4570 7 : R = sqrr(R); shiftr_inplace(R, -1);
4571 7 : R = mulrr(R, divri(addir(u1,rsqd),v));
4572 7 : break;
4573 : }
4574 63 : if (equalii(u,u1))
4575 : {
4576 14 : R = sqrr(R); shiftr_inplace(R, -1);
4577 14 : break;
4578 : }
4579 49 : R = mulrr(R, divri(addir(u1,rsqd),v));
4580 49 : Rexpo += expo(R); setexpo(R,0);
4581 49 : u = u1; v = v1;
4582 49 : if (Rexpo & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
4583 49 : if (gc_needed(av2,2))
4584 : {
4585 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
4586 0 : gerepileall(av2,3, &R,&u,&v);
4587 : }
4588 49 : }
4589 21 : R = logr_abs(divri(R,v));
4590 21 : if (Rexpo)
4591 : {
4592 21 : GEN t = mulsr(Rexpo, mplog2(prec));
4593 21 : shiftr_inplace(t, 1);
4594 21 : R = addrr(R,t);
4595 : }
4596 21 : return gerepileuptoleaf(av, R);
4597 : }
4598 :
4599 : /*************************************************************************/
4600 : /** **/
4601 : /** CLASS NUMBER **/
4602 : /** **/
4603 : /*************************************************************************/
4604 :
4605 : int
4606 10628904 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
4607 :
4608 15047085 : static GEN qfi_pow(void *E, GEN f, GEN n)
4609 15047085 : { return E? nupow(f,n,(GEN)E): powgi(f,n); }
4610 18997455 : static GEN qfi_comp(void *E, GEN f, GEN g)
4611 18997455 : { return E? nucomp(f,g,(GEN)E): qficomp(f,g); }
4612 : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
4613 : gidentical,qfb_equal1,NULL};
4614 :
4615 : GEN
4616 2418537 : qfi_order(GEN q, GEN o)
4617 2418537 : { return gen_order(q, o, NULL, &qfi_group); }
4618 :
4619 : GEN
4620 0 : qfi_log(GEN a, GEN g, GEN o)
4621 0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
4622 :
4623 : GEN
4624 521297 : qfi_Shanks(GEN a, GEN g, long n)
4625 : {
4626 521297 : pari_sp av = avma;
4627 : GEN T, X;
4628 : long rt_n, c;
4629 :
4630 521297 : a = redimag(a);
4631 521297 : g = redimag(g);
4632 :
4633 521297 : rt_n = sqrt((double)n);
4634 521297 : c = n / rt_n;
4635 521297 : c = (c * rt_n < n + 1) ? c + 1 : c;
4636 :
4637 521297 : T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
4638 521297 : X = gen_Shanks(T, a, c, NULL, &qfi_group);
4639 :
4640 521297 : if (!X) { avma = av; return X; }
4641 279692 : return gerepileuptoint(av, X);
4642 : }
4643 :
4644 : GEN
4645 140 : qfbclassno0(GEN x,long flag)
4646 : {
4647 140 : switch(flag)
4648 : {
4649 126 : case 0: return map_proto_G(classno,x);
4650 14 : case 1: return map_proto_G(classno2,x);
4651 0 : default: pari_err_FLAG("qfbclassno");
4652 : }
4653 0 : return NULL; /* not reached */
4654 : }
4655 :
4656 : /* f^h = 1, return order(f). Set *pfao to its factorization */
4657 : static GEN
4658 2270664 : find_order(void *E, GEN f, GEN h, GEN *pfao)
4659 : {
4660 2270664 : GEN v = gen_factored_order(f, h, E, &qfi_group);
4661 2270664 : *pfao = gel(v,2); return gel(v,1);
4662 : }
4663 :
4664 : static int
4665 5550 : ok_q(GEN q, GEN h, GEN d2, long r2)
4666 : {
4667 5550 : if (d2)
4668 : {
4669 7 : if (r2 <= 2 && !mpodd(q)) return 0;
4670 7 : return is_pm1(coprime_part(q,d2));
4671 : }
4672 : else
4673 : {
4674 5543 : if (r2 <= 1 && !mpodd(q)) return 0;
4675 5543 : return is_pm1(coprime_part(q,h));
4676 : }
4677 : }
4678 :
4679 : /* a,b given by their factorizations. Return factorization of lcm(a,b).
4680 : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
4681 : static GEN
4682 295900 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
4683 : {
4684 295900 : GEN P = ZV_union_shallow(gel(Fa,1), gel(Fb,1));
4685 295900 : GEN A = gen_1, B = gen_1;
4686 295900 : long i, l = lg(P);
4687 295900 : GEN E = cgetg(l, t_COL);
4688 874865 : for (i=1; i<l; i++)
4689 : {
4690 578965 : GEN p = gel(P,i);
4691 578965 : long va = Z_pval(a,p);
4692 578965 : long vb = Z_pval(b,p);
4693 578965 : if (va < vb)
4694 : {
4695 297621 : B = mulii(B,powiu(p,vb));
4696 297621 : gel(E,i) = utoi(vb);
4697 : }
4698 : else
4699 : {
4700 281344 : A = mulii(A,powiu(p,va));
4701 281344 : gel(E,i) = utoi(va);
4702 : }
4703 : }
4704 295900 : *pA = A;
4705 295900 : *pB = B; return mkmat2(P,E);
4706 : }
4707 :
4708 : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
4709 : static void
4710 295900 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
4711 : {
4712 295900 : GEN A,B, g1 = *pg1, d1 = *pd1;
4713 295900 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
4714 295900 : *pg1 = gmul(powgi(g1, diviiexact(d1,A)), powgi(f, diviiexact(o,B)));
4715 295900 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
4716 295900 : }
4717 :
4718 : /* Write x = Df^2, where D = fundamental discriminant,
4719 : * P^E = factorisation of conductor f, with E[i] >= 0 */
4720 : static void
4721 1720764 : corediscfact(GEN x, long xmod4, GEN *ptD, GEN *ptP, GEN *ptE)
4722 : {
4723 1720764 : long s = signe(x), l, i;
4724 1720764 : GEN fa = absZ_factor(x);
4725 1720764 : GEN d, P = gel(fa,1), E = gtovecsmall(gel(fa,2));
4726 :
4727 1720764 : l = lg(P); d = gen_1;
4728 4508027 : for (i=1; i<l; i++)
4729 : {
4730 2787263 : if (E[i] & 1) d = mulii(d, gel(P,i));
4731 2787263 : E[i] >>= 1;
4732 : }
4733 1720764 : if (!xmod4 && mod4(d) != ((s < 0)? 3: 1)) { d = shifti(d,2); E[1]--; }
4734 1720764 : *ptD = (s < 0)? negi(d): d;
4735 1720764 : *ptP = P;
4736 1720764 : *ptE = E;
4737 1720764 : }
4738 :
4739 : static GEN
4740 1678710 : conductor_part(GEN x, long xmod4, GEN *ptD, GEN *ptreg)
4741 : {
4742 1678710 : long l, i, s = signe(x);
4743 : GEN E, H, D, P, reg;
4744 :
4745 1678710 : corediscfact(x, xmod4, &D, &P, &E);
4746 1678710 : H = gen_1; l = lg(P);
4747 : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
4748 4360301 : for (i=1; i<l; i++)
4749 : {
4750 2681591 : long e = E[i];
4751 2681591 : if (e)
4752 : {
4753 7 : GEN p = gel(P,i);
4754 7 : H = mulii(H, subis(p, kronecker(D,p)));
4755 7 : if (e >= 2) H = mulii(H, powiu(p,e-1));
4756 : }
4757 : }
4758 :
4759 : /* divide by [ O_K^* : O^* ] */
4760 1678710 : if (s < 0)
4761 : {
4762 1678696 : reg = NULL;
4763 1678696 : switch(itou_or_0(D))
4764 : {
4765 0 : case 4: H = shifti(H,-1); break;
4766 0 : case 3: H = divis(H,3); break;
4767 : }
4768 : } else {
4769 14 : reg = quadregulator(D,DEFAULTPREC);
4770 14 : if (!equalii(x,D))
4771 0 : H = divii(H, roundr(divrr(quadregulator(x,DEFAULTPREC), reg)));
4772 : }
4773 1678710 : if (ptreg) *ptreg = reg;
4774 1678710 : *ptD = D; return H;
4775 : }
4776 :
4777 : static long
4778 1678689 : two_rank(GEN x)
4779 : {
4780 1678689 : GEN p = gel(absZ_factor(x),1);
4781 1678689 : long l = lg(p)-1;
4782 : #if 0 /* positive disc not needed */
4783 : if (signe(x) > 0)
4784 : {
4785 : long i;
4786 : for (i=1; i<=l; i++)
4787 : if (mod4(gel(p,i)) == 3) { l--; break; }
4788 : }
4789 : #endif
4790 1678689 : return l-1;
4791 : }
4792 :
4793 : static GEN
4794 31893332 : sqr_primeform(GEN x, ulong p) { return redimag(qfisqr(primeform_u(x, p))); }
4795 : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
4796 : static GEN
4797 1678689 : get_forms(GEN D, GEN *pL)
4798 : {
4799 1678689 : const long MAXFORM = 20;
4800 1678689 : GEN L, sqrtD = gsqrt(absi(D),DEFAULTPREC), forms = vectrunc_init(MAXFORM+1);
4801 1678689 : long s, nforms = 0;
4802 : ulong p;
4803 : forprime_t S;
4804 1678689 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
4805 1678689 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
4806 1678689 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
4807 1678689 : if (s < 10) s = 200;
4808 1554960 : else if (s < 20) s = 1000;
4809 1477 : else if (s < 5000) s = 5000;
4810 1678689 : u_forprime_init(&S, 2, s);
4811 284027476 : while ( (p = u_forprime_next(&S)) )
4812 : {
4813 280670098 : long d, k = kroiu(D,p);
4814 : pari_sp av2;
4815 280670098 : if (!k) continue;
4816 278845499 : if (k > 0)
4817 : {
4818 139862100 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
4819 139862100 : d = p-1;
4820 : }
4821 : else
4822 138983399 : d = p+1;
4823 278845499 : av2 = avma; affrr(divru(mulur(p,L),d), L); avma = av2;
4824 : }
4825 1678689 : *pL = L; return forms;
4826 : }
4827 :
4828 : /* h ~ #G, return o = order of f, set fao = its factorization */
4829 : static GEN
4830 1678738 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
4831 : {
4832 1678738 : long s = minss(itos(sqrti(h)), 10000);
4833 1678738 : GEN T = gen_Shanks_init(f, s, E, &qfi_group);
4834 1678738 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
4835 1678738 : return find_order(E, f, addiu(v,1), pfao);
4836 : }
4837 :
4838 : /* if g = 1 in G/<f> ? */
4839 : static int
4840 518 : equal1(void *E, GEN T, ulong N, GEN g)
4841 518 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
4842 :
4843 : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
4844 : * FIXME: should be gen_order, but equal1 has the wrong prototype */
4845 : static GEN
4846 112 : relative_order(void *E, GEN a, GEN o, ulong N, GEN T)
4847 : {
4848 112 : pari_sp av = avma;
4849 : long i, l;
4850 : GEN m;
4851 :
4852 112 : m = get_arith_ZZM(o);
4853 112 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
4854 112 : o = gel(m,1);
4855 112 : m = gel(m,2); l = lgcols(m);
4856 322 : for (i = l-1; i; i--)
4857 : {
4858 210 : GEN t, y, p = gcoeff(m,i,1);
4859 210 : long j, e = itos(gcoeff(m,i,2));
4860 210 : if (l == 2) {
4861 35 : t = gen_1;
4862 35 : y = a;
4863 : } else {
4864 175 : t = diviiexact(o, powiu(p,e));
4865 175 : y = powgi(a, t);
4866 : }
4867 210 : if (equal1(E, T,N,y)) o = t;
4868 : else {
4869 126 : for (j = 1; j < e; j++)
4870 : {
4871 28 : y = powgi(y, p);
4872 28 : if (equal1(E, T,N,y)) break;
4873 : }
4874 119 : if (j < e) {
4875 21 : if (j > 1) p = powiu(p, j);
4876 21 : o = mulii(t, p);
4877 : }
4878 : }
4879 : }
4880 112 : return gerepilecopy(av, o);
4881 : }
4882 :
4883 : /* h(x) for x<0 using Baby Step/Giant Step.
4884 : * Assumes G is not too far from being cyclic.
4885 : *
4886 : * Compute G^2 instead of G so as to kill most of the non-cyclicity */
4887 : GEN
4888 1680638 : classno(GEN x)
4889 : {
4890 1680638 : pari_sp av = avma;
4891 : long r2, k, s, i, l;
4892 : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
4893 : void *E;
4894 :
4895 1680638 : if (signe(x) >= 0) return classno2(x);
4896 :
4897 1680631 : check_quaddisc(x, &s, &k, "classno");
4898 1680631 : if (abscmpiu(x,12) <= 0) return gen_1;
4899 :
4900 1678689 : Hf = conductor_part(x, k, &D, NULL);
4901 1678689 : if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
4902 1678689 : forms = get_forms(D, &L);
4903 1678689 : r2 = two_rank(D);
4904 1678689 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
4905 :
4906 1678689 : l = lg(forms);
4907 1678689 : order_bound = const_vec(l-1, NULL);
4908 1678689 : E = expi(D) > 60? (void*)sqrtnint(shifti(absi(D),-2),4): NULL;
4909 1678689 : g1 = gel(forms,1);
4910 1678689 : gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
4911 1678689 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
4912 1678689 : d2 = NULL; /* not computed yet */
4913 1678689 : if (is_pm1(q)) goto END;
4914 417553 : for (i=2; i < l; i++)
4915 : {
4916 411940 : GEN o, fao, a, F, fd, f = gel(forms,i);
4917 411940 : fd = powgi(f, d1); if (is_pm1(gel(fd,1))) continue;
4918 295900 : F = powgi(fd, q);
4919 295900 : a = gel(F,1);
4920 295900 : o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
4921 : /* f^(d1 q) = 1 */
4922 295900 : fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
4923 295900 : o = find_order(E, f, fao, &fao);
4924 295900 : gel(order_bound,i) = o;
4925 : /* o = order of f, fao = factor(o) */
4926 295900 : update_g1(&g1,&d1,&fad1, f,o,fao);
4927 295900 : q = diviiround(hin, d1);
4928 295900 : if (is_pm1(q)) goto END;
4929 : }
4930 : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
4931 5613 : if (expi(q) > 3)
4932 : { /* q large: compute d2, 2nd elt divisor */
4933 70 : ulong N, n = 2*itou(sqrti(d1));
4934 70 : GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
4935 70 : d2 = gen_1;
4936 70 : N = itou( gceil(gdivgs(d1,n)) ); /* order(g1) <= n*N */
4937 287 : for (i = 1; i < l; i++)
4938 : {
4939 280 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
4940 280 : if (!B) B = find_order(E, f, fad1, /*junk*/&d);
4941 280 : f = powgi(f,d2);
4942 280 : if (equal1(E, T,N,f)) continue;
4943 112 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
4944 : /* f^B = 1 */
4945 112 : d = relative_order(E, f, B, N,T);
4946 112 : d2= mulii(d,d2);
4947 112 : D = mulii(d1,d2);
4948 112 : q = diviiround(hin,D);
4949 112 : if (is_pm1(q)) { d1 = D; goto END; }
4950 : }
4951 : /* very probably, d2 is the 2nd elementary divisor */
4952 7 : d1 = D; /* product of first two elt divisors */
4953 : }
4954 : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
4955 : * 2-rank */
4956 5550 : if (!ok_q(q,d1,d2,r2))
4957 : {
4958 0 : GEN q0 = q;
4959 : long d;
4960 0 : if (cmpii(mulii(q,d1), hin) < 0)
4961 : { /* try q = q0+1,-1,+2,-2 */
4962 0 : d = 1;
4963 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
4964 : }
4965 : else
4966 : { /* q0-1,+1,-2,+2 */
4967 0 : d = -1;
4968 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
4969 : }
4970 : }
4971 5550 : d1 = mulii(d1,q);
4972 :
4973 : END:
4974 1678689 : return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
4975 : }
4976 :
4977 : /* use Euler products */
4978 : GEN
4979 21 : classno2(GEN x)
4980 : {
4981 21 : pari_sp av = avma;
4982 21 : const long prec = DEFAULTPREC;
4983 : long n, i, r, s;
4984 : GEN p1, p2, S, p4, p5, p7, Hf, Pi, reg, logd, d, dr, D, half;
4985 :
4986 21 : check_quaddisc(x, &s, &r, "classno2");
4987 21 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
4988 :
4989 21 : Hf = conductor_part(x, r, &D, ®);
4990 21 : if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
4991 :
4992 21 : Pi = mppi(prec);
4993 21 : d = absi(D); dr = itor(d, prec);
4994 21 : logd = logr_abs(dr);
4995 21 : p1 = sqrtr(divrr(mulir(d,logd), gmul2n(Pi,1)));
4996 21 : if (s > 0)
4997 : {
4998 14 : GEN invlogd = invr(logd);
4999 14 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
5000 14 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
5001 : }
5002 21 : n = itos_or_0( mptrunc(p1) );
5003 21 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
5004 :
5005 21 : p4 = divri(Pi,d);
5006 21 : p7 = invr(sqrtr_abs(Pi));
5007 21 : half = real2n(-1, prec);
5008 21 : if (s > 0)
5009 : { /* i = 1, shortcut */
5010 14 : p1 = sqrtr_abs(dr);
5011 14 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5012 14 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
5013 546 : for (i=2; i<=n; i++)
5014 : {
5015 532 : long k = kroiu(D,i); if (!k) continue;
5016 434 : p2 = mulir(sqru(i), p4);
5017 434 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5018 434 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
5019 434 : S = (k>0)? addrr(S,p5): subrr(S,p5);
5020 : }
5021 14 : S = shiftr(divrr(S,reg),-1);
5022 : }
5023 : else
5024 : { /* i = 1, shortcut */
5025 7 : p1 = gdiv(sqrtr_abs(dr), Pi);
5026 7 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
5027 7 : S = addrr(p5, divrr(p1, mpexp(p4)));
5028 952 : for (i=2; i<=n; i++)
5029 : {
5030 945 : long k = kroiu(D,i); if (!k) continue;
5031 945 : p2 = mulir(sqru(i), p4);
5032 945 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
5033 945 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
5034 945 : S = (k>0)? addrr(S,p5): subrr(S,p5);
5035 : }
5036 : }
5037 21 : return gerepileuptoint(av, mulii(Hf, roundr(S)));
5038 : }
5039 :
5040 : static GEN
5041 5304 : hclassno2(GEN x)
5042 : {
5043 : long i, l, s, xmod4;
5044 : GEN Q, H, D, P, E;
5045 :
5046 5304 : x = negi(x);
5047 5304 : check_quaddisc(x, &s, &xmod4, "hclassno");
5048 5304 : corediscfact(x, xmod4, &D, &P, &E);
5049 :
5050 5304 : Q = quadclassunit0(D, 0, NULL, 0);
5051 5304 : H = gel(Q,1); l = lg(P);
5052 :
5053 : /* H \prod_{p^e||f} (1 + (p^e-1)/(p-1))[ p - (D/p) ] */
5054 24043 : for (i=1; i<l; i++)
5055 : {
5056 18739 : long e = E[i];
5057 18739 : if (e)
5058 : {
5059 4198 : GEN p = gel(P,i), t = subis(p, kronecker(D,p));
5060 4198 : if (e > 1) t = mulii(t, diviiexact(subis(powiu(p,e), 1), subis(p,1)));
5061 4198 : H = mulii(H, addsi(1, t));
5062 : }
5063 : }
5064 5304 : switch( itou_or_0(D) )
5065 : {
5066 0 : case 3: H = gdivgs(H, 3); break;
5067 0 : case 4: H = gdivgs(H, 2); break;
5068 : }
5069 5304 : return H;
5070 : }
5071 :
5072 : GEN
5073 76111 : hclassno(GEN x)
5074 : {
5075 : ulong a, b, b2, d, h;
5076 : long s;
5077 : int f;
5078 :
5079 76111 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
5080 76111 : s = signe(x);
5081 76111 : if (s < 0) return gen_0;
5082 76111 : if (!s) return gdivgs(gen_1, -12);
5083 :
5084 76111 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
5085 :
5086 76111 : d = itou_or_0(x);
5087 76111 : if (!d || d > 500000) return hclassno2(x);
5088 :
5089 70807 : h = 0; b = d&1; b2 = (1+d)>>2; f=0;
5090 70807 : if (!b)
5091 : {
5092 2074994 : for (a=1; a*a<b2; a++)
5093 2020913 : if (b2%a == 0) h++;
5094 54081 : f = (a*a==b2); b=2; b2=(4+d)>>2;
5095 : }
5096 2092598 : while (b2*3 < d)
5097 : {
5098 1950984 : if (b2%b == 0) h++;
5099 131085034 : for (a=b+1; a*a < b2; a++)
5100 129134050 : if (b2%a == 0) h += 2;
5101 1950984 : if (a*a == b2) h++;
5102 1950984 : b += 2; b2 = (b*b+d)>>2;
5103 : }
5104 70807 : if (b2*3 == d) retmkfrac(utoipos(3*h+1), utoipos(3));
5105 70702 : if (f) retmkfrac(utoipos(2*h+1), gen_2);
5106 66131 : return utoipos(h);
5107 : }
5108 :
5109 : /******************************************************************/
5110 : /* */
5111 : /* RAMANUJAN's TAU FUNCTION */
5112 : /* */
5113 : /******************************************************************/
5114 :
5115 : /* h(D) / (w(D)/2), D fundamental */
5116 : static GEN
5117 36750 : myh(GEN D)
5118 : {
5119 : GEN Q;
5120 36750 : if (equalis(D, -3)) return mkfrac(gen_1, utoipos(3));
5121 32690 : if (equalis(D, -4)) return ghalf;
5122 28161 : Q = quadclassunit0(D, 0, NULL, 0); return gel(Q,1);
5123 : }
5124 :
5125 : /* 1 + q + ... + q^v */
5126 : static GEN
5127 6083 : geomsum(GEN q, long v)
5128 : {
5129 : GEN u;
5130 6083 : if (v == 0) return gen_1;
5131 273 : u = addui(1, q);
5132 273 : for (; v > 1; v--) u = addui(1, mulii(q, u));
5133 273 : return u;
5134 : }
5135 :
5136 : /* 4|N, not fundamental at 2; Hurwitz class number in level 2,
5137 : * equal to H(N)+2H(N/4), H=qfbhclassno */
5138 : static GEN
5139 36750 : Hspec(GEN N)
5140 : {
5141 : GEN D0, P, E, H, s;
5142 : long j, lP;
5143 :
5144 36750 : corediscfact(negi(N), 0, &D0, &P, &E);
5145 36750 : H = myh(D0);
5146 36750 : lP = lg(P); /* E[1] > 0 since N not fundamental at */
5147 36750 : s = addui(3, muliu(subiu(int2n(E[1]+1), 3), 2 - kroiu(D0,2)));
5148 86933 : for (j = 2; j < lP; j++)
5149 : {
5150 50183 : long e = E[j];
5151 50183 : if (e)
5152 : {
5153 6083 : ulong p = itou(gel(P,j));
5154 6083 : GEN q = geomsum(utoipos(p), e-1);
5155 6083 : s = mulii(s, addui(1, muliu(q, p - kroiu(D0,p))));
5156 : }
5157 : }
5158 36750 : return gmul(H,s);
5159 : }
5160 :
5161 : /* Ramanujan tau function for p prime */
5162 : static GEN
5163 14903 : tauprime(GEN p)
5164 : {
5165 14903 : pari_sp av = avma, av2;
5166 : GEN s, p_27, p_9, T;
5167 : ulong lim, t, tin;
5168 :
5169 14903 : if (absequaliu(p, 2)) return utoineg(24);
5170 : /* p > 2 */
5171 11396 : p_27 = mulsi(7, sqri(p));
5172 11396 : p_9 = mulsi(9, p);
5173 11396 : av2 = avma;
5174 11396 : lim = itou(sqrtint(p));
5175 11396 : tin = mod4(p) == 3? 1: 0;
5176 11396 : s = gen_0;
5177 87178 : for (t = 1; t <= lim; ++t)
5178 : {
5179 75782 : GEN p2, t2 = sqru(t), t3 = shifti(subii(p, t2), 2); /* 4(p-t^2) */
5180 : /* t mod 2 != tin <=> t3 not fundamental at 2 */
5181 75782 : p2 = ((t&1) == tin) ? hclassno(t3) : Hspec(t3);
5182 75782 : s = gadd(s, gmul(mulii(powiu(t2,3), addii(p_27, mulii(t2, subii(mulsi(4, t2), p_9)))), p2));
5183 75782 : if (!(t & 255)) s = gerepileupto(av2, s);
5184 : }
5185 : /* 28p^3 - 28p^2 - 90p - 35 */
5186 11396 : T = addsi(-35, mulii(p, addsi(-90, mulii(p, addsi(-28, mulsi(28, p))))));
5187 11396 : return gerepileupto(av, subii(mulii(powiu(p,3),T), addsi(1, gmulsg(128, s))));
5188 : }
5189 :
5190 : /* Ramanujan tau function, return 0 for <= 0 */
5191 : GEN
5192 7035 : ramanujantau(GEN n)
5193 : {
5194 7035 : pari_sp ltop = avma;
5195 : GEN T, F, P, E;
5196 : long j, lP;
5197 :
5198 7035 : if (!(F = check_arith_all(n,"ramanujantau")))
5199 : {
5200 7014 : if (signe(n) <= 0) return gen_0;
5201 7007 : F = Z_factor(n);
5202 : }
5203 : else
5204 : {
5205 21 : P = gel(F,1);
5206 21 : if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
5207 : }
5208 :
5209 7014 : P = gel(F,1);
5210 7014 : E = gel(F,2); lP = lg(P);
5211 7014 : T = gen_1;
5212 21917 : for (j = 1; j < lP; j++)
5213 : {
5214 14903 : GEN p = gel(P,j), tp = tauprime(p), t1 = tp, t0 = gen_1;
5215 14903 : long k, e = itou(gel(E,j));
5216 20160 : for (k = 1; k < e; k++)
5217 : {
5218 5257 : GEN t2 = subii(mulii(tp, t1), mulii(powiu(p, 11), t0));
5219 5257 : t0 = t1; t1 = t2;
5220 : }
5221 14903 : T = mulii(T, t1);
5222 : }
5223 7014 : return gerepileuptoint(ltop, T);
5224 : }
|
