Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_arith
19 :
20 : /*********************************************************************/
21 : /** **/
22 : /** FUNDAMENTAL DISCRIMINANTS **/
23 : /** **/
24 : /*********************************************************************/
25 : static long
26 1407 : fa_isfundamental(GEN F)
27 : {
28 1407 : GEN P = gel(F,1), E = gel(F,2);
29 1407 : long i, s, l = lg(P);
30 :
31 1407 : if (l == 1) return 1;
32 1400 : s = signe(gel(P,1)); /* = signe(x) */
33 1400 : if (!s) return 0;
34 1393 : if (s < 0) { l--; P = vecslice(P,2,l); E = vecslice(E,2,l); }
35 1393 : if (l == 1) return 0;
36 1386 : if (!absequaliu(gel(P,1), 2))
37 686 : i = 1; /* need x = 1 mod 4 */
38 : else
39 : {
40 700 : i = 2;
41 700 : switch(itou(gel(E,1)))
42 : {
43 182 : case 2: s = -s; break; /* need x/4 = 3 mod 4 */
44 84 : case 3: s = 0; break; /* no condition mod 4 */
45 434 : default: return 0;
46 : }
47 : }
48 1974 : for(; i < l; i++)
49 : {
50 1190 : if (!equali1(gel(E,i))) return 0;
51 1022 : if (s && Mod4(gel(P,i)) == 3) s = -s;
52 : }
53 784 : return s >= 0;
54 : }
55 : long
56 20482 : isfundamental(GEN x)
57 : {
58 20482 : if (typ(x) != t_INT)
59 : {
60 1407 : pari_sp av = avma;
61 1407 : long v = fa_isfundamental(check_arith_all(x,"isfundamental"));
62 1407 : return gc_long(av,v);
63 : }
64 19075 : return Z_isfundamental(x);
65 : }
66 :
67 : /* x fundamental ? */
68 : long
69 16652 : uposisfundamental(ulong x)
70 : {
71 16652 : ulong r = x & 15; /* x mod 16 */
72 16652 : if (!r) return 0;
73 15840 : switch(r & 3)
74 : { /* x mod 4 */
75 3431 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
76 5958 : case 1: return uissquarefree(x);
77 6451 : default: return 0;
78 : }
79 : }
80 : /* -x fundamental ? */
81 : long
82 29337 : unegisfundamental(ulong x)
83 : {
84 29337 : ulong r = x & 15; /* x mod 16 */
85 29337 : if (!r) return 0;
86 27692 : switch(r & 3)
87 : { /* x mod 4 */
88 6938 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
89 10425 : case 3: return uissquarefree(x);
90 10329 : default: return 0;
91 : }
92 : }
93 : long
94 25053 : sisfundamental(long x)
95 25053 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
96 :
97 : long
98 19642 : Z_isfundamental(GEN x)
99 : {
100 : long r;
101 19642 : switch(lgefint(x))
102 : {
103 7 : case 2: return 0;
104 9233 : case 3: return signe(x) < 0? unegisfundamental(x[2])
105 26858 : : uposisfundamental(x[2]);
106 : }
107 2010 : r = mod16(x);
108 2010 : if (!r) return 0;
109 1884 : if ((r & 3) == 0)
110 : {
111 : pari_sp av;
112 376 : r >>= 2; /* |x|/4 mod 4 */
113 376 : if (signe(x) < 0) r = 4-r;
114 376 : if (r == 1) return 0;
115 250 : av = avma;
116 250 : r = Z_issquarefree( shifti(x,-2) );
117 250 : return gc_long(av, r);
118 : }
119 1508 : r &= 3; /* |x| mod 4 */
120 1508 : if (signe(x) < 0) r = 4-r;
121 1508 : return (r==1) ? Z_issquarefree(x) : 0;
122 : }
123 :
124 : static GEN
125 2058 : fa_quaddisc(GEN f)
126 : {
127 2058 : GEN P = gel(f,1), E = gel(f,2), s = gen_1;
128 2058 : long i, l = lg(P);
129 6783 : for (i = 1; i < l; i++) /* possibly including -1 */
130 4725 : if (mpodd(gel(E,i))) s = mulii(s, gel(P,i));
131 2058 : if (Mod4(s) > 1) s = shifti(s,2);
132 2058 : return s;
133 : }
134 :
135 : GEN
136 2912 : quaddisc(GEN x)
137 : {
138 2912 : const pari_sp av = avma;
139 2912 : long tx = typ(x);
140 : GEN F;
141 2912 : if (is_rational_t(tx)) F = factor(x);
142 : else
143 : {
144 1407 : F = check_arith_all(x,"quaddisc");
145 1407 : if (tx == t_VEC && typ(gel(x,1)) == t_INT
146 1407 : && Z_issquarefree_fact(gel(x,2)))
147 : {
148 854 : x = gel(x,1);
149 854 : return (Mod4(x) > 1)? shifti(x, 2): icopy(x);
150 : }
151 : }
152 2058 : return gerepileuptoint(av, fa_quaddisc(F));
153 : }
154 :
155 :
156 : /***********************************************************************/
157 : /** **/
158 : /** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
159 : /** **/
160 : /***********************************************************************/
161 : /* replace f by f * [u,1; 1,0] */
162 : static void
163 22156268 : update_f(GEN f, GEN u)
164 : {
165 22156268 : GEN a = gcoeff(f,1,1), b = gcoeff(f,1,2);
166 22156268 : GEN c = gcoeff(f,2,1), d = gcoeff(f,2,2);
167 22156268 : gcoeff(f,1,1) = addmulii(b, u,a); gcoeff(f,1,2) = a;
168 22156268 : gcoeff(f,2,1) = addmulii(d, u,c); gcoeff(f,2,2) = c;
169 22156268 : }
170 :
171 : /* f is a vector of matrices and i an index whose bits give the non-zero
172 : * entries; the product of the non zero entries is the actual result.
173 : * if i odd, f[1] is implicitely [f[1],1;1,0] */
174 :
175 : static long
176 22317344 : update_fm(GEN f, GEN a, long i)
177 : {
178 : #ifdef LONG_IS_64BIT
179 19129152 : const long LIM = 10;
180 : #else
181 3188192 : const long LIM = 18;
182 : #endif
183 22317344 : pari_sp av = avma;
184 : long k, v;
185 : GEN u;
186 22317344 : if (!odd(i)) { gel(f,1) = a; return i+1; }
187 22236806 : u = gel(f, 1);
188 22236806 : if (typ(u) == t_INT) /* [u,1;1,0] * [a,1;1,0] */
189 80538 : { gel(f,1) = mkmat22(addiu(mulii(a, u), 1), u, a, gen_1); return i; }
190 22156268 : update_f(u, a); if (lgefint(gcoeff(u,1,1)) < LIM) return i;
191 80538 : v = vals(i+1); gel(f,1) = gen_0;
192 161003 : for (k = 1; k < v; k++) { u = ZM2_mul(gel(f,k+1), u); gel(f,k+1) = gen_0; }
193 80538 : if (v != 1) u = gerepileupto(av, u);
194 80538 : gel(f,v+1) = u; return i+1;
195 : }
196 : /* \prod f[j]; if first only return column 1 */
197 : static GEN
198 69041 : prod_fm(GEN f, long i, long first)
199 : {
200 69041 : long k, v = vals(i) + 1;
201 69041 : GEN u = gel(f, v);
202 : /* i a power of 2: f[1] can't be a t_INT */
203 69041 : if (!(i >>= v)) return first? gel(u,1): u;
204 138 : for (k = v+1; i; i >>= 1, k++)
205 118 : if (odd(i))
206 : {
207 73 : GEN w = gel(f,k);
208 73 : switch(typ(u))
209 : {
210 0 : case t_INT: update_f(w, u);
211 0 : u = first? gel(w,1): w; break;
212 6 : case t_COL: /* implies 'first' */
213 6 : u = ZM_ZC_mul(w, u); break;
214 67 : default: /* t_MAT */
215 67 : u = first? ZM_ZC_mul(w, gel(u,1)): ZM2_mul(w, u); break;
216 : }
217 45 : }
218 20 : return u;
219 : }
220 :
221 : GEN
222 69048 : quadunit0(GEN x, long v)
223 : {
224 69048 : GEN y = quadunit(x);
225 69041 : if (v==-1) v = fetch_user_var("w");
226 69041 : setvarn(gel(y,1), v); return y;
227 : }
228 :
229 : struct uimod { GEN N, T; };
230 : static GEN
231 973 : ui_pow(void *E, GEN x, GEN n)
232 973 : { struct uimod *S = (struct uimod*)E; return FpXQ_pow(x, n, S->T, S->N); }
233 : static int
234 994 : ui_equal1(GEN x) { return degpol(x) < 1; }
235 : static const struct bb_group
236 : ui_group={ NULL,ui_pow,NULL,NULL,NULL,ui_equal1,NULL};
237 :
238 : static void
239 77 : quadunit_mod(GEN D, GEN d, GEN N, GEN *pu, GEN *pv)
240 : {
241 : GEN u1, u2, v1, v2, p, q, q1, u, v;
242 77 : int m = mpodd(D);
243 77 : pari_sp av = avma;
244 77 : p = (mpodd(d) == m)? d: subiu(d, 1);
245 77 : u1 = negi(p); u2 = gen_2;
246 77 : v1 = gen_1; v2 = gen_0; q = gen_2;
247 77 : q1 = shifti(subii(D, sqri(p)), -1);
248 : for(;;)
249 105 : {
250 182 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
251 182 : p = subii(d, r);
252 182 : if (equalii(p1, p) && signe(v2))
253 : { /* even period */
254 7 : u = addmulii(sqri(u2), D, sqri(v2));
255 7 : v = shifti(mulii(u2,v2), 1);
256 7 : break;
257 : }
258 175 : t = Fp_addmul(u1, A, u2, N); u1 = u2; u2 = t;
259 175 : t = Fp_addmul(v1, A, v2, N); v1 = v2; v2 = t;
260 175 : t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
261 175 : if (equalii(q, t))
262 : { /* odd period */
263 70 : u = addmulii(mulii(u1,u2), D, mulii(v1,v2));
264 70 : v = addmulii(mulii(u1,v2), u2, v1);
265 70 : break;
266 : }
267 105 : if (gc_needed(av, 2))
268 : {
269 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit_mod");
270 0 : gerepileall(av, 7, &p, &u1,&u2,&v1,&v2, &q,&q1);
271 : }
272 : }
273 77 : *pu = modii(u, N);
274 77 : *pv = modii(v, N); if (m) *pu = Fp_sub(*pu, *pv, N);
275 77 : }
276 : GEN
277 0 : quadunit_basecase(GEN D)
278 : {
279 0 : pari_sp av0 = avma;
280 : GEN d, u1, u2, v1, v2, p, q, q1, u, v, a, b, c;
281 : int m;
282 0 : long first = 1;
283 0 : check_quaddisc_real(D, NULL, "quadunit");
284 0 : m = mpodd(D); d = sqrtremi(D, &a);
285 0 : p = d; q1 = shifti(a, -1); q = gen_2;
286 0 : if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
287 0 : u1 = gen_2; u2 = p;
288 0 : v1 = gen_0; v2 = gen_1;
289 : for(;;)
290 0 : {
291 0 : GEN t = q;
292 0 : if (first) { first = 0; q = q1; }
293 : else
294 : {
295 0 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
296 0 : p = subii(d, r);
297 0 : if (equalii(p1, p)) /* even period */
298 0 : { a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break; }
299 0 : r = addmulii(u1, A, u2); u1 = u2; u2 = r;
300 0 : r = addmulii(v1, A, v2); v1 = v2; v2 = r;
301 0 : q = submulii(q1, A, subii(p, p1));
302 : }
303 0 : q1 = t;
304 0 : if (equalii(q, t))
305 : { /* odd period */
306 0 : a = mulii(u1, u2); b = mulii(v1, v2);
307 0 : c = mulii(addii(u1, v1), addii(u2, v2)); break;
308 : }
309 : }
310 0 : u = diviiexact(addmulii(a, D, b), q);
311 0 : v = diviiexact(subii(c, addii(a, b)), q);
312 0 : if (m == 1) u = subii(u, v);
313 0 : u = shifti(u, -1);
314 0 : return gerepilecopy(av0, mkquad(quadpoly_i(D), u, v));
315 : }
316 :
317 : GEN
318 69048 : quadunit(GEN D)
319 : {
320 69048 : pari_sp av0 = avma, av;
321 : GEN f, d, p, q, q1, u, v, a, b, c;
322 : int m;
323 69048 : long i = 0;
324 69048 : check_quaddisc_real(D, NULL, "quadunit");
325 69041 : m = mpodd(D); d = sqrtremi(D, &a); av = avma;
326 69041 : p = d; q1 = shifti(a, -1); q = gen_2;
327 69041 : if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
328 69041 : f = zerovec(2 + (expi(D)>>1));
329 69041 : gel(f,1) = mkmat22(p, gen_2, gen_1, gen_0);
330 : for(;;)
331 22369158 : {
332 22438199 : GEN t = q, u1,u2, v1,v2;
333 22438199 : if (!i) { i = 1; q = q1; }
334 : else
335 : {
336 22369158 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
337 22369158 : p = subii(d, r);
338 22369158 : if (equalii(p1, p))
339 : { /* even period */
340 51814 : f = prod_fm(f, i, 1); u2 = gel(f,1); v2 = gel(f,2);
341 51814 : a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break;
342 : }
343 22317344 : i = update_fm(f, A, i);
344 22317344 : q = submulii(q1, A, subii(p, p1));
345 : }
346 22386385 : q1 = t;
347 22386385 : if (equalii(q, t))
348 : { /* odd period */
349 17227 : f = prod_fm(f, i, 0);
350 17227 : u2 = gcoeff(f,1,1); u1 = gcoeff(f,1,2); a = mulii(u1, u2);
351 17227 : v2 = gcoeff(f,2,1); v1 = gcoeff(f,2,2); b = mulii(v1, v2);
352 17227 : c = mulii(addii(u1, v1), addii(u2, v2)); break;
353 : }
354 22369158 : if (gc_needed(av, 2))
355 : {
356 83 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit (%ld)", i);
357 83 : gerepileall(av, 4, &p, &f, &q,&q1);
358 : }
359 : }
360 69041 : u = diviiexact(addmulii(a, D, b), q);
361 69041 : v = diviiexact(subii(c, addii(a, b)), q);
362 69041 : if (m == 1) u = subii(u, v);
363 69041 : u = shifti(u, -1);
364 69041 : return gerepilecopy(av0, mkquad(quadpoly_i(D), u, v));
365 : }
366 : /* D > 0, d = sqrti(D) */
367 : static GEN
368 69097 : quadunit_q(GEN D, GEN d, long *pN)
369 : {
370 69097 : pari_sp av = avma;
371 : GEN p, q, q1;
372 69097 : long first = 1;
373 69097 : p = (Mod2(d) == Mod2(D))? d: subiu(d, 1);
374 69097 : q = gen_2;
375 69097 : q1 = shifti(subii(D, sqri(p)), -1);
376 : for(;;)
377 1000930 : {
378 1070027 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
379 1070027 : p = subii(d, r);
380 1087317 : if (!first && equalii(p1, p)) { *pN = 1; return q; } /* even period */
381 1018220 : first = 0;
382 1018220 : t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
383 1018220 : if (equalii(q, t)) { *pN = -1; return q; } /* odd period */
384 1000930 : if (gc_needed(av, 2))
385 : {
386 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunitnorm");
387 0 : gerepileall(av, 3, &p, &q, &q1);
388 : }
389 : }
390 : }
391 : long
392 69027 : quadunitnorm(GEN D)
393 : {
394 69027 : pari_sp av = avma;
395 : long s, r;
396 69027 : check_quaddisc(D, &s, &r, "quadunitnorm");
397 69020 : if (s < 0) return 1;
398 69020 : (void)quadunit_q(D, sqrti(D), &s); return gc_long(av, s);
399 : }
400 : GEN
401 119 : quadunitindex(GEN D, GEN N)
402 : {
403 119 : pari_sp av = avma, av2;
404 : GEN y, u, v, q, d, P, E, F, a, faH, H;
405 : struct uimod S;
406 : long r, s;
407 :
408 119 : check_quaddisc(D, &s, &r, "quadunitindex");
409 112 : if ((F = check_arith_pos(N,"quadunitindex")))
410 35 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
411 105 : if (equali1(N)) return gen_1;
412 98 : if (s < 0) switch(itos_or_0(D)) {
413 7 : case -3: return utoipos(3);
414 7 : case -4: return utoipos(2);
415 7 : default: return gen_1;
416 : }
417 77 : if (!F) F = Z_factor(N);
418 77 : d = sqrti(D); av2 = avma;
419 77 : q = gerepileuptoint(av2, quadunit_q(D, d, &s));
420 77 : if (mpodd(N) && equali1(gcdii(q, N)))
421 : {
422 7 : quadunit_mod(D, d, N, &u, &v);
423 7 : q = Fp_inv(shifti(q, 1), N);
424 7 : u = Fp_mul(u, q, N);
425 7 : v = Fp_mul(v, q, N); v = modii(shifti(v, 1), N);
426 : }
427 : else
428 : {
429 70 : GEN M = shifti(mulii(q, N), 1);
430 70 : quadunit_mod(D, d, M, &u, &v);
431 70 : u = diviiexact(u, q);
432 70 : v = diviiexact(v, q); u = shifti(u,-1);
433 : }
434 : /* fundamental unit = y mod N */
435 77 : S.N = N; S.T = quadpoly_i(D); y = deg1pol_shallow(v, u, 0);
436 77 : H = quadclassnoF_fact(D, gel(F,1), gel(F,2));
437 77 : a = Z_smoothen(H, gel(F,1), &P, &E);
438 77 : faH = mkmat2(P, E);
439 77 : if (a) faH = merge_factor(Z_factor(a), faH,(void*)&cmpii,cmp_nodata);
440 77 : y = gen_order(y, mkvec2(H, faH), (void*)&S, &ui_group);
441 77 : return gerepileupto(av, y);
442 : }
443 :
444 : GEN
445 42 : quadregulator(GEN x, long prec)
446 : {
447 42 : pari_sp av = avma, av2;
448 : GEN R, rsqd, u, v, sqd;
449 : long r, e;
450 :
451 42 : check_quaddisc_real(x, &r, "quadregulator");
452 42 : sqd = sqrti(x);
453 42 : rsqd = gsqrt(x,prec); av2 = avma;
454 42 : e = 0; R = real2n(1, prec); u = utoi(r); v = gen_2;
455 : for(;;)
456 49 : {
457 91 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
458 91 : GEN v1 = divii(subii(x,sqri(u1)),v);
459 91 : if (equalii(v,v1)) { R = mulrr(sqrr(R), divri(addir(u1,rsqd),v)); break; }
460 63 : if (equalii(u,u1)) { R = sqrr(R); break; }
461 49 : R = mulrr(R, divri(addir(u1,rsqd),v));
462 49 : e += expo(R); setexpo(R,0);
463 49 : u = u1; v = v1;
464 49 : if (e & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
465 49 : if (gc_needed(av2,2))
466 : {
467 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
468 0 : gerepileall(av2,3, &R,&u,&v);
469 : }
470 : }
471 42 : R = divri(R, v); e = 2*e - 1;
472 : /* avoid loss of accuracy */
473 42 : if (!((e + expo(R)) & ~EXPOBITS)) { setexpo(R, e + expo(R)); e = 0; }
474 42 : R = logr_abs(R);
475 42 : if (e) R = addrr(R, mulsr(e, mplog2(prec)));
476 42 : return gerepileuptoleaf(av, R);
477 : }
478 :
479 : /*************************************************************************/
480 : /** **/
481 : /** CLASS NUMBER **/
482 : /** **/
483 : /*************************************************************************/
484 :
485 : int
486 8243821 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
487 :
488 9198757 : static GEN qfi_pow(void *E, GEN f, GEN n)
489 9198757 : { return E? nupow(f,n,(GEN)E): qfbpow_i(f,n); }
490 6350443 : static GEN qfi_comp(void *E, GEN f, GEN g)
491 6350443 : { return E? nucomp(f,g,(GEN)E): qfbcomp_i(f,g); }
492 : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
493 : gidentical,qfb_equal1,NULL};
494 :
495 : GEN
496 3022086 : qfi_order(GEN q, GEN o)
497 3022086 : { return gen_order(q, o, NULL, &qfi_group); }
498 :
499 : GEN
500 0 : qfi_log(GEN a, GEN g, GEN o)
501 0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
502 :
503 : GEN
504 646322 : qfi_Shanks(GEN a, GEN g, long n)
505 : {
506 646322 : pari_sp av = avma;
507 : GEN T, X;
508 : long rt_n, c;
509 :
510 646322 : a = qfbred_i(a);
511 646322 : g = qfbred_i(g);
512 :
513 646322 : rt_n = sqrt((double)n);
514 646322 : c = n / rt_n;
515 646322 : c = (c * rt_n < n + 1) ? c + 1 : c;
516 :
517 646322 : T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
518 646322 : X = gen_Shanks(T, a, c, NULL, &qfi_group);
519 646322 : return X? gerepileuptoint(av, X): gc_NULL(av);
520 : }
521 :
522 : GEN
523 455 : qfbclassno0(GEN x,long flag)
524 : {
525 455 : switch(flag)
526 : {
527 441 : case 0: return map_proto_G(classno,x);
528 14 : case 1: return map_proto_G(classno2,x);
529 0 : default: pari_err_FLAG("qfbclassno");
530 : }
531 : return NULL; /* LCOV_EXCL_LINE */
532 : }
533 :
534 : /* f^h = 1, return order(f). Set *pfao to its factorization */
535 : static GEN
536 5159 : find_order(void *E, GEN f, GEN h, GEN *pfao)
537 : {
538 5159 : GEN v = gen_factored_order(f, h, E, &qfi_group);
539 5159 : *pfao = gel(v,2); return gel(v,1);
540 : }
541 :
542 : static int
543 301 : ok_q(GEN q, GEN h, GEN d2, long r2)
544 : {
545 301 : if (d2)
546 : {
547 245 : if (r2 <= 2 && !mpodd(q)) return 0;
548 245 : return is_pm1(Z_ppo(q,d2));
549 : }
550 : else
551 : {
552 56 : if (r2 <= 1 && !mpodd(q)) return 0;
553 56 : return is_pm1(Z_ppo(q,h));
554 : }
555 : }
556 :
557 : /* a,b given by their factorizations. Return factorization of lcm(a,b).
558 : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
559 : static GEN
560 182 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
561 : {
562 182 : GEN P = ZC_union_shallow(gel(Fa,1), gel(Fb,1));
563 182 : GEN A = gen_1, B = gen_1;
564 182 : long i, l = lg(P);
565 182 : GEN E = cgetg(l, t_COL);
566 707 : for (i=1; i<l; i++)
567 : {
568 525 : GEN p = gel(P,i);
569 525 : long va = Z_pval(a,p);
570 525 : long vb = Z_pval(b,p);
571 525 : if (va < vb)
572 : {
573 203 : B = mulii(B,powiu(p,vb));
574 203 : gel(E,i) = utoi(vb);
575 : }
576 : else
577 : {
578 322 : A = mulii(A,powiu(p,va));
579 322 : gel(E,i) = utoi(va);
580 : }
581 : }
582 182 : *pA = A;
583 182 : *pB = B; return mkmat2(P,E);
584 : }
585 :
586 : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
587 : static void
588 182 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
589 : {
590 182 : GEN A,B, g1 = *pg1, d1 = *pd1;
591 182 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
592 182 : *pg1 = gmul(qfbpow_i(g1, diviiexact(d1,A)), qfbpow_i(f, diviiexact(o,B)));
593 182 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
594 182 : }
595 :
596 : /* *pD = coredisc(x), *pR = regulator (x > 0) or NULL */
597 : GEN
598 455 : quadclassnoF(GEN x, GEN *pD)
599 : {
600 : GEN D, H, P, E;
601 455 : if (lgefint(x) == 3)
602 : {
603 : ulong d, h;
604 407 : if (signe(x) < 0)
605 : {
606 380 : h = unegquadclassnoF(x[2], &d);
607 380 : if (pD) *pD = utoineg(d);
608 : }
609 : else
610 : {
611 27 : h = uposquadclassnoF(x[2], &d);
612 27 : if (pD) *pD = utoipos(d);
613 : }
614 407 : return utoi(h);
615 : }
616 48 : D = coredisc2_fact(absZ_factor(x), signe(x), &P, &E);
617 48 : H = quadclassnoF_fact(D, P, E);
618 : /* divide by [ O_K^* : O^* ] */
619 48 : if (signe(x) < 0) switch(itou_or_0(D))
620 : { /* |x| >= 2^BIL, hence x != D */
621 0 : case 4: H = shifti(H,-1); break;
622 0 : case 3: H = divis(H,3); break;
623 : }
624 48 : else if (!equalii(x,D))
625 8 : H = diviiexact(H, quadunitindex(D, mkmat2(P, zc_to_ZC(E))));
626 48 : if (pD) *pD = D; return H;
627 : }
628 :
629 : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ];
630 : * s = 1 or -1; D = s * d; assume Df^2 fits in an ulong */
631 : ulong
632 4358 : uquadclassnoF_fact(ulong d, long s, GEN P, GEN E)
633 : {
634 4358 : long i, l = lg(P);
635 4358 : ulong H = 1;
636 4567 : for (i = 1; i < l; i++)
637 : {
638 209 : ulong p = P[i], e = E[i];
639 209 : long D = (long)(p == 2? d & 7: d % p), a;
640 209 : if (s < 0) D = -D;
641 209 : a = kross(D,p);
642 209 : if (!a)
643 56 : H *= upowuu(p, e);
644 : else
645 : {
646 153 : H *= p - a;
647 153 : if (e >= 2) H *= upowuu(p, e-1);
648 : }
649 : }
650 4358 : return H;
651 : }
652 : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
653 : GEN
654 125 : quadclassnoF_fact(GEN D, GEN P, GEN E)
655 : {
656 125 : long i, l = lg(P);
657 : GEN H;
658 125 : if (typ(E) != t_VECSMALL) E = vec_to_vecsmall(E);
659 245 : for (i = 1, H = gen_1; i < l; i++)
660 : {
661 120 : GEN p = gel(P,i);
662 120 : long e = E[i], s = kronecker(D,p);
663 120 : if (!s)
664 42 : H = mulii(H, e == 1? p: powiu(p, e));
665 : else
666 : {
667 78 : H = mulii(H, subis(p, s));
668 78 : if (e >= 2) H = mulii(H, e == 2? p: powiu(p,e-1));
669 : }
670 : }
671 125 : return H;
672 : }
673 : ulong
674 380 : unegquadclassnoF(ulong x, ulong *pD)
675 : {
676 380 : pari_sp av = avma;
677 : GEN E, P;
678 380 : ulong D = coredisc2u_fact(factoru(x), -1, &P, &E);
679 380 : long H = uquadclassnoF_fact(D, -1, P, E);
680 380 : switch(D)
681 : { /* divide by [ O_K^* : O^* ] */
682 0 : case 4: H >>= 1; break;
683 0 : case 3: H /= 3; break;
684 : }
685 380 : *pD = D; return gc_ulong(av, H);
686 : }
687 : ulong
688 27 : uposquadclassnoF(ulong x, ulong *pD)
689 : {
690 : GEN P, E;
691 27 : ulong D = coredisc2u_fact(factoru(x), 1, &P, &E);
692 27 : ulong H = uquadclassnoF_fact(D, 1, P, E);
693 27 : if (x != D)
694 : {
695 13 : GEN F = mkvec2(utoipos(usqrt(x / D)), mkmat2(zc_to_ZC(P), zc_to_ZC(E)));
696 13 : H /= itou(quadunitindex(utoipos(D), F));
697 : }
698 27 : *pD = D; return H;
699 : }
700 :
701 : static long
702 413 : two_rank(GEN x)
703 : {
704 413 : GEN p = gel(absZ_factor(x),1);
705 413 : long l = lg(p)-1;
706 : #if 0 /* positive disc not needed */
707 : if (signe(x) > 0)
708 : {
709 : long i;
710 : for (i=1; i<=l; i++)
711 : if (mod4(gel(p,i)) == 3) { l--; break; }
712 : }
713 : #endif
714 413 : return l-1;
715 : }
716 :
717 : static GEN
718 7847 : sqr_primeform(GEN x, ulong p) { return qfbsqr_i(primeform_u(x, p)); }
719 : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
720 : static GEN
721 413 : get_forms(GEN D, GEN *pL)
722 : {
723 413 : const long MAXFORM = 20;
724 413 : GEN L, sqrtD = gsqrt(absi_shallow(D),DEFAULTPREC);
725 413 : GEN forms = vectrunc_init(MAXFORM+1);
726 413 : long s, nforms = 0;
727 : ulong p;
728 : forprime_t S;
729 413 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
730 413 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
731 413 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
732 413 : if (s < 10) s = 200;
733 413 : else if (s < 20) s = 1000;
734 413 : else if (s < 5000) s = 5000;
735 413 : u_forprime_init(&S, 2, s);
736 13458732 : while ( (p = u_forprime_next(&S)) )
737 : {
738 13458319 : long d, k = kroiu(D,p);
739 : pari_sp av2;
740 13458319 : if (!k) continue;
741 13457605 : if (k > 0)
742 : {
743 6730451 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
744 6730451 : d = p-1;
745 : }
746 : else
747 6727154 : d = p+1;
748 13457605 : av2 = avma; affrr(divru(mulur(p,L),d), L); set_avma(av2);
749 : }
750 413 : *pL = L; return forms;
751 : }
752 :
753 : /* h ~ #G, return o = order of f, set fao = its factorization */
754 : static GEN
755 511 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
756 : {
757 511 : long s = minss(itos(sqrti(h)), 10000);
758 511 : GEN T = gen_Shanks_init(f, s, E, &qfi_group);
759 511 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
760 511 : return find_order(E, f, addiu(v,1), pfao);
761 : }
762 :
763 : /* if g = 1 in G/<f> ? */
764 : static int
765 5684 : equal1(void *E, GEN T, ulong N, GEN g)
766 5684 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
767 :
768 : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
769 : * FIXME: should be gen_order, but equal1 has the wrong prototype */
770 : static GEN
771 350 : relative_order(void *E, GEN a, GEN o, ulong N, GEN T)
772 : {
773 350 : pari_sp av = avma;
774 : long i, l;
775 : GEN m;
776 :
777 350 : m = get_arith_ZZM(o);
778 350 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
779 350 : o = gel(m,1);
780 350 : m = gel(m,2); l = lgcols(m);
781 1148 : for (i = l-1; i; i--)
782 : {
783 798 : GEN t, y, p = gcoeff(m,i,1);
784 798 : long j, e = itos(gcoeff(m,i,2));
785 798 : if (l == 2) {
786 49 : t = gen_1;
787 49 : y = a;
788 : } else {
789 749 : t = diviiexact(o, powiu(p,e));
790 749 : y = powgi(a, t);
791 : }
792 798 : if (equal1(E, T,N,y)) o = t;
793 : else {
794 364 : for (j = 1; j < e; j++)
795 : {
796 84 : y = powgi(y, p);
797 84 : if (equal1(E, T,N,y)) break;
798 : }
799 357 : if (j < e) {
800 77 : if (j > 1) p = powiu(p, j);
801 77 : o = mulii(t, p);
802 : }
803 : }
804 : }
805 350 : return gerepilecopy(av, o);
806 : }
807 :
808 : /* h(x) for x<0 using Baby Step/Giant Step.
809 : * Assumes G is not too far from being cyclic.
810 : *
811 : * Compute G^2 instead of G so as to kill most of the noncyclicity */
812 : GEN
813 441 : classno(GEN x)
814 : {
815 441 : pari_sp av = avma;
816 : long r2, k, s, i, l;
817 : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
818 : void *E;
819 :
820 441 : if (signe(x) >= 0) return classno2(x);
821 :
822 413 : check_quaddisc(x, &s, &k, "classno");
823 413 : if (abscmpiu(x,12) <= 0) return gen_1;
824 :
825 413 : Hf = quadclassnoF(x, &D);
826 413 : if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
827 413 : forms = get_forms(D, &L);
828 413 : r2 = two_rank(D);
829 413 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
830 :
831 413 : l = lg(forms);
832 413 : order_bound = const_vec(l-1, NULL);
833 413 : E = expi(D) > 60? (void*)sqrtnint(shifti(absi_shallow(D),-2),4): NULL;
834 413 : g1 = gel(forms,1);
835 413 : gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
836 413 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
837 413 : d2 = NULL; /* not computed yet */
838 413 : if (is_pm1(q)) goto END;
839 7315 : for (i=2; i < l; i++)
840 : {
841 6951 : GEN o, fao, a, F, fd, f = gel(forms,i);
842 6951 : fd = qfbpow_i(f, d1); if (is_pm1(gel(fd,1))) continue;
843 182 : F = qfbpow_i(fd, q);
844 182 : a = gel(F,1);
845 182 : o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
846 : /* f^(d1 q) = 1 */
847 182 : fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
848 182 : o = find_order(E, f, fao, &fao);
849 182 : gel(order_bound,i) = o;
850 : /* o = order of f, fao = factor(o) */
851 182 : update_g1(&g1,&d1,&fad1, f,o,fao);
852 182 : q = diviiround(hin, d1);
853 182 : if (is_pm1(q)) goto END;
854 : }
855 : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
856 364 : if (expi(q) > 3)
857 : { /* q large: compute d2, 2nd elt divisor */
858 308 : ulong N, n = 2*itou(sqrti(d1));
859 308 : GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
860 308 : d2 = gen_1;
861 308 : N = itou( gceil(gdivgu(d1,n)) ); /* order(g1) <= n*N */
862 5047 : for (i = 1; i < l; i++)
863 : {
864 4802 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
865 4802 : if (!B) B = find_order(E, f, fad1, /*junk*/&d);
866 4802 : f = qfbpow_i(f,d2);
867 4802 : if (equal1(E, T,N,f)) continue;
868 350 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
869 : /* f^B = 1 */
870 350 : d = relative_order(E, f, B, N,T);
871 350 : d2= mulii(d,d2);
872 350 : D = mulii(d1,d2);
873 350 : q = diviiround(hin,D);
874 350 : if (is_pm1(q)) { d1 = D; goto END; }
875 : }
876 : /* very probably, d2 is the 2nd elementary divisor */
877 245 : d1 = D; /* product of first two elt divisors */
878 : }
879 : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
880 : * 2-rank */
881 301 : if (!ok_q(q,d1,d2,r2))
882 : {
883 0 : GEN q0 = q;
884 : long d;
885 0 : if (cmpii(mulii(q,d1), hin) < 0)
886 : { /* try q = q0+1,-1,+2,-2 */
887 0 : d = 1;
888 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
889 : }
890 : else
891 : { /* q0-1,+1,-2,+2 */
892 0 : d = -1;
893 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
894 : }
895 : }
896 301 : d1 = mulii(d1,q);
897 :
898 413 : END:
899 413 : return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
900 : }
901 :
902 : /* use Euler products */
903 : GEN
904 42 : classno2(GEN x)
905 : {
906 42 : pari_sp av = avma;
907 42 : const long prec = DEFAULTPREC;
908 : long n, i, s;
909 42 : GEN p1, p2, S, p4, p5, p7, Hf, Pi, logd, sqrtd, D, half, reg = NULL;
910 :
911 42 : check_quaddisc(x, &s, NULL, "classno2");
912 42 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
913 :
914 42 : Hf = quadclassnoF(x, &D);
915 42 : if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
916 :
917 42 : Pi = mppi(prec);
918 42 : sqrtd = sqrtr_abs(itor(D, prec));
919 42 : logd = logr_abs(sqrtd); shiftr_inplace(logd,1);
920 42 : p1 = sqrtr_abs(divrr(mulir(D,logd), gmul2n(Pi,1)));
921 42 : if (s > 0)
922 : {
923 35 : GEN invlogd = invr(logd);
924 35 : reg = quadregulator(D, prec);
925 35 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
926 35 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
927 : }
928 42 : n = itos_or_0( mptrunc(p1) );
929 42 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
930 :
931 42 : p4 = divri(Pi, D); setsigne(p4, 1);
932 42 : p7 = invr(sqrtr_abs(Pi));
933 42 : half = real2n(-1, prec);
934 42 : if (s > 0)
935 : { /* i = 1, shortcut */
936 35 : p1 = sqrtd;
937 35 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
938 35 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
939 588 : for (i=2; i<=n; i++)
940 : {
941 553 : long k = kroiu(D,i); if (!k) continue;
942 455 : p2 = mulir(sqru(i), p4);
943 455 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
944 455 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
945 455 : S = (k>0)? addrr(S,p5): subrr(S,p5);
946 : }
947 35 : S = shiftr(divrr(S,reg),-1);
948 : }
949 : else
950 : { /* i = 1, shortcut */
951 7 : p1 = gdiv(sqrtd, Pi);
952 7 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
953 7 : S = addrr(p5, divrr(p1, mpexp(p4)));
954 952 : for (i=2; i<=n; i++)
955 : {
956 945 : long k = kroiu(D,i); if (!k) continue;
957 945 : p2 = mulir(sqru(i), p4);
958 945 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
959 945 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
960 945 : S = (k>0)? addrr(S,p5): subrr(S,p5);
961 : }
962 : }
963 42 : return gerepileuptoint(av, mulii(Hf, roundr(S)));
964 : }
965 :
966 : /* 1 + q + ... + q^v, v > 0 */
967 : static GEN
968 973 : geomsumu(ulong q, long v)
969 : {
970 973 : GEN u = utoipos(1+q);
971 1267 : for (; v > 1; v--) u = addui(1, mului(q, u));
972 973 : return u;
973 : }
974 : static GEN
975 973 : geomsum(GEN q, long v)
976 : {
977 : GEN u;
978 973 : if (lgefint(q) == 3) return geomsumu(q[2], v);
979 0 : u = addiu(q,1);
980 0 : for (; v > 1; v--) u = addui(1, mulii(q, u));
981 0 : return u;
982 : }
983 :
984 : /* 1+p+...+p^(e-1), e >= 1; assuming result fits in an ulong */
985 : static ulong
986 2433 : usumpow(ulong p, long e)
987 : {
988 : ulong q;
989 : long i;
990 2433 : if (p == 2) return (1UL << e) - 1; /* also OK if e = BITS_IN_LONG */
991 1523 : e--; for (i = 1, q = p + 1; i < e; i++) q = p*q + 1;
992 1488 : return q;
993 : }
994 : long
995 44025 : uhclassnoF_fact(GEN faF, long D)
996 : {
997 44025 : GEN P = gel(faF,1), E = gel(faF,2);
998 44025 : long i, t, l = lg(P);
999 101954 : for (i = t = 1; i < l; i++)
1000 : {
1001 57929 : long p = P[i], e = E[i], s = kross(D,p);
1002 57929 : if (e == 1) { t *= 1 + p - s; continue; }
1003 15137 : if (s == 1) { t *= upowuu(p,e); continue; }
1004 2433 : t *= 1 + usumpow(p, e) * (p - s);
1005 : }
1006 44025 : return t;
1007 : }
1008 : /* Hurwitz(D F^2)/ Hurwitz(D)
1009 : * = \sum_{f|F} f \prod_{p|f} (1-kro(D/p)/p)
1010 : * = \prod_{p^e || F} (1 + (p^e-1) / (p-1) * (p-kro(D/p))) */
1011 : GEN
1012 60723 : hclassnoF_fact(GEN P, GEN E, GEN D)
1013 : {
1014 : GEN H;
1015 60723 : long i, l = lg(P);
1016 60723 : if (l == 1) return gen_1;
1017 56169 : for (i = 1, H = NULL; i < l; i++)
1018 : {
1019 30242 : GEN t, p = gel(P,i);
1020 30242 : long e = E[i], s = kronecker(D,p);
1021 30243 : if (e == 1) t = addiu(p, 1-s);
1022 1555 : else if (s == 1) t = powiu(p, e);
1023 973 : else t = addui(1, mulii(subis(p, s), geomsum(p, e-1)));
1024 30244 : H = H? mulii(H,t): t;
1025 : }
1026 25927 : return H;
1027 : }
1028 : static GEN
1029 60726 : hclassno6_large(GEN x)
1030 : {
1031 60726 : GEN H = NULL, P, E, D = coredisc2_fact(absZ_factor(x), -1, &P, &E);
1032 60728 : long l = lg(P);
1033 :
1034 60728 : if (l > 1 && lgefint(x) == 3)
1035 : { /* F != 1, second chance */
1036 25928 : ulong h = hclassno6u_no_cache(x[2]);
1037 25928 : if (h) H = utoipos(h);
1038 : }
1039 60728 : if (!H)
1040 : {
1041 60728 : H = quadclassno(D);
1042 60727 : switch(itou_or_0(D))
1043 : {
1044 49 : case 3: H = shifti(H,1);break;
1045 7 : case 4: H = muliu(H,3); break;
1046 60671 : default:H = muliu(H,6); break;
1047 : }
1048 0 : }
1049 60723 : return mulii(H, hclassnoF_fact(P, E, D));
1050 : }
1051 :
1052 : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
1053 : GEN
1054 132477 : hclassno6(GEN x)
1055 : {
1056 132477 : ulong d = itou_or_0(x);
1057 132478 : if (d)
1058 : { /* create cache if d small */
1059 132477 : ulong h = d < 500000 ? hclassno6u(d): hclassno6u_no_cache(d);
1060 132479 : if (h) return utoipos(h);
1061 : }
1062 60726 : return hclassno6_large(x);
1063 : }
1064 :
1065 : GEN
1066 49 : hclassno(GEN x)
1067 : {
1068 : long a, s;
1069 49 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
1070 49 : s = signe(x);
1071 49 : if (s < 0) return gen_0;
1072 49 : if (!s) return gdivgs(gen_1, -12);
1073 49 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
1074 49 : return gdivgu(hclassno6(x), 6);
1075 : }
1076 :
1077 : /* return [N',v]; v contains all x mod N' s.t. x^2 + B x + C = 0 modulo N */
1078 : GEN
1079 2774818 : Zn_quad_roots(GEN N, GEN B, GEN C)
1080 : {
1081 2774818 : pari_sp av = avma;
1082 2774818 : GEN fa = NULL, D, w, v, P, E, F0, Q0, F, mF, A, Q, T, R, Np, N4;
1083 : long l, i, j, ct;
1084 :
1085 2774818 : if ((fa = check_arith_non0(N,"Zn_quad_roots")))
1086 : {
1087 35994 : N = typ(N) == t_VEC? gel(N,1): factorback(N);
1088 35994 : fa = clean_Z_factor(fa);
1089 : }
1090 2774827 : N = absi_shallow(N);
1091 2774826 : N4 = shifti(N,2);
1092 2774794 : D = modii(subii(sqri(B), shifti(C,2)), N4);
1093 2774776 : if (!signe(D))
1094 : { /* (x + B/2)^2 = 0 (mod N), D = B^2-4C = 0 (4N) => B even */
1095 812 : if (!fa) fa = Z_factor(N);
1096 812 : P = gel(fa,1);
1097 812 : E = ZV_to_zv(gel(fa,2));
1098 812 : l = lg(P);
1099 1757 : for (i = 1; i < l; i++) E[i] = (E[i]+1) >> 1;
1100 812 : Np = factorback2(P, E); /* x = -B mod N' */
1101 812 : B = shifti(B,-1);
1102 812 : return gerepilecopy(av, mkvec2(Np, mkvec(Fp_neg(B,Np))));
1103 : }
1104 2773964 : if (!fa)
1105 2738204 : fa = Z_factor(N4);
1106 : else /* convert to factorization of N4 = 4*N */
1107 35760 : fa = famat_reduce(famat_mulpows_shallow(fa, gen_2, 2));
1108 2774005 : P = gel(fa,1); l = lg(P);
1109 2774005 : E = ZV_to_zv(gel(fa,2));
1110 2774001 : F = cgetg(l, t_VEC);
1111 2773998 : mF= cgetg(l, t_VEC); F0 = gen_0;
1112 2773998 : Q = cgetg(l, t_VEC); Q0 = gen_1;
1113 6677279 : for (i = j = 1, ct = 0; i < l; i++)
1114 : {
1115 6039191 : GEN p = gel(P,i), q, f, mf, D0;
1116 6039191 : long t2, s = E[i], t = Z_pvalrem(D, p, &D0), d = s - t;
1117 6039207 : if (d <= 0)
1118 : {
1119 1377859 : q = powiu(p, (s+1)>>1);
1120 2385187 : Q0 = mulii(Q0, q); continue;
1121 : }
1122 : /* d > 0 */
1123 6611737 : if (odd(t)) return NULL;
1124 4475903 : t2 = t >> 1;
1125 4475903 : if (i > 1)
1126 : { /* p > 2 */
1127 2823211 : if (kronecker(D0, p) == -1) return NULL;
1128 1375067 : q = powiu(p, s - t2);
1129 1375069 : f = Zp_sqrt(D0, p, d);
1130 1375079 : if (!f) return NULL; /* p was not actually prime... */
1131 1375065 : if (t2) f = mulii(powiu(p,t2), f);
1132 1375065 : mf = Fp_neg(f, q);
1133 : }
1134 : else
1135 : { /* p = 2 */
1136 1652692 : if (d <= 3)
1137 : {
1138 1282770 : if (d == 3 && Mod8(D0) != 1) return NULL;
1139 1058231 : if (d == 2 && Mod4(D0) != 1) return NULL;
1140 1007327 : Q0 = int2n(1+t2); F0 = NULL; continue;
1141 : }
1142 369922 : if (Mod8(D0) != 1) return NULL;
1143 143143 : q = int2n(d - 1 + t2);
1144 143143 : f = Z2_sqrt(D0, d);
1145 143143 : if (t2) f = shifti(f, t2);
1146 143143 : mf = Fp_neg(f, q);
1147 : }
1148 1518167 : gel(Q,j) = q;
1149 1518167 : gel(F,j) = f;
1150 1518167 : gel(mF,j)= mf; j++;
1151 : }
1152 638088 : setlg(Q,j);
1153 638097 : setlg(F,j);
1154 638095 : setlg(mF,j);
1155 638096 : if (is_pm1(Q0)) A = leafcopy(F);
1156 : else
1157 : { /* append the fixed congruence (F0 mod Q0) */
1158 603181 : if (!F0) F0 = shifti(Q0,-1);
1159 603181 : A = shallowconcat(F, F0);
1160 603200 : Q = shallowconcat(Q, Q0);
1161 : }
1162 638124 : ct = 1 << (j-1);
1163 638124 : T = ZV_producttree(Q);
1164 638083 : R = ZV_chinesetree(Q,T);
1165 638093 : Np = gmael(T, lg(T)-1, 1);
1166 638093 : B = modii(B, Np);
1167 638104 : if (!signe(B)) B = NULL;
1168 638104 : Np = shifti(Np, -1); /* N' = (\prod_i Q[i]) / 2 */
1169 638091 : w = cgetg(3, t_VEC);
1170 638110 : gel(w,1) = icopy(Np);
1171 638131 : gel(w,2) = v = cgetg(ct+1, t_VEC);
1172 638133 : l = lg(F);
1173 2737040 : for (j = 1; j <= ct; j++)
1174 : {
1175 2098926 : pari_sp av2 = avma;
1176 2098926 : long m = j - 1;
1177 : GEN u;
1178 6101649 : for (i = 1; i < l; i++)
1179 : {
1180 4002723 : gel(A,i) = (m&1L)? gel(mF,i): gel(F,i);
1181 4002723 : m >>= 1;
1182 : }
1183 2098926 : u = ZV_chinese_tree(A,Q,T,R); /* u mod N' st u^2 = B^2-4C modulo 4N */
1184 2098886 : if (B) u = subii(u,B);
1185 2098886 : gel(v,j) = gerepileuptoint(av2, modii(shifti(u,-1), Np));
1186 : }
1187 638114 : return gerepileupto(av, w);
1188 : }
|
