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 1206 : fa_isfundamental(GEN F)
27 : {
28 1206 : GEN P = gel(F,1), E = gel(F,2);
29 1206 : long i, s, l = lg(P);
30 :
31 1206 : if (l == 1) return 1;
32 1200 : s = signe(gel(P,1)); /* = signe(x) */
33 1200 : if (!s) return 0;
34 1194 : if (s < 0) { l--; P = vecslice(P,2,l); E = vecslice(E,2,l); }
35 1194 : if (l == 1) return 0;
36 1188 : if (!absequaliu(gel(P,1), 2))
37 588 : i = 1; /* need x = 1 mod 4 */
38 : else
39 : {
40 600 : i = 2;
41 600 : switch(itou(gel(E,1)))
42 : {
43 156 : case 2: s = -s; break; /* need x/4 = 3 mod 4 */
44 72 : case 3: s = 0; break; /* no condition mod 4 */
45 372 : default: return 0;
46 : }
47 : }
48 1692 : for(; i < l; i++)
49 : {
50 1020 : if (!equali1(gel(E,i))) return 0;
51 876 : if (s && Mod4(gel(P,i)) == 3) s = -s;
52 : }
53 672 : return s >= 0;
54 : }
55 : long
56 17638 : isfundamental(GEN x)
57 : {
58 17638 : if (typ(x) != t_INT)
59 : {
60 1206 : pari_sp av = avma;
61 1206 : long v = fa_isfundamental(check_arith_all(x,"isfundamental"));
62 1206 : return gc_long(av,v);
63 : }
64 16432 : return Z_isfundamental(x);
65 : }
66 :
67 : /* x fundamental ? */
68 : long
69 13376 : uposisfundamental(ulong x)
70 : {
71 13376 : ulong r = x & 15; /* x mod 16 */
72 13376 : if (!r) return 0;
73 12752 : switch(r & 3)
74 : { /* x mod 4 */
75 2781 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
76 5013 : case 1: return uissquarefree(x);
77 4958 : default: return 0;
78 : }
79 : }
80 : /* -x fundamental ? */
81 : long
82 25806 : unegisfundamental(ulong x)
83 : {
84 25806 : ulong r = x & 15; /* x mod 16 */
85 25806 : if (!r) return 0;
86 24468 : switch(r & 3)
87 : { /* x mod 4 */
88 6495 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
89 9691 : case 3: return uissquarefree(x);
90 8282 : default: return 0;
91 : }
92 : }
93 : long
94 21522 : sisfundamental(long x)
95 21522 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
96 :
97 : long
98 16839 : Z_isfundamental(GEN x)
99 : {
100 : long r;
101 16839 : switch(lgefint(x))
102 : {
103 6 : case 2: return 0;
104 8574 : case 3: return signe(x) < 0? unegisfundamental(x[2])
105 23398 : : uposisfundamental(x[2]);
106 : }
107 2009 : r = mod16(x);
108 2009 : if (!r) return 0;
109 1883 : 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 1507 : r &= 3; /* |x| mod 4 */
120 1507 : if (signe(x) < 0) r = 4-r;
121 1507 : return (r==1) ? Z_issquarefree(x) : 0;
122 : }
123 :
124 : static GEN
125 1496 : fa_quaddisc(GEN f)
126 : {
127 1496 : GEN P = gel(f,1), E = gel(f,2), s = gen_1;
128 1496 : long i, l = lg(P);
129 4965 : for (i = 1; i < l; i++) /* possibly including -1 */
130 3469 : if (mpodd(gel(E,i))) s = mulii(s, gel(P,i));
131 1496 : if (Mod4(s) > 1) s = shifti(s,2);
132 1496 : return s;
133 : }
134 :
135 : GEN
136 2106 : quaddisc(GEN x)
137 : {
138 2106 : const pari_sp av = avma;
139 2106 : long tx = typ(x);
140 : GEN F;
141 2106 : if (is_rational_t(tx)) F = factor(x);
142 : else
143 : {
144 1005 : F = check_arith_all(x,"quaddisc");
145 1005 : if (tx == t_VEC && typ(gel(x,1)) == t_INT
146 1005 : && Z_issquarefree_fact(gel(x,2)))
147 : {
148 610 : x = gel(x,1);
149 610 : return (Mod4(x) > 1)? shifti(x, 2): icopy(x);
150 : }
151 : }
152 1496 : return gc_INT(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 15143094 : update_f(GEN f, GEN u)
164 : {
165 15143094 : GEN a = gcoeff(f,1,1), b = gcoeff(f,1,2);
166 15143094 : GEN c = gcoeff(f,2,1), d = gcoeff(f,2,2);
167 15143094 : gcoeff(f,1,1) = addmulii(b, u,a); gcoeff(f,1,2) = a;
168 15143094 : gcoeff(f,2,1) = addmulii(d, u,c); gcoeff(f,2,2) = c;
169 15143094 : }
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] may be an int: implicitely represent [f[1],1;1,0] */
174 : static long
175 15257670 : update_fm(GEN f, GEN a, long i)
176 : {
177 : #ifdef LONG_IS_64BIT
178 12206136 : const long LIM = 10;
179 : #else
180 3051534 : const long LIM = 18;
181 : #endif
182 15257670 : pari_sp av = avma;
183 : long k, v;
184 : GEN u;
185 15257670 : if (!odd(i)) { gel(f,1) = a; return i+1; }
186 15200382 : u = gel(f, 1);
187 15200382 : if (typ(u) == t_INT) /* [u,1;1,0] * [a,1;1,0] */
188 57288 : { gel(f,1) = mkmat22(addiu(mulii(a, u), 1), u, a, gen_1); return i; }
189 15143094 : update_f(u, a); if (lgefint(gcoeff(u,1,1)) < LIM) return i;
190 57288 : v = vals(i+1); gel(f,1) = gen_0;
191 114538 : for (k = 1; k < v; k++) { u = ZM2_mul(gel(f,k+1), u); gel(f,k+1) = gen_0; }
192 57288 : if (v != 1) u = gc_upto(av, u);
193 57288 : gel(f,v+1) = u; return i+1;
194 : }
195 : /* \prod f[j]; if first only return column 1 */
196 : static GEN
197 5 : prod_fm(GEN f, long i, long first)
198 : {
199 5 : long k, v = vals(i) + 1;
200 5 : GEN u = gel(f, v);
201 : /* i a power of 2: f[1] can't be a t_INT */
202 5 : if (!(i >>= v)) return first? gel(u,1): u;
203 75 : for (k = v+1; i; i >>= 1, k++)
204 70 : if (odd(i))
205 : {
206 38 : GEN w = gel(f,k);
207 38 : switch(typ(u))
208 : {
209 0 : case t_INT: update_f(w, u);
210 0 : u = first? gel(w,1): w; break;
211 0 : case t_COL: /* implies 'first' */
212 0 : u = ZM_ZC_mul(w, u); break;
213 38 : default: /* t_MAT */
214 38 : u = first? ZM_ZC_mul(w, gel(u,1)): ZM2_mul(w, u); break;
215 : }
216 : }
217 5 : return u;
218 : }
219 :
220 : GEN
221 49321 : quadunit0(GEN x, long v)
222 : {
223 49321 : GEN y = quadunit(x);
224 49316 : if (v==-1) v = fetch_user_var("w");
225 49316 : setvarn(gel(y,1), v); return y;
226 : }
227 :
228 : struct uimod { GEN N, T; };
229 : static GEN
230 14945 : ui_pow(void *E, GEN x, GEN n)
231 14945 : { struct uimod *S = (struct uimod*)E; return FpXQ_pow(x, n, S->T, S->N); }
232 : static int
233 31395 : ui_equal1(GEN x) { return degpol(x) < 1; }
234 : static const struct bb_group
235 : ui_group={ NULL,ui_pow,NULL,NULL,NULL,ui_equal1,NULL};
236 :
237 : static void
238 70 : quadunit_uvmod(GEN D, GEN d, GEN N, GEN *pu, GEN *pv)
239 : {
240 : GEN u1, u2, v1, v2, p, q, q1, u, v;
241 70 : int m = mpodd(D), first = 1;
242 70 : pari_sp av = avma;
243 70 : p = (mpodd(d) == m)? d: subiu(d, 1);
244 70 : u1 = negi(p); u2 = gen_2;
245 70 : v1 = gen_1; v2 = gen_0; q = gen_2;
246 70 : q1 = shifti(subii(D, sqri(p)), -1);
247 : for(;;)
248 220 : {
249 290 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
250 290 : p = subii(d, r);
251 290 : if (equalii(p1, p) && !first)
252 : { /* even period */
253 10 : u = addmulii(sqri(u2), D, sqri(v2));
254 10 : v = shifti(mulii(u2,v2), 1);
255 10 : break;
256 : }
257 280 : first = 0;
258 280 : t = Fp_addmul(u1, A, u2, N); u1 = u2; u2 = t;
259 280 : t = Fp_addmul(v1, A, v2, N); v1 = v2; v2 = t;
260 280 : t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
261 280 : if (equalii(q, t))
262 : { /* odd period */
263 60 : u = addmulii(mulii(u1,u2), D, mulii(v1,v2));
264 60 : v = addmulii(mulii(u1,v2), u2, v1);
265 60 : break;
266 : }
267 220 : if (gc_needed(av, 2))
268 : {
269 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit_uvmod");
270 0 : (void)gc_all(av, 7, &p, &u1,&u2,&v1,&v2, &q,&q1);
271 : }
272 : }
273 70 : *pu = modii(u, N);
274 70 : *pv = modii(v, N); if (m) *pu = Fp_sub(*pu, *pv, N);
275 70 : }
276 : /* fundamental unit is u + vx mod quadpoly(D); always called with D
277 : * fundamental and relatively small but would work in all cases. Should be
278 : * called whenever the fundamental unit is so "small" that asymptotically
279 : * fast multiplication is not used in the continued fraction loop */
280 : static void
281 49311 : quadunit_uv_basecase(GEN D, GEN *pu, GEN *pv)
282 : {
283 49311 : GEN u1, u2, v1, v2, p, q, q1, u, v, a, b, c, d = sqrtremi(D, &a);
284 49311 : int m = mpodd(D);
285 49311 : long first = 1;
286 :
287 49311 : p = d; q1 = shifti(a, -1); q = gen_2;
288 49311 : if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
289 49311 : u1 = gen_2; u2 = p;
290 49311 : v1 = gen_0; v2 = gen_1;
291 : for(;;)
292 622181 : {
293 671492 : GEN t = q;
294 671492 : if (first) { first = 0; q = q1; }
295 : else
296 : {
297 622181 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
298 622181 : p = subii(d, r);
299 622181 : if (equalii(p1, p)) /* even period */
300 32670 : { a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break; }
301 589511 : r = addmulii(u1, A, u2); u1 = u2; u2 = r;
302 589511 : r = addmulii(v1, A, v2); v1 = v2; v2 = r;
303 589511 : q = submulii(q1, A, subii(p, p1));
304 : }
305 638822 : q1 = t;
306 638822 : if (equalii(q, t))
307 : { /* odd period */
308 16641 : a = mulii(u1, u2); b = mulii(v1, v2);
309 16641 : c = mulii(addii(u1, v1), addii(u2, v2)); break;
310 : }
311 : }
312 49311 : u = diviiexact(addmulii(a, D, b), q);
313 49311 : v = diviiexact(subii(c, addii(a, b)), q);
314 49311 : if (m == 1) u = subii(u, v);
315 49311 : *pu = shifti(u, -1); *pv = v;
316 49311 : }
317 :
318 : /* D > 0, d = sqrti(D) */
319 : static GEN
320 10440 : quadunit_q(GEN D, GEN d, long *pN)
321 : {
322 10440 : pari_sp av = avma;
323 : GEN p, q, q1;
324 10440 : long first = 1;
325 10440 : p = (Mod2(d) == Mod2(D))? d: subiu(d, 1);
326 10440 : q = gen_2;
327 10440 : q1 = shifti(subii(D, sqri(p)), -1);
328 : for(;;)
329 109570 : {
330 120010 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
331 120010 : p = subii(d, r);
332 126745 : if (!first && equalii(p1, p)) { *pN = 1; return q; } /* even period */
333 116305 : first = 0;
334 116305 : t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
335 116305 : if (equalii(q, t)) { *pN = -1; return q; } /* odd period */
336 109570 : if (gc_needed(av, 2))
337 : {
338 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunitnorm");
339 0 : (void)gc_all(av, 3, &p, &q, &q1);
340 : }
341 : }
342 : }
343 : /* fundamental unit mod N */
344 : static GEN
345 70 : quadunit_mod(GEN D, GEN N)
346 : {
347 70 : GEN q, u, v, d = sqrti(D);
348 70 : pari_sp av = avma;
349 : long s;
350 70 : q = gc_INT(av, quadunit_q(D, d, &s));
351 70 : if (mpodd(N) && equali1(gcdii(q, N)))
352 : {
353 10 : quadunit_uvmod(D, d, N, &u, &v);
354 10 : q = Fp_inv(shifti(q, 1), N);
355 10 : u = Fp_mul(u, q, N);
356 10 : v = Fp_mul(v, q, N); v = modii(shifti(v, 1), N);
357 : }
358 : else
359 : {
360 60 : GEN M = shifti(mulii(q, N), 1);
361 60 : quadunit_uvmod(D, d, M, &u, &v);
362 60 : u = diviiexact(u, q);
363 60 : v = modii(diviiexact(v, q), N); u = shifti(u,-1);
364 : }
365 70 : return deg1pol_shallow(v, u, 0);
366 : }
367 :
368 : /* f \prod_{p|f} [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
369 : static GEN
370 19002 : quadclassnoEuler_fact(GEN D, GEN P, GEN E)
371 : {
372 19002 : long i, l = lg(P);
373 : GEN H;
374 19002 : if (typ(E) != t_VECSMALL) E = vec_to_vecsmall(E);
375 40383 : for (i = 1, H = gen_1; i < l; i++)
376 : {
377 21381 : GEN p = gel(P,i);
378 21381 : long e = E[i], s = kronecker(D,p);
379 21381 : if (!s)
380 6090 : H = mulii(H, e == 1? p: powiu(p, e));
381 : else
382 : {
383 15291 : H = mulii(H, subis(p, s));
384 15291 : if (e >= 2) H = mulii(H, e == 2? p: powiu(p,e-1));
385 : }
386 : }
387 19002 : return H;
388 : }
389 :
390 : /* D > 0; y mod (N,T) congruent to fundamental unit of maximal order and
391 : * disc D. Return unit index of order of conductor N */
392 : static GEN
393 18960 : quadunitindex_ii(GEN D, GEN N, GEN F, GEN y, GEN T)
394 : {
395 18960 : GEN H = quadclassnoEuler_fact(D, gel(F,1), gel(F,2));
396 18960 : GEN P, E, a = Z_smoothen(H, gel(F,1), &P, &E), faH = mkmat2(P, E);
397 : struct uimod S;
398 :
399 18960 : if (a) faH = ZM_merge_factor(Z_factor(a), faH);
400 : /* multiple of unit index, in [H, factor(H)] format */
401 18960 : S.N = N; S.T = FpX_red(T, N);
402 18960 : return gen_order(y, mkvec2(H,faH), (void*)&S, &ui_group);
403 : }
404 : static GEN
405 70 : quadunitindex_i(GEN D, GEN N, GEN F)
406 70 : { return quadunitindex_ii(D, N, F, quadunit_mod(D, N), quadpoly_i(D)); }
407 : GEN
408 80 : quadunitindex(GEN D, GEN N)
409 : {
410 80 : pari_sp av = avma;
411 : long r, s;
412 : GEN F;
413 80 : check_quaddisc(D, &s, &r, "quadunitindex");
414 75 : if ((F = check_arith_pos(N,"quadunitindex")))
415 10 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
416 70 : if (equali1(N)) return gen_1;
417 65 : if (s < 0) switch(itos_or_0(D)) {
418 5 : case -3: return utoipos(3);
419 5 : case -4: return utoipos(2);
420 5 : default: return gen_1;
421 : }
422 50 : return gc_INT(av, quadunitindex_i(D, N, F? F: Z_factor(N)));
423 : }
424 :
425 : /* fundamental unit is u + vx mod quadpoly(D); always called with D
426 : * fundamental but would work in all cases. Same algorithm as basecase,
427 : * except we compute the product of elementary matrices with a product tree */
428 : static void
429 5 : quadunit_uv(GEN D, GEN *pu, GEN *pv)
430 : {
431 5 : GEN a, b, c, u, v, p, q, q1, f, d = sqrtremi(D, &a);
432 5 : pari_sp av = avma;
433 5 : long i = 0;
434 5 : int m = mpodd(D);
435 :
436 5 : p = d; q1 = shifti(a, -1); q = gen_2;
437 5 : if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
438 5 : f = zerovec(2 + (expi(D)>>1));
439 5 : gel(f,1) = mkmat22(p, gen_2, gen_1, gen_0);
440 : for(;;)
441 15257670 : {
442 15257675 : GEN t = q, u1,u2, v1,v2;
443 15257675 : if (!i) { i = 1; q = q1; }
444 : else
445 : {
446 15257670 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
447 15257670 : p = subii(d, r);
448 15257670 : if (equalii(p1, p))
449 : { /* even period */
450 0 : f = prod_fm(f, i, 1); u2 = gel(f,1); v2 = gel(f,2);
451 0 : a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break;
452 : }
453 15257670 : i = update_fm(f, A, i);
454 15257670 : q = submulii(q1, A, subii(p, p1));
455 : }
456 15257675 : q1 = t;
457 15257675 : if (equalii(q, t))
458 : { /* odd period */
459 5 : f = prod_fm(f, i, 0);
460 5 : u2 = gcoeff(f,1,1); u1 = gcoeff(f,1,2); a = mulii(u1, u2);
461 5 : v2 = gcoeff(f,2,1); v1 = gcoeff(f,2,2); b = mulii(v1, v2);
462 5 : c = mulii(addii(u1, v1), addii(u2, v2)); break;
463 : }
464 15257670 : if (gc_needed(av, 2))
465 : {
466 67 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit (%ld)", i);
467 67 : (void)gc_all(av, 4, &p, &f, &q,&q1);
468 : }
469 : }
470 5 : u = diviiexact(addmulii(a, D, b), q);
471 5 : v = diviiexact(subii(c, addii(a, b)), q);
472 5 : if (m == 1) u = subii(u, v);
473 5 : *pu = shifti(u, -1); *pv = v;
474 5 : }
475 : GEN
476 49321 : quadunit(GEN D0)
477 : {
478 49321 : pari_sp av = avma;
479 : GEN P, E, D, u, v;
480 49321 : long s = signe(D0);
481 : /* check_quaddisc_real omitting test for squares */
482 49321 : if (typ(D0) != t_INT) pari_err_TYPE("quadunit", D0);
483 49316 : if (s <= 0) pari_err_DOMAIN("quadunit", "disc","<=",gen_0,D0);
484 49316 : if (mod4(D0) > 1) pari_err_DOMAIN("quadunit","disc % 4",">", gen_1,D0);
485 49316 : D = coredisc2_fact(Z_factor(D0), s, &P, &E);
486 : /* test for squares done here for free */
487 49316 : if (equali1(D)) pari_err_DOMAIN("quadunit","issquare(disc)","=", gen_1,D0);
488 49316 : if (cmpiu(D, 2000000) < 0)
489 49311 : quadunit_uv_basecase(D, &u, &v);
490 : else
491 5 : quadunit_uv(D, &u, &v);
492 49316 : if (lg(P) != 1)
493 : { /* non-trivial conductor N > 1 */
494 18890 : GEN N = factorback2(P,E), qD = quadpoly_i(D);
495 18890 : GEN n, y = deg1pol_shallow(v, u, 0); /* maximal order fund unit */
496 18890 : n = quadunitindex_ii(D, N, mkvec2(P,E), FpX_red(y,N), qD); /* unit index */
497 18890 : y = ZXQ_powu(y, itou(n), qD); /* fund unit of order of conductor N */
498 18890 : v = gel(y,3); u = gel(y,2); /* u + v w_D */
499 18890 : if (mpodd(D))
500 : { /* w_D = (1+sqrt(D))/2 */
501 12395 : if (mpodd(D0))
502 : { /* w_D0 = (1 + N sqrt(D)) / 2 */
503 4405 : GEN v0 = v;
504 4405 : v = diviiexact(v, N);
505 4405 : u = addii(u, shifti(subii(v0, v), -1));
506 : }
507 : else
508 : { /* w_D0 = N sqrt(D)/2, N is even */
509 7990 : v = shifti(v, -1);
510 7990 : u = addii(u, v);
511 7990 : v = diviiexact(v, shifti(N,-1));
512 : }
513 : }
514 : else /* w_D = sqrt(D), w_D0 = N sqrt(D) */
515 6495 : v = diviiexact(v, N);
516 : }
517 49316 : return gc_GEN(av, mkquad(quadpoly_i(D0), u, v));
518 : }
519 : long
520 49310 : quadunitnorm(GEN D)
521 : {
522 49310 : pari_sp av = avma;
523 : long s, r;
524 49310 : check_quaddisc(D, &s, &r, "quadunitnorm");
525 49305 : if (s < 0) return 1;
526 49305 : if (mod4(D) == 1 && BPSW_psp(D)) return -1;
527 43680 : if (Z_has_prime_3mod4(D)) return 1;
528 10370 : (void)quadunit_q(D, sqrti(D), &s); return gc_long(av, s);
529 : }
530 :
531 : GEN
532 38 : quadregulator(GEN x, long prec)
533 : {
534 38 : pari_sp av = avma, av2;
535 : GEN R, rsqd, u, v, sqd;
536 : long r, e;
537 :
538 38 : check_quaddisc_real(x, &r, "quadregulator");
539 38 : sqd = sqrti(x);
540 38 : rsqd = gsqrt(x,prec); av2 = avma;
541 38 : e = 0; R = real2n(1, prec); u = utoi(r); v = gen_2;
542 : for(;;)
543 107 : {
544 145 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
545 145 : GEN v1 = divii(subii(x,sqri(u1)),v);
546 145 : if (equalii(v,v1)) { R = mulrr(sqrr(R), divri(addir(u1,rsqd),v)); break; }
547 119 : if (equalii(u,u1)) { R = sqrr(R); break; }
548 107 : R = mulrr(R, divri(addir(u1,rsqd),v));
549 107 : e += expo(R); setexpo(R,0);
550 107 : u = u1; v = v1;
551 107 : if (e & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
552 107 : if (gc_needed(av2,2))
553 : {
554 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
555 0 : (void)gc_all(av2,3, &R,&u,&v);
556 : }
557 : }
558 38 : R = divri(R, v); e = 2*e - 1;
559 : /* avoid loss of accuracy */
560 38 : if (!((e + expo(R)) & ~EXPOBITS)) { setexpo(R, e + expo(R)); e = 0; }
561 38 : R = logr_abs(R);
562 38 : if (e) R = addrr(R, mulsr(e, mplog2(prec)));
563 38 : return gc_leaf(av, R);
564 : }
565 :
566 : /*************************************************************************/
567 : /** **/
568 : /** CLASS NUMBER **/
569 : /** **/
570 : /*************************************************************************/
571 :
572 : int
573 9227345 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
574 :
575 7995951 : static GEN qfi_pow(void *E, GEN f, GEN n)
576 7995951 : { return E? nupow(f,n,(GEN)E): qfbpow_i(f,n); }
577 5459561 : static GEN qfi_comp(void *E, GEN f, GEN g)
578 5459561 : { return E? nucomp(f,g,(GEN)E): qfbcomp_i(f,g); }
579 : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
580 : gidentical,qfb_equal1,NULL};
581 :
582 : GEN
583 2600370 : qfi_order(GEN q, GEN o)
584 2600370 : { return gen_order(q, o, NULL, &qfi_group); }
585 :
586 : GEN
587 0 : qfi_log(GEN a, GEN g, GEN o)
588 0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
589 :
590 : GEN
591 560712 : qfi_Shanks(GEN a, GEN g, long n)
592 : {
593 560712 : pari_sp av = avma;
594 : GEN T, X;
595 : long rt_n, c;
596 :
597 560712 : a = qfi_red(a);
598 560712 : g = qfi_red(g);
599 :
600 560712 : rt_n = sqrt((double)n);
601 560712 : c = n / rt_n;
602 560712 : c = (c * rt_n < n + 1) ? c + 1 : c;
603 :
604 560712 : T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
605 560712 : X = gen_Shanks(T, a, c, NULL, &qfi_group);
606 560712 : return X? gc_INT(av, X): gc_NULL(av);
607 : }
608 :
609 : GEN
610 334 : qfbclassno0(GEN x,long flag)
611 : {
612 334 : switch(flag)
613 : {
614 322 : case 0: return map_proto_G(classno,x);
615 12 : case 1: return map_proto_G(classno2,x);
616 0 : default: pari_err_FLAG("qfbclassno");
617 : }
618 : return NULL; /* LCOV_EXCL_LINE */
619 : }
620 :
621 : /* f^h = 1, return order(f). Set *pfao to its factorization */
622 : static GEN
623 3876 : find_order(void *E, GEN f, GEN h, GEN *pfao)
624 : {
625 3876 : GEN v = gen_factored_order(f, h, E, &qfi_group);
626 3876 : *pfao = gel(v,2); return gel(v,1);
627 : }
628 :
629 : static int
630 215 : ok_q(GEN q, GEN h, GEN d2, long r2)
631 : {
632 215 : if (d2)
633 : {
634 175 : if (r2 <= 2 && !mpodd(q)) return 0;
635 175 : return is_pm1(Z_ppo(q,d2));
636 : }
637 : else
638 : {
639 40 : if (r2 <= 1 && !mpodd(q)) return 0;
640 40 : return is_pm1(Z_ppo(q,h));
641 : }
642 : }
643 :
644 : /* a,b given by their factorizations. Return factorization of lcm(a,b).
645 : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
646 : static GEN
647 130 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
648 : {
649 130 : GEN P = ZC_union_shallow(gel(Fa,1), gel(Fb,1));
650 130 : GEN A = gen_1, B = gen_1;
651 130 : long i, l = lg(P);
652 130 : GEN E = cgetg(l, t_COL);
653 505 : for (i=1; i<l; i++)
654 : {
655 375 : GEN p = gel(P,i);
656 375 : long va = Z_pval(a,p);
657 375 : long vb = Z_pval(b,p);
658 375 : if (va < vb)
659 : {
660 145 : B = mulii(B,powiu(p,vb));
661 145 : gel(E,i) = utoi(vb);
662 : }
663 : else
664 : {
665 230 : A = mulii(A,powiu(p,va));
666 230 : gel(E,i) = utoi(va);
667 : }
668 : }
669 130 : *pA = A;
670 130 : *pB = B; return mkmat2(P,E);
671 : }
672 :
673 : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
674 : static void
675 130 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
676 : {
677 130 : GEN A,B, g1 = *pg1, d1 = *pd1;
678 130 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
679 130 : *pg1 = gmul(qfbpow_i(g1, diviiexact(d1,A)), qfbpow_i(f, diviiexact(o,B)));
680 130 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
681 130 : }
682 :
683 : /* Let s = 1 or -1; D = s * d; assume Df^2 fits in an ulong
684 : * Return f / [O_{Df^2}^*:O_D^*] * \prod_{p|f} [ 1 - (D/p) p^-1 ]
685 : * The Euler product is \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
686 : ulong
687 5441 : uquadclassnoF_fact(ulong d, long s, GEN P, GEN E)
688 : {
689 5441 : long i, l = lg(P);
690 5441 : ulong H = 1;
691 6940 : for (i = 1; i < l; i++)
692 : {
693 1499 : ulong p = P[i], e = E[i];
694 1499 : long D = (long)(p == 2? d & 7: d % p), a;
695 1499 : if (s < 0) D = -D;
696 1499 : a = kross(D,p);
697 1499 : if (!a)
698 393 : H *= upowuu(p, e);
699 : else
700 : {
701 1106 : H *= p - a;
702 1106 : if (e >= 2) H *= upowuu(p, e-1);
703 : }
704 : }
705 5441 : if (l == 1) return H;
706 1133 : if (s < 0)
707 : {
708 1119 : switch(d)
709 : { /* divide by [ O_K^* : O^* ] */
710 80 : case 4: H >>= 1; break;
711 215 : case 3: H /= 3; break;
712 : }
713 : }
714 : else
715 : {
716 14 : GEN fa = mkmat2(zc_to_ZC(P), zc_to_ZC(E));
717 14 : H /= itou(quadunitindex_i(utoipos(d), factorback(fa), fa));
718 : }
719 1133 : return H;
720 : }
721 : GEN
722 42 : quadclassnoF_fact(GEN D, GEN P, GEN E)
723 : {
724 42 : GEN H = quadclassnoEuler_fact(D, P, E);
725 42 : if (lg(P) == 1) return H;
726 11 : if (signe(D) < 0)
727 : {
728 5 : switch(itou_or_0(D))
729 : { /* divide by [ O_K^* : O^* ] */
730 0 : case 4: H = shifti(H,-1); break;
731 0 : case 3: H = diviuexact(H,3); break;
732 : }
733 : }
734 : else
735 : {
736 6 : GEN fa = mkvec2(P, E);
737 6 : H = diviiexact(H, quadunitindex_i(D, factorback2(P, E), fa));
738 : }
739 11 : return H;
740 : }
741 :
742 : static ulong
743 482 : quadclassnoF_u(ulong x, long s, ulong *pD)
744 : {
745 482 : pari_sp av = avma;
746 : GEN P, E;
747 482 : ulong D = coredisc2u_fact(factoru(x), s, &P, &E);
748 482 : long H = uquadclassnoF_fact(D, s, P, E);
749 482 : *pD = D; return gc_ulong(av, H);
750 : }
751 : ulong
752 0 : unegquadclassnoF(ulong x, ulong *pD) { return quadclassnoF_u(x, -1, pD); }
753 : ulong
754 0 : uposquadclassnoF(ulong x, ulong *pD) { return quadclassnoF_u(x, 1, pD); }
755 :
756 : /* *pD = coredisc(x), *pR = regulator (x > 0) or NULL */
757 : GEN
758 524 : quadclassnoF(GEN x, GEN *pD)
759 : {
760 : GEN D, P, E;
761 524 : if (lgefint(x) == 3)
762 : {
763 482 : long s = signe(x);
764 482 : ulong d, h = quadclassnoF_u(x[2], s, &d);
765 482 : if (pD) *pD = s > 0? utoipos(d): utoineg(d);
766 482 : return utoipos(h);
767 : }
768 42 : D = coredisc2_fact(absZ_factor(x), signe(x), &P, &E);
769 42 : if (pD) *pD = D;
770 42 : return quadclassnoF_fact(D, P, E);
771 : }
772 :
773 : static long
774 486 : two_rank(GEN x)
775 : {
776 486 : GEN p = gel(absZ_factor(x),1);
777 486 : long l = lg(p)-1;
778 : #if 0 /* positive disc not needed */
779 : if (signe(x) > 0)
780 : {
781 : long i;
782 : for (i=1; i<=l; i++)
783 : if (mod4(gel(p,i)) == 3) { l--; break; }
784 : }
785 : #endif
786 486 : return l-1;
787 : }
788 :
789 : static GEN
790 9209 : sqr_primeform(GEN x, ulong p) { return qfbsqr_i(primeform_u(x, p)); }
791 : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
792 : static GEN
793 486 : get_forms(GEN D, GEN *pL)
794 : {
795 486 : const long MAXFORM = 20;
796 486 : GEN L, sqrtD = gsqrt(absi_shallow(D),DEFAULTPREC);
797 486 : GEN forms = vectrunc_init(MAXFORM+1);
798 486 : long s, nforms = 0;
799 : ulong p;
800 : forprime_t S;
801 486 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
802 486 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
803 486 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
804 486 : if (s < 10) s = 200;
805 318 : else if (s < 20) s = 1000;
806 296 : else if (s < 5000) s = 5000;
807 486 : u_forprime_init(&S, 2, s);
808 9625664 : while ( (p = u_forprime_next(&S)) )
809 : {
810 9625178 : long d, k = kroiu(D,p);
811 : pari_sp av2;
812 9625178 : if (!k) continue;
813 9624329 : if (k > 0)
814 : {
815 4813059 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
816 4813059 : d = p-1;
817 : }
818 : else
819 4811270 : d = p+1;
820 9624329 : av2 = avma; affrr(divru(mulur(p,L),d), L); set_avma(av2);
821 : }
822 486 : *pL = L; return forms;
823 : }
824 :
825 : /* h ~ #G, return o = order of f, set fao = its factorization */
826 : static GEN
827 556 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
828 : {
829 556 : long s = minss(itos(sqrti(h)), 10000);
830 556 : GEN T = gen_Shanks_init(f, s, E, &qfi_group);
831 556 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
832 556 : return find_order(E, f, addiu(v,1), pfao);
833 : }
834 :
835 : /* if g = 1 in G/<f> ? */
836 : static int
837 4060 : equal1(void *E, GEN T, ulong N, GEN g)
838 4060 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
839 :
840 : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
841 : * FIXME: should be gen_order, but equal1 has the wrong prototype */
842 : static GEN
843 250 : relative_order(void *E, GEN a, GEN o, ulong N, GEN T)
844 : {
845 250 : pari_sp av = avma;
846 : long i, l;
847 : GEN m;
848 :
849 250 : m = get_arith_ZZM(o);
850 250 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
851 250 : o = gel(m,1);
852 250 : m = gel(m,2); l = lgcols(m);
853 820 : for (i = l-1; i; i--)
854 : {
855 570 : GEN t, y, p = gcoeff(m,i,1);
856 570 : long j, e = itos(gcoeff(m,i,2));
857 570 : if (l == 2) {
858 35 : t = gen_1;
859 35 : y = a;
860 : } else {
861 535 : t = diviiexact(o, powiu(p,e));
862 535 : y = powgi(a, t);
863 : }
864 570 : if (equal1(E, T,N,y)) o = t;
865 : else {
866 260 : for (j = 1; j < e; j++)
867 : {
868 60 : y = powgi(y, p);
869 60 : if (equal1(E, T,N,y)) break;
870 : }
871 255 : if (j < e) {
872 55 : if (j > 1) p = powiu(p, j);
873 55 : o = mulii(t, p);
874 : }
875 : }
876 : }
877 250 : return gc_GEN(av, o);
878 : }
879 :
880 : /* h(x) for x<0 using Baby Step/Giant Step.
881 : * Assumes G is not too far from being cyclic.
882 : *
883 : * Compute G^2 instead of G so as to kill most of the noncyclicity */
884 : GEN
885 598 : classno(GEN x)
886 : {
887 598 : pari_sp av = avma;
888 : long r2, k, s, i, l;
889 : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
890 : void *E;
891 :
892 598 : if (signe(x) >= 0) return classno2(x);
893 :
894 572 : check_quaddisc(x, &s, &k, "classno");
895 572 : if (abscmpiu(x,12) <= 0) return gen_1;
896 :
897 486 : Hf = quadclassnoF(x, &D);
898 486 : if (abscmpiu(D,12) <= 0) return gc_GEN(av, Hf);
899 486 : forms = get_forms(D, &L);
900 486 : r2 = two_rank(D);
901 486 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
902 :
903 486 : l = lg(forms);
904 486 : order_bound = const_vec(l-1, NULL);
905 486 : E = expi(D) > 60? (void*)sqrtnint(shifti(absi_shallow(D),-2),4): NULL;
906 486 : g1 = gel(forms,1);
907 486 : gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
908 486 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
909 486 : d2 = NULL; /* not computed yet */
910 486 : if (is_pm1(q)) goto END;
911 5225 : for (i=2; i < l; i++)
912 : {
913 4965 : GEN o, fao, a, F, fd, f = gel(forms,i);
914 4965 : fd = qfbpow_i(f, d1); if (is_pm1(gel(fd,1))) continue;
915 130 : F = qfbpow_i(fd, q);
916 130 : a = gel(F,1);
917 130 : o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
918 : /* f^(d1 q) = 1 */
919 130 : fao = ZM_merge_factor(fad1,fao);
920 130 : o = find_order(E, f, fao, &fao);
921 130 : gel(order_bound,i) = o;
922 : /* o = order of f, fao = factor(o) */
923 130 : update_g1(&g1,&d1,&fad1, f,o,fao);
924 130 : q = diviiround(hin, d1);
925 130 : if (is_pm1(q)) goto END;
926 : }
927 : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
928 260 : if (expi(q) > 3)
929 : { /* q large: compute d2, 2nd elt divisor */
930 220 : ulong N, n = 2*itou(sqrti(d1));
931 220 : GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
932 220 : d2 = gen_1;
933 220 : N = itou( gceil(gdivgu(d1,n)) ); /* order(g1) <= n*N */
934 3605 : for (i = 1; i < l; i++)
935 : {
936 3430 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
937 3430 : if (!B) B = find_order(E, f, fad1, /*junk*/&d);
938 3430 : f = qfbpow_i(f,d2);
939 3430 : if (equal1(E, T,N,f)) continue;
940 250 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
941 : /* f^B = 1 */
942 250 : d = relative_order(E, f, B, N,T);
943 250 : d2= mulii(d,d2);
944 250 : D = mulii(d1,d2);
945 250 : q = diviiround(hin,D);
946 250 : if (is_pm1(q)) { d1 = D; goto END; }
947 : }
948 : /* very probably, d2 is the 2nd elementary divisor */
949 175 : d1 = D; /* product of first two elt divisors */
950 : }
951 : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
952 : * 2-rank */
953 215 : if (!ok_q(q,d1,d2,r2))
954 : {
955 0 : GEN q0 = q;
956 : long d;
957 0 : if (cmpii(mulii(q,d1), hin) < 0)
958 : { /* try q = q0+1,-1,+2,-2 */
959 0 : d = 1;
960 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
961 : }
962 : else
963 : { /* q0-1,+1,-2,+2 */
964 0 : d = -1;
965 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
966 : }
967 : }
968 215 : d1 = mulii(d1,q);
969 :
970 486 : END:
971 486 : return gc_INT(av, shifti(mulii(d1,Hf), r2));
972 : }
973 :
974 : /* use Euler products */
975 : GEN
976 38 : classno2(GEN x)
977 : {
978 38 : pari_sp av = avma;
979 38 : const long prec = DEFAULTPREC;
980 : long n, i, s;
981 38 : GEN p1, p2, S, p4, p5, p7, Hf, Pi, logd, sqrtd, D, half, reg = NULL;
982 :
983 38 : check_quaddisc(x, &s, NULL, "classno2");
984 38 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
985 :
986 38 : Hf = quadclassnoF(x, &D);
987 38 : if (s < 0 && abscmpiu(D,12) <= 0) return gc_GEN(av, Hf); /* |D| < 12*/
988 :
989 38 : Pi = mppi(prec);
990 38 : sqrtd = sqrtr_abs(itor(D, prec));
991 38 : logd = logr_abs(sqrtd); shiftr_inplace(logd,1);
992 38 : p1 = sqrtr_abs(divrr(mulir(D,logd), gmul2n(Pi,1)));
993 38 : if (s > 0)
994 : {
995 32 : GEN invlogd = invr(logd);
996 32 : reg = quadregulator(D, prec);
997 32 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
998 32 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
999 : }
1000 38 : n = itos_or_0( mptrunc(p1) );
1001 38 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
1002 :
1003 38 : p4 = divri(Pi, D); setsigne(p4, 1);
1004 38 : p7 = invr(sqrtr_abs(Pi));
1005 38 : half = real2n(-1, prec);
1006 38 : if (s > 0)
1007 : { /* i = 1, shortcut */
1008 32 : p1 = sqrtd;
1009 32 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
1010 32 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
1011 1048 : for (i=2; i<=n; i++)
1012 : {
1013 1016 : long k = kroiu(D,i); if (!k) continue;
1014 932 : p2 = mulir(sqru(i), p4);
1015 932 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
1016 932 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
1017 932 : S = (k>0)? addrr(S,p5): subrr(S,p5);
1018 : }
1019 32 : S = shiftr(divrr(S,reg),-1);
1020 : }
1021 : else
1022 : { /* i = 1, shortcut */
1023 6 : p1 = gdiv(sqrtd, Pi);
1024 6 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
1025 6 : S = addrr(p5, divrr(p1, mpexp(p4)));
1026 816 : for (i=2; i<=n; i++)
1027 : {
1028 810 : long k = kroiu(D,i); if (!k) continue;
1029 810 : p2 = mulir(sqru(i), p4);
1030 810 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
1031 810 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
1032 810 : S = (k>0)? addrr(S,p5): subrr(S,p5);
1033 : }
1034 : }
1035 38 : return gc_INT(av, mulii(Hf, roundr(S)));
1036 : }
1037 :
1038 : /* 1 + q + ... + q^v, v > 0 */
1039 : static GEN
1040 695 : geomsumu(ulong q, long v)
1041 : {
1042 695 : GEN u = utoipos(1+q);
1043 905 : for (; v > 1; v--) u = addui(1, mului(q, u));
1044 695 : return u;
1045 : }
1046 : static GEN
1047 695 : geomsum(GEN q, long v)
1048 : {
1049 : GEN u;
1050 695 : if (lgefint(q) == 3) return geomsumu(q[2], v);
1051 0 : u = addiu(q,1);
1052 0 : for (; v > 1; v--) u = addui(1, mulii(q, u));
1053 0 : return u;
1054 : }
1055 :
1056 : /* 1+p+...+p^(e-1), e >= 1; assuming result fits in an ulong */
1057 : static ulong
1058 7811 : usumpow(ulong p, long e)
1059 : {
1060 : ulong q;
1061 : long i;
1062 7811 : if (p == 2) return (1UL << e) - 1; /* also OK if e = BITS_IN_LONG */
1063 5521 : e--; for (i = 1, q = p + 1; i < e; i++) q = p*q + 1;
1064 5480 : return q;
1065 : }
1066 : long
1067 128720 : uhclassnoF_fact(GEN faF, long D)
1068 : {
1069 128720 : GEN P = gel(faF,1), E = gel(faF,2);
1070 128720 : long i, t, l = lg(P);
1071 281144 : for (i = t = 1; i < l; i++)
1072 : {
1073 152424 : long p = P[i], e = E[i], s = kross(D,p);
1074 152424 : if (e == 1) { t *= 1 + p - s; continue; }
1075 40627 : if (s == 1) { t *= upowuu(p,e); continue; }
1076 7811 : t *= 1 + usumpow(p, e) * (p - s);
1077 : }
1078 128720 : return t;
1079 : }
1080 : /* Hurwitz(D F^2)/ Hurwitz(D)
1081 : * = \sum_{f|F} f \prod_{p|f} (1-kro(D/p)/p)
1082 : * = \prod_{p^e || F} (1 + (p^e-1) / (p-1) * (p-kro(D/p))) */
1083 : GEN
1084 43346 : hclassnoF_fact(GEN P, GEN E, GEN D)
1085 : {
1086 : GEN H;
1087 43346 : long i, l = lg(P);
1088 43346 : if (l == 1) return gen_1;
1089 40090 : for (i = 1, H = NULL; i < l; i++)
1090 : {
1091 21585 : GEN t, p = gel(P,i);
1092 21585 : long e = E[i], s = kronecker(D,p);
1093 21585 : if (e == 1) t = addiu(p, 1-s);
1094 1110 : else if (s == 1) t = powiu(p, e);
1095 695 : else t = addui(1, mulii(subis(p, s), geomsum(p, e-1)));
1096 21585 : H = H? mulii(H,t): t;
1097 : }
1098 18505 : return H;
1099 : }
1100 : static GEN
1101 43346 : hclassno6_large(GEN x)
1102 : {
1103 43346 : GEN H = NULL, P, E, D = coredisc2_fact(absZ_factor(x), -1, &P, &E);
1104 43346 : long l = lg(P);
1105 :
1106 43346 : if (l > 1 && lgefint(x) == 3)
1107 : { /* F != 1, second chance */
1108 18504 : ulong h = hclassno6u_no_cache(x[2]);
1109 18504 : if (h) H = utoipos(h);
1110 : }
1111 43346 : if (!H)
1112 : {
1113 43346 : H = quadclassno(D);
1114 43346 : switch(itou_or_0(D))
1115 : {
1116 35 : case 3: H = shifti(H,1);break;
1117 5 : case 4: H = muliu(H,3); break;
1118 43306 : default:H = muliu(H,6); break;
1119 : }
1120 : }
1121 43346 : return mulii(H, hclassnoF_fact(P, E, D));
1122 : }
1123 :
1124 : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
1125 : GEN
1126 94636 : hclassno6(GEN x)
1127 : {
1128 94636 : ulong d = itou_or_0(x);
1129 94636 : if (d)
1130 : { /* create cache if d small */
1131 94635 : ulong h = d < 500000 ? hclassno6u(d): hclassno6u_no_cache(d);
1132 94635 : if (h) return utoipos(h);
1133 : }
1134 43346 : return hclassno6_large(x);
1135 : }
1136 :
1137 : GEN
1138 41 : hclassno(GEN x)
1139 : {
1140 : long a, s;
1141 41 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
1142 41 : s = signe(x);
1143 41 : if (s < 0) return gen_0;
1144 41 : if (!s) return gdivgs(gen_1, -12);
1145 41 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
1146 41 : return gdivgu(hclassno6(x), 6);
1147 : }
1148 :
1149 : /* return [N',v]; v contains all x mod N' s.t. x^2 + B x + C = 0 modulo N */
1150 : GEN
1151 1827403 : Zn_quad_roots(GEN N, GEN B, GEN C)
1152 : {
1153 1827403 : pari_sp av = avma;
1154 : GEN fa, D, w, v, P, E, F0, Q0, F, mF, A, Q, T, R, Np, N4;
1155 : long l, i, j, ct;
1156 :
1157 1827403 : if ((fa = check_arith_non0(N,"Zn_quad_roots")))
1158 : {
1159 5628 : N = typ(N) == t_VEC? gel(N,1): factorback(N);
1160 5628 : fa = clean_Z_factor(fa);
1161 : }
1162 1827403 : N = absi_shallow(N);
1163 1827403 : N4 = shifti(N,2);
1164 1827403 : D = modii(subii(sqri(B), shifti(C,2)), N4);
1165 1827403 : if (!signe(D))
1166 : { /* (x + B/2)^2 = 0 (mod N), D = B^2-4C = 0 (4N) => B even */
1167 495 : if (!fa) fa = Z_factor(N);
1168 495 : P = gel(fa,1);
1169 495 : E = ZV_to_zv(gel(fa,2));
1170 495 : l = lg(P);
1171 1085 : for (i = 1; i < l; i++) E[i] = (E[i]+1) >> 1;
1172 495 : Np = factorback2(P, E); /* x = -B mod N' */
1173 495 : B = shifti(B,-1);
1174 495 : return gc_GEN(av, mkvec2(Np, mkvec(Fp_neg(B,Np))));
1175 : }
1176 1826908 : if (!fa)
1177 1821400 : fa = Z_factor(N4);
1178 : else /* convert to factorization of N4 = 4*N */
1179 5508 : fa = famat_reduce(famat_mulpows_shallow(fa, gen_2, 2));
1180 1826908 : P = gel(fa,1); l = lg(P);
1181 1826908 : E = ZV_to_zv(gel(fa,2));
1182 1826908 : F = cgetg(l, t_VEC);
1183 1826908 : mF= cgetg(l, t_VEC); F0 = gen_0;
1184 1826908 : Q = cgetg(l, t_VEC); Q0 = gen_1;
1185 4383375 : for (i = j = 1, ct = 0; i < l; i++)
1186 : {
1187 4014607 : GEN p = gel(P,i), q, f, mf, D0;
1188 4014607 : long t2, s = E[i], t = Z_pvalrem(D, p, &D0), d = s - t;
1189 4014607 : if (d <= 0)
1190 : {
1191 894942 : q = powiu(p, (s+1)>>1);
1192 1545300 : Q0 = mulii(Q0, q); continue;
1193 : }
1194 : /* d > 0 */
1195 4445287 : if (odd(t)) return NULL;
1196 2987147 : t2 = t >> 1;
1197 2987147 : if (i > 1)
1198 : { /* p > 2 */
1199 1871837 : if (kronecker(D0, p) == -1) return NULL;
1200 904842 : q = powiu(p, s - t2);
1201 904842 : f = Zp_sqrt(D0, p, d);
1202 904842 : if (!f) return NULL; /* p was not actually prime... */
1203 904842 : if (t2) f = mulii(powiu(p,t2), f);
1204 904842 : mf = Fp_neg(f, q);
1205 : }
1206 : else
1207 : { /* p = 2 */
1208 1115310 : if (d <= 3)
1209 : {
1210 846989 : if (d == 3 && Mod8(D0) != 1) return NULL;
1211 686724 : if (d == 2 && Mod4(D0) != 1) return NULL;
1212 650358 : Q0 = int2n(1+t2); F0 = NULL; continue;
1213 : }
1214 268321 : if (Mod8(D0) != 1) return NULL;
1215 106325 : q = int2n(d - 1 + t2);
1216 106325 : f = Z2_sqrt(D0, d);
1217 106325 : if (t2) f = shifti(f, t2);
1218 106325 : mf = Fp_neg(f, q);
1219 : }
1220 1011167 : gel(Q,j) = q;
1221 1011167 : gel(F,j) = f;
1222 1011167 : gel(mF,j)= mf; j++;
1223 : }
1224 368768 : setlg(Q,j);
1225 368768 : setlg(F,j);
1226 368768 : setlg(mF,j);
1227 368768 : if (is_pm1(Q0)) A = leafcopy(F);
1228 : else
1229 : { /* append the fixed congruence (F0 mod Q0) */
1230 340937 : if (!F0) F0 = shifti(Q0,-1);
1231 340937 : A = shallowconcat(F, F0);
1232 340937 : Q = shallowconcat(Q, Q0);
1233 : }
1234 368768 : if (j-1 >= LGnumBITS) pari_err_OVERFLOW("Zn_quad_roots");
1235 368768 : ct = 1L << (j-1);
1236 368768 : T = ZV_producttree(Q);
1237 368768 : R = ZV_chinesetree(Q,T);
1238 368768 : Np = gmael(T, lg(T)-1, 1);
1239 368768 : B = modii(B, Np);
1240 368768 : if (!signe(B)) B = NULL;
1241 368768 : Np = shifti(Np, -1); /* N' = (\prod_i Q[i]) / 2 */
1242 368768 : w = cgetg(3, t_VEC);
1243 368768 : gel(w,1) = icopy(Np);
1244 368768 : gel(w,2) = v = cgetg(ct+1, t_VEC);
1245 368768 : l = lg(F);
1246 1737613 : for (j = 1; j <= ct; j++)
1247 : {
1248 1368845 : pari_sp av2 = avma;
1249 1368845 : long m = j - 1;
1250 : GEN u;
1251 4248709 : for (i = 1; i < l; i++)
1252 : {
1253 2879864 : gel(A,i) = (m&1L)? gel(mF,i): gel(F,i);
1254 2879864 : m >>= 1;
1255 : }
1256 1368845 : u = ZV_chinese_tree(A,Q,T,R); /* u mod N' st u^2 = B^2-4C modulo 4N */
1257 1368845 : if (B) u = subii(u,B);
1258 1368845 : gel(v,j) = gc_INT(av2, modii(shifti(u,-1), Np));
1259 : }
1260 368768 : return gc_upto(av, w);
1261 : }
1262 :
1263 : GEN
1264 0 : Zn_sqrtall(GEN D, GEN fa)
1265 : {
1266 0 : pari_sp av = avma;
1267 : long i, j, k, lB, aN;
1268 : GEN a, L, V, B, N;
1269 0 : a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
1270 0 : V = Zn_quad_roots(fa, gen_0, negi(D));
1271 0 : if (!V) return NULL;
1272 0 : N = gel(V,1); B = gel(V,2); lB = lg(B);
1273 0 : aN = itou(diviiexact(a, N)); /* |a|/N */
1274 0 : L = cgetg((lB-1)*aN+1, t_VEC);
1275 0 : for (k = 1, i = 1; i < lB; i++)
1276 : {
1277 0 : GEN b = gel(B,i);
1278 0 : for (j = 0;; b = addii(b, N))
1279 : {
1280 0 : gel(L, k++) = b;
1281 0 : if (++j == aN) break;
1282 : }
1283 : }
1284 0 : return gc_upto(av, ZV_sort(L));
1285 : }
|
