Branch data 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 [ + - ]: 34 : remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
30 : : static ulong
31 : 18206 : u_remove2(ulong q) { return q >> vals(q); }
32 : : GEN
33 : 34 : odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
34 : : static GEN
35 : 18206 : 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 : 27 : is_gener_expo(GEN p, GEN L0)
40 : : {
41 : 27 : GEN L, q = shifti(p,-1);
42 : : long i, l;
43 [ + + ]: 27 : if (L0) {
44 : 8 : l = lg(L0);
45 : 8 : L = cgetg(l, t_VEC);
46 : : } else {
47 : 19 : L0 = L = odd_prime_divisors(q);
48 : 19 : l = lg(L);
49 : : }
50 [ + + ]: 87 : for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
51 : 27 : return L;
52 : : }
53 : : static GEN
54 : 18022 : u_is_gener_expo(ulong p, GEN L0)
55 : : {
56 : 18022 : const ulong q = p >> 1;
57 : : long i, l;
58 : : GEN L;
59 [ + + ]: 18022 : if (L0) {
60 : 352 : l = lg(L0);
61 : 352 : L = cgetg(l, t_VECSMALL);
62 : : } else {
63 : 17670 : L0 = L = u_odd_prime_divisors(q);
64 : 17670 : l = lg(L);
65 : : }
66 [ + + ]: 51271 : for (i=1; i<l; i++) L[i] = q / (ulong)L0[i];
67 : 18022 : return L;
68 : : }
69 : :
70 : : int
71 : 39499 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
72 : : {
73 : : long i;
74 [ + + ]: 39499 : if (krouu(x, p) >= 0) return 0;
75 [ + + ]: 56914 : for (i=lg(L)-1; i; i--)
76 : : {
77 : 37790 : ulong t = Fl_powu(x, (ulong)L[i], p);
78 [ + + ][ + + ]: 37790 : if (t == p_1 || t == 1) return 0;
79 : : }
80 : 39499 : return 1;
81 : : }
82 : : /* assume p prime */
83 : : ulong
84 : 43325 : pgener_Fl_local(ulong p, GEN L0)
85 : : {
86 : 43325 : const pari_sp av = avma;
87 : 43325 : const ulong p_1 = p-1;
88 : : long x;
89 : : GEN L;
90 [ + + ]: 43325 : if (p <= 19) switch(p)
[ - + + ]
91 : : { /* quick trivial cases */
92 : 0 : case 2: return 1;
93 : : case 7:
94 : 2988 : case 17: return 3;
95 : 22339 : default: return 2;
96 : : }
97 : 17998 : L = u_is_gener_expo(p,L0);
98 [ + + ]: 37815 : for (x=2;;x++) { if (is_gener_Fl(x,p,p_1,L)) break; }
99 : 43325 : avma = av; return x;
100 : : }
101 : : ulong
102 : 35943 : 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 : 119 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
108 : : {
109 [ + + ]: 119 : long i, t = lgefint(x)==3? krosi(x[2], p): kronecker(x, p);
110 [ + + ]: 119 : if (t >= 0) return 0;
111 [ + + ]: 261 : for (i = lg(L)-1; i; i--)
112 : : {
113 : 190 : GEN t = Fp_pow(x, gel(L,i), p);
114 [ + + ][ + + ]: 190 : if (equalii(t, p_1) || equali1(t)) return 0;
115 : : }
116 : 119 : return 1;
117 : : }
118 : :
119 : : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
120 : : GEN
121 : 9302 : pgener_Fp_local(GEN p, GEN L0)
122 : : {
123 : 9302 : pari_sp av0 = avma;
124 : : GEN x, p_1, L;
125 [ + + ]: 9302 : if (lgefint(p) == 3)
126 : : {
127 : : ulong z;
128 [ + + ]: 9279 : if (p[2] == 2) return gen_1;
129 [ + + ]: 6984 : if (L0) L0 = ZV_to_nv(L0);
130 : 6984 : z = pgener_Fl_local((ulong)p[2], L0);
131 : 6984 : avma = av0; return utoipos(z);
132 : : }
133 : 23 : p_1 = subis(p,1); L = is_gener_expo(p, L0);
134 : 23 : x = utoipos(2);
135 [ + + ]: 26 : for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
136 : 9302 : avma = av0; return utoipos((ulong)x[2]);
137 : : }
138 : :
139 : : GEN
140 : 8867 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
141 : :
142 : : ulong
143 : 52 : pgener_Zl(ulong p)
144 : : {
145 [ - + ]: 52 : 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 [ - + ]: 52 : if (p == 40487) return 10;
148 : : #ifndef LONG_IS_64BIT
149 : 4 : return pgener_Fl(p);
150 : : #else
151 [ + + ]: 48 : if (p < (1UL<<32)) return pgener_Fl(p);
152 : : else
153 : : {
154 : 24 : const pari_sp av = avma;
155 : 24 : const ulong p_1 = p-1;
156 : : long x ;
157 : 24 : GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
158 : 24 : for (x=2;;x++)
159 [ + + ][ + - ]: 96 : if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2))) break;
160 : 120 : avma = av; return x;
161 : : }
162 : : #endif
163 : : }
164 : :
165 : : /* p prime. Return a primitive root modulo p^e, e > 1 */
166 : : GEN
167 : 56 : pgener_Zp(GEN p)
168 : : {
169 [ + + ]: 56 : if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
170 : : else
171 : : {
172 : 4 : const pari_sp av = avma;
173 : 4 : GEN p_1 = subis(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
174 : 4 : GEN x = utoipos(2);
175 : 12 : for (;; x[2]++)
176 [ + + ][ + - ]: 16 : if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
177 : 68 : avma = av; return utoipos((ulong)x[2]);
178 : : }
179 : : }
180 : :
181 : : static GEN
182 : 357 : gener_Zp(GEN q, GEN F)
183 : : {
184 : 357 : GEN p = NULL;
185 : 357 : long e = 0;
186 [ + + ]: 357 : 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 (equaliu(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 : 343 : e = Z_isanypower(q, &p);
201 [ + + ]: 357 : return e > 1? pgener_Zp(p): pgener_Fp(q);
202 : : }
203 : :
204 : : GEN
205 : 420 : znprimroot(GEN N)
206 : : {
207 : 420 : pari_sp av = avma;
208 : : GEN x, n, F;
209 : :
210 [ + + ]: 420 : 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 [ - + ]: 413 : if (signe(N) < 0) N = absi(N);
216 [ + + ]: 413 : if (cmpiu(N, 4) <= 0) { avma = av; return mkintmodu(N[2]-1,N[2]); }
217 [ + + + ]: 371 : 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 : 336 : x = gener_Zp(N,F);
228 : 336 : break;
229 : : }
230 : 399 : 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 : 305 : rootsof1u_Fp(ulong n, GEN p)
246 : : {
247 : 305 : pari_sp av = avma;
248 : 305 : GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
249 : 305 : z = pgener_Fp_local(p, Flv_to_ZV(L));
250 : 305 : z = Fp_pow(z, diviuexact(subis(p,1), n), p); /* prim. n-th root of 1 */
251 : 305 : return gerepileuptoint(av, z);
252 : : }
253 : :
254 : : ulong
255 : 231 : rootsof1_Fl(ulong n, ulong p)
256 : : {
257 : 231 : pari_sp av = avma;
258 : 231 : GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
259 : 231 : ulong z = pgener_Fl_local(p, L);
260 : 231 : z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
261 : 231 : avma = av; return z;
262 : : }
263 : :
264 : : GEN
265 : 2870 : znstar(GEN N)
266 : : {
267 : 2870 : GEN F = NULL, P, E, cyc, gen, mod;
268 : : long i, j, nbp, sizeh;
269 : 2870 : pari_sp av = avma;
270 : :
271 [ + + ]: 2870 : if ((F = check_arith_all(N,"znstar")))
272 : : {
273 : 14 : F = clean_Z_factor(F);
274 [ - + ]: 14 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
275 : : }
276 [ + + ]: 2870 : if (!signe(N))
277 : : {
278 : 14 : avma = av;
279 : 14 : retmkvec3(gen_2, mkvec(gen_2), mkvec(gen_m1));
280 : : }
281 [ + + ]: 2856 : if (signe(N) < 0) N = absi(N);
282 [ + + ]: 2856 : if (cmpiu(N,2) <= 0)
283 : : {
284 : 7 : avma = av;
285 : 7 : retmkvec3(gen_1, cgetg(1,t_VEC), cgetg(1,t_VEC));
286 : : }
287 [ + + ]: 2849 : if (!F) F = Z_factor(N);
288 : 2849 : P = gel(F,1);
289 : 2849 : E = gel(F,2); nbp = lg(P)-1;
290 : 2849 : cyc = cgetg(nbp+2,t_VEC);
291 : 2849 : gen = cgetg(nbp+2,t_VEC);
292 : 2849 : mod = cgetg(nbp+2,t_VEC);
293 [ + + + + ]: 2849 : switch(mod8(N))
294 : : {
295 : : case 0: {
296 : 147 : long v2 = itos(gel(E,1));
297 : 147 : gel(cyc,1) = int2n(v2-2);
298 : 147 : gel(cyc,2) = gen_2;
299 : 147 : gel(gen,1) = utoipos(5);
300 : 147 : gel(gen,2) = addis(int2n(v2-1), -1);
301 : 147 : gel(mod,1) = gel(mod,2) = int2n(v2);
302 : 147 : sizeh = nbp+1; i = 3; j = 2; break;
303 : : }
304 : : case 4:
305 : 217 : gel(cyc,1) = gen_2;
306 : 217 : gel(gen,1) = utoipos(3);
307 : 217 : gel(mod,1) = utoipos(4);
308 : 217 : sizeh = nbp; i = j = 2; break;
309 : : case 2: case 6:
310 : 21 : sizeh = nbp-1; i=1; j=2; break;
311 : : default: /* 1, 3, 5, 7 */
312 : 2464 : sizeh = nbp; i = j = 1;
313 : : }
314 [ + + ]: 5432 : for ( ; j<=nbp; i++,j++)
315 : : {
316 : 2583 : long e = itos(gel(E,j));
317 : 2583 : GEN p = gel(P,j), q = powiu(p, e-1), Q = mulii(p, q);
318 : 2583 : gel(cyc,i) = subii(Q, q); /* phi(p^e) */
319 [ + + ]: 2583 : gel(gen,i) = e > 1? pgener_Zp(p): pgener_Fp(p);
320 : 2583 : gel(mod,i) = Q;
321 : : }
322 : 2849 : setlg(gen, sizeh+1);
323 : 2849 : setlg(cyc, sizeh+1);
324 [ + + ]: 2849 : if (nbp > 1)
325 [ + + ]: 280 : for (i=1; i<=sizeh; i++)
326 : : {
327 : 203 : GEN Q = gel(mod,i), g = gel(gen,i), qinv = Fp_inv(Q, diviiexact(N,Q));
328 : 203 : g = addii(g, mulii(mulii(subsi(1,g),qinv),Q));
329 : 203 : gel(gen,i) = modii(g, N);
330 : : }
331 : :
332 : : /*The cyc[i] are > 1. They remain so in the loop*/
333 [ + + ]: 3094 : for (i=sizeh; i>=2; i--)
334 : : {
335 : 245 : GEN ci = gel(cyc,i), gi = gel(gen,i);
336 [ + + ]: 595 : for (j=i-1; j>=1; j--) /* we want cyc[i] | cyc[j] */
337 : : {
338 : 350 : GEN cj = gel(cyc,j), gj, qj, v, d;
339 : :
340 : 350 : d = bezout(ci,cj,NULL,&v); /* > 1 */
341 [ + + ]: 350 : if (absi_equal(ci, d)) continue; /* ci | cj */
342 [ + + ]: 147 : if (absi_equal(cj, d)) { /* cj | ci */
343 : 126 : swap(gel(gen,j),gel(gen,i)); gi = gel(gen,i);
344 : 126 : swap(gel(cyc,j),gel(cyc,i)); ci = gel(cyc,i); continue;
345 : : }
346 : :
347 : 21 : gj = gel(gen,j);
348 : 21 : qj = diviiexact(cj,d);
349 : 21 : gel(cyc,j) = mulii(ci,qj);
350 : 21 : gel(cyc,i) = d;
351 : : /* [1,v*cj/d; 0,1]*[1,0;-1,1]*diag(cj,ci)*[ci/d,-v; cj/d,u]
352 : : * = diag(lcm,gcd), with u ci + v cj = d */
353 : :
354 : : /* (gj, gi) *= [1,0; -1,1]^-1 */
355 : 21 : gj = Fp_mul(gj, gi, N); /* order ci*qj = lcm(ci,cj) */
356 : : /* (gj,gi) *= [1,v*qj; 0,1]^-1 */
357 : 21 : togglesign_safe(&v);
358 [ - + ]: 21 : if (signe(v) < 0) v = modii(v,ci); /* >= 0 to avoid inversions */
359 : 21 : gi = Fp_mul(gi, Fp_pow(gj, mulii(qj, v), N), N);
360 : 21 : gel(gen,i) = gi;
361 : 21 : gel(gen,j) = gj;
362 [ + - ]: 350 : ci = d; if (equaliu(ci, 2)) break;
363 : : }
364 : : }
365 : 2870 : return gerepilecopy(av, mkvec3(ZV_prod(cyc), cyc, FpV_to_mod(gen,N)));
366 : : }
367 : :
368 : : /*********************************************************************/
369 : : /** **/
370 : : /** INVERSE TOTIENT FUNCTION **/
371 : : /** **/
372 : : /*********************************************************************/
373 : : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
374 : : * primes. Return factor(N) */
375 : : GEN
376 : 350651 : Z_factor_listP(GEN N, GEN L)
377 : : {
378 : 350651 : long i, k, l = lg(L);
379 : 350651 : GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
380 [ + + ]: 1346688 : for (i = k = 1; i < l; i++)
381 : : {
382 : 996037 : GEN p = gel(L,i);
383 : 996037 : long v = Z_pvalrem(N, p, &N);
384 [ + + ]: 996037 : if (v)
385 : : {
386 : 792176 : gel(P,k) = p;
387 : 792176 : gel(E,k) = utoipos(v);
388 : 792176 : k++;
389 : : }
390 : : }
391 : 350651 : setlg(P, k);
392 : 350651 : setlg(E, k); return mkmat2(P,E);
393 : : }
394 : :
395 : : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
396 : : * L is a list of primes containing all prime divisors of n. */
397 : : static long
398 : 621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
399 : : {
400 : 621565 : pari_sp av = avma, av2;
401 : : GEN k, D;
402 : : long i, v;
403 [ + + ][ + + ]: 621565 : if (m && mod2(n))
404 : : {
405 [ + + ]: 270914 : if (!equali1(n)) return 0;
406 [ + - ]: 69986 : if (px) *px = gen_1;
407 : 69986 : return 1;
408 : : }
409 : 350651 : D = divisors(Z_factor_listP(shifti(n, -1), L));
410 : : /* loop through primes p > m, d = p-1 | n */
411 : 350651 : av2 = avma;
412 [ + + ]: 350651 : if (!m)
413 : : { /* special case p = 2, d = 1 */
414 : 69986 : k = n;
415 : 69986 : for (v = 1;; v++) {
416 [ + - ]: 69986 : if (istotient_i(k, gen_2, L, px)) {
417 [ + - ]: 69986 : if (px) *px = shifti(*px, v);
418 : 69986 : return 1;
419 : : }
420 [ # # ]: 0 : if (mod2(k)) break;
421 : 0 : k = shifti(k,-1);
422 : 0 : }
423 : 0 : avma = av2;
424 : : }
425 [ + + ]: 1099462 : for (i = 1; i < lg(D); ++i)
426 : : {
427 : 1001588 : GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
428 [ + - ][ + + ]: 1001588 : if (m && cmpii(d, m) < 0) continue;
429 : 677782 : p = addiu(d, 1);
430 [ + + ]: 677782 : if (!isprime(p)) continue;
431 : 442064 : k = diviiexact(n, d);
432 : 442064 : for (v = 1;; v++) {
433 : : GEN r;
434 [ + + ]: 481593 : if (istotient_i(k, p, L, px)) {
435 [ + - ]: 182791 : if (px) *px = mulii(*px, powiu(p, v));
436 : 182791 : return 1;
437 : : }
438 : 298802 : k = dvmdii(k, p, &r);
439 [ + + ]: 298802 : if (r != gen_0) break;
440 : 39529 : }
441 : 259273 : avma = av2;
442 : : }
443 : 621565 : avma = av; return 0;
444 : : }
445 : :
446 : : /* find x such that phi(x) = n */
447 : : long
448 : 70000 : istotient(GEN n, GEN *px)
449 : : {
450 : 70000 : pari_sp av = avma;
451 [ - + ]: 70000 : if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
452 [ - + ]: 70000 : if (signe(n) < 1) return 0;
453 [ + + ]: 70000 : if (mod2(n))
454 : : {
455 [ - + ]: 14 : if (!equali1(n)) return 0;
456 [ + - ]: 14 : if (px) *px = gen_1;
457 : 14 : return 1;
458 : : }
459 [ + - ]: 69986 : if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
460 : : {
461 [ - + ]: 69986 : if (!px) avma = av;
462 : : else
463 : 69986 : *px = gerepileuptoint(av, *px);
464 : 69986 : return 1;
465 : : }
466 : 70000 : avma = av; return 0;
467 : : }
468 : :
469 : : /*********************************************************************/
470 : : /** **/
471 : : /** INTEGRAL LOGARITHM **/
472 : : /** **/
473 : : /*********************************************************************/
474 : :
475 : : /* y > 1 and B > 0 integers. Return e such that y^(e-1) <= B < y^e, i.e
476 : : * e = 1 + floor(log_y B). Set *ptq = y^e if non-NULL */
477 : : long
478 : 56550 : logint(GEN B, GEN y, GEN *ptq)
479 : : {
480 : 56550 : pari_sp av = avma, av2;
481 : : long eB, ey, e, i, fl;
482 : 56550 : GEN q,pow2, r = y;
483 : :
484 [ - + ]: 56550 : if (typ(B) != t_INT) B = ceil_safe(B);
485 : 56550 : eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
486 : 56550 : ey = expi(y); /* result e satisfies e > eB / (ey+1) */
487 [ + + ]: 56550 : if (eB <= 13 * ey) /* e small, be naive */
488 : : {
489 : 35734 : for (e=1;; e++)
490 : : { /* here, r = y^e */
491 : 250908 : fl = cmpii(r, B);
492 [ + + ]: 250908 : if (fl > 0) goto END;
493 : 215174 : r = mulii(r,y);
494 : 215174 : }
495 : : }
496 : : /* e > 13 ey / (ey + 1) >= 6.5 */
497 : :
498 : : /* binary splitting: compute bits of e one by one */
499 : : /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
500 : 20816 : pow2 = new_chunk((long)log2(eB)+2);
501 : 20816 : gel(pow2,0) = y;
502 : 20816 : for (i=0,q=r;; )
503 : : { /* r = y^2^i */
504 : 112376 : fl = cmpii(r,B);
505 [ - + ]: 112376 : if (!fl) { e = 1L<<i; e++; r = mulii(r,y); goto END; }
506 [ + + ]: 112376 : if (fl > 0) break;
507 : 91560 : q = r; r = sqri(q);
508 : 91560 : i++; gel(pow2,i) = r;
509 : 91560 : }
510 : :
511 : 20816 : av2 = avma;
512 : 20816 : for (i--, e=1L<<i;;)
513 : : { /* y^e = q < B <= r = q * y^(2^i) */
514 [ + + ]: 91546 : if (!fl) break; /* B = r */
515 : : /* q < B < r */
516 [ + + ][ + + ]: 91539 : if (--i < 0) { if (fl > 0) e++; break; }
517 : 70730 : r = mulii(q, gel(pow2,i));
518 : 70730 : fl = cmpii(r, B);
519 [ + + ]: 70730 : if (fl <= 0) { e += (1L<<i); q = r = gerepileuptoint(av2, r); }
520 : 70730 : }
521 [ + + ]: 20816 : if (fl <= 0) { e++; r = mulii(r,y); }
522 : : END:
523 [ + + ]: 56550 : if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
524 : 56550 : return e;
525 : : }
526 : :
527 : : long
528 : 42 : logint0(GEN B, GEN y, GEN *ptq)
529 : : {
530 : : long e;
531 [ - + ]: 42 : if (typ(B) != t_INT) pari_err_TYPE("logint",B);
532 [ - + ]: 42 : if (signe(B)<=0)
533 : 0 : pari_err_DOMAIN("logint", "x" ,"<=", gen_0, B);
534 [ - + ]: 42 : if (typ(y) != t_INT) pari_err_TYPE("logint",y);
535 [ + - ][ - + ]: 42 : if (signe(y)<=0 || equali1(y))
536 : 0 : pari_err_DOMAIN("logint", "b" ,"<=", gen_1, y);
537 [ + + ]: 42 : if (equaliu(y, 2))
538 : : {
539 : 21 : e = expi(B);
540 [ + + ]: 21 : if (ptq) *ptq = int2n(e);
541 : 21 : return e;
542 : : }
543 : 21 : e = logint(B,y,ptq)-1;
544 [ - + ]: 21 : if (ptq)
545 : : {
546 : 0 : pari_sp av = avma;
547 : 0 : *ptq = gerepileuptoint(av, diviiexact(*ptq, y));
548 : : }
549 : 42 : return e;
550 : : }
551 : :
552 : : /*********************************************************************/
553 : : /** **/
554 : : /** INTEGRAL SQUARE ROOT **/
555 : : /** **/
556 : : /*********************************************************************/
557 : : GEN
558 : 2183 : sqrtint(GEN a)
559 : : {
560 [ - + ]: 2183 : if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
561 [ + + + ]: 2183 : switch (signe(a))
562 : : {
563 : 2169 : case 1: return sqrti(a);
564 : 7 : case 0: return gen_0;
565 : 7 : default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
566 : : }
567 : 2176 : return NULL; /* not reached */
568 : : }
569 : :
570 : : /*********************************************************************/
571 : : /** **/
572 : : /** PERFECT SQUARE **/
573 : : /** **/
574 : : /*********************************************************************/
575 : : static int
576 : 5126695 : carremod(ulong A)
577 : : {
578 : 5126695 : const int carresmod64[]={
579 : : 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,
580 : : 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};
581 : 5126695 : const int carresmod63[]={
582 : : 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,
583 : : 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};
584 : 5126695 : const int carresmod65[]={
585 : : 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,
586 : : 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};
587 : 5126695 : const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
588 : 5126695 : return (carresmod64[A & 0x3fUL]
589 [ + + ]: 2273035 : && carresmod63[A % 63UL]
590 [ + + ]: 1579357 : && carresmod65[A % 65UL]
591 [ + + ][ + + ]: 7399730 : && carresmod11[A % 11UL]);
592 : : }
593 : :
594 : : /* emulate Z_issquareall on single-word integers */
595 : : long
596 : 4961492 : uissquareall(ulong A, ulong *sqrtA)
597 : : {
598 [ - + ]: 4961492 : if (!A) { *sqrtA = 0; return 1; }
599 [ + + ]: 4961492 : if (carremod(A))
600 : : {
601 : 1339261 : ulong a = usqrt(A);
602 [ + + ]: 1339261 : if (a * a == A) { *sqrtA = a; return 1; }
603 : : }
604 : 4961492 : return 0;
605 : : }
606 : : long
607 : 115541 : uissquare(ulong A)
608 : : {
609 [ - + ]: 115541 : if (!A) return 1;
610 [ + + ]: 115541 : if (carremod(A))
611 : : {
612 : 2330 : ulong a = usqrt(A);
613 [ + + ]: 2330 : if (a * a == A) return 1;
614 : : }
615 : 115541 : return 0;
616 : : }
617 : :
618 : : long
619 : 2778876 : Z_issquareall(GEN x, GEN *pt)
620 : : {
621 : : pari_sp av;
622 : : GEN y, r;
623 : :
624 [ + + + ]: 2778876 : switch(signe(x))
625 : : {
626 : 111358 : case -1: return 0;
627 [ - + ]: 35 : case 0: if (pt) *pt=gen_0; return 1;
628 : : }
629 [ + + ]: 2667483 : if (lgefint(x) == 3)
630 : : {
631 : 2617821 : ulong u = (ulong)x[2], a;
632 [ + + ]: 2617821 : if (!pt) return uissquare(u);
633 [ + + ]: 2502280 : if (!uissquareall(u, &a)) return 0;
634 : 2617821 : *pt = utoipos(a); return 1;
635 : : }
636 [ + + ]: 49662 : if (!carremod(umodiu(x, 64*63*65*11))) return 0;
637 : 5782 : av = avma; y = sqrtremi(x, &r);
638 [ + + ]: 5782 : if (r != gen_0) { avma = av; return 0; }
639 [ + + ]: 5202 : if (pt) { *pt = y; avma = (pari_sp)y; } else avma = av;
640 : 2778876 : return 1;
641 : : }
642 : :
643 : : /* a t_INT, p prime */
644 : : long
645 : 0 : Zp_issquare(GEN a, GEN p)
646 : : {
647 : : long v;
648 : : GEN ap;
649 : :
650 [ # # ][ # # ]: 0 : if (!signe(a) || gequal1(a)) return 1;
651 : 0 : v = Z_pvalrem(a, p, &ap);
652 [ # # ]: 0 : if (v&1) return 0;
653 [ # # ]: 0 : return equaliu(p, 2)? umodiu(ap, 8) == 1
654 : 0 : : kronecker(ap,p) == 1;
655 : : }
656 : :
657 : : static int
658 : 535234 : is_char_2(GEN a)
659 : : {
660 : : long j;
661 : : GEN b;
662 [ + + - + : 535234 : switch(typ(a))
+ ]
663 : : {
664 : : case t_INTMOD:
665 : 56 : b = gel(a,1);
666 [ + + ]: 56 : if (!mod2(b))
667 : : {
668 [ - + ]: 21 : if (!equaliu(b, 2)) pari_err_IMPL( "issquare for this input");
669 : 21 : return 1;
670 : : }
671 : 35 : return 0;
672 : : case t_FFELT:
673 [ + + ]: 91 : if (equaliu(FF_p_i(a), 2)) return 1;
674 : 84 : return 0;
675 : : case t_POLMOD:
676 [ # # ][ # # ]: 0 : if (is_char_2(gel(a,1)) || is_char_2(gel(a,2))) return 1;
677 : 0 : return 0;
678 : : case t_POL:
679 [ + + ]: 535206 : for (j = 2; j < lg(a); j++)
680 [ + + ]: 534422 : if (is_char_2(gel(a,j))) return 1;
681 : 784 : return 0;
682 : : }
683 : 535234 : return 0;
684 : : }
685 : :
686 : : static long
687 : 1106 : polissquareall(GEN x, GEN *pt)
688 : : {
689 : : pari_sp av;
690 : 1106 : long v, l = degpol(x);
691 : : GEN y, a, b;
692 : :
693 [ + + ]: 1106 : if (!signe(x))
694 : : {
695 [ + - ]: 7 : if (pt) *pt = gcopy(x);
696 : 7 : return 1;
697 : : }
698 [ + + ]: 1099 : if (pt) *pt = gen_0;
699 [ + + ]: 1099 : if (l&1) return 0; /* odd degree */
700 : 875 : av = avma;
701 : 875 : v = RgX_valrem(x, &x);
702 [ + + ]: 875 : if (v) {
703 : 35 : l = degpol(x);
704 [ + + ]: 35 : if (l&1) return 0;
705 : : }
706 : 868 : a = gel(x,2);
707 [ + + + ]: 868 : switch (typ(a))
708 : : {
709 : : case t_INT:
710 [ - + ]: 805 : if (!Z_issquareall(a,&b)) { avma = av; return 0; }
711 : 805 : break;
712 : : case t_POL:
713 [ - + ]: 14 : if (!polissquareall(a,&b)) { avma = av; return 0; }
714 : 14 : break;
715 : : default:
716 [ - + ]: 49 : if (!issquare(a)) { avma = av; return 0; }
717 : 49 : b = NULL; break;
718 : : }
719 [ + + ]: 868 : if (!l) {
720 [ + + ]: 56 : if (!pt) { avma = av; return 1; }
721 [ - + ]: 28 : if (!b) b = gsqrt(a,DEFAULTPREC);
722 : 28 : y = scalarpol(b, varn(x)); goto END;
723 : : }
724 [ + + ]: 812 : if (is_char_2(x))
725 : : {
726 : : long i, lx;
727 : 28 : x = gmul(x, mkintmod(gen_1, gen_2));
728 : 28 : lx = lg(x);
729 [ + + ]: 28 : if ((lx-3) & 1) { avma = av; return 0; }
730 [ + + ]: 49 : for (i = 3; i < lx; i+=2)
731 [ - + ]: 28 : if (!gequal0(gel(x,i))) { avma = av; return 0; }
732 [ + + ]: 21 : if (pt) {
733 : 14 : y = cgetg((lx+3) / 2, t_POL);
734 [ + + ]: 49 : for (i = 2; i < lx; i+=2)
735 [ - + ]: 35 : if (!issquareall(gel(x,i), &gel(y, (i+2)>>1))) { avma = av; return 0; }
736 : 14 : y[1] = evalsigne(1) | evalvarn(varn(x));
737 : 14 : goto END;
738 : : } else {
739 [ + + ]: 21 : for (i = 2; i < lx; i+=2)
740 [ - + ]: 14 : if (!issquare(gel(x,i))) { avma = av; return 0; }
741 : 7 : avma = av; return 1;
742 : : }
743 : : }
744 : : else
745 : : {
746 : 784 : x = RgX_Rg_div(x,a);
747 : 784 : y = gtrunc(gsqrt(RgX_to_ser(x,2+l),0));
748 [ + + ]: 784 : if (!RgX_equal(gsqr(y), x)) { avma = av; return 0; }
749 [ - + ]: 560 : if (!pt) { avma = av; return 1; }
750 [ + + ]: 560 : if (!gequal1(a))
751 : : {
752 [ + - ]: 14 : if (!b) b = gsqrt(a,DEFAULTPREC);
753 : 14 : y = gmul(b, y);
754 : : }
755 : : }
756 : : END:
757 [ + + ]: 602 : *pt = v? gerepilecopy(av, RgX_shift_shallow(y, v >> 1)): gerepileupto(av, y);
758 : 1106 : return 1;
759 : : }
760 : :
761 : : /* b unit mod p */
762 : : static int
763 : 238 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
764 : : {
765 [ + + ]: 238 : if (d == 1)
766 : : { /* mod p: faster */
767 [ - + ]: 168 : if (!Fp_ispower(b, K, p)) return 0;
768 [ + + ]: 168 : if (pt) *pt = Fp_sqrtn(b, K, p, NULL);
769 : : }
770 : : else
771 : : { /* mod p^{2 +} */
772 [ + + ]: 70 : if (!ispower(cvtop(b, p, d), K, pt)) return 0;
773 [ + + ]: 49 : if (pt) *pt = gtrunc(*pt);
774 : : }
775 : 238 : return 1;
776 : : }
777 : :
778 : : /* We're studying whether a mod (q*p^e) is a K-th power, (q,p) = 1.
779 : : * Decide mod p^e, then reduce a mod q unless q = NULL. */
780 : : static int
781 : 378 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
782 : : {
783 : : GEN t, A;
784 : 378 : long v = Z_pvalrem(*pa, p, &A), d = e - v;
785 [ + + ]: 378 : if (d <= 0) t = gen_0;
786 : : else
787 : : {
788 : : ulong r;
789 : 322 : v = udivui_rem(v, K, &r);
790 [ + + ][ + + ]: 322 : if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
[ + + ]
791 [ + + ][ + + ]: 217 : if (L && v) t = mulii(t, powiu(p, v));
792 : : }
793 [ + - ]: 273 : if (q) *pa = modii(*pa, q);
794 [ + + ]: 273 : if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
795 : 378 : return 1;
796 : : }
797 : : long
798 : 266 : Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
799 : : {
800 : : GEN L, N;
801 : : pari_sp av;
802 : : long e, i, l;
803 : : ulong pp;
804 : : forprime_t S;
805 : :
806 [ + + ]: 266 : if (!signe(a))
807 : : {
808 [ + - ]: 7 : if (pt) {
809 : 7 : GEN t = cgetg(3, t_INTMOD);
810 : 7 : gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
811 : : }
812 : 7 : return 1;
813 : : }
814 : : /* a != 0 */
815 : 259 : av = avma;
816 : :
817 [ - + ]: 259 : if (typ(q) != t_INT) /* integer factorization */
818 : : {
819 : 0 : GEN P = gel(q,1), E = gel(q,2);
820 : 0 : l = lg(P);
821 [ # # ]: 0 : L = pt? vectrunc_init(l): NULL;
822 [ # # ]: 0 : for (i = 1; i < l; i++)
823 : : {
824 : 0 : GEN p = gel(P,i);
825 : 0 : long e = itos(gel(E,i));
826 [ # # ]: 0 : if (!handle_pe(&a, NULL, L, K, p, e)) { avma = av; return 0; }
827 : : }
828 : 0 : goto END;
829 : : }
830 [ + + ]: 259 : if (!mod2(K)
831 [ + + ]: 168 : && kronecker(a, shifti(q,-vali(q))) == -1) { avma = av; return 0; }
832 [ + + ]: 252 : L = pt? vectrunc_init(expi(q)+1): NULL;
833 : 252 : u_forprime_init(&S, 2, tridiv_bound(q));
834 [ + + ]: 883302 : while ((pp = u_forprime_next(&S)))
835 : : {
836 : : int stop;
837 : 883148 : e = Z_lvalrem_stop(&q, pp, &stop);
838 [ + + ]: 883148 : if (!e) continue;
839 [ + + ]: 133 : if (!handle_pe(&a, q, L, K, utoipos(pp), e)) { avma = av; return 0; }
840 [ + + ]: 112 : if (stop)
841 : : {
842 [ + + ][ - + ]: 133 : if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
843 : : goto END;
844 : : }
845 : : }
846 : 154 : l = lg(primetab);
847 [ - + ]: 154 : for (i = 1; i < l; i++)
848 : : {
849 : 0 : GEN p = gel(primetab,i);
850 : 0 : e = Z_pvalrem(q, p, &q);
851 [ # # ]: 0 : if (!e) continue;
852 [ # # ]: 0 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
853 [ # # ]: 0 : if (is_pm1(q)) goto END;
854 : : }
855 : 154 : N = gcdii(a,q);
856 [ + + ]: 154 : if (!is_pm1(N))
857 : : {
858 [ + + ]: 112 : if (ifac_isprime(N))
859 : : {
860 : 70 : e = Z_pvalrem(q, N, &q);
861 [ + - ]: 70 : if (!handle_pe(&a, q, L, K, N, e)) { avma = av; return 0; }
862 : : }
863 : : else
864 : : {
865 : 42 : GEN part = ifac_start(N, 0);
866 : : for(;;)
867 : : {
868 : : long e;
869 : : GEN p;
870 [ + + ]: 84 : if (!ifac_next(&part, &p, &e)) break;
871 : 42 : e = Z_pvalrem(q, p, &q);
872 [ - + ]: 42 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
873 : 84 : }
874 : : }
875 : : }
876 [ + - ]: 84 : if (!is_pm1(q))
877 : : {
878 [ + + ]: 84 : if (ifac_isprime(q))
879 : : {
880 [ - + ]: 28 : if (!handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
881 : : }
882 : : else
883 : : {
884 : 56 : GEN part = ifac_start(q, 0);
885 : : for(;;)
886 : : {
887 : : long e;
888 : : GEN p;
889 [ + + ]: 140 : if (!ifac_next(&part, &p, &e)) break;
890 [ + + ]: 98 : if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
891 : 126 : }
892 : : }
893 : : }
894 : : END:
895 [ + + ]: 147 : if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
896 : 266 : return 1;
897 : : }
898 : :
899 : : long
900 : 144145 : issquareall(GEN x, GEN *pt)
901 : : {
902 : 144145 : long tx = typ(x);
903 : : GEN F;
904 : : pari_sp av;
905 : :
906 [ + + ]: 144145 : if (!pt) return issquare(x);
907 [ + + + + : 1428 : switch(tx)
+ + + - ]
908 : : {
909 : 224 : case t_INT: return Z_issquareall(x, pt);
910 : 77 : case t_FRAC: av = avma;
911 : 77 : F = cgetg(3, t_FRAC);
912 [ + + ]: 77 : if ( !Z_issquareall(gel(x,1), &gel(F,1))
913 [ - + ]: 77 : || !Z_issquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
914 : 28 : *pt = F; return 1;
915 : :
916 : 1036 : case t_POL: return polissquareall(x,pt);
917 : 7 : case t_RFRAC: av = avma;
918 : 7 : F = cgetg(3, t_RFRAC);
919 [ + - ]: 7 : if ( !issquareall(gel(x,1), &gel(F,1))
920 [ - + ]: 7 : || !polissquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
921 : 7 : *pt = F; return 1;
922 : :
923 : : case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
924 [ - + ]: 7 : if (!issquare(x)) return 0;
925 : 7 : *pt = gsqrt(x, DEFAULTPREC); return 1;
926 : :
927 : : case t_INTMOD:
928 : 49 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
929 : :
930 : 28 : case t_FFELT: return FF_issquareall(x, pt);
931 : :
932 : : }
933 : 0 : pari_err_TYPE("issquareall",x);
934 : 144145 : return 0; /* not reached */
935 : : }
936 : :
937 : : long
938 : 142892 : issquare(GEN x)
939 : : {
940 : : pari_sp av;
941 : : GEN a, p;
942 : : long i, v;
943 : :
944 [ + + + + : 142892 : switch(typ(x))
+ + + + +
+ - ]
945 : : {
946 : : case t_INT:
947 : 142619 : return Z_issquare(x);
948 : :
949 : : case t_REAL:
950 : 14 : return (signe(x)>=0);
951 : :
952 : : case t_INTMOD:
953 : 70 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, NULL);
954 : :
955 : : case t_FRAC:
956 [ + - ][ + + ]: 21 : return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
957 : :
958 : 21 : case t_FFELT: return FF_issquareall(x, NULL);
959 : :
960 : : case t_COMPLEX:
961 : 7 : return 1;
962 : :
963 : : case t_PADIC:
964 [ - + ]: 63 : a = gel(x,4); if (!signe(a)) return 1;
965 [ + + ]: 63 : if (valp(x)&1) return 0;
966 : 49 : p = gel(x,2);
967 [ + + ]: 49 : if (!equaliu(p, 2)) return (kronecker(a,p) != -1);
968 : :
969 : 14 : v = precp(x); /* here p=2, a is odd */
970 [ + - ][ + + ]: 14 : if ((v>=3 && mod8(a) != 1 ) ||
[ - + ]
971 [ # # ]: 7 : (v==2 && mod4(a) != 1)) return 0;
972 : 7 : return 1;
973 : :
974 : : case t_POL:
975 : 49 : return polissquareall(x,NULL);
976 : :
977 : : case t_SER:
978 [ + + ]: 21 : if (!signe(x)) return 1;
979 [ + + ]: 14 : if (valp(x)&1) return 0;
980 : 7 : return issquare(gel(x,2));
981 : :
982 : : case t_RFRAC:
983 : 7 : av = avma; i = issquare(gmul(gel(x,1),gel(x,2)));
984 : 7 : avma = av; return i;
985 : : }
986 : 0 : pari_err_TYPE("issquare",x);
987 : 142892 : return 0; /* not reached */
988 : : }
989 [ # # ]: 0 : GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
990 [ # # ]: 0 : GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
991 : :
992 : : long
993 : 1386 : ispolygonal(GEN x, GEN S, GEN *N)
994 : : {
995 : 1386 : pari_sp av = avma;
996 : : GEN D, d, n;
997 [ - + ]: 1386 : if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
998 [ - + ]: 1386 : if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
999 [ - + ]: 1386 : if (cmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
1000 [ - + ]: 1386 : if (signe(x) < 0) return 0;
1001 [ + + ][ + - ]: 1386 : if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
1002 [ + + ][ + - ]: 1260 : if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
1003 : : /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
1004 [ + + ]: 1134 : if (cmpiu(S, 1<<16) < 0) /* common case ! */
1005 : : {
1006 : 441 : ulong s = S[2], r;
1007 [ + + ]: 441 : if (s == 4) return Z_issquareall(x, N);
1008 [ - + ]: 378 : if (s == 3)
1009 : 0 : D = addiu(shifti(x, 3), 1);
1010 : : else
1011 : 378 : D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
1012 [ - + ]: 378 : if (!Z_issquareall(D, &d)) { avma = av; return 0; }
1013 [ - + ]: 378 : if (s == 3)
1014 : 0 : d = subiu(d, 1);
1015 : : else
1016 : 378 : d = addiu(d, s - 4);
1017 : 378 : n = diviu_rem(d, 2*s - 4, &r);
1018 [ - + ]: 441 : if (r) { avma = av; return 0; }
1019 : : }
1020 : : else
1021 : : {
1022 : 693 : GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
1023 : 693 : D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
1024 [ - + ]: 693 : if (!Z_issquareall(D, &d)) { avma = av; return 0; }
1025 : 693 : d = addii(d, S_4);
1026 : 693 : n = dvmdii(shifti(d,-1), S_2, &r);
1027 [ - + ]: 693 : if (r != gen_0) { avma = av; return 0; }
1028 : : }
1029 [ + - ]: 1071 : if (N) *N = gerepileuptoint(av, n); else avma = av;
1030 : 1386 : return 1;
1031 : : }
1032 : :
1033 : : /*********************************************************************/
1034 : : /** **/
1035 : : /** PERFECT POWER **/
1036 : : /** **/
1037 : : /*********************************************************************/
1038 : : static long
1039 : 56 : polispower(GEN x, GEN K, GEN *pt)
1040 : : {
1041 : : pari_sp av;
1042 : 56 : long v, l = degpol(x), k = itos(K);
1043 : : GEN y, a, b;
1044 : :
1045 [ + + ]: 56 : if (!signe(x)) return 1;
1046 [ + + ]: 49 : if (l % k) return 0; /* degree not multiple of k */
1047 : 42 : v = RgX_valrem(x, &x);
1048 [ + + ]: 42 : if (v % k) return 0;
1049 : 35 : av = avma; a = gel(x,2); b = NULL;
1050 [ + + ]: 35 : if (!ispower(a, K, &b)) { avma = av; return 0; }
1051 : 21 : av = avma;
1052 [ + + ]: 21 : if (degpol(x))
1053 : : {
1054 : 7 : x = RgX_Rg_div(x,a);
1055 : 7 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
1056 [ - + ]: 7 : if (!RgX_equal(powgi(y, K), x)) { avma = av; return 0; }
1057 : : }
1058 : 14 : else y = pol_1(varn(x));
1059 [ + - ]: 21 : if (pt)
1060 : : {
1061 [ + - ]: 21 : if (!gequal1(a))
1062 : : {
1063 [ - + ]: 21 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
1064 : 21 : y = gmul(b,y);
1065 : : }
1066 [ + + ]: 21 : *pt = v? gerepilecopy(av, RgX_shift_shallow(y, v/k)): gerepileupto(av, y);
1067 : : }
1068 : 0 : else avma = av;
1069 : 56 : return 1;
1070 : : }
1071 : :
1072 : : long
1073 : 693 : Z_ispowerall(GEN x, ulong k, GEN *pt)
1074 : : {
1075 : 693 : long s = signe(x);
1076 : : ulong mask;
1077 [ - + ][ # # ]: 693 : if (!s) { if (pt) *pt = gen_0; return 1; }
1078 [ + + ]: 693 : if (s > 0) {
1079 [ + + ]: 658 : if (k == 2) return Z_issquareall(x, pt);
1080 [ + + ]: 364 : if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
1081 [ + + ]: 91 : if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
1082 [ + + ]: 77 : if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
1083 : 70 : return is_kth_power(x, k, pt);
1084 : : }
1085 [ + + ]: 35 : if (!odd(k)) return 0;
1086 [ + - ]: 28 : if (Z_ispowerall(absi(x), k, pt))
1087 : : {
1088 [ + + ]: 28 : if (pt) *pt = negi(*pt);
1089 : 28 : return 1;
1090 : : };
1091 : 693 : return 0;
1092 : : }
1093 : :
1094 : : /* is x a K-th power mod p ? Assume p prime. */
1095 : : int
1096 : 168 : Fp_ispower(GEN x, GEN K, GEN p)
1097 : : {
1098 : 168 : pari_sp av = avma;
1099 : : GEN p_1;
1100 : : long r;
1101 : 168 : x = modii(x, p);
1102 [ + - ][ + + ]: 168 : if (!signe(x) || equali1(x)) { avma = av; return 1; }
1103 : : /* implies p > 2 */
1104 : 91 : p_1 = subiu(p,1);
1105 : 91 : K = gcdii(K, p_1);
1106 [ + + ]: 91 : if (equaliu(K, 2)) { r = kronecker(x,p); avma = av; return (r > 0); }
1107 : 35 : x = Fp_pow(x, diviiexact(p_1,K), p);
1108 : 168 : avma = av; return equali1(x);
1109 : : }
1110 : :
1111 : : /* x unit defined modulo 2^e, e > 0, p prime */
1112 : : static int
1113 : 2373 : U2_issquare(GEN x, long e)
1114 : : {
1115 [ - + ]: 2373 : long r = signe(x)>=0?mod8(x):8-mod8(x);
1116 [ - + ]: 2373 : if (e==1) return 1;
1117 [ + + ]: 2373 : if (e==2) return (r&3L) == 1;
1118 : 2373 : return r == 1;
1119 : : }
1120 : : /* x unit defined modulo p^e, e > 0, p prime */
1121 : : static int
1122 : 4690 : Up_issquare(GEN x, GEN p, long e)
1123 [ + + ]: 4690 : { return (equaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
1124 : :
1125 : : long
1126 : 2548 : Zn_issquare(GEN d, GEN fn)
1127 : : {
1128 : : long j, np;
1129 [ - + ]: 2548 : if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
1130 [ - + ]: 2548 : if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
1131 : : /* integer factorization */
1132 : 2548 : np = nbrows(fn);
1133 [ + + ]: 5320 : for (j = 1; j <= np; ++j)
1134 : : {
1135 : 4970 : GEN r, p = gcoeff(fn, j, 1);
1136 : 4970 : long e = itos(gcoeff(fn, j, 2));
1137 : 4970 : long v = Z_pvalrem(d,p,&r);
1138 [ + + ][ + - ]: 4970 : if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
[ + + ]
1139 : : }
1140 : 2548 : return 1;
1141 : : }
1142 : :
1143 : : long
1144 : 7001603 : ispower(GEN x, GEN K, GEN *pt)
1145 : : {
1146 : : GEN z;
1147 : :
1148 [ + + ]: 7001603 : if (!K) return gisanypower(x, pt);
1149 [ - + ]: 1449 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
1150 [ - + ]: 1449 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
1151 [ - + ][ # # ]: 1449 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
1152 [ + + + + : 1449 : switch(typ(x)) {
+ + + + +
+ - ]
1153 : : case t_INT:
1154 : 77 : return Z_ispowerall(x, itou(K), pt);
1155 : : case t_FRAC:
1156 : : {
1157 : 21 : GEN a = gel(x,1), b = gel(x,2);
1158 : 21 : ulong k = itou(K);
1159 [ + + ]: 21 : if (pt) {
1160 : 14 : z = cgetg(3, t_FRAC);
1161 [ + - ][ + - ]: 14 : if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
1162 : 14 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1163 : : }
1164 : 0 : avma = (pari_sp)(z + 3); return 0;
1165 : : }
1166 [ + - ][ + - ]: 21 : return Z_ispower(a, k) && Z_ispower(b, k);
1167 : : }
1168 : : case t_INTMOD:
1169 : 147 : return Zn_ispower(gel(x,2), gel(x,1), K, pt);
1170 : : case t_FFELT:
1171 : 28 : return FF_ispower(x, K, pt);
1172 : :
1173 : : case t_PADIC:
1174 : 1099 : z = Qp_sqrtn(x, K, NULL);
1175 [ + + ]: 1099 : if (!z) return 0;
1176 [ + + ]: 805 : if (pt) *pt = z;
1177 : 805 : return 1;
1178 : :
1179 : : case t_POL:
1180 : 49 : return polispower(x, K, pt);
1181 : : case t_RFRAC: {
1182 : 7 : GEN a = gel(x,1), b = gel(x,2);
1183 [ + - ]: 7 : if (pt) {
1184 : 7 : z = cgetg(3, t_RFRAC);
1185 [ + - ][ + - ]: 7 : if (ispower(a, K, &a) && polispower(b, K, &b)) {
1186 : 7 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
1187 : : }
1188 : 0 : avma = (pari_sp)(z + 3); return 0;
1189 : : }
1190 [ # # ][ # # ]: 7 : return (ispower(a, K, NULL) && polispower(b, K, NULL));
1191 : : }
1192 : : case t_REAL:
1193 [ + - ][ - + ]: 7 : if (signe(x) < 0 && !mpodd(K)) return 0;
1194 : : case t_COMPLEX:
1195 [ - + ]: 14 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1196 : 14 : return 1;
1197 : :
1198 : : case t_SER:
1199 [ + - ][ + - ]: 7 : if (signe(x) && (!dvdsi(valp(x), K) || !ispower(gel(x,2), K, NULL)))
[ - + ]
1200 : 0 : return 0;
1201 [ + - ]: 7 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
1202 : 7 : return 1;
1203 : :
1204 : 0 : default: pari_err_TYPE("ispower",x);
1205 : 7001603 : return 0; /* not reached */
1206 : : }
1207 : : }
1208 : :
1209 : : long
1210 : 7000154 : gisanypower(GEN x, GEN *pty)
1211 : : {
1212 : 7000154 : long tx = typ(x);
1213 : : ulong k, h;
1214 [ + + ]: 7000154 : if (tx == t_INT) return Z_isanypower(x, pty);
1215 [ + - ]: 14 : if (tx == t_FRAC)
1216 : : {
1217 : 14 : pari_sp av = avma;
1218 : 14 : GEN fa, P, E, a = gel(x,1), b = gel(x,2);
1219 : : long i, j, p, e;
1220 : 14 : int sw = (absi_cmp(a, b) > 0);
1221 : :
1222 [ - + ]: 14 : if (sw) swap(a, b);
1223 [ + + ]: 14 : k = Z_isanypower(a, pty? &a: NULL);
1224 [ + + ]: 14 : if (!k)
1225 : : { /* a = -1,1 or not a pure power */
1226 [ - + ]: 7 : if (!is_pm1(a)) { avma = av; return 0; }
1227 [ + - ]: 7 : if (signe(a) < 0) b = negi(b);
1228 [ + - ]: 7 : k = Z_isanypower(b, pty? &b: NULL);
1229 [ + - ][ - + ]: 7 : if (!k || !pty) { avma = av; return k; }
1230 : 7 : *pty = gerepileupto(av, ginv(b));
1231 : 7 : return k;
1232 : : }
1233 : 7 : fa = factoru(k);
1234 : 7 : P = gel(fa,1);
1235 : 7 : E = gel(fa,2); h = k;
1236 [ + + ]: 14 : for (i = lg(P) - 1; i > 0; i--)
1237 : : {
1238 : 7 : p = P[i];
1239 : 7 : e = E[i];
1240 [ + + ]: 21 : for (j = 0; j < e; j++)
1241 [ - + ]: 14 : if (!is_kth_power(b, p, &b)) break;
1242 [ - + ]: 7 : if (j < e) k /= upowuu(p, e - j);
1243 : : }
1244 [ - + ]: 7 : if (k == 1) { avma = av; return 0; }
1245 [ + - ]: 7 : if (!pty) { avma = av; return k; }
1246 [ # # ]: 0 : if (k != h) a = powiu(a, h/k);
1247 : 0 : *pty = gerepilecopy(av, mkfrac(a, b));
1248 : 14 : return k;
1249 : : }
1250 : 0 : pari_err_TYPE("gisanypower", x);
1251 : 7000154 : return 0; /* not reached */
1252 : : }
1253 : :
1254 : : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
1255 : : * No need to optimize for 2,3,5,7 powers (done before) */
1256 : : static long
1257 : 505708 : split_exponent(ulong e, GEN *x)
1258 : : {
1259 : : GEN fa, P, E;
1260 : 505708 : long i, j, l, k = 1;
1261 [ + + ]: 505708 : if (e == 1) return 1;
1262 : 7 : fa = factoru(e);
1263 : 7 : P = gel(fa,1);
1264 : 7 : E = gel(fa,2); l = lg(P);
1265 [ + + ]: 14 : for (i = 1; i < l; i++)
1266 : : {
1267 : 7 : ulong p = P[i];
1268 [ + + ]: 14 : for (j = 0; j < E[i]; j++)
1269 : : {
1270 : : GEN y;
1271 [ + - ]: 7 : if (!is_kth_power(*x, p, &y)) break;
1272 : 7 : k *= p; *x = y;
1273 : : }
1274 : : }
1275 : 505708 : return k;
1276 : : }
1277 : :
1278 : : static long
1279 : 864264 : Z_isanypower_nosmalldiv(GEN *px)
1280 : : { /* any prime divisor of x is > 102 */
1281 : 864264 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
1282 : 864264 : const double LOG103 = 4.6347; /* lower bound for log(103) */
1283 : : forprime_t T;
1284 : 864264 : ulong mask = 7, e2;
1285 : : long k, ex;
1286 : 864264 : GEN y, x = *px;
1287 : :
1288 : 864264 : k = 1;
1289 [ + + ]: 865419 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
1290 [ + + ]: 864355 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
1291 : 864264 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
1292 [ + + ]: 864264 : if (u_forprime_init(&T, 11, e2))
1293 : : {
1294 : 16786 : GEN logx = NULL;
1295 : 16786 : const ulong Q = 30011; /* prime */
1296 : : ulong p, xmodQ;
1297 : 16786 : double dlogx = 0;
1298 : : /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
1299 : : * for large p the modular checks are no longer competitively fast */
1300 [ + + ]: 16807 : while ( (ex = is_pth_power(x, &y, &T, 30)) )
1301 : : {
1302 : 21 : k *= ex; x = y;
1303 : 21 : e2 = (ulong)((expi(x) + 1) / LOG2_103);
1304 : 21 : u_forprime_restrict(&T, e2);
1305 : : }
1306 [ - + ]: 16786 : if (DEBUGLEVEL>4) err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x));
1307 : 16786 : xmodQ = umodiu(x, Q);
1308 : : /* test Q | x, just in case */
1309 [ + + ]: 16786 : if (!xmodQ) return k * split_exponent(Z_lval(x,Q), px);
1310 : : /* x^(1/p) < 2^31 */
1311 : 16779 : p = T.p;
1312 [ + + ]: 16779 : if (p <= e2)
1313 : : {
1314 : 16765 : logx = logr_abs( itor(x, DEFAULTPREC) );
1315 : 16765 : dlogx = rtodbl(logx);
1316 : 16765 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1317 : : }
1318 [ + + ][ + + ]: 134617 : while (p && p <= e2)
1319 : : { /* is x a p-th power ? By computing y = round(x^(1/p)).
1320 : : * Check whether y^p = x, first mod Q, then exactly. */
1321 : 117838 : pari_sp av = avma;
1322 : : long e;
1323 : 117838 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
1324 : 117838 : ulong ymodQ = umodiu(y,Q);
1325 [ + + ][ + + ]: 117838 : if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
1326 [ - + ]: 21 : || !equalii(powiu(y, p), x)) avma = av;
1327 : : else
1328 : : {
1329 : 21 : k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
1330 : 21 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1331 : 21 : u_forprime_restrict(&T, e2);
1332 : 21 : continue; /* if success, retry same p */
1333 : : }
1334 : 117817 : p = u_forprime_next(&T);
1335 : : }
1336 : : }
1337 : 864264 : *px = x; return k;
1338 : : }
1339 : :
1340 : : static ulong tinyprimes[] = {
1341 : : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
1342 : : 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
1343 : : 157, 163, 167, 173, 179, 181, 191, 193, 197, 199
1344 : : };
1345 : :
1346 : : /* disregard the sign of x, caller will take care of x < 0 */
1347 : : static long
1348 : 7001729 : Z_isanypower_aux(GEN x, GEN *pty)
1349 : : {
1350 : : long ex, v, i, l, k;
1351 : : GEN y, P, E;
1352 : 7001729 : ulong mask, e = 0;
1353 : :
1354 [ + + ]: 7001729 : if (absi_cmp(x, gen_2) < 0) return 0; /* -1,0,1 */
1355 : :
1356 : 7001715 : k = l = 1;
1357 : 7001715 : P = cgetg(26 + 1, t_VECSMALL);
1358 : 7001715 : E = cgetg(26 + 1, t_VECSMALL);
1359 : : /* trial division */
1360 [ + + ]: 67137931 : for(i = 0; i < 26; i++)
1361 : : {
1362 : 60136216 : ulong p = tinyprimes[i];
1363 : : int stop;
1364 : 60136216 : v = Z_lvalrem_stop(&x, p, &stop);
1365 [ + + ]: 60136216 : if (v)
1366 : : {
1367 : 7923174 : P[l] = p;
1368 : 7923174 : E[l] = v; l++;
1369 [ + + ]: 7923174 : e = ugcd(e, v); if (e == 1) goto END;
1370 : : }
1371 [ + + ]: 54717887 : if (stop) {
1372 [ + + ]: 249032 : if (is_pm1(x)) k = e;
1373 : : goto END;
1374 : : }
1375 : : }
1376 : :
1377 [ + + ]: 1334354 : if (e)
1378 : : { /* Bingo. Result divides e */
1379 : : long v3, v5, v7;
1380 : 505701 : ulong e2 = e;
1381 : 505701 : v = u_lvalrem(e2, 2, &e2);
1382 [ + + ]: 505701 : if (v)
1383 : : {
1384 [ + + ]: 375249 : for (i = 0; i < v; i++)
1385 : : {
1386 [ + + ]: 374171 : if (!Z_issquareall(x, &y)) break;
1387 : 1288 : k <<= 1; x = y;
1388 : : }
1389 : : }
1390 : 505701 : mask = 0;
1391 [ + + ]: 505701 : v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
1392 [ + + ]: 505701 : v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
1393 [ + + ]: 505701 : v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
1394 [ + + ]: 505778 : while ( (ex = is_357_power(x, &y, &mask)) ) {
1395 : 77 : x = y;
1396 [ + + + - ]: 77 : switch(ex)
1397 : : {
1398 [ + - ]: 28 : case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
1399 [ + - ]: 28 : case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
1400 [ + - ]: 21 : case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
1401 : : }
1402 : : }
1403 : 505701 : k *= split_exponent(e2, &x);
1404 : : }
1405 : : else
1406 : 828653 : k = Z_isanypower_nosmalldiv(&x);
1407 : : END:
1408 [ + + ][ + + ]: 7001715 : if (pty && k != 1)
1409 : : {
1410 [ + + ]: 8869 : if (e)
1411 : : { /* add missing small factors */
1412 : 7616 : y = powuu(P[1], E[1] / k);
1413 [ + + ]: 14770 : for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
1414 [ + + ]: 7616 : x = equali1(x)? y: mulii(x,y);
1415 : : }
1416 : 8869 : *pty = x;
1417 : : }
1418 [ + + ]: 7001729 : return k == 1? 0: k;
1419 : : }
1420 : :
1421 : : long
1422 : 7001729 : Z_isanypower(GEN x, GEN *pty)
1423 : : {
1424 : 7001729 : pari_sp av = avma;
1425 : 7001729 : long k = Z_isanypower_aux(x, pty);
1426 [ + + ]: 7001729 : if (!k) { avma = av; return 0; }
1427 [ + + ]: 8904 : if (signe(x) < 0)
1428 : : {
1429 : 35 : long v = vals(k);
1430 [ + + ]: 35 : if (v)
1431 : : {
1432 : : GEN y;
1433 : 21 : k >>= v;
1434 [ + + ]: 21 : if (k == 1) { avma = av; return 0; }
1435 [ - + ]: 14 : if (!pty) { avma = av; return k; }
1436 : 14 : y = *pty;
1437 : 14 : y = powiu(y, 1<<v);
1438 : 14 : togglesign(y);
1439 : 14 : *pty = gerepileuptoint(av, y);
1440 : 14 : return k;
1441 : : }
1442 [ + - ]: 14 : if (pty) togglesign_safe(pty);
1443 : : }
1444 [ + + ]: 8883 : if (pty) *pty = gerepilecopy(av, *pty); else avma = av;
1445 : 7001729 : return k;
1446 : : }
1447 : :
1448 : : /* Faster than */
1449 : : /* !cmpii(n, int2n(vali(n))) */
1450 : : /* !cmpis(shifti(n, -vali(n)), 1) */
1451 : : /* expi(n) == vali(n) */
1452 : : /* hamming(n) == 1 */
1453 : : /* even for single-word values, and much faster for multiword values. */
1454 : : /* If all you have is a word, you can just use n & !(n & (n-1)). */
1455 : : long
1456 : 119 : Z_ispow2(GEN n)
1457 : : {
1458 : : GEN xp;
1459 : : long i, lx;
1460 : : ulong u;
1461 [ - + ]: 119 : if (signe(n) != 1) return 0;
1462 : 119 : xp = int_LSW(n);
1463 : 119 : lx = lgefint(n);
1464 : 119 : u = *xp;
1465 [ + + ]: 135 : for (i = 3; i < lx; ++i)
1466 : : {
1467 [ + + ]: 24 : if (u) return 0;
1468 : 16 : xp = int_nextW(xp);
1469 : 16 : u = *xp;
1470 : : }
1471 : 119 : return !(u & (u-1)); /* faster than hamming_word(u) == 1 */
1472 : : }
1473 : :
1474 : : long
1475 : 840091 : isprimepower(GEN n, GEN *pt)
1476 : : {
1477 : 840091 : pari_sp av = avma;
1478 : : long i, v;
1479 : :
1480 [ - + ]: 840091 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1481 [ - + ]: 840091 : if (signe(n) <= 0) return 0;
1482 : :
1483 [ + + ]: 840091 : if (lgefint(n) == 3)
1484 : : {
1485 : : ulong p;
1486 : 540054 : v = uisprimepower(n[2], &p);
1487 [ + + ]: 540054 : if (v)
1488 : : {
1489 [ + + ]: 53842 : if (pt) *pt = utoipos(p);
1490 : 53842 : return v;
1491 : : }
1492 : 540054 : return 0;
1493 : : }
1494 [ + + ]: 1655397 : for (i = 0; i < 26; i++)
1495 : : {
1496 : 1619786 : ulong p = tinyprimes[i];
1497 : 1619786 : v = Z_lvalrem(n, p, &n);
1498 [ + + ]: 1619786 : if (v)
1499 : : {
1500 : 264426 : avma = av;
1501 [ + + ]: 264426 : if (!is_pm1(n)) return 0;
1502 [ + + ]: 30 : if (pt) *pt = utoipos(p);
1503 : 30 : return v;
1504 : : }
1505 : : }
1506 : : /* p | n => p >= 103 */
1507 : 35611 : v = Z_isanypower_nosmalldiv(&n); /* expensive */
1508 [ + + ]: 35611 : if (!isprime(n)) { avma = av; return 0; }
1509 [ + + ]: 5299 : if (pt) *pt = gerepilecopy(av, n); else avma = av;
1510 : 840091 : return v;
1511 : : }
1512 : :
1513 : : long
1514 : 540061 : uisprimepower(ulong n, ulong *pp)
1515 : : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1516 : : * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1517 : 540061 : const ulong CUTOFF = 200UL;
1518 : 540061 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1519 : 540061 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1520 : : #ifdef LONG_IS_64BIT
1521 : : /* primes preceeding the appropriate root of ULONG_MAX. */
1522 : 480054 : const ulong ROOT9 = 137;
1523 : 480054 : const ulong ROOT8 = 251;
1524 : 480054 : const ulong ROOT7 = 563;
1525 : 480054 : const ulong ROOT5 = 7129;
1526 : 480054 : const ulong ROOT4 = 65521;
1527 : : #else
1528 : 60007 : const ulong ROOT9 = 11;
1529 : 60007 : const ulong ROOT8 = 13;
1530 : 60007 : const ulong ROOT7 = 23;
1531 : 60007 : const ulong ROOT5 = 83;
1532 : 60007 : const ulong ROOT4 = 251;
1533 : : #endif
1534 : : ulong mask;
1535 : : long v, i;
1536 : : int e;
1537 [ + + ]: 540061 : if (n < 2) return 0;
1538 [ + + ]: 540047 : if (!odd(n)) {
1539 [ + + ]: 270047 : if (n & (n-1)) return 0;
1540 : 264 : *pp = 2; return vals(n);
1541 : : }
1542 [ + + ]: 270000 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1543 [ + + ]: 3647666 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1544 : : {
1545 : 3588672 : ulong p = tinyprimes[i];
1546 [ + + ]: 3588672 : if (n % p == 0)
1547 : : {
1548 : 210964 : v = u_lvalrem(n, p, &n);
1549 [ + + ]: 210964 : if (n == 1) { *pp = p; return v; }
1550 : 209438 : return 0;
1551 : : }
1552 : : }
1553 : : /* p | n => p >= CUTOFF */
1554 : :
1555 [ + + ]: 58994 : if (n < CUTOFF3)
1556 : : {
1557 [ + + ][ + - ]: 46298 : if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1558 [ # # ]: 0 : if (uissquareall(n, &n)) { *pp = n; return 2; }
1559 : 0 : return 0;
1560 : : }
1561 : :
1562 : : /* Check for squares, fourth powers, and eighth powers as appropriate. */
1563 : 12696 : v = 1;
1564 [ - + ]: 12696 : if (uissquareall(n, &n)) {
1565 : 0 : v <<= 1;
1566 [ # # ][ # # ]: 0 : if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1567 : 0 : v <<= 1;
1568 [ # # ][ # # ]: 0 : if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1569 : : }
1570 : : }
1571 : :
1572 [ - + ]: 12696 : if (CUTOFF > ROOT5) mask = 1;
1573 : : else
1574 : : {
1575 : 12696 : const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1576 [ + - ]: 12696 : if (n < CUTOFF5) mask = 1; else mask = 3;
1577 [ + - ]: 12696 : if (CUTOFF <= ROOT7)
1578 : : {
1579 : 12696 : const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1580 [ - + ]: 12696 : if (n >= CUTOFF7) mask = 7;
1581 : : }
1582 : : }
1583 : :
1584 [ - + ][ # # ]: 12696 : if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1585 [ - + ]: 12696 : if ((e = uis_357_power(n, &n, &mask))) v *= e;
1586 : :
1587 [ + + ]: 12696 : if (uisprime_101(n)) { *pp = n; return v; }
1588 : 540061 : return 0;
1589 : : }
1590 : :
1591 : : /*********************************************************************/
1592 : : /** **/
1593 : : /** KRONECKER SYMBOL **/
1594 : : /** **/
1595 : : /*********************************************************************/
1596 : : /* t = 3,5 mod 8 ? (= 2 not a square mod t) */
1597 : : static int
1598 : 19487896 : ome(long t)
1599 : : {
1600 [ + + ]: 19487896 : switch(t & 7)
1601 : : {
1602 : : case 3:
1603 : 11523992 : case 5: return 1;
1604 : 19487896 : default: return 0;
1605 : : }
1606 : : }
1607 : : /* t a t_INT, is t = 3,5 mod 8 ? */
1608 : : static int
1609 : 2309886 : gome(GEN t)
1610 [ + - ]: 2309886 : { return signe(t)? ome( mod2BIL(t) ): 0; }
1611 : :
1612 : : /* t a t_INT, return 1 if t = 3 (mod 4), 0 otherwise */
1613 : : static int
1614 : 56 : eps(GEN t)
1615 : : {
1616 [ - + - ]: 56 : switch(signe(t))
1617 : : {
1618 : 0 : case -1: return mod4(t) == 1;
1619 : 56 : case 1: return mod4(t) == 3;
1620 : 56 : default: return 0;
1621 : : }
1622 : : }
1623 : :
1624 : : /* assume y odd, return kronecker(x,y) * s */
1625 : : static long
1626 : 20106508 : krouu_s(ulong x, ulong y, long s)
1627 : : {
1628 : 20106508 : ulong x1 = x, y1 = y, z;
1629 [ + + ]: 65081918 : while (x1)
1630 : : {
1631 : 44975410 : long r = vals(x1);
1632 [ + + ]: 44975410 : if (r)
1633 : : {
1634 [ + + ][ + + ]: 24243300 : if (odd(r) && ome(y1)) s = -s;
1635 : 24243300 : x1 >>= r;
1636 : : }
1637 [ + + ]: 44975410 : if (x1 & y1 & 2) s = -s;
1638 : 44975410 : z = y1 % x1; y1 = x1; x1 = z;
1639 : : }
1640 [ + + ]: 20106508 : return (y1 == 1)? s: 0;
1641 : : }
1642 : :
1643 : : long
1644 : 4614366 : kronecker(GEN x, GEN y)
1645 : : {
1646 : 4614366 : pari_sp av = avma, lim;
1647 : 4614366 : long s = 1, r;
1648 : : ulong xu, yu;
1649 : :
1650 [ - + ]: 4614366 : if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
1651 [ - + ]: 4614366 : if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
1652 [ - - + ]: 4614366 : switch (signe(y))
1653 : : {
1654 [ # # ]: 0 : case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
1655 : 0 : case 0: return is_pm1(x);
1656 : : }
1657 : 4614366 : r = vali(y);
1658 [ + + ]: 4614366 : if (r)
1659 : : {
1660 [ + + ]: 448 : if (!mpodd(x)) { avma = av; return 0; }
1661 [ + - ][ + + ]: 254 : if (odd(r) && gome(x)) s = -s;
1662 : 254 : y = shifti(y,-r);
1663 : : }
1664 : 4614172 : lim = stack_lim(av,2);
1665 : 4614172 : x = modii(x,y);
1666 [ + + ]: 4670317 : while (lgefint(x) > 3) /* x < y */
1667 : : {
1668 : : GEN z;
1669 : 56145 : r = vali(x);
1670 [ + + ]: 56145 : if (r)
1671 : : {
1672 [ + + ][ + + ]: 30664 : if (odd(r) && gome(y)) s = -s;
1673 : 30664 : x = shifti(x,-r);
1674 : : }
1675 : : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1676 [ + + ]: 56145 : if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
1677 : 56145 : z = remii(y,x); y = x; x = z;
1678 [ - + ]: 56145 : if (low_stack(lim, stack_lim(av,2)))
1679 : : {
1680 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
1681 : 0 : gerepileall(av, 2, &x, &y);
1682 : : }
1683 : : }
1684 : 4614172 : xu = itou(x);
1685 [ + + ][ + + ]: 4614172 : if (!xu) return is_pm1(y)? s: 0;
1686 : 4604833 : r = vals(xu);
1687 [ + + ]: 4604833 : if (r)
1688 : : {
1689 [ + + ][ + + ]: 3216263 : if (odd(r) && gome(y)) s = -s;
1690 : 3216263 : xu >>= r;
1691 : : }
1692 : : /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1693 [ + + ]: 4604833 : if (xu & mod2BIL(y) & 2) s = -s;
1694 : 4604833 : yu = umodiu(y, xu);
1695 : 4614366 : avma = av; return krouu_s(yu, xu, s);
1696 : : }
1697 : :
1698 : : long
1699 : 633 : krois(GEN x, long y)
1700 : : {
1701 : : ulong yu;
1702 : 633 : long s = 1, r;
1703 : :
1704 [ - + ]: 633 : if (y <= 0)
1705 : : {
1706 [ # # ]: 0 : if (y == 0) return is_pm1(x);
1707 [ # # ]: 0 : yu = (ulong)-y; if (signe(x) < 0) s = -1;
1708 : : }
1709 : : else
1710 : 633 : yu = (ulong)y;
1711 : 633 : r = vals(yu);
1712 [ - + ]: 633 : if (r)
1713 : : {
1714 [ # # ]: 0 : if (!mpodd(x)) return 0;
1715 [ # # ][ # # ]: 0 : if (odd(r) && gome(x)) s = -s;
1716 : 0 : yu >>= r;
1717 : : }
1718 : 633 : return krouu_s(umodiu(x, yu), yu, s);
1719 : : }
1720 : : /* assume y != 0 */
1721 : : long
1722 : 1128312 : kroiu(GEN x, ulong y)
1723 : : {
1724 : 1128312 : long s = 1, r = vals(y);
1725 [ + + ]: 1128312 : if (r)
1726 : : {
1727 [ + + ]: 9877 : if (!mpodd(x)) return 0;
1728 [ + + ][ + + ]: 6825 : if (odd(r) && gome(x)) s = -s;
1729 : 6825 : y >>= r;
1730 : : }
1731 : 1128312 : return krouu_s(umodiu(x, y), y, s);
1732 : : }
1733 : :
1734 : : long
1735 : 25109 : krosi(long x, GEN y)
1736 : : {
1737 : 25109 : const pari_sp av = avma;
1738 : 25109 : long s = 1, r;
1739 : : ulong u, xu;
1740 : :
1741 [ - - + ]: 25109 : switch (signe(y))
1742 : : {
1743 [ # # ]: 0 : case -1: y = negi(y); if (x < 0) s = -1; break;
1744 [ # # ][ # # ]: 0 : case 0: return (x==1 || x==-1);
1745 : : }
1746 : 25109 : r = vali(y);
1747 [ - + ]: 25109 : if (r)
1748 : : {
1749 [ # # ]: 0 : if (!odd(x)) { avma = av; return 0; }
1750 [ # # ][ # # ]: 0 : if (odd(r) && ome(x)) s = -s;
1751 : 0 : y = shifti(y,-r);
1752 : : }
1753 [ + + ][ + + ]: 25109 : if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
1754 : 25109 : xu = (ulong)x;
1755 [ + + ]: 25109 : if (lgefint(y) == 3)
1756 : 17922 : return krouu_s(xu, itou(y), s);
1757 [ - + ]: 7187 : if (!xu) return 0; /* y != 1 */
1758 : 7187 : r = vals(xu);
1759 [ + + ]: 7187 : if (r)
1760 : : {
1761 [ + + ][ + + ]: 4995 : if (odd(r) && gome(y)) s = -s;
1762 : 4995 : xu >>= r;
1763 : : }
1764 : : /* xu=3 mod 4 && y=3 mod 4 ? (both are odd here) */
1765 [ + + ]: 7187 : if (xu & mod2BIL(y) & 2) s = -s;
1766 : 7187 : u = umodiu(y, xu);
1767 : 25109 : avma = av; return krouu_s(u, xu, s);
1768 : : }
1769 : :
1770 : : long
1771 : 64701 : kross(long x, long y)
1772 : : {
1773 : : ulong yu;
1774 : 64701 : long s = 1, r;
1775 : :
1776 [ + + ]: 64701 : if (y <= 0)
1777 : : {
1778 [ - + ]: 378 : if (y == 0) return (labs(x)==1);
1779 [ - + ]: 378 : yu = (ulong)-y; if (x < 0) s = -1;
1780 : : }
1781 : : else
1782 : 64323 : yu = (ulong)y;
1783 : 64701 : r = vals(yu);
1784 [ - + ]: 64701 : if (r)
1785 : : {
1786 [ # # ]: 0 : if (!odd(x)) return 0;
1787 [ # # ][ # # ]: 0 : if (odd(r) && ome(x)) s = -s;
1788 : 0 : yu >>= r;
1789 : : }
1790 [ + + ]: 64701 : x %= (long)yu; if (x < 0) x += yu;
1791 : 64701 : return krouu_s((ulong)x, yu, s);
1792 : : }
1793 : :
1794 : : long
1795 : 14285972 : krouu(ulong x, ulong y)
1796 : : {
1797 : : long r;
1798 [ + + ]: 14285972 : if (y & 1) return krouu_s(x, y, 1);
1799 [ - + ]: 960 : if (!odd(x)) return 0;
1800 : 960 : r = vals(y);
1801 [ + - ][ - + ]: 14285972 : return krouu_s(x, y >> r, (odd(r) && ome(x))? -1: 1);
1802 : : }
1803 : :
1804 : : /*********************************************************************/
1805 : : /** **/
1806 : : /** HILBERT SYMBOL **/
1807 : : /** **/
1808 : : /*********************************************************************/
1809 : : /* x,y are t_INT or t_REAL */
1810 : : static long
1811 : 49 : mphilbertoo(GEN x, GEN y)
1812 : : {
1813 : 49 : long sx = signe(x), sy = signe(y);
1814 [ + - ][ - + ]: 49 : if (!sx || !sy) return 0;
1815 [ + + ][ + - ]: 49 : return (sx < 0 && sy < 0)? -1: 1;
1816 : : }
1817 : :
1818 : : long
1819 : 182 : hilbertii(GEN x, GEN y, GEN p)
1820 : : {
1821 : : pari_sp av;
1822 : : long oddvx, oddvy, z;
1823 : :
1824 [ + + ]: 182 : if (!p) return mphilbertoo(x,y);
1825 [ + - ][ - + ]: 154 : if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
1826 [ + + ][ + + ]: 154 : if (!signe(x) || !signe(y)) return 0;
1827 : 133 : av = avma;
1828 : 133 : oddvx = odd(Z_pvalrem(x,p,&x));
1829 : 133 : oddvy = odd(Z_pvalrem(y,p,&y));
1830 : : /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
1831 [ + + ]: 133 : if (equaliu(p, 2))
1832 : : {
1833 [ + + ][ + + ]: 42 : z = (eps(x) && eps(y))? -1: 1;
1834 [ + + ][ - + ]: 42 : if (oddvx && gome(y)) z = -z;
1835 [ + + ][ - + ]: 42 : if (oddvy && gome(x)) z = -z;
1836 : : }
1837 : : else
1838 : : {
1839 [ + + ][ - + ]: 91 : z = (oddvx && oddvy && eps(p))? -1: 1;
[ # # ]
1840 [ + + ][ - + ]: 91 : if (oddvx && kronecker(y,p) < 0) z = -z;
1841 [ - + ][ # # ]: 91 : if (oddvy && kronecker(x,p) < 0) z = -z;
1842 : : }
1843 : 182 : avma = av; return z;
1844 : : }
1845 : :
1846 : : static void
1847 : 196 : err_prec(void) { pari_err_PREC("hilbert"); }
1848 : : static void
1849 : 161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
1850 : : static void
1851 : 56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
1852 : :
1853 : : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
1854 : : * Return lift(x) provided it's p-adic accuracy is large enough to decide
1855 : : * hilbert()'s value [ problem at p = 2 ] */
1856 : : static GEN
1857 : 434 : lift_intmod(GEN x, GEN *pp)
1858 : : {
1859 : 434 : GEN p = *pp, N = gel(x,1);
1860 : 434 : x = gel(x,2);
1861 [ + + ]: 434 : if (!p)
1862 : : {
1863 : 273 : *pp = p = N;
1864 [ + + ]: 273 : switch(itos_or_0(p))
1865 : : {
1866 : : case 2:
1867 : 126 : case 4: err_prec();
1868 : : }
1869 : 147 : return x;
1870 : : }
1871 [ + + ]: 161 : if (!signe(p)) err_oo(N);
1872 [ + + ]: 119 : if (equaliu(p,2))
1873 [ + + ]: 42 : { if (vali(N) <= 2) err_prec(); }
1874 : : else
1875 [ + + ]: 77 : { if (!dvdii(N,p)) err_p(N,p); }
1876 [ + + ]: 35 : if (!signe(x)) err_prec();
1877 : 175 : return x;
1878 : : }
1879 : : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
1880 : : * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
1881 : : * to decide hilbert()'s value [ problem at p = 2 ]*/
1882 : : static GEN
1883 : 210 : lift_padic(GEN x, GEN *pp)
1884 : : {
1885 : 210 : GEN p = *pp, q = gel(x,2), y = gel(x,4);
1886 [ + + ]: 210 : if (!p) *pp = p = q;
1887 [ + + ]: 147 : else if (!equalii(p,q)) err_p(p, q);
1888 [ + + ][ + + ]: 105 : if (equaliu(p,2) && precp(x) <= 2) err_prec();
1889 [ - + ]: 70 : if (!signe(y)) err_prec();
1890 [ - + ]: 70 : return odd(valp(x))? mulii(p,y): y;
1891 : : }
1892 : :
1893 : : long
1894 : 672 : hilbert(GEN x, GEN y, GEN p)
1895 : : {
1896 : 672 : pari_sp av = avma;
1897 : 672 : long tx = typ(x), ty = typ(y), z;
1898 : :
1899 [ + + ][ - + ]: 672 : if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
1900 [ + + ]: 672 : if (tx == t_REAL)
1901 : : {
1902 [ + + ][ + + ]: 77 : if (p && signe(p)) err_oo(p);
1903 [ + - + ]: 63 : switch (ty)
1904 : : {
1905 : : case t_INT:
1906 : 7 : case t_REAL: return mphilbertoo(x,y);
1907 : 0 : case t_FRAC: return mphilbertoo(x,gel(y,1));
1908 : 56 : default: pari_err_TYPE2("hilbert",x,y);
1909 : : }
1910 : : }
1911 [ + + ]: 595 : if (ty == t_REAL)
1912 : : {
1913 [ - + ][ # # ]: 14 : if (p && signe(p)) err_oo(p);
1914 [ + - - ]: 14 : switch (tx)
1915 : : {
1916 : : case t_INT:
1917 : 14 : case t_REAL: return mphilbertoo(x,y);
1918 : 0 : case t_FRAC: return mphilbertoo(gel(x,1),y);
1919 : 0 : default: pari_err_TYPE2("hilbert",x,y);
1920 : : }
1921 : : }
1922 [ + + ]: 581 : if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
1923 [ + + ]: 378 : if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
1924 : :
1925 [ + + ]: 322 : if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
1926 [ + + ]: 259 : if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
1927 : :
1928 [ + + ][ + - ]: 182 : if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
1929 [ + + ][ + - ]: 182 : if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
1930 : :
1931 [ + - ][ - + ]: 182 : if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
1932 [ + + ][ + + ]: 182 : if (p && !signe(p)) p = NULL;
1933 : 203 : z = hilbertii(x,y,p); avma = av; return z;
1934 : : }
1935 : :
1936 : : /*******************************************************************/
1937 : : /* */
1938 : : /* SQUARE ROOT MODULO p */
1939 : : /* */
1940 : : /*******************************************************************/
1941 : :
1942 : : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
1943 : : ulong
1944 : 4998994 : Fl_sqrt(ulong a, ulong p)
1945 : : {
1946 : : long i, e, k;
1947 : 4998994 : ulong pi = get_Fl_red(p);
1948 : : ulong p1, q, v, y, w, m;
1949 : :
1950 [ + + ]: 4998994 : if (!a) return 0;
1951 : 4936284 : p1 = p - 1; e = vals(p1);
1952 [ + + ]: 4936284 : if (e == 0) /* p = 2 */
1953 : : {
1954 [ + + ]: 28 : if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1955 : 21 : return ((a & 1) == 0)? 0: 1;
1956 : : }
1957 : 4936256 : q = p1 >> e; /* q = (p-1)/2^oo is odd */
1958 [ + + ]: 4936256 : if (e == 1) y = p1;
1959 : : else /* look for an odd power of a primitive root */
1960 : 3065478 : for (k=2; ; k++)
1961 : : { /* loop terminates for k < p (even if p composite) */
1962 : 7537016 : i = krouu(k, p);
1963 [ + + ]: 7537016 : if (i >= 0)
1964 : : {
1965 [ + + ]: 4471545 : if (i) continue;
1966 : 7 : pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
1967 : : }
1968 : 3065471 : y = m = Fl_powu_pre(k, q, p, pi);
1969 [ + + ]: 9265879 : for (i=1; i<e; i++)
1970 [ - + ]: 6200408 : if ((m = Fl_sqr_pre(m, p, pi)) == 1) break;
1971 [ + - ]: 3065471 : if (i == e) break; /* success */
1972 : 4471538 : }
1973 : :
1974 : 4936249 : p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
1975 [ - + ]: 4936249 : if (!p1) return 0;
1976 : 4936249 : v = Fl_mul_pre(a, p1, p, pi);
1977 : 4936249 : w = Fl_mul_pre(v, p1, p, pi);
1978 [ + + ]: 7425527 : while (w != 1)
1979 : : { /* a*w = v^2, y primitive 2^e-th root of 1
1980 : : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
1981 : 2489331 : p1 = Fl_sqr_pre(w,p,pi);
1982 [ + + ][ + + ]: 4986173 : for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
1983 [ + + ]: 2489331 : if (k == e) return ~0UL;
1984 : : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
1985 : 2489278 : p1 = y;
1986 [ + + ]: 3748610 : for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
1987 : 2489278 : y = Fl_sqr_pre(p1, p, pi); e = k;
1988 : 2489278 : w = Fl_mul_pre(y, w, p, pi);
1989 : 2489278 : v = Fl_mul_pre(v, p1, p, pi);
1990 : : }
1991 [ + + ]: 4936196 : p1 = p - v; if (v > p1) v = p1;
1992 : 4998980 : return v;
1993 : : }
1994 : :
1995 : : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
1996 : : * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
1997 : : * and in average is about 2 or 2.5 times worse. But try both algorithms for
1998 : : * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
1999 : : *
2000 : : * If X^2 := t^2 - a is not a square in F_p (so X is in F_p^2), then
2001 : : * (t+X)^(p+1) = (t-X)(t+X) = a, hence sqrt(a) = (t+X)^((p+1)/2) in F_p^2.
2002 : : * If (a|p)=1, then sqrt(a) is in F_p.
2003 : : * cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */
2004 : :
2005 : : /* compute y^2, y = y[1] + y[2] X */
2006 : : static GEN
2007 : 449 : sqrt_Cipolla_sqr(void *data, GEN y)
2008 : : {
2009 : 449 : GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
2010 : 449 : GEN u2 = sqri(u), v2 = sqri(v);
2011 : 449 : v = subii(sqri(addii(v,u)), addii(u2,v2));
2012 : 449 : u = addii(u2, mulii(v2,n));
2013 : : /* NOT mkvec2: must be gerepileupto-able */
2014 : 449 : retmkvec2(modii(u,p), modii(v,p));
2015 : : }
2016 : : /* compute (t+X) y^2 */
2017 : : static GEN
2018 : 23 : sqrt_Cipolla_msqr(void *data, GEN y)
2019 : : {
2020 : 23 : GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
2021 : 23 : ulong t = gt[2];
2022 : 23 : GEN d = addii(u, mului(t,v)), d2= sqri(d);
2023 : 23 : GEN b = remii(mulii(a,v), p);
2024 : 23 : u = subii(mului(t,d2), mulii(b,addii(u,d)));
2025 : 23 : v = subii(d2, mulii(b,v));
2026 : : /* NOT mkvec2: must be gerepileupto-able */
2027 : 23 : retmkvec2(modii(u,p), modii(v,p));
2028 : : }
2029 : : /* assume a reduced mod p [ otherwise correct but inefficient ] */
2030 : : static GEN
2031 : 8 : sqrt_Cipolla(GEN a, GEN p)
2032 : : {
2033 : : pari_sp av1;
2034 : : GEN u, v, n, y, pov2;
2035 : : ulong t;
2036 : :
2037 [ - + ]: 8 : if (kronecker(a, p) < 0) return NULL;
2038 : 8 : pov2 = shifti(p,-1);
2039 [ + - ]: 8 : if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
2040 : :
2041 : 8 : av1 = avma;
2042 : 8 : for(t=1; ; t++)
2043 : : {
2044 : 41 : n = subsi((long)(t*t), a);
2045 [ + + ]: 41 : if (kronecker(n, p) < 0) break;
2046 : 33 : avma = av1;
2047 : 33 : }
2048 : :
2049 : : /* compute (t+X)^((p-1)/2) =: u+vX */
2050 : 8 : u = utoipos(t);
2051 : 8 : y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
2052 : : sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
2053 : : /* Now u+vX = (t+X)^((p-1)/2); thus
2054 : : * (u+vX)(t+X) = sqrt(a) + 0 X
2055 : : * Whence,
2056 : : * sqrt(a) = (u+vt)t - v*a
2057 : : * 0 = (u+vt)
2058 : : * Thus a square root is v*a */
2059 : :
2060 : 8 : v = Fp_mul(gel(y, 2), a, p);
2061 [ + + ]: 8 : if (cmpii(v,pov2) > 0) v = subii(p,v);
2062 : 8 : return v;
2063 : : }
2064 : :
2065 : : #define sqrmod(x,p) (remii(sqri(x),p))
2066 : :
2067 : : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
2068 : : GEN
2069 : 1516028 : Fp_sqrt(GEN a, GEN p)
2070 : : {
2071 : 1516028 : pari_sp av = avma, av1,lim;
2072 : : long i, k, e;
2073 : : GEN p1, q, v, y, w, m;
2074 : :
2075 [ - + ]: 1516028 : if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
2076 [ - + ]: 1516028 : if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
2077 [ + - ][ - + ]: 1516028 : if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
2078 [ + + ]: 1516028 : if (lgefint(p) == 3)
2079 : : {
2080 : 1510920 : ulong u = (ulong)p[2]; u = Fl_sqrt(umodiu(a, u), u);
2081 [ + + ]: 1510906 : if (u == ~0UL) return NULL;
2082 : 1510871 : return utoi(u);
2083 : : }
2084 : :
2085 : 5108 : p1 = addsi(-1,p); e = vali(p1);
2086 : 5108 : a = modii(a, p);
2087 : :
2088 : : /* On average, the algorithm of Cipolla is better than the algorithm of
2089 : : * Tonelli and Shanks if and only if e(e-1)>8*log2(n)+20
2090 : : * see LNCS 2286 pp 430 [GTL] */
2091 [ + + ]: 5108 : if (e*(e-1) > 20 + 8 * expi(p))
2092 : : {
2093 : 8 : v = sqrt_Cipolla(a,p);
2094 [ - + ]: 8 : if (!v) { avma = av; return NULL; }
2095 : 8 : return gerepileuptoint(av,v);
2096 : : }
2097 : :
2098 [ - + ]: 5100 : if (e == 0) /* p = 2 */
2099 : : {
2100 : 0 : avma = av;
2101 [ # # ]: 0 : if (!equaliu(p,2)) pari_err_PRIME("Fp_sqrt [modulus]",p);
2102 [ # # ][ # # ]: 0 : if (!signe(a) || !mod2(a)) return gen_0;
2103 : 0 : return gen_1;
2104 : : }
2105 : 5100 : q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
2106 [ + + ]: 5100 : if (e == 1) y = p1;
2107 : : else /* look for an odd power of a primitive root */
2108 : 2774 : for (k=2; ; k++)
2109 : : { /* loop terminates for k < p (even if p composite) */
2110 : :
2111 : 3746 : i = krosi(k,p);
2112 [ + + ]: 3746 : if (i >= 0)
2113 : : {
2114 [ + - ]: 972 : if (i) continue;
2115 : 0 : pari_err_PRIME("Fp_sqrt [modulus]",p);
2116 : : }
2117 : 2774 : av1 = avma;
2118 : 2774 : y = m = Fp_pow(utoipos((ulong)k),q,p);
2119 [ + + ]: 6503 : for (i=1; i<e; i++)
2120 [ - + ]: 3729 : if (gequal1(m = sqrmod(m,p))) break;
2121 [ + - ]: 2774 : if (i == e) break; /* success */
2122 : 0 : avma = av1;
2123 : 972 : }
2124 : :
2125 : 5100 : p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ [(q-1)/2] */
2126 [ - + ]: 5100 : if (!signe(p1)) { avma=av; return gen_0; }
2127 : 5100 : v = Fp_mul(a, p1, p);
2128 : 5100 : w = Fp_mul(v, p1, p);
2129 : 5100 : lim = stack_lim(av,1);
2130 [ + + ]: 7021 : while (!equali1(w))
2131 : : { /* a*w = v^2, y primitive 2^e-th root of 1
2132 : : a square --> w even power of y, hence w^(2^(e-1)) = 1 */
2133 : 1924 : p1 = sqrmod(w,p);
2134 [ + + ][ + + ]: 2821 : for (k=1; !equali1(p1) && k < e; k++) p1 = sqrmod(p1,p);
2135 [ + + ]: 1924 : if (k == e) { avma=av; return NULL; } /* p composite or (a/p) != 1 */
2136 : : /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
2137 : 1921 : p1 = y;
2138 [ + + ]: 2266 : for (i=1; i < e-k; i++) p1 = sqrmod(p1,p);
2139 : 1921 : y = sqrmod(p1, p); e = k;
2140 : 1921 : w = Fp_mul(y, w, p);
2141 : 1921 : v = Fp_mul(v, p1, p);
2142 [ - + ]: 1921 : if (low_stack(lim, stack_lim(av,1)))
2143 : : {
2144 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
2145 : 0 : gerepileall(av,3, &y,&w,&v);
2146 : : }
2147 : : }
2148 : 5097 : av1 = avma;
2149 [ + + ]: 5097 : p1 = subii(p,v); if (cmpii(v,p1) > 0) v = p1; else avma = av1;
2150 : 1516014 : return gerepileuptoint(av, v);
2151 : : }
2152 : :
2153 : : /*********************************************************************/
2154 : : /** **/
2155 : : /** GCD & BEZOUT **/
2156 : : /** **/
2157 : : /*********************************************************************/
2158 : :
2159 : : GEN
2160 : 2265806 : lcmii(GEN x, GEN y)
2161 : : {
2162 : : pari_sp av;
2163 : : GEN a, b;
2164 [ + - ][ - + ]: 2265806 : if (!signe(x) || !signe(y)) return gen_0;
2165 : 2265806 : av = avma;
2166 [ + + ]: 2265806 : a = gcdii(x,y); if (!equali1(a)) y = diviiexact(y,a);
2167 : 2265806 : b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
2168 : : }
2169 : :
2170 : : /*********************************************************************/
2171 : : /** **/
2172 : : /** CHINESE REMAINDERS **/
2173 : : /** **/
2174 : : /*********************************************************************/
2175 : :
2176 : : /* P.M. & M.H.
2177 : : *
2178 : : * Chinese Remainder Theorem. x and y must have the same type (integermod,
2179 : : * polymod, or polynomial/vector/matrix recursively constructed with these
2180 : : * as coefficients). Creates (with the same type) a z in the same residue
2181 : : * class as x and the same residue class as y, if it is possible.
2182 : : *
2183 : : * We also allow (during recursion) two identical objects even if they are
2184 : : * not integermod or polymod. For example, if
2185 : : *
2186 : : * x = [1. mod(5, 11), mod(X + mod(2, 7), X^2 + 1)]
2187 : : * y = [1, mod(7, 17), mod(X + mod(0, 3), X^2 + 1)],
2188 : : *
2189 : : * then chinese(x, y) returns
2190 : : *
2191 : : * [1, mod(16, 187), mod(X + mod(9, 21), X^2 + 1)]
2192 : : *
2193 : : * Someone else may want to allow power series, complex numbers, and
2194 : : * quadratic numbers.
2195 : : */
2196 : :
2197 : : GEN
2198 : 35 : chinese1(GEN x) { return gassoc_proto(chinese,x,NULL); }
2199 : :
2200 : : GEN
2201 : 1554 : chinese(GEN x, GEN y)
2202 : : {
2203 : : pari_sp av,tetpil;
2204 : 1554 : long tx = typ(x);
2205 : : GEN z,p1,p2,d,u,v;
2206 : :
2207 [ + + ]: 1554 : if (!y) return chinese1(x);
2208 [ + + ]: 1519 : if (gequal(x,y)) return gcopy(x);
2209 [ + - ][ + + : 1512 : if (tx == typ(y)) switch(tx)
+ + - ]
2210 : : {
2211 : : case t_POLMOD:
2212 : : {
2213 : 7 : GEN A = gel(x,1), B = gel(y,1);
2214 : 7 : GEN a = gel(x,2), b = gel(y,2);
2215 : 7 : z = cgetg(3, t_POLMOD);
2216 [ - + ]: 7 : if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
2217 [ - + ]: 7 : if (RgX_equal(A,B)) /* same modulus */
2218 : : {
2219 : 0 : gel(z,1) = gcopy(A);
2220 : 0 : gel(z,2) = chinese(a,b);
2221 : 0 : return z;
2222 : : }
2223 : 7 : av = avma;
2224 : 7 : d = RgX_extgcd(A,B,&u,&v);
2225 : 7 : p2 = gsub(b, a);
2226 [ - + ]: 7 : if (!gequal0(gmod(p2, d))) break;
2227 : 7 : p1 = gdiv(A,d);
2228 : 7 : p2 = gadd(a, gmul(gmul(u,p1), p2));
2229 : :
2230 : 7 : tetpil = avma;
2231 : 7 : gel(z,1) = gmul(p1,B);
2232 : 7 : gel(z,2) = gmod(p2,gel(z,1));
2233 : 7 : gerepilecoeffssp(av,tetpil,z+1,2); return z;
2234 : : }
2235 : : case t_INTMOD:
2236 : : {
2237 : 1491 : GEN A = gel(x,1), B = gel(y,1);
2238 : 1491 : GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
2239 : 1491 : z = cgetg(3,t_INTMOD);
2240 : 1491 : Z_chinese_pre(A, B, &C, &U, &d);
2241 : 1491 : c = Z_chinese_post(a, b, C, U, d);
2242 [ - + ]: 1491 : if (!c) pari_err_OP("chinese", x,y);
2243 : 1491 : gel(z,1) = icopy_avma(C, (pari_sp)z);
2244 : 1491 : gel(z,2) = icopy_avma(c, (pari_sp)gel(z,1));
2245 : 1491 : avma = (pari_sp)gel(z,2); return z;
2246 : : }
2247 : : case t_POL:
2248 : : {
2249 : 7 : long i, lx = lg(x), ly = lg(y);
2250 [ - + ]: 7 : if (varn(x) != varn(y)) break;
2251 [ + - ]: 7 : if (lx < ly) { swap(x,y); lswap(lx,ly); }
2252 : 7 : z = cgetg(lx, t_POL); z[1] = x[1];
2253 [ + + ]: 21 : for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2254 [ + + ]: 14 : for ( ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
2255 : 7 : return z;
2256 : : }
2257 : :
2258 : : case t_VEC: case t_COL: case t_MAT:
2259 : : {
2260 : : long i, lx;
2261 [ + - ]: 7 : z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
2262 [ + + ]: 21 : for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
2263 : 7 : return z;
2264 : : }
2265 : : }
2266 : 0 : pari_err_OP("chinese",x,y);
2267 : 1554 : return NULL; /* not reached */
2268 : : }
2269 : :
2270 : : /* init chinese(Mod(.,A), Mod(.,B)) */
2271 : : void
2272 : 195470 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
2273 : : {
2274 : 195470 : GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
2275 : 195470 : GEN t = diviiexact(A,d);
2276 : 195470 : *pU = mulii(u, t);
2277 : 195470 : *pC = mulii(t, B);
2278 [ + + ]: 195470 : if (pd) *pd = d;
2279 : 195470 : }
2280 : : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
2281 : : * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
2282 : : * If d not NULL, check wether a = b mod d. */
2283 : : GEN
2284 : 268168 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
2285 : : {
2286 : : GEN b_a;
2287 [ + + ]: 268168 : if (!signe(a))
2288 : : {
2289 [ + + ][ - + ]: 212694 : if (d && remii(b, d) != gen_0) return NULL;
2290 : 212694 : return Fp_mul(b, U, C);
2291 : : }
2292 : 55474 : b_a = subii(b,a);
2293 [ + + ][ - + ]: 55474 : if (d && remii(b_a, d) != gen_0) return NULL;
2294 : 268168 : return modii(addii(a, mulii(U, b_a)), C);
2295 : : }
2296 : : GEN
2297 : 1988 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
2298 : : {
2299 : 1988 : pari_sp av = avma;
2300 : 1988 : GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
2301 : 1988 : return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
2302 : : }
2303 : : GEN
2304 : 191984 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
2305 : : {
2306 : 191984 : GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
2307 : 191984 : return Z_chinese_post(a,b, *pC, U, NULL);
2308 : : }
2309 : :
2310 : : /* return lift(chinese(a mod A, b mod B))
2311 : : * assume(A,B)=1, a,b,A,B integers and C = A*B */
2312 : : GEN
2313 : 1050 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
2314 : : {
2315 : 1050 : pari_sp av = avma;
2316 : 1050 : GEN U = mulii(Fp_inv(A,B), A);
2317 : 1050 : return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2318 : : }
2319 : :
2320 : : /* chinese1 for coprime moduli in Z */
2321 : : static GEN
2322 : 71641 : chinese1_coprime_Z_aux(GEN x, GEN y)
2323 : : {
2324 : 71641 : GEN z = cgetg(3, t_INTMOD);
2325 : 71641 : GEN A = gel(x,1), a = gel(x, 2);
2326 : 71641 : GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
2327 : 71641 : pari_sp av = avma;
2328 : 71641 : GEN U = mulii(Fp_inv(A,B), A);
2329 : 71641 : gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
2330 : 71641 : gel(z,1) = C; return z;
2331 : : }
2332 : : GEN
2333 : 87489 : chinese1_coprime_Z(GEN x) {return gassoc_proto(chinese1_coprime_Z_aux,x,NULL);}
2334 : :
2335 : : /*********************************************************************/
2336 : : /** **/
2337 : : /** MODULAR EXPONENTIATION **/
2338 : : /** **/
2339 : : /*********************************************************************/
2340 : :
2341 : : /* modified Barrett reduction with one fold */
2342 : : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
2343 : :
2344 : : static GEN
2345 : 11262 : Fp_invmBarrett(GEN p, long s)
2346 : : {
2347 : 11262 : GEN R, Q = dvmdii(int2n(3*s),p,&R);
2348 : 11262 : return mkvec2(Q,R);
2349 : : }
2350 : :
2351 : : static GEN
2352 : 338508 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN p)
2353 : : {
2354 : 338508 : pari_sp av = avma;
2355 : 338508 : GEN Q = gel(B, 1), R = gel(B, 2);
2356 : 338508 : long sQ = expi(Q);
2357 : 338508 : GEN A = addii(remi2n(a, 3*s), mulii(R,shifti(a, -3*s)));
2358 : 338508 : GEN q = shifti(mulii(shifti(A, sQ-3*s), Q), -sQ);
2359 : 338508 : GEN r = subii(A, mulii(q, p));
2360 : 338508 : GEN sr= subii(r,p); /* Now 0 <= r < 4*p */
2361 [ + + ]: 338508 : if (signe(sr)<0) return gerepileuptoint(av, r);
2362 : 323673 : r=sr; sr = subii(r,p); /* Now 0 <= r < 3*p */
2363 [ + + ]: 323673 : if (signe(sr)<0) return gerepileuptoint(av, r);
2364 : 164007 : r=sr; sr = subii(r,p); /* Now 0 <= r < 2*p */
2365 [ + + ]: 338508 : return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
2366 : : }
2367 : :
2368 : : /* Montgomery reduction */
2369 : :
2370 : : INLINE ulong
2371 : 247659 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
2372 : :
2373 : : typedef struct muldata {
2374 : : GEN N;
2375 : : GEN iM;
2376 : : ulong inv, s;
2377 : : GEN (*res)(struct muldata *,GEN);
2378 : : GEN (*mul2)(struct muldata *,GEN);
2379 : : } muldata;
2380 : :
2381 : : /* Montgomery reduction */
2382 : : static GEN
2383 : 11006523 : _montred(muldata *D, GEN x)
2384 : : {
2385 : 11006523 : return red_montgomery(x, D->N, D->inv);
2386 : : }
2387 : :
2388 : : static GEN
2389 : 2731852 : _remii(muldata *D, GEN x) { return remii(x, D->N); }
2390 : :
2391 : : static GEN
2392 : 338508 : _remiibar(muldata *D, GEN x) { return Fp_rem_mBarrett(x, D->iM, D->s, D->N); }
2393 : :
2394 : : /* 2x mod N */
2395 : : static GEN
2396 : 1893202 : _muli2red(muldata *D, GEN x)
2397 : : {
2398 : 1893202 : GEN z = shifti(x,1);
2399 [ + + ]: 1893202 : return (cmpii(z,D->N) >= 0)? subii(z,D->N): z;
2400 : : }
2401 : : static GEN
2402 : 1583412 : _muli2montred(muldata *D, GEN x)
2403 : : {
2404 : 1583412 : GEN z = _muli2red(D,x);
2405 : 1583412 : long l = lgefint(D->N);
2406 [ + + ]: 1585155 : while (lgefint(z) > l) z = subii(z,D->N);
2407 : 1583412 : return z;
2408 : : }
2409 : : static GEN
2410 : 755526 : _mul(void *data, GEN x, GEN y)
2411 : : {
2412 : 755526 : muldata *D = (muldata *)data;
2413 : 755526 : return D->res(D, mulii(x,y));
2414 : : }
2415 : : static GEN
2416 : 11180496 : _sqr(void *data, GEN x)
2417 : : {
2418 : 11180496 : muldata *D = (muldata *)data;
2419 : 11180496 : return D->res(D, sqri(x));
2420 : : }
2421 : : static GEN
2422 : 1893202 : _m2sqr(void *data, GEN x)
2423 : : {
2424 : 1893202 : muldata *D = (muldata *)data;
2425 : 1893202 : return D->mul2(D, D->res(D, sqri(x)));
2426 : : }
2427 : :
2428 : : ulong
2429 : 9205409 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
2430 : : {
2431 : : ulong y, z, n;
2432 [ + + ]: 9205409 : if (n0 <= 1)
2433 : : { /* frequent special cases */
2434 [ + + ]: 225989 : if (n0 == 1) return x;
2435 [ + - ]: 88675 : if (n0 == 0) return 1;
2436 : : }
2437 [ + + ]: 8979420 : if (x <= 1) return x; /* 0 or 1 */
2438 : 8951963 : y = 1; z = x; n = n0;
2439 : : for(;;)
2440 : : {
2441 [ + + ]: 103635628 : if (n&1) y = Fl_mul_pre(y,z,p,pi);
2442 [ + + ]: 103635628 : n>>=1; if (!n) return y;
2443 : 94683665 : z = Fl_sqr_pre(z,p,pi);
2444 : 103889074 : }
2445 : : }
2446 : :
2447 : : ulong
2448 : 59667625 : Fl_powu(ulong x, ulong n0, ulong p)
2449 : : {
2450 : : ulong y, z, n;
2451 [ + + ]: 59667625 : if (n0 <= 2)
2452 : : { /* frequent special cases */
2453 [ + + ]: 48101121 : if (n0 == 2) return Fl_sqr(x,p);
2454 [ + + ]: 4477552 : if (n0 == 1) return x;
2455 [ + - ]: 2142 : if (n0 == 0) return 1;
2456 : : }
2457 [ + + ]: 11566504 : if (x <= 1) return x; /* 0 or 1 */
2458 [ + + ]: 11493344 : if (!SMALL_ULONG(p))
2459 : 1203689 : return Fl_powu_pre(x, n0, p, get_Fl_red(p));
2460 : 10289655 : y = 1; z = x; n = n0;
2461 : : for(;;)
2462 : : {
2463 [ + + ]: 123470484 : if (n&1) y = Fl_mul(y,z,p);
2464 [ + + ]: 123470484 : n>>=1; if (!n) return y;
2465 : 113180829 : z = Fl_sqr(z,p);
2466 : 172848454 : }
2467 : : }
2468 : :
2469 : : static long
2470 : 310007 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D)
2471 : : {
2472 : 310007 : D->N = N;
2473 [ + + ][ + + ]: 310007 : if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
[ + + ]
2474 : : {
2475 : 11262 : D->mul2 = &_muli2red;
2476 : 11262 : D->res = &_remiibar;
2477 : 11262 : D->s = 1+(expi(N)>>1);
2478 : 11262 : D->iM = Fp_invmBarrett(N, D->s);
2479 : 11262 : return 0;
2480 : : }
2481 [ + + ][ + + ]: 298745 : else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
2482 : : {
2483 : 247659 : *y = remii(shifti(*y, bit_accuracy(lN)), N);
2484 : 247659 : D->mul2 = &_muli2montred;
2485 : 247659 : D->res = &_montred;
2486 : 247659 : D->inv = init_montdata(N);
2487 : 247659 : return 1;
2488 : : }
2489 : : else
2490 : : {
2491 : 51086 : D->mul2 = &_muli2red;
2492 : 51086 : D->res = &_remii;
2493 : 310007 : return 0;
2494 : : }
2495 : : }
2496 : :
2497 : : GEN
2498 : 6002555 : Fp_powu(GEN A, ulong k, GEN N)
2499 : : {
2500 : 6002555 : long lN = lgefint(N), sA;
2501 : : int base_is_2, use_montgomery;
2502 : : muldata D;
2503 : : pari_sp av;
2504 : :
2505 [ + + ]: 6002555 : if (lN == 3) {
2506 : 3817597 : ulong n = (ulong)N[2];
2507 : 3817597 : return utoi( Fl_powu(umodiu(A, n), k, n) );
2508 : : }
2509 [ + + ]: 2184958 : if (k <= 2)
2510 : : { /* frequent special cases */
2511 [ + + ]: 1944748 : if (k == 2) return Fp_sqr(A,N);
2512 [ + - ]: 85144 : if (k == 1) return A;
2513 [ # # ]: 0 : if (k == 0) return gen_1;
2514 : : }
2515 [ + + ][ - + ]: 240210 : sA = signe(A)==-1 && odd(k);
2516 : 240210 : base_is_2 = 0;
2517 [ + + ]: 240210 : if (lgefint(A) == 3) switch(A[2])
[ + + + ]
2518 : : {
2519 [ - + ]: 2402 : case 1: return sA ? gen_m1 : gen_1;
2520 : 36719 : case 2: base_is_2 = 1; break;
2521 : : }
2522 : :
2523 : : /* TODO: Move this out of here and use for general modular computations */
2524 : 237808 : av = avma;
2525 : 237808 : use_montgomery = Fp_select_red(&A, k, N, lN, &D);
2526 [ + + ]: 237808 : if (base_is_2)
2527 : 36719 : A = gen_powu_fold_i(A, k, (void*)&D, &_sqr, &_m2sqr);
2528 : : else
2529 : 201089 : A = gen_powu_i(A, k, (void*)&D, &_sqr, &_mul);
2530 [ + + ]: 237808 : if (use_montgomery)
2531 : : {
2532 : 194703 : A = _montred(&D, A);
2533 [ - + ]: 194703 : if (cmpii(A,N) >= 0) A = subii(A,N);
2534 [ - + ]: 194703 : if (sA) A = subii(N, A);
2535 : : }
2536 : 6002555 : return gerepileuptoint(av, A);
2537 : : }
2538 : :
2539 : : GEN
2540 : 14336 : Fp_pows(GEN A, long k, GEN N)
2541 : : {
2542 [ + + ]: 14336 : if (lgefint(N) == 3) {
2543 : 309 : ulong n = N[2];
2544 : 309 : ulong a = umodiu(A, n);
2545 [ - + ]: 309 : if (k < 0) {
2546 : 0 : a = Fl_inv(a, n);
2547 : 0 : k = -k;
2548 : : }
2549 : 309 : return utoi( Fl_powu(a, (ulong)k, n) );
2550 : : }
2551 [ - + ]: 14027 : if (k < 0) { A = Fp_inv(A, N); k = -k; };
2552 : 14336 : return Fp_powu(A, (ulong)k, N);
2553 : : }
2554 : :
2555 : : /* A^K mod N */
2556 : : GEN
2557 : 2291644 : Fp_pow(GEN A, GEN K, GEN N)
2558 : : {
2559 : 2291644 : pari_sp av = avma;
2560 : 2291644 : long t,s, lN = lgefint(N), sA;
2561 : : int base_is_2, use_montgomery;
2562 : : GEN y;
2563 : : muldata D;
2564 : :
2565 : 2291644 : s = signe(K);
2566 [ + + ]: 2291644 : if (!s)
2567 : : {
2568 : 5674 : t = signe(remii(A,N)); avma = av;
2569 [ + - ]: 5674 : return t? gen_1: gen_0;
2570 : : }
2571 [ + + ]: 2285970 : if (lN == 3)
2572 : : {
2573 : 2109008 : ulong k, n = N[2], a = umodiu(A, n);
2574 [ + + ]: 2109008 : if (s < 0) a = Fl_inv(a, n);
2575 [ + + ]: 2108952 : if (a <= 1) return utoi(a); /* 0 or 1 */
2576 [ + + ]: 2010772 : if (lgefint(K) > 3)
2577 : : { /* silly case : huge exponent, small modulus */
2578 : 7 : pari_warn(warner, "Mod(a,b)^n with n >> b : wasteful");
2579 [ + - ]: 7 : if (s > 0)
2580 : : {
2581 : 7 : ulong d = ugcd(a, n);
2582 [ + - ]: 7 : if (d != 1)
2583 : : { /* write n = n1 n2, with n2 maximal such that (n1,a) = 1 */
2584 : 7 : ulong n1 = ucoprime_part(n, d), n2 = n/n1;
2585 : :
2586 : 7 : k = umodiu(K, eulerphiu(n1));
2587 : : /* CRT: = a^K (mod n1), = 0 (mod n2)*/
2588 : 7 : return utoi( Fl_mul(Fl_powu(a, k, n1), n2 * Fl_inv(n2,n1), n) );
2589 : : }
2590 : : }
2591 : : /* gcd(a,n) = 1 */
2592 : 0 : k = umodiu(K, eulerphiu(n));
2593 : : }
2594 : : else
2595 : 2010765 : k = (ulong)K[2];
2596 : 2010765 : return utoi(Fl_powu(a, k, n));
2597 : : }
2598 : :
2599 [ + + ]: 176962 : if (s < 0) y = Fp_inv(A,N);
2600 : : else
2601 : : {
2602 : 176630 : y = modii(A,N);
2603 [ - + ]: 176630 : if (!signe(y)) { avma = av; return gen_0; }
2604 : : }
2605 [ + + ]: 176962 : if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
2606 : :
2607 : 72240 : base_is_2 = 0;
2608 [ - + ][ # # ]: 72240 : sA = signe(y)==-1 && mod2(K);
2609 [ + + ]: 72240 : if (lgefint(y) == 3) switch(y[2])
[ + + + ]
2610 : : {
2611 [ - + ]: 41 : case 1: return sA ? gen_m1 : gen_1;
2612 : 53178 : case 2: base_is_2 = 1; break;
2613 : : }
2614 : :
2615 : : /* TODO: Move this out of here and use for general modular computations */
2616 : 72199 : use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D);
2617 [ + + ]: 72199 : if (base_is_2)
2618 : 53178 : y = gen_pow_fold_i(y, K, (void*)&D, &_sqr, &_m2sqr);
2619 : : else
2620 : 19021 : y = gen_pow_i(y, K, (void*)&D, &_sqr, &_mul);
2621 [ + + ]: 72199 : if (use_montgomery)
2622 : : {
2623 : 52956 : y = _montred(&D,y);
2624 [ - + ]: 52956 : if (cmpii(y,N) >= 0) y = subii(y,N);
2625 [ - + ]: 52956 : if (sA) y = subii(N, y);
2626 : : }
2627 : 2291588 : return gerepileuptoint(av,y);
2628 : : }
2629 : :
2630 : : static GEN
2631 : 109512 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
2632 : :
2633 : : static GEN
2634 : 530600 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
2635 : :
2636 : : static GEN
2637 : 330 : _Fp_rand(void *E) { return addis(randomi(subis((GEN)E,1)),1); }
2638 : :
2639 : : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
2640 : :
2641 : : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
2642 : : equalii,equali1,Fp_easylog};
2643 : :
2644 : : static GEN
2645 : 8700992 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
2646 : :
2647 : : static GEN
2648 : 36175801 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
2649 : :
2650 : : static GEN
2651 : 718811 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
2652 : :
2653 : : static GEN
2654 : 37742109 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
2655 : :
2656 : : static GEN
2657 : 123153 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
2658 : :
2659 : : static int
2660 : 2170862 : _Fp_equal0(GEN x) { return signe(x)==0; }
2661 : :
2662 : : static GEN
2663 : 669979 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
2664 : :
2665 : : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
2666 : : _Fp_inv,_Fp_equal0,_Fp_s};
2667 : :
2668 : 12934 : const struct bb_field *get_Fp_field(void **E, GEN p)
2669 : : {
2670 : 12934 : *E = (void*)p; return &Fp_field;
2671 : : }
2672 : :
2673 : : /*********************************************************************/
2674 : : /** **/
2675 : : /** ORDER of INTEGERMOD x in (Z/nZ)* **/
2676 : : /** **/
2677 : : /*********************************************************************/
2678 : : ulong
2679 : 2947 : Fl_order(ulong a, ulong o, ulong p)
2680 : : {
2681 : 2947 : pari_sp av = avma;
2682 : : GEN m, P, E;
2683 : : long i;
2684 [ - + ]: 2947 : if (!o) o = p-1;
2685 : 2947 : m = factoru(o);
2686 : 2947 : P = gel(m,1);
2687 : 2947 : E = gel(m,2);
2688 [ + + ]: 8085 : for (i = lg(P)-1; i; i--)
2689 : : {
2690 : 5138 : ulong j, l=P[i], e=E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
2691 [ + + ]: 5138 : if (y == 1) o = t;
2692 : : else {
2693 [ + + ][ + + ]: 5292 : for (j = 1; j < e; j++) { y = Fl_powu(y, l, p); if (y == 1) break; }
2694 : 3829 : o = t * upowuu(l, j);
2695 : : }
2696 : : }
2697 : 2947 : avma = av; return o;
2698 : : }
2699 : :
2700 : : /*Find the exact order of a assuming a^o==1*/
2701 : : GEN
2702 : 6151 : Fp_order(GEN a, GEN o, GEN p) {
2703 [ + + ][ + + ]: 6151 : if (lgefint(p) == 3 && typ(o) == t_INT && lgefint(o)==3)
[ + - ]
2704 : : {
2705 : 28 : ulong pp = p[2], oo = o[2];
2706 : 28 : return utoi( Fl_order(umodiu(a, pp), oo, pp) );
2707 : : }
2708 : 6151 : return gen_order(a, o, (void*)p, &Fp_star);
2709 : : }
2710 : : GEN
2711 : 49 : Fp_factored_order(GEN a, GEN o, GEN p)
2712 : 49 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
2713 : :
2714 : : /* return order of a mod p^e, e > 0, pe = p^e */
2715 : : static GEN
2716 : 70 : Zp_order(GEN a, GEN p, long e, GEN pe)
2717 : : {
2718 : : GEN ap, op;
2719 [ + + ]: 70 : if (equaliu(p, 2))
2720 : : {
2721 [ - + ]: 49 : if (e == 1) return gen_1;
2722 [ - + ][ # # ]: 49 : if (e == 2) return mod4(a) == 1? gen_1: gen_2;
2723 [ + + ]: 49 : if (mod4(a) == 1)
2724 : 14 : op = gen_1;
2725 : : else {
2726 : 35 : op = gen_2;
2727 : 35 : a = Fp_sqr(a, pe);
2728 : : }
2729 : : } else {
2730 [ - + ]: 21 : ap = (e == 1)? a: remii(a,p);
2731 : 21 : op = Fp_order(ap, subis(p,1), p);
2732 [ + - ]: 21 : if (e == 1) return op;
2733 : 0 : a = Fp_pow(a, op, pe); /* 1 mod p */
2734 : : }
2735 [ + + ]: 49 : if (equali1(a)) return op;
2736 : 70 : return mulii(op, powiu(p, e - Z_pval(subis(a,1), p)));
2737 : : }
2738 : :
2739 : : GEN
2740 : 70 : znorder(GEN x, GEN o)
2741 : : {
2742 : 70 : pari_sp av = avma;
2743 : : GEN b, a;
2744 : :
2745 [ + + ]: 70 : if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
2746 : 63 : b = gel(x,1); a = gel(x,2);
2747 [ + + ]: 63 : if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
2748 [ + + ]: 56 : if (!o)
2749 : : {
2750 : 42 : GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
2751 : 42 : long i, l = lg(P);
2752 : 42 : o = gen_1;
2753 [ + + ]: 84 : for (i = 1; i < l; i++)
2754 : : {
2755 : 42 : GEN p = gel(P,i);
2756 : 42 : long e = itos(gel(E,i));
2757 : :
2758 [ + - ]: 42 : if (l == 2)
2759 : 42 : o = Zp_order(a, p, e, b);
2760 : : else {
2761 : 0 : GEN pe = powiu(p,e);
2762 : 0 : o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
2763 : : }
2764 : : }
2765 : 42 : return gerepileuptoint(av, o);
2766 : : }
2767 : 56 : return Fp_order(a, o, b);
2768 : : }
2769 : : GEN
2770 : 7 : order(GEN x) { return znorder(x, NULL); }
2771 : :
2772 : : /*********************************************************************/
2773 : : /** **/
2774 : : /** DISCRETE LOGARITHM in (Z/nZ)* **/
2775 : : /** **/
2776 : : /*********************************************************************/
2777 : : static GEN
2778 : 55599 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
2779 : : {
2780 : 55599 : pari_sp av = avma;
2781 : : GEN h1, h2, F, G;
2782 [ + + ]: 55599 : if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return NULL;
2783 [ + + ][ + + ]: 33440 : if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
2784 : : {
2785 : 84 : GEN M = cgetg(3, t_MAT);
2786 : 84 : gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
2787 : 84 : gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
2788 : 84 : return gerepileupto(av, M);
2789 : : }
2790 : 55599 : avma = av; return NULL;
2791 : : }
2792 : :
2793 : : static GEN
2794 : 55599 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
2795 : : {
2796 : : GEN rel;
2797 : : do
2798 : : {
2799 : 55599 : (*e)++; *g = Fp_mul(*g, b, p);
2800 : 55599 : rel = Fp_log_halfgcd(bnd, C, *g, p);
2801 [ + + ]: 55599 : } while (!rel);
2802 : 84 : return rel;
2803 : : }
2804 : :
2805 : : struct Fp_log_rel
2806 : : {
2807 : : GEN rel;
2808 : : long *sieve;
2809 : : ulong prmax;
2810 : : long nbrel, nbmax;
2811 : : };
2812 : :
2813 : : /* add u^e */
2814 : : static long
2815 : 22197 : addifsmooth1(struct Fp_log_rel *r, GEN h, long u, long e)
2816 : : {
2817 : 22197 : pari_sp av = avma;
2818 : : GEN z;
2819 [ + + ]: 22197 : if ((z = Z_issmooth_fact(h, r->prmax)))
2820 : : {
2821 : 924 : long off = r->prmax+1;
2822 : 924 : GEN F = cgetg(3, t_MAT);
2823 : 924 : gel(F,1) = vecsmall_append(gel(z,1), off+u);
2824 : 924 : gel(F,2) = vecsmall_append(gel(z,2), e);
2825 : 924 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
2826 : : }
2827 : 22197 : return r->nbrel==r->nbmax;
2828 : : }
2829 : :
2830 : : /* add u^-1 v^-1 */
2831 : : static long
2832 : 153230 : addifsmooth2(struct Fp_log_rel *r, GEN h, long u, long v)
2833 : : {
2834 : 153230 : pari_sp av = avma;
2835 : : GEN z;
2836 [ + + ]: 153230 : if ((z = Z_issmooth_fact(h, r->prmax)))
2837 : : {
2838 : 56308 : long off = r->prmax+1;
2839 : 56308 : GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
2840 : 56308 : GEN F = cgetg(3, t_MAT);
2841 : 56308 : gel(F,1) = vecsmall_concat(gel(z,1), P);
2842 : 56308 : gel(F,2) = vecsmall_concat(gel(z,2), E);
2843 : 56308 : gel(r->rel,++r->nbrel) = gerepileupto(av, F);
2844 : : }
2845 : 153230 : return r->nbrel==r->nbmax;
2846 : : }
2847 : :
2848 : : /*
2849 : : Let p=C^2+c
2850 : : Solve h = (C+x)*(C+a)-p = 0 [mod l]
2851 : : h= -c+x*(C+a)+C*a = 0 [mod l]
2852 : : x = (c-C*a)/(C+a) [mod l]
2853 : : h = -c+C*(x+a)+a*x
2854 : : */
2855 : :
2856 : : static void
2857 : 21994 : Fp_log_sieve_h(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, long a, GEN pr, GEN sz)
2858 : : {
2859 : 21994 : long th = expi(C), n = lg(pr)-1;
2860 : : long i,j;
2861 [ - + ]: 21994 : if (addifsmooth1(r, addis(C,a), a, -1)) return;
2862 [ + + ]: 19854541 : for(j=0; j<=a; j++)
2863 : 19832547 : r->sieve[j]=0;
2864 [ + + ]: 12839582 : for(i=1; i<=n; i++)
2865 : : {
2866 : 12817588 : ulong li = pr[i], s = sz[i], al = a % li;
2867 : 12817588 : ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
2868 [ + + ]: 12817588 : if (!iv) continue;
2869 : 12765361 : u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
2870 [ + + ]: 51272431 : for(j = u; j<=a; j+=li)
2871 : 38507070 : r->sieve[j] += s;
2872 : : }
2873 : 21994 : th = th - expu(th)-1;
2874 [ + + ]: 19814347 : for(j=0; j<a; j++)
2875 [ + + ]: 19792367 : if (r->sieve[j]>=th)
2876 : : {
2877 : 153230 : GEN h = addiu(subii(muliu(C,a+j),c), a*j);
2878 [ + + ]: 153230 : if (addifsmooth2(r, h, a, j)) return;
2879 : : }
2880 : : /* j = a */
2881 [ + + ]: 21980 : if (r->sieve[a]>=th)
2882 : : {
2883 : 203 : GEN h = addiu(subii(muliu(C,2*a),c), a*a);
2884 [ - + ]: 21994 : if (addifsmooth1(r, h, a, -2)) return;
2885 : : }
2886 : : }
2887 : :
2888 : : static GEN
2889 : 511 : _psi(void*E, GEN y)
2890 : : {
2891 : 511 : GEN lx = (GEN) E;
2892 : 511 : long prec = lg(lx);
2893 : 511 : GEN ly = glog(y, prec);
2894 : 511 : GEN u = gdiv(lx, ly);
2895 : 511 : return gsub(gdiv(y ,ly), gpow(u, u, prec));
2896 : : }
2897 : :
2898 : : static GEN
2899 : 14 : opt_param(GEN x, long prec)
2900 : : {
2901 : 14 : return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
2902 : : }
2903 : :
2904 : : static GEN
2905 : 14 : check_kernel(long N, long prmax, GEN C, GEN M, GEN p, GEN m)
2906 : : {
2907 : 14 : pari_sp av = avma;
2908 : 14 : GEN K = FpMs_leftkernel_elt(M, N, m);
2909 : 14 : long i, f=0;
2910 : 14 : long l = lg(K), lm = lgefint(m);
2911 : 14 : GEN idx = diviiexact(subis(p,1),m), g;
2912 : : pari_timer ti;
2913 [ - + ]: 14 : if (DEBUGLEVEL) timer_start(&ti);
2914 [ + - ]: 28 : for(i=1; i<l; i++)
2915 [ + + ]: 28 : if (signe(gel(K,i)))
2916 : 14 : break;
2917 : 14 : g = utoi(i);
2918 : 14 : K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
2919 [ + + ]: 73836 : for(i=1; i<l; i++)
2920 : : {
2921 : 73822 : GEN k = gel(K,i);
2922 [ + + ]: 73822 : GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
2923 [ + + ][ - + ]: 73822 : if (signe(k)==0 || !equalii(Fp_pow(g, mulii(k,idx), p),
2924 : : Fp_pow(j, idx, p)))
2925 : 44814 : gel(K,i) = cgetineg(lm);
2926 : : else
2927 : 29008 : f++;
2928 : : }
2929 [ - + ]: 14 : if (DEBUGLEVEL) timer_printf(&ti,"found %ld logs", f);
2930 : 14 : return gerepileupto(av, K);
2931 : : }
2932 : :
2933 : : static GEN
2934 : 28 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
2935 : : {
2936 : 28 : pari_sp av=avma;
2937 : 28 : GEN aa = gen_1;
2938 : 28 : long AV = 0;
2939 : : for(;;)
2940 : : {
2941 : 84 : GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
2942 : 84 : GEN F = gel(A,1), E = gel(A,2);
2943 : 84 : GEN Ao = gen_0;
2944 : 84 : long i, l = lg(F);
2945 [ + + ]: 418 : for(i=1; i<l; i++)
2946 : : {
2947 : 390 : GEN Ki = gel(K,F[i]);
2948 [ + + ]: 390 : if (signe(Ki)<0) break;
2949 : 334 : Ao = addii(Ao, mulis(Ki, E[i]));
2950 : : }
2951 [ + + ]: 84 : if (i==l) return Fp_div(Ao, utoi(AV), m);
2952 : 56 : aa = gerepileuptoint(av, aa);
2953 : 56 : }
2954 : : }
2955 : :
2956 : : static GEN
2957 : 14 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
2958 : : {
2959 : 14 : pari_sp av = avma, av2;
2960 : : long i, nbi, nbrow;
2961 : : GEN C, c, Ci, ci, pr, sz, l, Ao, Bo, K, d, p_1;
2962 : : pari_timer ti;
2963 : : struct Fp_log_rel r;
2964 : 14 : ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
2965 : 14 : ulong bnd = 4*bnds;
2966 [ + - ][ - + ]: 14 : if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
2967 : :
2968 : 14 : p_1 = subiu(p,1);
2969 [ - + ]: 14 : if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
2970 : 0 : m = diviiexact(p_1, coprime_part(p_1, m));
2971 : 14 : pr = primes_upto_zv(bnd);
2972 : 14 : nbi = lg(pr)-1;
2973 [ - + ]: 14 : if (DEBUGLEVEL)
2974 : : {
2975 : 0 : err_printf("bnd=%lu Size FB=%ld\n", bnd, nbi);
2976 : 0 : timer_start(&ti);
2977 : : }
2978 : 14 : C = sqrtremi(p, &c);
2979 : 14 : av2 = avma;
2980 : 14 : Ci = cgetg(nbi+1,t_VECSMALL);
2981 : 14 : ci = cgetg(nbi+1,t_VECSMALL);
2982 : 14 : sz = cgetg(nbi+1,t_VECSMALL);
2983 [ + + ]: 7168 : for (i = 1; i <= nbi; ++i)
2984 : : {
2985 : 7154 : ulong lp = pr[i];
2986 : 7154 : Ci[i] = umodiu(C, lp);
2987 : 7154 : ci[i] = umodiu(c, lp);
2988 : 7154 : sz[i] = expu(lp);
2989 : : }
2990 : 14 : r.nbrel = 0;
2991 : 14 : r.nbmax = 8*nbi;
2992 : 14 : r.rel = cgetg(r.nbmax+1,t_VEC);
2993 : 14 : r.sieve = cgetg(r.nbmax+2,t_VECSMALL)+1;
2994 : 14 : r.prmax = pr[nbi];
2995 [ + + ]: 22008 : for(i=0; r.nbrel < r.nbmax; i++)
2996 : : {
2997 : 21994 : Fp_log_sieve_h(&r, C, c, Ci, ci, i, pr, sz);
2998 [ - + ][ # # ]: 21994 : if (DEBUGLEVEL && (i&127)==0)
2999 : 0 : err_printf("%ld%% ",100*r.nbrel/(r.nbmax));
3000 : : }
3001 : 14 : nbrow = r.prmax+i;
3002 [ - + ]: 14 : if (DEBUGLEVEL)
3003 : : {
3004 : 0 : err_printf("\n");
3005 : 0 : timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+i);
3006 : : }
3007 : 14 : setlg(r.rel,r.nbrel+1);
3008 : 14 : r.rel = gerepileupto(av2, r.rel);
3009 : 14 : K = check_kernel(nbrow, r.prmax, C, r.rel, p, m);
3010 [ - + ]: 14 : if (DEBUGLEVEL) timer_start(&ti);
3011 : 14 : Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
3012 [ - + ]: 14 : if (DEBUGLEVEL) timer_printf(&ti," log element");
3013 : 14 : Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
3014 [ - + ]: 14 : if (DEBUGLEVEL) timer_printf(&ti," log generator");
3015 : 14 : d = gcdii(Ao,Bo);
3016 : 14 : l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
3017 [ - + ]: 14 : if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
3018 : 14 : return gerepileuptoint(av, l);
3019 : : }
3020 : :
3021 : : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
3022 : : static GEN
3023 : 213205 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
3024 : : {
3025 : 213205 : pari_sp av = avma;
3026 : 213205 : GEN p = (GEN)E;
3027 : : /* assume a reduced mod p, p not necessarily prime */
3028 [ + + ]: 213205 : if (equali1(a)) return gen_0;
3029 : : /* p > 2 */
3030 [ + + ]: 148854 : if (equalii(subis(p,1), a)) /* -1 */
3031 : : {
3032 : : pari_sp av2;
3033 : : GEN t;
3034 : 72682 : ord = dlog_get_ord(ord);
3035 [ + + ]: 72682 : if (mpodd(ord)) { avma = av; return cgetg(1, t_VEC); } /* no solution */
3036 : 72668 : t = shifti(ord,-1); /* only possible solution */
3037 : 72668 : av2 = avma;
3038 [ + + ]: 72668 : if (!equalii(Fp_pow(g, t, p), a)) { avma = av; return cgetg(1, t_VEC); }
3039 : 72640 : avma = av2; return gerepileuptoint(av, t);
3040 : : }
3041 [ + + ][ + + ]: 76172 : if (typ(ord)==t_INT && expi(ord)>=27 && BPSW_psp(p))
[ + + ]
3042 : 14 : return Fp_log_index(a, g, ord, p);
3043 : 213205 : avma = av; return NULL; /* not easy */
3044 : : }
3045 : :
3046 : : GEN
3047 : 238709 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
3048 : : {
3049 : 238709 : GEN v = dlog_get_ordfa(ord);
3050 : 238681 : ord = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
3051 : 238681 : return gen_PH_log(a,g,ord,(void*)p,&Fp_star);
3052 : : }
3053 : :
3054 : : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
3055 : : * PHI[l] = eulerphi(N / P[l]^E[l]). Destroys P/E */
3056 : : static GEN
3057 : 91 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
3058 : : {
3059 : 91 : long l = lg(P) - 1, e = E[l];
3060 [ + + ]: 91 : GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
3061 : : GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
3062 : :
3063 [ + + ]: 91 : if (l == 1) {
3064 : 63 : hpe = h;
3065 : 63 : gpe = g;
3066 : : } else {
3067 : 28 : hpe = modii(h, pe);
3068 : 28 : gpe = modii(g, pe);
3069 : : }
3070 [ + + ]: 91 : if (e == 1) {
3071 : 21 : hp = hpe;
3072 : 21 : gp = gpe;
3073 : : } else {
3074 : 70 : hp = remii(hpe, p);
3075 : 70 : gp = remii(gpe, p);
3076 : : }
3077 [ + + ][ - + ]: 91 : if (hp == gen_0 || gp == gen_0) return NULL;
3078 [ + + ]: 77 : if (equaliu(p, 2))
3079 : : {
3080 : 28 : GEN N = int2n(e);
3081 : 28 : ogpe = Zp_order(gpe, gen_2, e, N);
3082 : 28 : a = Fp_log(hpe, gpe, ogpe, N);
3083 [ + + ]: 28 : if (typ(a) != t_INT) return NULL;
3084 : : }
3085 : : else
3086 : : { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
3087 : : is trivial */
3088 : : /* [order(gp), factor(order(gp))] */
3089 : 49 : GEN v = Fp_factored_order(gp, subis(p,1), p);
3090 : 49 : GEN ogp = gel(v,1);
3091 [ - + ]: 49 : if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
3092 : 49 : a = Fp_log(hp, gp, v, p);
3093 [ - + ]: 49 : if (typ(a) != t_INT) return NULL;
3094 [ + + ]: 49 : if (e == 1) ogpe = ogp;
3095 : : else
3096 : : { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
3097 : : /* use p-adic log: O(log p + e) mul*/
3098 : : long vpogpe, vpohpe;
3099 : :
3100 : 28 : hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
3101 : 28 : gpe = Fp_pow(gpe, ogp, pe);
3102 : : /* g,h = 1 mod p; compute b s.t. h = g^b */
3103 : :
3104 : : /* v_p(order g mod pe) */
3105 [ + - ]: 28 : vpogpe = equali1(gpe)? 0: e - Z_pval(subis(gpe,1), p);
3106 : : /* v_p(order h mod pe) */
3107 [ + - ]: 28 : vpohpe = equali1(gpe)? 0: e - Z_pval(subis(hpe,1), p);
3108 [ - + ]: 28 : if (vpohpe > vpogpe) return NULL;
3109 : :
3110 : 28 : ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
3111 [ - + ][ # # ]: 28 : if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
3112 : 28 : b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
3113 : 28 : a = addii(a, mulii(ogp, gtrunc(b)));
3114 : : }
3115 : : }
3116 : : /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
3117 [ + + ]: 63 : if (l == 1) return a;
3118 : :
3119 : 28 : N = diviiexact(N, pe); /* make N coprime to p */
3120 : 28 : h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
3121 : 28 : g = Fp_pow(g, modii(ogpe, phi), N);
3122 : 28 : setlg(P, l); /* remove last element */
3123 : 28 : setlg(E, l);
3124 : 28 : b = znlog_rec(h, g, N, P, E, PHI);
3125 [ - + ]: 28 : if (!b) return NULL;
3126 : 91 : return addmulii(a, b, ogpe);
3127 : : }
3128 : :
3129 : : static GEN
3130 : 63 : get_PHI(GEN P, GEN E)
3131 : : {
3132 : 63 : long i, l = lg(P);
3133 : 63 : GEN PHI = cgetg(l, t_VEC);
3134 : 63 : gel(PHI,1) = gen_1;
3135 [ + + ]: 91 : for (i=1; i<l-1; i++)
3136 : : {
3137 : 28 : GEN t, p = gel(P,i);
3138 : 28 : long e = E[i];
3139 : 28 : t = mulii(powiu(p, e-1), subis(p,1));
3140 [ + + ]: 28 : if (i > 1) t = mulii(t, gel(PHI,i));
3141 : 28 : gel(PHI,i+1) = t;
3142 : : }
3143 : 63 : return PHI;
3144 : : }
3145 : :
3146 : : GEN
3147 : 168 : znlog(GEN h, GEN g, GEN o)
3148 : : {
3149 : 168 : pari_sp av = avma;
3150 : : GEN N, fa, P, E, x;
3151 [ + + - ]: 168 : switch (typ(g))
3152 : : {
3153 : : case t_PADIC:
3154 : : {
3155 : 28 : GEN p = gel(g,2);
3156 : 28 : long v = valp(g);
3157 [ - + ]: 28 : if (v < 0) pari_err_DIM("znlog");
3158 [ - + ]: 28 : if (v > 0) {
3159 : 0 : long k = gvaluation(h, p);
3160 [ # # ]: 0 : if (k % v) return cgetg(1,t_VEC);
3161 : 0 : k /= v;
3162 [ # # ]: 0 : if (!gequal(h, gpowgs(g,k))) { avma = av; return cgetg(1,t_VEC); }
3163 : 0 : avma = av; return stoi(k);
3164 : : }
3165 : 28 : N = gel(g,3);
3166 : 28 : g = Rg_to_Fp(g, N);
3167 : 28 : break;
3168 : : }
3169 : : case t_INTMOD:
3170 : 140 : N = gel(g,1);
3171 : 140 : g = gel(g,2); break;
3172 : 0 : default: pari_err_TYPE("znlog", g);
3173 : 0 : return NULL; /* not reached */
3174 : : }
3175 [ - + ]: 168 : if (equali1(N)) { avma = av; return gen_0; }
3176 : 168 : h = Rg_to_Fp(h, N);
3177 [ + + ]: 168 : if (o) return gerepileupto(av, Fp_log(h, g, o, N));
3178 : 63 : fa = Z_factor(N);
3179 : 63 : P = gel(fa,1);
3180 : 63 : E = vec_to_vecsmall(gel(fa,2));
3181 : 63 : x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
3182 [ + + ]: 63 : if (!x) { avma = av; return cgetg(1,t_VEC); }
3183 : 140 : return gerepileuptoint(av, x);
3184 : : }
3185 : :
3186 : : GEN
3187 : 1306 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
3188 : : {
3189 : 1306 : a = modii(a,p);
3190 [ + + ]: 1306 : if (!signe(a))
3191 : : {
3192 [ - + ]: 7 : if (zeta) *zeta = gen_1;
3193 [ + - ]: 7 : if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
3194 : 0 : return gen_0;
3195 : : }
3196 [ + + ]: 1299 : if (equaliu(n,2))
3197 : : {
3198 [ - + ]: 480 : if (zeta) *zeta = addis(p,-1);
3199 : 480 : return Fp_sqrt(a,p);
3200 : : }
3201 : 1299 : return gen_Shanks_sqrtn(a,n,addis(p,-1),zeta,(void*)p,&Fp_star);
3202 : : }
3203 : :
3204 : : /*********************************************************************/
3205 : : /** **/
3206 : : /** FUNDAMENTAL DISCRIMINANTS **/
3207 : : /** **/
3208 : : /*********************************************************************/
3209 : : long
3210 : 14028 : isfundamental(GEN x) {
3211 [ - + ]: 14028 : if (typ(x) != t_INT) pari_err_TYPE("isfundamental",x);
3212 : 14028 : return Z_isfundamental(x);
3213 : : }
3214 : :
3215 : : /* x fundamental ? */
3216 : : long
3217 : 6153 : uposisfundamental(ulong x)
3218 : : {
3219 : 6153 : ulong r = x & 15; /* x mod 16 */
3220 [ + + ]: 6153 : if (!r) return 0;
3221 [ + + + ]: 5775 : switch(r & 3)
3222 : : { /* x mod 4 */
3223 [ + + ]: 1163 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
3224 : 1612 : case 1: return uissquarefree(x);
3225 : 6153 : default: return 0;
3226 : : }
3227 : : }
3228 : : /* -x fundamental ? */
3229 : : long
3230 : 9387 : unegisfundamental(ulong x)
3231 : : {
3232 : 9387 : ulong r = x & 15; /* x mod 16 */
3233 [ + + ]: 9387 : if (!r) return 0;
3234 [ + + + ]: 8904 : switch(r & 3)
3235 : : { /* x mod 4 */
3236 [ + + ]: 1674 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
3237 : 4230 : case 3: return uissquarefree(x);
3238 : 9387 : default: return 0;
3239 : : }
3240 : : }
3241 : : long
3242 : 14581 : Z_isfundamental(GEN x)
3243 : : {
3244 : : long r;
3245 [ - + + ]: 14581 : switch(lgefint(x))
3246 : : {
3247 : 0 : case 2: return 0;
3248 : 12579 : case 3: return signe(x) < 0? unegisfundamental(x[2])
3249 [ + + ]: 12579 : : uposisfundamental(x[2]);
3250 : : }
3251 : 2002 : r = mod16(x);
3252 [ + + ]: 2002 : if (!r) return 0;
3253 [ + + ]: 1876 : if ((r & 3) == 0)
3254 : : {
3255 : : pari_sp av;
3256 : 376 : r >>= 2; /* |x|/4 mod 4 */
3257 [ + + ]: 376 : if (signe(x) < 0) r = 4-r;
3258 [ + + ]: 376 : if (r == 1) return 0;
3259 : 250 : av = avma;
3260 : 250 : r = Z_issquarefree( shifti(x,-2) );
3261 : 250 : avma = av; return r;
3262 : : }
3263 : 1500 : r &= 3; /* |x| mod 4 */
3264 [ + + ]: 1500 : if (signe(x) < 0) r = 4-r;
3265 [ + + ]: 14581 : return (r==1) ? Z_issquarefree(x) : 0;
3266 : : }
3267 : :
3268 : : GEN
3269 : 14 : quaddisc(GEN x)
3270 : : {
3271 : 14 : const pari_sp av = avma;
3272 : 14 : long i,r,tx=typ(x);
3273 : : GEN P,E,f,s;
3274 : :
3275 [ - + ]: 14 : if (!is_rational_t(tx)) pari_err_TYPE("quaddisc",x);
3276 : 14 : f = factor(x);
3277 : 14 : P = gel(f,1);
3278 : 14 : E = gel(f,2); s = gen_1;
3279 [ + + ]: 70 : for (i=1; i<lg(P); i++)
3280 [ + + ]: 56 : if (odd(mael(E,i,2))) s = mulii(s,gel(P,i));
3281 [ + - ]: 14 : r = mod4(s); if (gsigne(x) < 0) r = 4-r;
3282 [ - + ]: 14 : if (r>1) s = shifti(s,2);
3283 : 14 : return gerepileuptoint(av, s);
3284 : : }
3285 : :
3286 : : /*********************************************************************/
3287 : : /** **/
3288 : : /** FACTORIAL **/
3289 : : /** **/
3290 : : /*********************************************************************/
3291 : : /* return a * (a+1) * ... * b. Assume a <= b [ note: factoring out powers of 2
3292 : : * first is slower ... ] */
3293 : : GEN
3294 : 99396 : mulu_interval(ulong a, ulong b)
3295 : : {
3296 : 99396 : pari_sp av = avma;
3297 : 99396 : ulong k, l, N, n = b - a + 1;
3298 : : long lx;
3299 : : GEN x;
3300 : :
3301 [ + + ]: 99396 : if (n < 61)
3302 : : {
3303 : 94749 : x = utoi(a);
3304 [ + + ]: 975675 : for (k=a+1; k<=b; k++) x = mului(k,x);
3305 : 94749 : return gerepileuptoint(av, x);
3306 : : }
3307 : 4647 : lx = 1; x = cgetg(2 + n/2, t_VEC);
3308 : 4647 : N = b + a;
3309 : 4647 : for (k = a;; k++)
3310 : : {
3311 [ + + ]: 670919 : l = N - k; if (l <= k) break;
3312 : 666272 : gel(x,lx++) = muluu(k,l);
3313 : 666272 : }
3314 [ + + ]: 4647 : if (l == k) gel(x,lx++) = utoipos(k);
3315 : 4647 : setlg(x, lx);
3316 : 99396 : return gerepileuptoint(av, divide_conquer_prod(x, mulii));
3317 : : }
3318 : :
3319 : : GEN
3320 : 56423 : mpfact(long n)
3321 : : {
3322 [ + + ]: 56423 : if (n < 2)
3323 : : {
3324 [ - + ]: 1855 : if (n < 0) pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
3325 : 1855 : return gen_1;
3326 : : }
3327 : 56423 : return mulu_interval(2UL, (ulong)n);
3328 : : }
3329 : :
3330 : : /*******************************************************************/
3331 : : /** **/
3332 : : /** LUCAS & FIBONACCI **/
3333 : : /** **/
3334 : : /*******************************************************************/
3335 : : static void
3336 : 112 : lucas(ulong n, GEN *a, GEN *b)
3337 : : {
3338 : : GEN z, t, zt;
3339 [ + + ]: 210 : if (!n) { *a = gen_2; *b = gen_1; return; }
3340 : 98 : lucas(n >> 1, &z, &t); zt = mulii(z, t);
3341 [ + + + + : 98 : switch(n & 3) {
- ]
3342 : 28 : case 0: *a = addsi(-2,sqri(z)); *b = addsi(-1,zt); break;
3343 : 28 : case 1: *a = addsi(-1,zt); *b = addsi(2,sqri(t)); break;
3344 : 14 : case 2: *a = addsi(2,sqri(z)); *b = addsi(1,zt); break;
3345 : 98 : case 3: *a = addsi(1,zt); *b = addsi(-2,sqri(t));
3346 : : }
3347 : : }
3348 : :
3349 : : GEN
3350 : 14 : fibo(long n)
3351 : : {
3352 : 14 : pari_sp av = avma;
3353 : : GEN a, b;
3354 [ - + ]: 14 : if (!n) return gen_0;
3355 : 14 : lucas((ulong)(labs(n)-1), &a, &b);
3356 : 14 : a = diviuexact(addii(shifti(a,1),b), 5);
3357 [ - + ][ # # ]: 14 : if (n < 0 && !odd(n)) setsigne(a, -1);
3358 : 14 : return gerepileuptoint(av, a);
3359 : : }
3360 : :
3361 : : /*******************************************************************/
3362 : : /* */
3363 : : /* CONTINUED FRACTIONS */
3364 : : /* */
3365 : : /*******************************************************************/
3366 : : static GEN
3367 : 203 : icopy_lg(GEN x, long l)
3368 : : {
3369 : 203 : long lx = lgefint(x);
3370 : : GEN y;
3371 : :
3372 [ + + ]: 203 : if (lx >= l) return icopy(x);
3373 : 203 : y = cgeti(l); affii(x, y); return y;
3374 : : }
3375 : :
3376 : : /* continued fraction of a/b. If y != NULL, stop when partial quotients
3377 : : * differ from y */
3378 : : static GEN
3379 : 406 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
3380 : : {
3381 : : GEN z, c;
3382 : 406 : ulong i, l, ly = lgefint(b);
3383 : :
3384 : : /* times log(2) / log2( (1+sqrt(5)) / 2 ) */
3385 : 406 : l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
3386 [ + + ][ + - ]: 406 : if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
[ + - ]
3387 [ - + ]: 406 : if (l > LGBITS) l = LGBITS;
3388 : :
3389 : 406 : z = cgetg(l,t_VEC);
3390 : 406 : l--;
3391 [ + + ]: 406 : if (y) {
3392 : 203 : pari_sp av = avma;
3393 [ + + ]: 203 : if (l >= (ulong)lg(y)) l = lg(y)-1;
3394 [ + + ]: 8732 : for (i = 1; i <= l; i++)
3395 : : {
3396 : 8529 : GEN q = gel(y,i);
3397 : 8529 : gel(z,i) = q;
3398 [ + + ]: 8529 : c = b; if (!gequal1(q)) c = mulii(q, b);
3399 : 8529 : c = subii(a, c);
3400 [ + + ]: 8529 : if (signe(c) < 0)
3401 : : { /* partial quotient too large */
3402 : 50 : c = addii(c, b);
3403 [ + + ]: 50 : if (signe(c) >= 0) i++; /* by 1 */
3404 : 50 : break;
3405 : : }
3406 [ - + ]: 8479 : if (cmpii(c, b) >= 0)
3407 : : { /* partial quotient too small */
3408 : 0 : c = subii(c, b);
3409 [ # # ]: 0 : if (cmpii(c, b) < 0) {
3410 : : /* by 1. If next quotient is 1 in y, add 1 */
3411 [ # # ][ # # ]: 0 : if (i < l && equali1(gel(y,i+1))) gel(z,i) = addis(q,1);
3412 : 0 : i++;
3413 : : }
3414 : 0 : break;
3415 : : }
3416 [ - + ]: 8479 : if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
3417 : 8479 : a = b; b = c;
3418 : : }
3419 : : } else {
3420 : 203 : a = icopy_lg(a, ly);
3421 : 203 : b = icopy(b);
3422 [ + + ]: 9258 : for (i = 1; i <= l; i++)
3423 : : {
3424 : 9069 : gel(z,i) = truedvmdii(a,b,&c);
3425 [ + + ]: 9069 : if (c == gen_0) { i++; break; }
3426 : 9055 : affii(c, a); cgiv(c); c = a;
3427 : 9055 : a = b; b = c;
3428 : : }
3429 : : }
3430 : 406 : i--;
3431 [ + - ][ + + ]: 406 : if (i > 1 && gequal1(gel(z,i)))
3432 : : {
3433 : 36 : cgiv(gel(z,i)); --i;
3434 : 36 : gel(z,i) = addsi(1, gel(z,i)); /* unclean: leave old z[i] on stack */
3435 : : }
3436 : 406 : setlg(z,i+1); return z;
3437 : : }
3438 : :
3439 : : static GEN
3440 : 0 : sersfcont(GEN a, GEN b, long k)
3441 : : {
3442 [ # # ]: 0 : long i, l = typ(a) == t_POL? lg(a): 3;
3443 : : GEN y, c;
3444 [ # # ]: 0 : if (lg(b) > l) l = lg(b);
3445 [ # # ][ # # ]: 0 : if (k > 0 && l > k+1) l = k+1;
3446 : 0 : y = cgetg(l,t_VEC);
3447 [ # # ]: 0 : for (i=1; i<l; i++)
3448 : : {
3449 : 0 : gel(y,i) = poldivrem(a,b,&c);
3450 [ # # ]: 0 : if (gequal0(c)) { i++; break; }
3451 : 0 : a = b; b = c;
3452 : : }
3453 : 0 : setlg(y, i); return y;
3454 : : }
3455 : :
3456 : : GEN
3457 : 217 : gboundcf(GEN x, long k)
3458 : : {
3459 : : pari_sp av;
3460 : 217 : long tx = typ(x), e;
3461 : : GEN y, a, b, c;
3462 : :
3463 [ + + ]: 217 : if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
3464 [ + - ]: 210 : if (is_scalar_t(tx))
3465 : : {
3466 [ - + ]: 210 : if (gequal0(x)) return mkvec(gen_0);
3467 [ - + - - ]: 210 : switch(tx)
3468 : : {
3469 : 0 : case t_INT: return mkveccopy(x);
3470 : : case t_REAL:
3471 : 210 : av = avma;
3472 : 210 : c = mantissa_real(x,&e);
3473 [ + + ]: 210 : if (e < 0) pari_err_PREC("gboundcf");
3474 : 203 : y = int2n(e);
3475 : 203 : a = Qsfcont(c,y, NULL, k);
3476 : 203 : b = addsi(signe(x), c);
3477 : 203 : return gerepilecopy(av, Qsfcont(b,y, a, k));
3478 : :
3479 : : case t_FRAC:
3480 : 0 : av = avma;
3481 : 0 : return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
3482 : : }
3483 : 0 : pari_err_TYPE("gboundcf",x);
3484 : : }
3485 : :
3486 [ # # # # ]: 0 : switch(tx)
3487 : : {
3488 : 0 : case t_POL: return mkveccopy(x);
3489 : : case t_SER:
3490 : 0 : av = avma;
3491 : 0 : return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
3492 : : case t_RFRAC:
3493 : 0 : av = avma;
3494 : 0 : return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
3495 : : }
3496 : 0 : pari_err_TYPE("gboundcf",x);
3497 : 203 : return NULL; /* not reached */
3498 : : }
3499 : :
3500 : : static GEN
3501 : 21 : sfcont2(GEN b, GEN x, long k)
3502 : : {
3503 : 21 : pari_sp av = avma;
3504 : 21 : long lb = lg(b), tx = typ(x), i;
3505 : : GEN y,p1;
3506 : :
3507 [ + + ]: 21 : if (k)
3508 : : {
3509 [ + - ]: 7 : if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
3510 : 0 : lb = k+1;
3511 : : }
3512 : 14 : y = cgetg(lb,t_VEC);
3513 [ - + ]: 14 : if (lb==1) return y;
3514 [ + - ]: 14 : if (is_scalar_t(tx))
3515 : : {
3516 [ - + ][ # # ]: 14 : if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
3517 : : }
3518 [ # # ]: 0 : else if (tx == t_SER) x = ser2rfrac_i(x);
3519 : :
3520 [ - + ]: 14 : if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
3521 : 14 : for (i = 1;;)
3522 : : {
3523 [ + - ]: 70 : if (tx == t_REAL)
3524 : : {
3525 : 70 : long e = expo(x);
3526 [ + + ][ - + ]: 70 : if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
3527 : 70 : gel(y,i) = floorr(x);
3528 : 70 : p1 = subri(x, gel(y,i));
3529 : : }
3530 : : else
3531 : : {
3532 : 0 : gel(y,i) = gfloor(x);
3533 : 0 : p1 = gsub(x, gel(y,i));
3534 : : }
3535 [ + + ]: 70 : if (++i >= lb) break;
3536 [ - + ]: 56 : if (gequal0(p1)) break;
3537 : 56 : x = gdiv(gel(b,i),p1);
3538 : 56 : }
3539 : 14 : setlg(y,i);
3540 : 14 : return gerepilecopy(av,y);
3541 : : }
3542 : :
3543 : :
3544 : : GEN
3545 : 7 : gcf(GEN x) { return gboundcf(x,0); }
3546 : : GEN
3547 : 7 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
3548 : : GEN
3549 : 231 : contfrac0(GEN x, GEN b, long nmax)
3550 : : {
3551 : : long tb;
3552 : :
3553 [ + + ]: 231 : if (!b) return gboundcf(x,nmax);
3554 : 35 : tb = typ(b);
3555 [ + + ]: 35 : if (tb == t_INT) return gboundcf(x,itos(b));
3556 [ - + ]: 28 : if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
3557 [ + + ]: 28 : if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
3558 : 210 : return sfcont2(b,x,nmax);
3559 : : }
3560 : :
3561 : : GEN
3562 : 77 : contfracpnqn(GEN x, long n)
3563 : : {
3564 : 77 : pari_sp av = avma;
3565 : 77 : long i, lx = lg(x);
3566 : : GEN M,A,B, p0,p1, q0,q1;
3567 : :
3568 : 77 : p0 = gen_1; q0 = gen_0;
3569 [ - + ]: 77 : if (lx == 1)
3570 : : {
3571 [ # # ]: 0 : if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
3572 [ # # ]: 0 : if (n == 0) retmkmat(mkcol2(p0, q0));
3573 : 0 : return matid(2);
3574 : : }
3575 [ + + - ]: 77 : switch(typ(x))
3576 : : {
3577 : 35 : case t_VEC: case t_COL: A = x; B = NULL; break;
3578 : : case t_MAT:
3579 [ - + + ]: 42 : switch(lgcols(x))
3580 : : {
3581 : 0 : case 2: A = row(x,1); B = NULL; break;
3582 : 35 : case 3: A = row(x,2); B = row(x,1); break;
3583 : 7 : default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
3584 : 0 : return NULL; /*not reached*/
3585 : : }
3586 : 35 : break;
3587 : 0 : default: pari_err_TYPE("pnqn",x);
3588 : 0 : return NULL; /*not reached*/
3589 : : }
3590 [ + + ]: 70 : if (n >= 0)
3591 : : {
3592 [ + + ]: 28 : if (n == 0) { avma = av; retmkmat(mkcol2(p0, q0)); }
3593 : 14 : lx = minss(lx, n+1);
3594 : : }
3595 : 56 : M = cgetg(lx+1, t_MAT);
3596 : 56 : p1 = gel(A,1);
3597 [ + + ]: 56 : q1 = B? gel(B,1): gen_1;
3598 : 56 : gel(M,1) = mkcol2(p0,q0);
3599 : 56 : gel(M,2) = mkcol2(p1,q1);
3600 [ + + ]: 280 : for (i=2; i<lx; i++)
3601 : : {
3602 : 224 : GEN a = gel(A,i), p2,q2;
3603 [ + + ]: 224 : if (B) {
3604 : 119 : GEN b = gel(B,i);
3605 : 119 : p0 = gmul(b,p0);
3606 : 119 : q0 = gmul(b,q0);
3607 : : }
3608 : 224 : p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
3609 : 224 : q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
3610 : 224 : gel(M,i+1) = mkcol2(p1,q1);
3611 : : }
3612 [ + + ]: 56 : if (n < 0) M = mkmat2(gel(M,lx), gel(M,lx-1));
3613 : 70 : return gerepilecopy(av, M);
3614 : : }
3615 : : GEN
3616 : 14 : pnqn(GEN x) { return contfracpnqn(x,-1); }
3617 : :
3618 : : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
3619 : : static GEN
3620 : 35 : mod_to_frac(GEN x, GEN N, GEN B)
3621 : : {
3622 : : GEN a, b, A;
3623 [ + + ]: 35 : if (B)
3624 : : {
3625 : 21 : A = divii(shifti(N, -1), B);
3626 : : /* denominator bound useless, don't use it */
3627 [ - + ]: 21 : if (cmpii(A, B) < 0) B = NULL;
3628 : : }
3629 [ + + ]: 35 : if (!B)
3630 : : {
3631 : 14 : A = sqrti(shifti(N, -1));
3632 : 14 : B = A;
3633 : : }
3634 [ + + ][ - + ]: 35 : if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
3635 [ + - ]: 35 : return equali1(b)? a: mkfrac(a,b);
3636 : : }
3637 : :
3638 : : static GEN
3639 : 56 : mod_to_rfrac(GEN x, GEN N, long B)
3640 : : {
3641 : : GEN a, b;
3642 : 56 : long A, d = degpol(N);
3643 [ + + ]: 56 : if (B >= 0)
3644 : : {
3645 : 21 : A = d-1 - B;
3646 : : /* denominator bound useless, don't use it */
3647 [ + + ]: 21 : if (A < B) B = -1;
3648 : : }
3649 [ + + ]: 56 : if (B < 0)
3650 : : {
3651 : 42 : B = d >> 1;
3652 [ + + ]: 42 : A = odd(d)? B : B-1;
3653 : : }
3654 [ - + ]: 56 : if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
3655 [ - + ]: 56 : if (! RgXQ_ratlift(x, N, A, B, &a,&b)) return NULL;
3656 [ + + ]: 56 : if (degpol(RgX_gcd(a,b)) > 0) return NULL;
3657 : 56 : return gdiv(a,b);
3658 : : }
3659 : :
3660 : : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
3661 : : * of the continued fraction of x with b <= k maximal */
3662 : : static GEN
3663 : 0 : bestappr_frac(GEN x, GEN k)
3664 : : {
3665 : : pari_sp av;
3666 : : GEN p0, p1, p, q0, q1, q, a, y;
3667 : :
3668 [ # # ]: 0 : if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
3669 : 0 : av = avma; y = x;
3670 : 0 : p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
3671 : 0 : q1 = gen_0; q0 = gen_1;
3672 : 0 : x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
3673 : : for(;;)
3674 : : {
3675 : 0 : x = ginv(x); /* > 1 */
3676 [ # # ]: 0 : a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
3677 [ # # ]: 0 : if (cmpii(a,k) > 0)
3678 : : { /* next partial quotient will overflow limits */
3679 : : GEN n, d;
3680 : 0 : a = divii(subii(k, q1), q0);
3681 : 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3682 : 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3683 : : /* compare |y-p0/q0|, |y-p1/q1| */
3684 : 0 : n = gel(y,1);
3685 : 0 : d = gel(y,2);
3686 [ # # ]: 0 : if (absi_cmp(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
3687 : : mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
3688 : 0 : { p1 = p0; q1 = q0; }
3689 : 0 : break;
3690 : : }
3691 : 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3692 : 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3693 : :
3694 [ # # ]: 0 : if (cmpii(q0,k) > 0) break;
3695 : 0 : x = gsub(x,a); /* 0 <= x < 1 */
3696 [ # # ]: 0 : if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
3697 : :
3698 : 0 : }
3699 : 0 : return gerepileupto(av, gdiv(p1,q1));
3700 : : }
3701 : : /* bestappr(t_REAL != 0), to maximal accuracy */
3702 : : static GEN
3703 : 0 : bestappr_real_max(GEN x)
3704 : : {
3705 : 0 : pari_sp av = avma;
3706 : : GEN p0, p1, p, q0, q1, q, a;
3707 : : long e;
3708 : 0 : p1 = gen_1; a = p0 = floorr(x); q1 = gen_0; q0 = gen_1;
3709 : 0 : x = subri(x,a); /* 0 <= x < 1 */
3710 : 0 : e = bit_prec(x) - expo(x);
3711 : : for(;;)
3712 : : {
3713 : : long d;
3714 [ # # ][ # # ]: 0 : if (!signe(x) || expi(q0) > e) { p1 = p0; q1 = q0; break; }
3715 : 0 : x = invr(x); /* > 1 */
3716 : 0 : d = nbits2prec(expo(x) + 1);
3717 [ # # ]: 0 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
3718 : :
3719 : 0 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
3720 : 0 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3721 : 0 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3722 : 0 : x = subri(x,a); /* 0 <= x < 1 */
3723 : 0 : }
3724 : 0 : return gerepileupto(av, gdiv(p1,q1));
3725 : : }
3726 : : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
3727 : : * of the continued fraction of x with b <= k maximal */
3728 : : static GEN
3729 : 184976 : bestappr_real(GEN x, GEN k)
3730 : : {
3731 : 184976 : pari_sp av = avma;
3732 : : GEN kr, p0, p1, p, q0, q1, q, a, y;
3733 : :
3734 : 184976 : y = x;
3735 : 184976 : p1 = gen_1; a = p0 = floorr(x);
3736 : 184976 : q1 = gen_0; q0 = gen_1;
3737 : 184976 : x = subri(x,a); /* 0 <= x < 1 */
3738 [ + + ]: 184976 : if (!signe(x)) { cgiv(x); return a; }
3739 : 181449 : kr = itor(k, realprec(x));
3740 : : for(;;)
3741 : : {
3742 : : long d;
3743 : 478806 : x = invr(x); /* > 1 */
3744 [ + + ]: 478806 : if (cmprr(x,kr) > 0)
3745 : : { /* next partial quotient will overflow limits */
3746 : 180380 : a = divii(subii(k, q1), q0);
3747 : 180380 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3748 : 180380 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3749 : : /* compare |y-p0/q0|, |y-p1/q1| */
3750 [ + + ]: 180380 : if (absr_cmp(mulir(q1, subri(mulir(q0,y), p0)),
3751 : : mulir(q0, subri(mulir(q1,y), p1))) < 0)
3752 : 1563 : { p1 = p0; q1 = q0; }
3753 : 180380 : break;
3754 : : }
3755 : 298426 : d = nbits2prec(expo(x) + 1);
3756 [ - + ]: 298426 : if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
3757 : :
3758 : 298426 : a = truncr(x); /* truncr(x) will NOT raise e_PREC */
3759 : 298426 : p = addii(mulii(a,p0), p1); p1=p0; p0=p;
3760 : 298426 : q = addii(mulii(a,q0), q1); q1=q0; q0=q;
3761 : :
3762 [ + + ]: 298426 : if (cmpii(q0,k) > 0) break;
3763 : 297776 : x = subri(x,a); /* 0 <= x < 1 */
3764 [ + + ]: 297776 : if (!signe(x)) { p1 = p0; q1 = q0; break; }
3765 : 297357 : }
3766 : 184976 : return gerepileupto(av, gdiv(p1,q1));
3767 : : }
3768 : :
3769 : : /* k t_INT or NULL */
3770 : : static GEN
3771 : 280801 : bestappr_Q(GEN x, GEN k)
3772 : : {
3773 : 280801 : long lx, tx = typ(x), i;
3774 : : GEN a, y;
3775 : :
3776 [ - - + + : 280801 : switch(tx)
+ + - ]
3777 : : {
3778 : 0 : case t_INT: return icopy(x);
3779 [ # # ]: 0 : case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
3780 : : case t_REAL:
3781 [ + + ]: 215183 : if (!signe(x)) return gen_0;
3782 [ + - ]: 184976 : return k? bestappr_real(x, k): bestappr_real_max(x);
3783 : :
3784 : : case t_INTMOD: {
3785 : 21 : pari_sp av = avma;
3786 [ + + ]: 21 : a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
3787 : 14 : return gerepilecopy(av, a);
3788 : : }
3789 : : case t_PADIC: {
3790 : 14 : pari_sp av = avma;
3791 : 14 : long v = valp(x);
3792 [ + + ]: 14 : a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
3793 [ - + ]: 7 : if (v) a = gmul(a, powis(gel(x,2), v));
3794 : 7 : return gerepilecopy(av, a);
3795 : : }
3796 : :
3797 : : case t_COMPLEX: case t_POLMOD: case t_POL: case t_SER: case t_RFRAC:
3798 : : case t_VEC: case t_COL: case t_MAT:
3799 : 65583 : y = cgetg_copy(x, &lx);
3800 [ + - ]: 65583 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
3801 [ + + ]: 343114 : for (; i<lx; i++)
3802 : : {
3803 [ - + ]: 277531 : a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
3804 : 277531 : gel(y,i) = a;
3805 : : }
3806 [ - + ]: 65583 : if (tx == t_POL) return normalizepol(y);
3807 [ - + ]: 65583 : if (tx == t_SER) return normalize(y);
3808 : 65583 : return y;
3809 : : }
3810 : 0 : pari_err_TYPE("bestappr_Q",x);
3811 : 280801 : return NULL; /* not reached */
3812 : : }
3813 : :
3814 : : static GEN
3815 : 49 : bestappr_ser(GEN x, long B)
3816 : : {
3817 : 49 : long v = valp(x), lx = lg(x);
3818 : : GEN N, t;
3819 : 49 : x = normalizepol(ser2pol_i(x, lx));
3820 : 49 : N = monomial(gen_1, lx-2, varn(x));
3821 [ + + ]: 49 : t = mod_to_rfrac(x, N, B); if (!t) return NULL;
3822 [ + + ]: 42 : if (v)
3823 : : {
3824 : : GEN a, b;
3825 : : long vx;
3826 [ - + ]: 14 : if (typ(t) == t_POL) return RgX_mulXn(t, v);
3827 : : /* t_RFRAC */
3828 : 14 : vx = varn(x);
3829 : 14 : a = gel(t,1);
3830 : 14 : b = gel(t,2);
3831 : 14 : v -= RgX_valrem(b, &b);
3832 [ - + ][ # # ]: 14 : if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
3833 [ + + ]: 14 : if (v < 0) b = RgX_shift(b, -v);
3834 [ + - ]: 7 : else if (v > 0) {
3835 [ - + ][ # # ]: 7 : if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
3836 : 7 : a = RgX_shift(a, v);
3837 : : }
3838 : 14 : t = mkrfraccopy(a, b);
3839 : : }
3840 : 49 : return t;
3841 : : }
3842 : : static GEN bestappr_RgX(GEN x, long B);
3843 : : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
3844 : : * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
3845 : : static GEN
3846 : 63 : bestappr_RgX(GEN x, long B)
3847 : : {
3848 : 63 : long i, lx, tx = typ(x);
3849 : : GEN y, t;
3850 [ - + + + : 63 : switch(tx)
- - ]
3851 : : {
3852 : : case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
3853 : : case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
3854 : 0 : return gcopy(x);
3855 : :
3856 : : case t_RFRAC: {
3857 : 14 : pari_sp av = avma;
3858 [ + + ][ - + ]: 14 : if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
3859 : 7 : x = rfractoser(x, varn(gel(x,2)), 2*B+1);
3860 [ - + ]: 7 : t = bestappr_ser(x, B); if (!t) return NULL;
3861 : 7 : return gerepileupto(av, t);
3862 : : }
3863 : : case t_POLMOD: {
3864 : 7 : pari_sp av = avma;
3865 [ - + ]: 7 : t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
3866 : 7 : return gerepileupto(av, t);
3867 : : }
3868 : : case t_SER: {
3869 : 42 : pari_sp av = avma;
3870 [ + + ]: 42 : t = bestappr_ser(x, B); if (!t) return NULL;
3871 : 35 : return gerepileupto(av, t);
3872 : : }
3873 : :
3874 : : case t_VEC: case t_COL: case t_MAT:
3875 : 0 : y = cgetg_copy(x, &lx);
3876 [ # # ]: 0 : if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
3877 [ # # ]: 0 : for (; i<lx; i++)
3878 : : {
3879 [ # # ]: 0 : t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
3880 : 0 : gel(y,i) = t;
3881 : : }
3882 : 0 : return y;
3883 : : }
3884 : 0 : pari_err_TYPE("bestappr_RgX",x);
3885 : 63 : return NULL; /* not reached */
3886 : : }
3887 : :
3888 : : /* allow k = NULL: maximal accuracy */
3889 : : GEN
3890 : 3270 : bestappr(GEN x, GEN k)
3891 : : {
3892 : 3270 : pari_sp av = avma;
3893 [ + + ]: 3270 : if (k) { /* replace by floor(k) */
3894 [ + + - ]: 3256 : switch(typ(k))
3895 : : {
3896 : : case t_INT:
3897 : 49 : break;
3898 : : case t_REAL: case t_FRAC:
3899 : 3207 : k = floor_safe(k); /* left on stack for efficiency */
3900 [ - + ]: 3207 : if (!signe(k)) k = gen_1;
3901 : 3207 : break;
3902 : : default:
3903 : 0 : pari_err_TYPE("bestappr [bound type]", k);
3904 : 0 : break;
3905 : : }
3906 : : }
3907 : 3270 : x = bestappr_Q(x, k);
3908 [ + + ]: 3270 : if (!x) { avma = av; return cgetg(1,t_VEC); }
3909 : 3270 : return x;
3910 : : }
3911 : : GEN
3912 : 63 : bestapprPade(GEN x, long B)
3913 : : {
3914 : 63 : pari_sp av = avma;
3915 : 63 : GEN t = bestappr_RgX(x, B);
3916 [ + + ]: 63 : if (!t) { avma = av; return cgetg(1,t_VEC); }
3917 : 63 : return t;
3918 : : }
3919 : :
3920 : : /***********************************************************************/
3921 : : /** **/
3922 : : /** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
3923 : : /** **/
3924 : : /***********************************************************************/
3925 : :
3926 : : static GEN
3927 : 28 : get_quad(GEN f, GEN pol, long r)
3928 : : {
3929 : 28 : GEN p1 = gcoeff(f,1,2), q1 = gcoeff(f,2,2);
3930 [ + - ]: 28 : return mkquad(pol, r? subii(p1,q1): p1, q1);
3931 : : }
3932 : :
3933 : : /* replace f by f * [a,1; 1,0] */
3934 : : static void
3935 : 28 : update_f(GEN f, GEN a)
3936 : : {
3937 : : GEN p1;
3938 : 28 : p1 = gcoeff(f,1,1);
3939 : 28 : gcoeff(f,1,1) = addii(mulii(a,p1), gcoeff(f,1,2));
3940 : 28 : gcoeff(f,1,2) = p1;
3941 : :
3942 : 28 : p1 = gcoeff(f,2,1);
3943 : 28 : gcoeff(f,2,1) = addii(mulii(a,p1), gcoeff(f,2,2));
3944 : 28 : gcoeff(f,2,2) = p1;
3945 : 28 : }
3946 : :
3947 : : GEN
3948 : 14 : quadunit(GEN x)
3949 : : {
3950 : 14 : pari_sp av = avma, av2, lim;
3951 : : GEN pol, y, a, u, v, sqd, f;
3952 : : long r;
3953 : :
3954 : 14 : check_quaddisc_real(x, &r, "quadunit");
3955 : 14 : pol = quadpoly(x);
3956 : 14 : sqd = sqrti(x); av2 = avma; lim = stack_lim(av2,2);
3957 : 14 : a = shifti(addsi(r,sqd),-1);
3958 : 14 : f = mkmat2(mkcol2(a, gen_1), mkcol2(gen_1, gen_0)); /* [a,0; 1,0] */
3959 : 14 : u = stoi(r); v = gen_2;
3960 : : for(;;)
3961 : : {
3962 : : GEN u1, v1;
3963 : 28 : u1 = subii(mulii(a,v),u);
3964 : 28 : v1 = divii(subii(x,sqri(u1)),v);
3965 [ + + ]: 28 : if ( equalii(v,v1) ) {
3966 : 14 : y = get_quad(f,pol,r);
3967 : 14 : update_f(f,a);
3968 : 14 : y = gdiv(get_quad(f,pol,r), gconj(y));
3969 : 14 : break;
3970 : : }
3971 : 14 : a = divii(addii(sqd,u1), v1);
3972 [ - + ]: 14 : if ( equalii(u,u1) ) {
3973 : 0 : y = get_quad(f,pol,r);
3974 : 0 : y = gdiv(y, gconj(y));
3975 : 0 : break;
3976 : : }
3977 : 14 : update_f(f,a);
3978 : 14 : u = u1; v = v1;
3979 [ - + ]: 14 : if (low_stack(lim, stack_lim(av2,2)))
3980 : : {
3981 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit");
3982 : 0 : gerepileall(av2,4, &a,&f,&u,&v);
3983 : : }
3984 : 14 : }
3985 [ + - ]: 14 : if (signe(gel(y,3)) < 0) y = gneg(y);
3986 : 14 : return gerepileupto(av, y);
3987 : : }
3988 : :
3989 : : GEN
3990 : 42 : quadregulator(GEN x, long prec)
3991 : : {
3992 : 42 : pari_sp av = avma, av2, lim;
3993 : : GEN R, rsqd, u, v, sqd;
3994 : : long r, Rexpo;
3995 : :
3996 : 42 : check_quaddisc_real(x, &r, "quadregulator");
3997 : 42 : sqd = sqrti(x);
3998 : 42 : rsqd = gsqrt(x,prec);
3999 : 42 : Rexpo = 0; R = real2n(1, prec); /* = 2 */
4000 : 42 : av2 = avma; lim = stack_lim(av2,2);
4001 : 42 : u = stoi(r); v = gen_2;
4002 : : for(;;)
4003 : : {
4004 : 140 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
4005 : 140 : GEN v1 = divii(subii(x,sqri(u1)),v);
4006 [ + + ]: 140 : if (equalii(v,v1))
4007 : : {
4008 : 14 : R = sqrr(R); shiftr_inplace(R, -1);
4009 : 14 : R = mulrr(R, divri(addir(u1,rsqd),v));
4010 : 14 : break;
4011 : : }
4012 [ + + ]: 126 : if (equalii(u,u1))
4013 : : {
4014 : 28 : R = sqrr(R); shiftr_inplace(R, -1);
4015 : 28 : break;
4016 : : }
4017 : 98 : R = mulrr(R, divri(addir(u1,rsqd),v));
4018 : 98 : Rexpo += expo(R); setexpo(R,0);
4019 : 98 : u = u1; v = v1;
4020 [ - + ]: 98 : if (Rexpo & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
4021 [ - + ]: 98 : if (low_stack(lim, stack_lim(av2,2)))
4022 : : {
4023 [ # # ]: 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
4024 : 0 : gerepileall(av2,3, &R,&u,&v);
4025 : : }
4026 : 98 : }
4027 : 42 : R = logr_abs(divri(R,v));
4028 [ + - ]: 42 : if (Rexpo)
4029 : : {
4030 : 42 : GEN t = mulsr(Rexpo, mplog2(prec));
4031 : 42 : shiftr_inplace(t, 1);
4032 : 42 : R = addrr(R,t);
4033 : : }
4034 : 42 : return gerepileuptoleaf(av, R);
4035 : : }
4036 : :
4037 : : /*************************************************************************/
4038 : : /** **/
4039 : : /** CLASS NUMBER **/
4040 : : /** **/
4041 : : /*************************************************************************/
4042 : 42 : static int qfb_is_1(GEN f) { return equali1(gel(f,1)); }
4043 : 42 : static GEN qfb_pow(void *E, GEN f, GEN n) { (void)E; return powgi(f,n); }
4044 : : static const struct bb_group qfb_group={ NULL,qfb_pow,NULL,NULL,
4045 : : NULL,qfb_is_1,NULL};
4046 : :
4047 : : GEN
4048 : 28 : qfbclassno0(GEN x,long flag)
4049 : : {
4050 [ + + - ]: 28 : switch(flag)
4051 : : {
4052 : 14 : case 0: return map_proto_G(classno,x);
4053 : 14 : case 1: return map_proto_G(classno2,x);
4054 : 0 : default: pari_err_FLAG("qfbclassno");
4055 : : }
4056 : 28 : return NULL; /* not reached */
4057 : : }
4058 : :
4059 : : /* f^h = 1, return order(f) */
4060 : : static GEN
4061 : 14 : find_order(GEN f, GEN h) { return gen_order(f, h, NULL, &qfb_group); }
4062 : :
4063 : : static GEN
4064 : 14 : end_classno(GEN h, GEN hin, GEN forms)
4065 : : {
4066 : 14 : GEN a, b, p1, q, fh, fg, f = gel(forms,1);
4067 : 14 : long i, com, l = lg(forms);
4068 : :
4069 : 14 : h = find_order(f,h); /* H = <f> */
4070 : 14 : q = diviiround(hin, h); /* approximate order of G/H */
4071 [ + + ]: 154 : for (i=2; i < l; i++)
4072 : : {
4073 : 140 : pari_sp av = avma;
4074 : 140 : fg = powgi(gel(forms,i), h);
4075 : 140 : fh = powgi(fg, q);
4076 : 140 : a = gel(fh,1);
4077 [ + - ]: 140 : if (equali1(a)) continue;
4078 : 0 : b = gel(fh,2); p1 = fg;
4079 : 0 : for (com=1; ; com++, p1 = gmul(p1,fg))
4080 [ # # ][ # # ]: 0 : if (equalii(gel(p1,1), a) && absi_equal(gel(p1,2), b)) break;
4081 [ # # ]: 0 : if (signe(gel(p1,2)) == signe(b)) com = -com;
4082 : : /* f_i ^ h(q+com) = 1 */
4083 : 0 : q = addsi(com,q);
4084 [ # # ]: 0 : if (gequal0(q))
4085 : : { /* f^(ih) != 1 for all 0 < i <= oldq. Happen if the original upper bound
4086 : : for h was wrong */
4087 : : ulong c;
4088 : 0 : p1 = fh;
4089 : 0 : for (c=1; ; c++, p1 = gmul(p1,fh))
4090 [ # # ]: 0 : if (qfb_is_1(p1)) break;
4091 : 0 : q = mulsi(-com, find_order(fh, utoipos(c)));
4092 : : }
4093 : 0 : q = gerepileuptoint(av, q);
4094 : : }
4095 : 14 : return mulii(q,h);
4096 : : }
4097 : :
4098 : : /* Write x = Df^2, where D = fundamental discriminant,
4099 : : * P^E = factorisation of conductor f, with E[i] >= 0 */
4100 : : static void
4101 : 70 : corediscfact(GEN x, long xmod4, GEN *ptD, GEN *ptP, GEN *ptE)
4102 : : {
4103 : 70 : long s = signe(x), l, i;
4104 : 70 : GEN fa = absi_factor(x);
4105 : 70 : GEN d, P = gel(fa,1), E = gtovecsmall(gel(fa,2));
4106 : :
4107 : 70 : l = lg(P); d = gen_1;
4108 [ + + ]: 168 : for (i=1; i<l; i++)
4109 : : {
4110 [ + - ]: 98 : if (E[i] & 1) d = mulii(d, gel(P,i));
4111 : 98 : E[i] >>= 1;
4112 : : }
4113 [ - + ][ # # ]: 70 : if (!xmod4 && mod4(d) != ((s < 0)? 3: 1)) { d = shifti(d,2); E[1]--; }
[ # # ]
4114 [ + + ]: 70 : *ptD = (s < 0)? negi(d): d;
4115 : 70 : *ptP = P;
4116 : 70 : *ptE = E;
4117 : 70 : }
4118 : :
4119 : : static GEN
4120 : 56 : conductor_part(GEN x, long xmod4, GEN *ptD, GEN *ptreg)
4121 : : {
4122 : 56 : long l, i, s = signe(x);
4123 : : GEN E, H, D, P, reg;
4124 : :
4125 : 56 : corediscfact(x, xmod4, &D, &P, &E);
4126 : 56 : H = gen_1; l = lg(P);
4127 : : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
4128 [ + + ]: 140 : for (i=1; i<l; i++)
4129 : : {
4130 : 84 : long e = E[i];
4131 [ - + ]: 84 : if (e)
4132 : : {
4133 : 0 : GEN p = gel(P,i);
4134 : 0 : H = mulii(H, subis(p, kronecker(D,p)));
4135 [ # # ]: 0 : if (e >= 2) H = mulii(H, powiu(p,e-1));
4136 : : }
4137 : : }
4138 : :
4139 : : /* divide by [ O_K^* : O^* ] */
4140 [ + + ]: 56 : if (s < 0)
4141 : : {
4142 : 28 : reg = NULL;
4143 [ - - + ]: 28 : switch(itou_or_0(D))
4144 : : {
4145 : 0 : case 4: H = shifti(H,-1); break;
4146 : 28 : case 3: H = divis(H,3); break;
4147 : : }
4148 : : } else {
4149 : 28 : reg = quadregulator(D,DEFAULTPREC);
4150 [ - + ]: 28 : if (!equalii(x,D))
4151 : 0 : H = divii(H, roundr(divrr(quadregulator(x,DEFAULTPREC), reg)));
4152 : : }
4153 [ + + ]: 56 : if (ptreg) *ptreg = reg;
4154 : 56 : *ptD = D; return H;
4155 : : }
4156 : :
4157 : : static long
4158 : 14 : two_rank(GEN x)
4159 : : {
4160 : 14 : GEN p = gel(absi_factor(x),1);
4161 : 14 : long l = lg(p)-1;
4162 : : #if 0 /* positive disc not needed */
4163 : : if (signe(x) > 0)
4164 : : {
4165 : : long i;
4166 : : for (i=1; i<=l; i++)
4167 : : if (mod4(gel(p,i)) == 3) { l--; break; }
4168 : : }
4169 : : #endif
4170 : 14 : return l-1;
4171 : : }
4172 : :
4173 : : static GEN
4174 : 154 : sqr_primeform(GEN x, ulong p) { return redimag(qfisqr(primeform_u(x, p))); }
4175 : :
4176 : : static ulong
4177 [ + - ]: 224 : _low(GEN x) { return signe(x)? mod2BIL(x): 0; }
4178 : :
4179 : : /* h(x) for x<0 using Baby Step/Giant Step.
4180 : : * Assumes G is not too far from being cyclic.
4181 : : *
4182 : : * Compute G^2 instead of G so as to kill most of the non-cyclicity */
4183 : : GEN
4184 : 28 : classno(GEN x)
4185 : : {
4186 : 28 : const long MAXFORM = 12;
4187 : 28 : pari_sp av = avma, av2, lim;
4188 : : long r2, p, nforms, k, i, j, com, s;
4189 : : GEN forms, count, index, tabla, tablb, hash, p1, p2;
4190 : : GEN hin, h, f, fh, fg, ftest, Hf, D;
4191 : : forprime_t S;
4192 : :
4193 [ + + ]: 28 : if (signe(x) >= 0) return classno2(x);
4194 : :
4195 : 14 : check_quaddisc(x, &s, &k, "classno");
4196 [ - + ]: 14 : if (cmpiu(x,12) <= 0) return gen_1;
4197 : :
4198 : 14 : Hf = conductor_part(x, k, &D, NULL);
4199 [ - + ]: 14 : if (cmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
4200 : :
4201 : 14 : p2 = gsqrt(absi(D),DEFAULTPREC);
4202 : 14 : p1 = mulrr(divrr(p2,mppi(DEFAULTPREC)), dbltor(1.005)); /*overshoot by 0.5%*/
4203 : 14 : s = itos_or_0( truncr(shiftr(sqrtr(p2), 1)) );
4204 [ - + ]: 14 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
4205 [ - + ]: 14 : if (s < 10) s = 200;
4206 [ - + ]: 14 : else if (s < 20) s = 1000;
4207 [ + - ]: 14 : else if (s < 5000) s = 5000;
4208 : 14 : forms = vectrunc_init(MAXFORM+1);
4209 : :
4210 : 14 : u_forprime_init(&S, 2, s);
4211 : 14 : nforms = 0;
4212 [ + + ]: 9380 : while ( (p = u_forprime_next(&S)) )
4213 : : {
4214 : 9366 : long d, k = kroiu(D,p);
4215 : : pari_sp av3;
4216 [ - + ]: 9366 : if (!k) continue;
4217 [ + + ]: 9366 : if (k > 0)
4218 : : {
4219 [ + + ]: 4760 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
4220 : 4760 : d = p - 1;
4221 : : }
4222 : : else
4223 : 4606 : d = p + 1;
4224 : 9366 : av3 = avma; affrr(divru(mulur(p,p1),d), p1);
4225 : 9366 : avma = av3;
4226 : : }
4227 : 14 : r2 = two_rank(D);
4228 : 14 : h = hin = roundr(shiftr(p1, -r2));
4229 : 14 : s = 2*itos(gceil(sqrtnr(p1, 4)));
4230 [ - + ]: 14 : if (s > 10000) s = 10000;
4231 : :
4232 [ + + ]: 3598 : count = new_chunk(256); for (i=0; i<=255; i++) count[i]=0;
4233 : 14 : index = new_chunk(257);
4234 : 14 : tabla = new_chunk(10000);
4235 : 14 : tablb = new_chunk(10000);
4236 : 14 : hash = new_chunk(10000);
4237 : 14 : f = gel(forms,1);
4238 : 14 : p1 = fh = powgi(f, h);
4239 [ + + ]: 98 : for (i=0; i<s; i++, p1 = qficomp(p1,f))
4240 : : {
4241 : 84 : tabla[i] = _low(gel(p1,1));
4242 : 84 : tablb[i] = _low(gel(p1,2)); count[tabla[i]&255]++;
4243 : : }
4244 : : /* follow the idea of counting sort to avoid maintaining linked lists in
4245 : : * hashtable */
4246 [ + + ]: 3584 : index[0]=0; for (i=0; i< 255; i++) index[i+1] = index[i]+count[i];
4247 : : /* index[i] = # of forms hashing to <= i */
4248 [ + + ]: 98 : for (i=0; i<s; i++) hash[ index[tabla[i]&255]++ ] = i;
4249 [ + + ]: 3598 : index[0]=0; for (i=0; i<=255; i++) index[i+1] = index[i]+count[i];
4250 : : /* hash[index[i-1]..index[i]-1] = forms hashing to i */
4251 : :
4252 : 14 : fg = gpowgs(f,s); av2 = avma; lim = stack_lim(av2,2);
4253 : 14 : ftest = gpowgs(p1,0);
4254 : 14 : for (com=0; ; com++)
4255 : : {
4256 : : long j1, k, l;
4257 : : GEN a, b;
4258 : 28 : a = gel(ftest,1); k = _low(a);
4259 : 28 : b = gel(ftest,2); l = _low(b); j = k&255;
4260 [ + + ]: 28 : for (j1=index[j]; j1 < index[j+1]; j1++)
4261 : : {
4262 : 14 : long j2 = hash[j1];
4263 [ + - ][ + - ]: 14 : if (tabla[j2] == k && tablb[j2] == l)
4264 : : {
4265 : 14 : p1 = gmul(gpowgs(f,j2),fh);
4266 [ + - ][ + - ]: 14 : if (equalii(gel(p1,1), a) && absi_equal(gel(p1,2), b))
4267 : : { /* p1 = ftest or ftest^(-1), we are done */
4268 [ + - ]: 14 : if (signe(gel(p1,2)) == signe(b)) com = -com;
4269 : 14 : h = addii(addis(h,j2), mulss(s,com));
4270 : 14 : h = end_classno(h, hin, forms);
4271 : 14 : return gerepileuptoint(av, shifti(mulii(h,Hf), r2));
4272 : : }
4273 : : }
4274 : : }
4275 : 14 : ftest = gmul(ftest,fg);
4276 [ - + ]: 14 : if (equali1(gel(ftest,1))) pari_err_IMPL("classno with too small order");
4277 [ - + ]: 14 : if (low_stack(lim, stack_lim(av2,2))) ftest = gerepileupto(av2,ftest);
4278 : 42 : }
4279 : : }
4280 : :
4281 : : /* use Euler products */
4282 : : GEN
4283 : 42 : classno2(GEN x)
4284 : : {
4285 : 42 : pari_sp av = avma;
4286 : 42 : const long prec = DEFAULTPREC;
4287 : : long n, i, r, s;
4288 : : GEN p1, p2, S, p4, p5, p7, Hf, Pi, reg, logd, d, dr, D, half;
4289 : :
4290 : 42 : check_quaddisc(x, &s, &r, "classno2");
4291 [ + + ][ - + ]: 42 : if (s < 0 && cmpiu(x,12) <= 0) return gen_1;
4292 : :
4293 : 42 : Hf = conductor_part(x, r, &D, ®);
4294 [ + + ][ - + ]: 42 : if (s < 0 && cmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
4295 : :
4296 : 42 : Pi = mppi(prec);
4297 : 42 : d = absi(D); dr = itor(d, prec);
4298 : 42 : logd = logr_abs(dr);
4299 : 42 : p1 = sqrtr(divrr(mulir(d,logd), gmul2n(Pi,1)));
4300 [ + + ]: 42 : if (s > 0)
4301 : : {
4302 : 28 : GEN invlogd = invr(logd);
4303 : 28 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
4304 [ - + ]: 28 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
4305 : : }
4306 : 42 : n = itos_or_0( mptrunc(p1) );
4307 [ - + ]: 42 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
4308 : :
4309 : 42 : p4 = divri(Pi,d);
4310 : 42 : p7 = invr(sqrtr_abs(Pi));
4311 : 42 : half = real2n(-1, prec);
4312 [ + + ]: 42 : if (s > 0)
4313 : : { /* i = 1, shortcut */
4314 : 28 : p1 = sqrtr_abs(dr);
4315 : 28 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
4316 : 28 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
4317 [ + + ]: 1092 : for (i=2; i<=n; i++)
4318 : : {
4319 [ + + ]: 1064 : long k = kroiu(D,i); if (!k) continue;
4320 : 868 : p2 = mulir(sqru(i), p4);
4321 : 868 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
4322 : 868 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
4323 [ + + ]: 868 : S = (k>0)? addrr(S,p5): subrr(S,p5);
4324 : : }
4325 : 28 : S = shiftr(divrr(S,reg),-1);
4326 : : }
4327 : : else
4328 : : { /* i = 1, shortcut */
4329 : 14 : p1 = gdiv(sqrtr_abs(dr), Pi);
4330 : 14 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
4331 : 14 : S = addrr(p5, divrr(p1, mpexp(p4)));
4332 [ + + ]: 1904 : for (i=2; i<=n; i++)
4333 : : {
4334 [ - + ]: 1890 : long k = kroiu(D,i); if (!k) continue;
4335 : 1890 : p2 = mulir(sqru(i), p4);
4336 : 1890 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
4337 : 1890 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
4338 [ + + ]: 1890 : S = (k>0)? addrr(S,p5): subrr(S,p5);
4339 : : }
4340 : : }
4341 : 42 : return gerepileuptoint(av, mulii(Hf, roundr(S)));
4342 : : }
4343 : :
4344 : : static GEN
4345 : 14 : hclassno2(GEN x)
4346 : : {
4347 : : long i, l, s, xmod4;
4348 : : GEN Q, H, D, P, E;
4349 : :
4350 : 14 : x = negi(x);
4351 : 14 : check_quaddisc(x, &s, &xmod4, "hclassno");
4352 : 14 : corediscfact(x, xmod4, &D, &P, &E);
4353 : :
4354 : 14 : Q = quadclassunit0(D, 0, NULL, 0);
4355 : 14 : H = gel(Q,1); l = lg(P);
4356 : :
4357 : : /* H \prod_{p^e||f} (1 + (p^e-1)/(p-1))[ p - (D/p) ] */
4358 [ + + ]: 28 : for (i=1; i<l; i++)
4359 : : {
4360 : 14 : long e = E[i];
4361 [ - + ]: 14 : if (e)
4362 : : {
4363 : 0 : GEN p = gel(P,i), t = subis(p, kronecker(D,p));
4364 [ # # ]: 0 : if (e > 1) t = mulii(t, diviiexact(subis(powiu(p,e), 1), subis(p,1)));
4365 : 0 : H = mulii(H, addsi(1, t));
4366 : : }
4367 : : }
4368 [ - - + ]: 14 : switch( itou_or_0(D) )
4369 : : {
4370 : 0 : case 3: H = gdivgs(H, 3); break;
4371 : 0 : case 4: H = gdivgs(H, 2); break;
4372 : : }
4373 : 14 : return H;
4374 : : }
4375 : :
4376 : : GEN
4377 : 14 : hclassno(GEN x)
4378 : : {
4379 : : ulong a, b, b2, d, h;
4380 : : long s;
4381 : : int f;
4382 : :
4383 [ - + ]: 14 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
4384 : 14 : s = signe(x);
4385 [ - + ]: 14 : if (s < 0) return gen_0;
4386 [ - + ]: 14 : if (!s) return gdivgs(gen_1, -12);
4387 : :
4388 [ + - ][ - + ]: 14 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
4389 : :
4390 : 14 : d = itou_or_0(x);
4391 [ + - ][ + - ]: 14 : if (!d || d > 500000) return hclassno2(x);
4392 : :
4393 : 0 : h = 0; b = d&1; b2 = (1+d)>>2; f=0;
4394 [ # # ]: 0 : if (!b)
4395 : : {
4396 [ # # ]: 0 : for (a=1; a*a<b2; a++)
4397 [ # # ]: 0 : if (b2%a == 0) h++;
4398 : 0 : f = (a*a==b2); b=2; b2=(4+d)>>2;
4399 : : }
4400 [ # # ]: 0 : while (b2*3 < d)
4401 : : {
4402 [ # # ]: 0 : if (b2%b == 0) h++;
4403 [ # # ]: 0 : for (a=b+1; a*a < b2; a++)
4404 [ # # ]: 0 : if (b2%a == 0) h += 2;
4405 [ # # ]: 0 : if (a*a == b2) h++;
4406 : 0 : b += 2; b2 = (b*b+d)>>2;
4407 : : }
4408 [ # # ]: 0 : if (b2*3 == d)
4409 : : {
4410 : 0 : GEN y = cgetg(3,t_FRAC);
4411 : 0 : gel(y,1) = utoipos(3*h+1);
4412 : 0 : gel(y,2) = utoipos(3); return y;
4413 : : }
4414 [ # # ]: 0 : if (f)
4415 : : {
4416 : 0 : GEN y = cgetg(3,t_FRAC);
4417 : 0 : gel(y,1) = utoipos(2*h+1);
4418 : 0 : gel(y,2) = gen_2; return y;
4419 : : }
4420 : 14 : return utoipos(h);
4421 : : }
4422 : :
|
