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 22890 : isfundamental(GEN x)
57 : {
58 22890 : 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 21483 : return Z_isfundamental(x);
65 : }
66 :
67 : /* x fundamental ? */
68 : long
69 16953 : uposisfundamental(ulong x)
70 : {
71 16953 : ulong r = x & 15; /* x mod 16 */
72 16953 : if (!r) return 0;
73 16141 : switch(r & 3)
74 : { /* x mod 4 */
75 3487 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
76 6203 : case 1: return uissquarefree(x);
77 6451 : default: return 0;
78 : }
79 : }
80 : /* -x fundamental ? */
81 : long
82 31500 : unegisfundamental(ulong x)
83 : {
84 31500 : ulong r = x & 15; /* x mod 16 */
85 31500 : if (!r) return 0;
86 29855 : switch(r & 3)
87 : { /* x mod 4 */
88 7855 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
89 11671 : case 3: return uissquarefree(x);
90 10329 : default: return 0;
91 : }
92 : }
93 : long
94 25109 : sisfundamental(long x)
95 25109 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
96 :
97 : long
98 22050 : Z_isfundamental(GEN x)
99 : {
100 : long r;
101 22050 : switch(lgefint(x))
102 : {
103 7 : case 2: return 0;
104 11396 : case 3: return signe(x) < 0? unegisfundamental(x[2])
105 31429 : : 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 2079 : fa_quaddisc(GEN f)
126 : {
127 2079 : GEN P = gel(f,1), E = gel(f,2), s = gen_1;
128 2079 : long i, l = lg(P);
129 6881 : for (i = 1; i < l; i++) /* possibly including -1 */
130 4802 : if (mpodd(gel(E,i))) s = mulii(s, gel(P,i));
131 2079 : if (Mod4(s) > 1) s = shifti(s,2);
132 2079 : return s;
133 : }
134 :
135 : GEN
136 2933 : quaddisc(GEN x)
137 : {
138 2933 : const pari_sp av = avma;
139 2933 : long tx = typ(x);
140 : GEN F;
141 2933 : 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 2079 : 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 21199714 : update_f(GEN f, GEN u)
164 : {
165 21199714 : GEN a = gcoeff(f,1,1), b = gcoeff(f,1,2);
166 21199714 : GEN c = gcoeff(f,2,1), d = gcoeff(f,2,2);
167 21199714 : gcoeff(f,1,1) = addmulii(b, u,a); gcoeff(f,1,2) = a;
168 21199714 : gcoeff(f,2,1) = addmulii(d, u,c); gcoeff(f,2,2) = c;
169 21199714 : }
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 21360738 : update_fm(GEN f, GEN a, long i)
176 : {
177 : #ifdef LONG_IS_64BIT
178 18309204 : const long LIM = 10;
179 : #else
180 3051534 : const long LIM = 18;
181 : #endif
182 21360738 : pari_sp av = avma;
183 : long k, v;
184 : GEN u;
185 21360738 : if (!odd(i)) { gel(f,1) = a; return i+1; }
186 21280226 : u = gel(f, 1);
187 21280226 : if (typ(u) == t_INT) /* [u,1;1,0] * [a,1;1,0] */
188 80512 : { gel(f,1) = mkmat22(addiu(mulii(a, u), 1), u, a, gen_1); return i; }
189 21199714 : update_f(u, a); if (lgefint(gcoeff(u,1,1)) < LIM) return i;
190 80512 : v = vals(i+1); gel(f,1) = gen_0;
191 160970 : for (k = 1; k < v; k++) { u = ZM2_mul(gel(f,k+1), u); gel(f,k+1) = gen_0; }
192 80512 : if (v != 1) u = gerepileupto(av, u);
193 80512 : gel(f,v+1) = u; return i+1;
194 : }
195 : /* \prod f[j]; if first only return column 1 */
196 : static GEN
197 7 : prod_fm(GEN f, long i, long first)
198 : {
199 7 : long k, v = vals(i) + 1;
200 7 : GEN u = gel(f, v);
201 : /* i a power of 2: f[1] can't be a t_INT */
202 7 : if (!(i >>= v)) return first? gel(u,1): u;
203 105 : for (k = v+1; i; i >>= 1, k++)
204 98 : if (odd(i))
205 : {
206 54 : GEN w = gel(f,k);
207 54 : 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 54 : default: /* t_MAT */
214 54 : u = first? ZM_ZC_mul(w, gel(u,1)): ZM2_mul(w, u); break;
215 : }
216 : }
217 7 : return u;
218 : }
219 :
220 : GEN
221 69048 : quadunit0(GEN x, long v)
222 : {
223 69048 : GEN y = quadunit(x);
224 69041 : if (v==-1) v = fetch_user_var("w");
225 69041 : setvarn(gel(y,1), v); return y;
226 : }
227 :
228 : struct uimod { GEN N, T; };
229 : static GEN
230 21378 : ui_pow(void *E, GEN x, GEN n)
231 21378 : { struct uimod *S = (struct uimod*)E; return FpXQ_pow(x, n, S->T, S->N); }
232 : static int
233 38612 : 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 98 : 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 98 : int m = mpodd(D), first = 1;
242 98 : pari_sp av = avma;
243 98 : p = (mpodd(d) == m)? d: subiu(d, 1);
244 98 : u1 = negi(p); u2 = gen_2;
245 98 : v1 = gen_1; v2 = gen_0; q = gen_2;
246 98 : q1 = shifti(subii(D, sqri(p)), -1);
247 : for(;;)
248 308 : {
249 406 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
250 406 : p = subii(d, r);
251 406 : if (equalii(p1, p) && !first)
252 : { /* even period */
253 14 : u = addmulii(sqri(u2), D, sqri(v2));
254 14 : v = shifti(mulii(u2,v2), 1);
255 14 : break;
256 : }
257 392 : first = 0;
258 392 : t = Fp_addmul(u1, A, u2, N); u1 = u2; u2 = t;
259 392 : t = Fp_addmul(v1, A, v2, N); v1 = v2; v2 = t;
260 392 : t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
261 392 : if (equalii(q, t))
262 : { /* odd period */
263 84 : u = addmulii(mulii(u1,u2), D, mulii(v1,v2));
264 84 : v = addmulii(mulii(u1,v2), u2, v1);
265 84 : break;
266 : }
267 308 : if (gc_needed(av, 2))
268 : {
269 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit_uvmod");
270 0 : gerepileall(av, 7, &p, &u1,&u2,&v1,&v2, &q,&q1);
271 : }
272 : }
273 98 : *pu = modii(u, N);
274 98 : *pv = modii(v, N); if (m) *pu = Fp_sub(*pu, *pv, N);
275 98 : }
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 69034 : quadunit_uv_basecase(GEN D, GEN *pu, GEN *pv)
282 : {
283 69034 : GEN u1, u2, v1, v2, p, q, q1, u, v, a, b, c, d = sqrtremi(D, &a);
284 69034 : int m = mpodd(D);
285 69034 : long first = 1;
286 :
287 69034 : p = d; q1 = shifti(a, -1); q = gen_2;
288 69034 : if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
289 69034 : u1 = gen_2; u2 = p;
290 69034 : v1 = gen_0; v2 = gen_1;
291 : for(;;)
292 871052 : {
293 940086 : GEN t = q;
294 940086 : if (first) { first = 0; q = q1; }
295 : else
296 : {
297 871052 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
298 871052 : p = subii(d, r);
299 871052 : if (equalii(p1, p)) /* even period */
300 45738 : { a = sqri(u2); b = sqri(v2); c = sqri(addii(u2, v2)); break; }
301 825314 : r = addmulii(u1, A, u2); u1 = u2; u2 = r;
302 825314 : r = addmulii(v1, A, v2); v1 = v2; v2 = r;
303 825314 : q = submulii(q1, A, subii(p, p1));
304 : }
305 894348 : q1 = t;
306 894348 : if (equalii(q, t))
307 : { /* odd period */
308 23296 : a = mulii(u1, u2); b = mulii(v1, v2);
309 23296 : c = mulii(addii(u1, v1), addii(u2, v2)); break;
310 : }
311 : }
312 69034 : u = diviiexact(addmulii(a, D, b), q);
313 69034 : v = diviiexact(subii(c, addii(a, b)), q);
314 69034 : if (m == 1) u = subii(u, v);
315 69034 : *pu = shifti(u, -1); *pv = v;
316 69034 : }
317 :
318 : /* D > 0, d = sqrti(D) */
319 : static GEN
320 69118 : quadunit_q(GEN D, GEN d, long *pN)
321 : {
322 69118 : pari_sp av = avma;
323 : GEN p, q, q1;
324 69118 : long first = 1;
325 69118 : p = (Mod2(d) == Mod2(D))? d: subiu(d, 1);
326 69118 : q = gen_2;
327 69118 : q1 = shifti(subii(D, sqri(p)), -1);
328 : for(;;)
329 1001133 : {
330 1070251 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
331 1070251 : p = subii(d, r);
332 1087555 : if (!first && equalii(p1, p)) { *pN = 1; return q; } /* even period */
333 1018437 : first = 0;
334 1018437 : t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
335 1018437 : if (equalii(q, t)) { *pN = -1; return q; } /* odd period */
336 1001133 : if (gc_needed(av, 2))
337 : {
338 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunitnorm");
339 0 : gerepileall(av, 3, &p, &q, &q1);
340 : }
341 : }
342 : }
343 : /* fundamental unit mod N */
344 : static GEN
345 98 : quadunit_mod(GEN D, GEN N)
346 : {
347 98 : GEN q, u, v, d = sqrti(D);
348 98 : pari_sp av = avma;
349 : long s;
350 98 : q = gerepileuptoint(av, quadunit_q(D, d, &s));
351 98 : if (mpodd(N) && equali1(gcdii(q, N)))
352 : {
353 14 : quadunit_uvmod(D, d, N, &u, &v);
354 14 : q = Fp_inv(shifti(q, 1), N);
355 14 : u = Fp_mul(u, q, N);
356 14 : v = Fp_mul(v, q, N); v = modii(shifti(v, 1), N);
357 : }
358 : else
359 : {
360 84 : GEN M = shifti(mulii(q, N), 1);
361 84 : quadunit_uvmod(D, d, M, &u, &v);
362 84 : u = diviiexact(u, q);
363 84 : v = modii(diviiexact(v, q), N); u = shifti(u,-1);
364 : }
365 98 : 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 26592 : quadclassnoEuler_fact(GEN D, GEN P, GEN E)
371 : {
372 26592 : long i, l = lg(P);
373 : GEN H;
374 26592 : if (typ(E) != t_VECSMALL) E = vec_to_vecsmall(E);
375 56525 : for (i = 1, H = gen_1; i < l; i++)
376 : {
377 29933 : GEN p = gel(P,i);
378 29933 : long e = E[i], s = kronecker(D,p);
379 29933 : if (!s)
380 8526 : H = mulii(H, e == 1? p: powiu(p, e));
381 : else
382 : {
383 21407 : H = mulii(H, subis(p, s));
384 21407 : if (e >= 2) H = mulii(H, e == 2? p: powiu(p,e-1));
385 : }
386 : }
387 26592 : 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 26544 : quadunitindex_ii(GEN D, GEN N, GEN F, GEN y, GEN T)
394 : {
395 26544 : GEN H = quadclassnoEuler_fact(D, gel(F,1), gel(F,2));
396 26544 : GEN P, E, a = Z_smoothen(H, gel(F,1), &P, &E), faH = mkmat2(P, E);
397 : struct uimod S;
398 :
399 26544 : if (a) faH = ZM_merge_factor(Z_factor(a), faH);
400 : /* multiple of unit index, in [H, factor(H)] format */
401 26544 : S.N = N; S.T = FpX_red(T, N);
402 26544 : return gen_order(y, mkvec2(H,faH), (void*)&S, &ui_group);
403 : }
404 : static GEN
405 98 : quadunitindex_i(GEN D, GEN N, GEN F)
406 98 : { return quadunitindex_ii(D, N, F, quadunit_mod(D, N), quadpoly_i(D)); }
407 : GEN
408 112 : quadunitindex(GEN D, GEN N)
409 : {
410 112 : pari_sp av = avma;
411 : long r, s;
412 : GEN F;
413 112 : check_quaddisc(D, &s, &r, "quadunitindex");
414 105 : if ((F = check_arith_pos(N,"quadunitindex")))
415 14 : N = typ(N) == t_VEC? gel(N,1): factorback(F);
416 98 : if (equali1(N)) return gen_1;
417 91 : if (s < 0) switch(itos_or_0(D)) {
418 7 : case -3: return utoipos(3);
419 7 : case -4: return utoipos(2);
420 7 : default: return gen_1;
421 : }
422 70 : return gerepileuptoint(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 7 : quadunit_uv(GEN D, GEN *pu, GEN *pv)
430 : {
431 7 : GEN a, b, c, u, v, p, q, q1, f, d = sqrtremi(D, &a);
432 7 : pari_sp av = avma;
433 7 : long i = 0;
434 7 : int m = mpodd(D);
435 :
436 7 : p = d; q1 = shifti(a, -1); q = gen_2;
437 7 : if (mpodd(d) != m) { p = subiu(d,1); q1 = addii(q1,d); } /* q1 = (D-p^2)/2 */
438 7 : f = zerovec(2 + (expi(D)>>1));
439 7 : gel(f,1) = mkmat22(p, gen_2, gen_1, gen_0);
440 : for(;;)
441 21360738 : {
442 21360745 : GEN t = q, u1,u2, v1,v2;
443 21360745 : if (!i) { i = 1; q = q1; }
444 : else
445 : {
446 21360738 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p;
447 21360738 : p = subii(d, r);
448 21360738 : 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 21360738 : i = update_fm(f, A, i);
454 21360738 : q = submulii(q1, A, subii(p, p1));
455 : }
456 21360745 : q1 = t;
457 21360745 : if (equalii(q, t))
458 : { /* odd period */
459 7 : f = prod_fm(f, i, 0);
460 7 : u2 = gcoeff(f,1,1); u1 = gcoeff(f,1,2); a = mulii(u1, u2);
461 7 : v2 = gcoeff(f,2,1); v1 = gcoeff(f,2,2); b = mulii(v1, v2);
462 7 : c = mulii(addii(u1, v1), addii(u2, v2)); break;
463 : }
464 21360738 : if (gc_needed(av, 2))
465 : {
466 96 : if(DEBUGMEM>1) pari_warn(warnmem,"quadunit (%ld)", i);
467 96 : gerepileall(av, 4, &p, &f, &q,&q1);
468 : }
469 : }
470 7 : u = diviiexact(addmulii(a, D, b), q);
471 7 : v = diviiexact(subii(c, addii(a, b)), q);
472 7 : if (m == 1) u = subii(u, v);
473 7 : *pu = shifti(u, -1); *pv = v;
474 7 : }
475 : GEN
476 69048 : quadunit(GEN D0)
477 : {
478 69048 : pari_sp av = avma;
479 : GEN P, E, D, u, v;
480 69048 : long s = signe(D0);
481 : /* check_quaddisc_real omitting test for squares */
482 69048 : if (typ(D0) != t_INT) pari_err_TYPE("quadunit", D0);
483 69041 : if (s <= 0) pari_err_DOMAIN("quadunit", "disc","<=",gen_0,D0);
484 69041 : if (mod4(D0) > 1) pari_err_DOMAIN("quadunit","disc % 4",">", gen_1,D0);
485 69041 : D = coredisc2_fact(Z_factor(D0), s, &P, &E);
486 : /* test for squares done here for free */
487 69041 : if (equali1(D)) pari_err_DOMAIN("quadunit","issquare(disc)","=", gen_1,D0);
488 69041 : if (cmpiu(D, 2000000) < 0)
489 69034 : quadunit_uv_basecase(D, &u, &v);
490 : else
491 7 : quadunit_uv(D, &u, &v);
492 69041 : if (lg(P) != 1)
493 : { /* non-trivial conductor N > 1 */
494 26446 : GEN N = factorback2(P,E), qD = quadpoly_i(D);
495 26446 : GEN n, y = deg1pol_shallow(v, u, 0); /* maximal order fund unit */
496 26446 : n = quadunitindex_ii(D, N, mkvec2(P,E), FpX_red(y,N), qD); /* unit index */
497 26446 : y = ZXQ_powu(y, itou(n), qD); /* fund unit of order of conductor N */
498 26446 : v = gel(y,3); u = gel(y,2); /* u + v w_D */
499 26446 : if (mpodd(D))
500 : { /* w_D = (1+sqrt(D))/2 */
501 17353 : if (mpodd(D0))
502 : { /* w_D0 = (1 + N sqrt(D)) / 2 */
503 6167 : GEN v0 = v;
504 6167 : v = diviiexact(v, N);
505 6167 : u = addii(u, shifti(subii(v0, v), -1));
506 : }
507 : else
508 : { /* w_D0 = N sqrt(D)/2, N is even */
509 11186 : v = shifti(v, -1);
510 11186 : u = addii(u, v);
511 11186 : v = diviiexact(v, shifti(N,-1));
512 : }
513 : }
514 : else /* w_D = sqrt(D), w_D0 = N sqrt(D) */
515 9093 : v = diviiexact(v, N);
516 : }
517 69041 : return gerepilecopy(av, mkquad(quadpoly_i(D0), u, v));
518 : }
519 : long
520 69027 : quadunitnorm(GEN D)
521 : {
522 69027 : pari_sp av = avma;
523 : long s, r;
524 69027 : check_quaddisc(D, &s, &r, "quadunitnorm");
525 69020 : if (s < 0) return 1;
526 69020 : (void)quadunit_q(D, sqrti(D), &s); return gc_long(av, s);
527 : }
528 :
529 : GEN
530 49 : quadregulator(GEN x, long prec)
531 : {
532 49 : pari_sp av = avma, av2;
533 : GEN R, rsqd, u, v, sqd;
534 : long r, e;
535 :
536 49 : check_quaddisc_real(x, &r, "quadregulator");
537 49 : sqd = sqrti(x);
538 49 : rsqd = gsqrt(x,prec); av2 = avma;
539 49 : e = 0; R = real2n(1, prec); u = utoi(r); v = gen_2;
540 : for(;;)
541 140 : {
542 189 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
543 189 : GEN v1 = divii(subii(x,sqri(u1)),v);
544 189 : if (equalii(v,v1)) { R = mulrr(sqrr(R), divri(addir(u1,rsqd),v)); break; }
545 154 : if (equalii(u,u1)) { R = sqrr(R); break; }
546 140 : R = mulrr(R, divri(addir(u1,rsqd),v));
547 140 : e += expo(R); setexpo(R,0);
548 140 : u = u1; v = v1;
549 140 : if (e & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
550 140 : if (gc_needed(av2,2))
551 : {
552 0 : if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
553 0 : gerepileall(av2,3, &R,&u,&v);
554 : }
555 : }
556 49 : R = divri(R, v); e = 2*e - 1;
557 : /* avoid loss of accuracy */
558 49 : if (!((e + expo(R)) & ~EXPOBITS)) { setexpo(R, e + expo(R)); e = 0; }
559 49 : R = logr_abs(R);
560 49 : if (e) R = addrr(R, mulsr(e, mplog2(prec)));
561 49 : return gerepileuptoleaf(av, R);
562 : }
563 :
564 : /*************************************************************************/
565 : /** **/
566 : /** CLASS NUMBER **/
567 : /** **/
568 : /*************************************************************************/
569 :
570 : int
571 9204105 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
572 :
573 10288260 : static GEN qfi_pow(void *E, GEN f, GEN n)
574 10288260 : { return E? nupow(f,n,(GEN)E): qfbpow_i(f,n); }
575 7015563 : static GEN qfi_comp(void *E, GEN f, GEN g)
576 7015563 : { return E? nucomp(f,g,(GEN)E): qfbcomp_i(f,g); }
577 : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
578 : gidentical,qfb_equal1,NULL};
579 :
580 : GEN
581 3370157 : qfi_order(GEN q, GEN o)
582 3370157 : { return gen_order(q, o, NULL, &qfi_group); }
583 :
584 : GEN
585 0 : qfi_log(GEN a, GEN g, GEN o)
586 0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
587 :
588 : GEN
589 722080 : qfi_Shanks(GEN a, GEN g, long n)
590 : {
591 722080 : pari_sp av = avma;
592 : GEN T, X;
593 : long rt_n, c;
594 :
595 722080 : a = qfbred_i(a);
596 722080 : g = qfbred_i(g);
597 :
598 722080 : rt_n = sqrt((double)n);
599 722080 : c = n / rt_n;
600 722080 : c = (c * rt_n < n + 1) ? c + 1 : c;
601 :
602 722080 : T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
603 722080 : X = gen_Shanks(T, a, c, NULL, &qfi_group);
604 722080 : return X? gerepileuptoint(av, X): gc_NULL(av);
605 : }
606 :
607 : GEN
608 462 : qfbclassno0(GEN x,long flag)
609 : {
610 462 : switch(flag)
611 : {
612 448 : case 0: return map_proto_G(classno,x);
613 14 : case 1: return map_proto_G(classno2,x);
614 0 : default: pari_err_FLAG("qfbclassno");
615 : }
616 : return NULL; /* LCOV_EXCL_LINE */
617 : }
618 :
619 : /* f^h = 1, return order(f). Set *pfao to its factorization */
620 : static GEN
621 5397 : find_order(void *E, GEN f, GEN h, GEN *pfao)
622 : {
623 5397 : GEN v = gen_factored_order(f, h, E, &qfi_group);
624 5397 : *pfao = gel(v,2); return gel(v,1);
625 : }
626 :
627 : static int
628 301 : ok_q(GEN q, GEN h, GEN d2, long r2)
629 : {
630 301 : if (d2)
631 : {
632 245 : if (r2 <= 2 && !mpodd(q)) return 0;
633 245 : return is_pm1(Z_ppo(q,d2));
634 : }
635 : else
636 : {
637 56 : if (r2 <= 1 && !mpodd(q)) return 0;
638 56 : return is_pm1(Z_ppo(q,h));
639 : }
640 : }
641 :
642 : /* a,b given by their factorizations. Return factorization of lcm(a,b).
643 : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
644 : static GEN
645 182 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
646 : {
647 182 : GEN P = ZC_union_shallow(gel(Fa,1), gel(Fb,1));
648 182 : GEN A = gen_1, B = gen_1;
649 182 : long i, l = lg(P);
650 182 : GEN E = cgetg(l, t_COL);
651 707 : for (i=1; i<l; i++)
652 : {
653 525 : GEN p = gel(P,i);
654 525 : long va = Z_pval(a,p);
655 525 : long vb = Z_pval(b,p);
656 525 : if (va < vb)
657 : {
658 203 : B = mulii(B,powiu(p,vb));
659 203 : gel(E,i) = utoi(vb);
660 : }
661 : else
662 : {
663 322 : A = mulii(A,powiu(p,va));
664 322 : gel(E,i) = utoi(va);
665 : }
666 : }
667 182 : *pA = A;
668 182 : *pB = B; return mkmat2(P,E);
669 : }
670 :
671 : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
672 : static void
673 182 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
674 : {
675 182 : GEN A,B, g1 = *pg1, d1 = *pd1;
676 182 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
677 182 : *pg1 = gmul(qfbpow_i(g1, diviiexact(d1,A)), qfbpow_i(f, diviiexact(o,B)));
678 182 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
679 182 : }
680 :
681 : /* Let s = 1 or -1; D = s * d; assume Df^2 fits in an ulong
682 : * Return f / [O_{Df^2}^*:O_D^*] * \prod_{p|f} [ 1 - (D/p) p^-1 ]
683 : * The Euler product is \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
684 : ulong
685 7081 : uquadclassnoF_fact(ulong d, long s, GEN P, GEN E)
686 : {
687 7081 : long i, l = lg(P);
688 7081 : ulong H = 1;
689 9054 : for (i = 1; i < l; i++)
690 : {
691 1973 : ulong p = P[i], e = E[i];
692 1973 : long D = (long)(p == 2? d & 7: d % p), a;
693 1973 : if (s < 0) D = -D;
694 1973 : a = kross(D,p);
695 1973 : if (!a)
696 504 : H *= upowuu(p, e);
697 : else
698 : {
699 1469 : H *= p - a;
700 1469 : if (e >= 2) H *= upowuu(p, e-1);
701 : }
702 : }
703 7081 : if (l == 1) return H;
704 1490 : if (s < 0)
705 : {
706 1470 : switch(d)
707 : { /* divide by [ O_K^* : O^* ] */
708 105 : case 4: H >>= 1; break;
709 287 : case 3: H /= 3; break;
710 : }
711 : }
712 : else
713 : {
714 20 : GEN fa = mkmat2(zc_to_ZC(P), zc_to_ZC(E));
715 20 : H /= itou(quadunitindex_i(utoipos(d), factorback(fa), fa));
716 : }
717 1490 : return H;
718 : }
719 : GEN
720 48 : quadclassnoF_fact(GEN D, GEN P, GEN E)
721 : {
722 48 : GEN H = quadclassnoEuler_fact(D, P, E);
723 48 : if (lg(P) == 1) return H;
724 15 : if (signe(D) < 0)
725 : {
726 7 : switch(itou_or_0(D))
727 : { /* divide by [ O_K^* : O^* ] */
728 0 : case 4: H = shifti(H,-1); break;
729 0 : case 3: H = diviuexact(H,3); break;
730 : }
731 : }
732 : else
733 : {
734 8 : GEN fa = mkvec2(P, E);
735 8 : H = diviiexact(H, quadunitindex_i(D, factorback2(P, E), fa));
736 : }
737 15 : return H;
738 : }
739 :
740 : static ulong
741 652 : quadclassnoF_u(ulong x, long s, ulong *pD)
742 : {
743 652 : pari_sp av = avma;
744 : GEN P, E;
745 652 : ulong D = coredisc2u_fact(factoru(x), s, &P, &E);
746 652 : long H = uquadclassnoF_fact(D, s, P, E);
747 652 : *pD = D; return gc_ulong(av, H);
748 : }
749 : ulong
750 0 : unegquadclassnoF(ulong x, ulong *pD) { return quadclassnoF_u(x, -1, pD); }
751 : ulong
752 0 : uposquadclassnoF(ulong x, ulong *pD) { return quadclassnoF_u(x, 1, pD); }
753 :
754 : /* *pD = coredisc(x), *pR = regulator (x > 0) or NULL */
755 : GEN
756 700 : quadclassnoF(GEN x, GEN *pD)
757 : {
758 : GEN D, P, E;
759 700 : if (lgefint(x) == 3)
760 : {
761 652 : long s = signe(x);
762 652 : ulong d, h = quadclassnoF_u(x[2], s, &d);
763 652 : if (pD) *pD = s > 0? utoipos(d): utoineg(d);
764 652 : return utoipos(h);
765 : }
766 48 : D = coredisc2_fact(absZ_factor(x), signe(x), &P, &E);
767 48 : if (pD) *pD = D;
768 48 : return quadclassnoF_fact(D, P, E);
769 : }
770 :
771 : static long
772 651 : two_rank(GEN x)
773 : {
774 651 : GEN p = gel(absZ_factor(x),1);
775 651 : long l = lg(p)-1;
776 : #if 0 /* positive disc not needed */
777 : if (signe(x) > 0)
778 : {
779 : long i;
780 : for (i=1; i<=l; i++)
781 : if (mod4(gel(p,i)) == 3) { l--; break; }
782 : }
783 : #endif
784 651 : return l-1;
785 : }
786 :
787 : static GEN
788 12334 : sqr_primeform(GEN x, ulong p) { return qfbsqr_i(primeform_u(x, p)); }
789 : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
790 : static GEN
791 651 : get_forms(GEN D, GEN *pL)
792 : {
793 651 : const long MAXFORM = 20;
794 651 : GEN L, sqrtD = gsqrt(absi_shallow(D),DEFAULTPREC);
795 651 : GEN forms = vectrunc_init(MAXFORM+1);
796 651 : long s, nforms = 0;
797 : ulong p;
798 : forprime_t S;
799 651 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
800 651 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
801 651 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
802 651 : if (s < 10) s = 200;
803 427 : else if (s < 20) s = 1000;
804 413 : else if (s < 5000) s = 5000;
805 651 : u_forprime_init(&S, 2, s);
806 13471626 : while ( (p = u_forprime_next(&S)) )
807 : {
808 13470975 : long d, k = kroiu(D,p);
809 : pari_sp av2;
810 13470975 : if (!k) continue;
811 13469841 : if (k > 0)
812 : {
813 6736219 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
814 6736219 : d = p-1;
815 : }
816 : else
817 6733622 : d = p+1;
818 13469841 : av2 = avma; affrr(divru(mulur(p,L),d), L); set_avma(av2);
819 : }
820 651 : *pL = L; return forms;
821 : }
822 :
823 : /* h ~ #G, return o = order of f, set fao = its factorization */
824 : static GEN
825 749 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
826 : {
827 749 : long s = minss(itos(sqrti(h)), 10000);
828 749 : GEN T = gen_Shanks_init(f, s, E, &qfi_group);
829 749 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
830 749 : return find_order(E, f, addiu(v,1), pfao);
831 : }
832 :
833 : /* if g = 1 in G/<f> ? */
834 : static int
835 5684 : equal1(void *E, GEN T, ulong N, GEN g)
836 5684 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
837 :
838 : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
839 : * FIXME: should be gen_order, but equal1 has the wrong prototype */
840 : static GEN
841 350 : relative_order(void *E, GEN a, GEN o, ulong N, GEN T)
842 : {
843 350 : pari_sp av = avma;
844 : long i, l;
845 : GEN m;
846 :
847 350 : m = get_arith_ZZM(o);
848 350 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
849 350 : o = gel(m,1);
850 350 : m = gel(m,2); l = lgcols(m);
851 1148 : for (i = l-1; i; i--)
852 : {
853 798 : GEN t, y, p = gcoeff(m,i,1);
854 798 : long j, e = itos(gcoeff(m,i,2));
855 798 : if (l == 2) {
856 49 : t = gen_1;
857 49 : y = a;
858 : } else {
859 749 : t = diviiexact(o, powiu(p,e));
860 749 : y = powgi(a, t);
861 : }
862 798 : if (equal1(E, T,N,y)) o = t;
863 : else {
864 364 : for (j = 1; j < e; j++)
865 : {
866 84 : y = powgi(y, p);
867 84 : if (equal1(E, T,N,y)) break;
868 : }
869 357 : if (j < e) {
870 77 : if (j > 1) p = powiu(p, j);
871 77 : o = mulii(t, p);
872 : }
873 : }
874 : }
875 350 : return gerepilecopy(av, o);
876 : }
877 :
878 : /* h(x) for x<0 using Baby Step/Giant Step.
879 : * Assumes G is not too far from being cyclic.
880 : *
881 : * Compute G^2 instead of G so as to kill most of the noncyclicity */
882 : GEN
883 798 : classno(GEN x)
884 : {
885 798 : pari_sp av = avma;
886 : long r2, k, s, i, l;
887 : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
888 : void *E;
889 :
890 798 : if (signe(x) >= 0) return classno2(x);
891 :
892 763 : check_quaddisc(x, &s, &k, "classno");
893 763 : if (abscmpiu(x,12) <= 0) return gen_1;
894 :
895 651 : Hf = quadclassnoF(x, &D);
896 651 : if (abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
897 651 : forms = get_forms(D, &L);
898 651 : r2 = two_rank(D);
899 651 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
900 :
901 651 : l = lg(forms);
902 651 : order_bound = const_vec(l-1, NULL);
903 651 : E = expi(D) > 60? (void*)sqrtnint(shifti(absi_shallow(D),-2),4): NULL;
904 651 : g1 = gel(forms,1);
905 651 : gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
906 651 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
907 651 : d2 = NULL; /* not computed yet */
908 651 : if (is_pm1(q)) goto END;
909 7315 : for (i=2; i < l; i++)
910 : {
911 6951 : GEN o, fao, a, F, fd, f = gel(forms,i);
912 6951 : fd = qfbpow_i(f, d1); if (is_pm1(gel(fd,1))) continue;
913 182 : F = qfbpow_i(fd, q);
914 182 : a = gel(F,1);
915 182 : o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
916 : /* f^(d1 q) = 1 */
917 182 : fao = ZM_merge_factor(fad1,fao);
918 182 : o = find_order(E, f, fao, &fao);
919 182 : gel(order_bound,i) = o;
920 : /* o = order of f, fao = factor(o) */
921 182 : update_g1(&g1,&d1,&fad1, f,o,fao);
922 182 : q = diviiround(hin, d1);
923 182 : if (is_pm1(q)) goto END;
924 : }
925 : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
926 364 : if (expi(q) > 3)
927 : { /* q large: compute d2, 2nd elt divisor */
928 308 : ulong N, n = 2*itou(sqrti(d1));
929 308 : GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
930 308 : d2 = gen_1;
931 308 : N = itou( gceil(gdivgu(d1,n)) ); /* order(g1) <= n*N */
932 5047 : for (i = 1; i < l; i++)
933 : {
934 4802 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
935 4802 : if (!B) B = find_order(E, f, fad1, /*junk*/&d);
936 4802 : f = qfbpow_i(f,d2);
937 4802 : if (equal1(E, T,N,f)) continue;
938 350 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
939 : /* f^B = 1 */
940 350 : d = relative_order(E, f, B, N,T);
941 350 : d2= mulii(d,d2);
942 350 : D = mulii(d1,d2);
943 350 : q = diviiround(hin,D);
944 350 : if (is_pm1(q)) { d1 = D; goto END; }
945 : }
946 : /* very probably, d2 is the 2nd elementary divisor */
947 245 : d1 = D; /* product of first two elt divisors */
948 : }
949 : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
950 : * 2-rank */
951 301 : if (!ok_q(q,d1,d2,r2))
952 : {
953 0 : GEN q0 = q;
954 : long d;
955 0 : if (cmpii(mulii(q,d1), hin) < 0)
956 : { /* try q = q0+1,-1,+2,-2 */
957 0 : d = 1;
958 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
959 : }
960 : else
961 : { /* q0-1,+1,-2,+2 */
962 0 : d = -1;
963 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
964 : }
965 : }
966 301 : d1 = mulii(d1,q);
967 :
968 651 : END:
969 651 : return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
970 : }
971 :
972 : /* use Euler products */
973 : GEN
974 49 : classno2(GEN x)
975 : {
976 49 : pari_sp av = avma;
977 49 : const long prec = DEFAULTPREC;
978 : long n, i, s;
979 49 : GEN p1, p2, S, p4, p5, p7, Hf, Pi, logd, sqrtd, D, half, reg = NULL;
980 :
981 49 : check_quaddisc(x, &s, NULL, "classno2");
982 49 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
983 :
984 49 : Hf = quadclassnoF(x, &D);
985 49 : if (s < 0 && abscmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
986 :
987 49 : Pi = mppi(prec);
988 49 : sqrtd = sqrtr_abs(itor(D, prec));
989 49 : logd = logr_abs(sqrtd); shiftr_inplace(logd,1);
990 49 : p1 = sqrtr_abs(divrr(mulir(D,logd), gmul2n(Pi,1)));
991 49 : if (s > 0)
992 : {
993 42 : GEN invlogd = invr(logd);
994 42 : reg = quadregulator(D, prec);
995 42 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
996 42 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
997 : }
998 49 : n = itos_or_0( mptrunc(p1) );
999 49 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
1000 :
1001 49 : p4 = divri(Pi, D); setsigne(p4, 1);
1002 49 : p7 = invr(sqrtr_abs(Pi));
1003 49 : half = real2n(-1, prec);
1004 49 : if (s > 0)
1005 : { /* i = 1, shortcut */
1006 42 : p1 = sqrtd;
1007 42 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
1008 42 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
1009 1358 : for (i=2; i<=n; i++)
1010 : {
1011 1316 : long k = kroiu(D,i); if (!k) continue;
1012 1218 : p2 = mulir(sqru(i), p4);
1013 1218 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
1014 1218 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
1015 1218 : S = (k>0)? addrr(S,p5): subrr(S,p5);
1016 : }
1017 42 : S = shiftr(divrr(S,reg),-1);
1018 : }
1019 : else
1020 : { /* i = 1, shortcut */
1021 7 : p1 = gdiv(sqrtd, Pi);
1022 7 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
1023 7 : S = addrr(p5, divrr(p1, mpexp(p4)));
1024 952 : for (i=2; i<=n; i++)
1025 : {
1026 945 : long k = kroiu(D,i); if (!k) continue;
1027 945 : p2 = mulir(sqru(i), p4);
1028 945 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
1029 945 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
1030 945 : S = (k>0)? addrr(S,p5): subrr(S,p5);
1031 : }
1032 : }
1033 49 : return gerepileuptoint(av, mulii(Hf, roundr(S)));
1034 : }
1035 :
1036 : /* 1 + q + ... + q^v, v > 0 */
1037 : static GEN
1038 973 : geomsumu(ulong q, long v)
1039 : {
1040 973 : GEN u = utoipos(1+q);
1041 1267 : for (; v > 1; v--) u = addui(1, mului(q, u));
1042 973 : return u;
1043 : }
1044 : static GEN
1045 973 : geomsum(GEN q, long v)
1046 : {
1047 : GEN u;
1048 973 : if (lgefint(q) == 3) return geomsumu(q[2], v);
1049 0 : u = addiu(q,1);
1050 0 : for (; v > 1; v--) u = addui(1, mulii(q, u));
1051 0 : return u;
1052 : }
1053 :
1054 : /* 1+p+...+p^(e-1), e >= 1; assuming result fits in an ulong */
1055 : static ulong
1056 9426 : usumpow(ulong p, long e)
1057 : {
1058 : ulong q;
1059 : long i;
1060 9426 : if (p == 2) return (1UL << e) - 1; /* also OK if e = BITS_IN_LONG */
1061 6643 : e--; for (i = 1, q = p + 1; i < e; i++) q = p*q + 1;
1062 6587 : return q;
1063 : }
1064 : long
1065 163465 : uhclassnoF_fact(GEN faF, long D)
1066 : {
1067 163465 : GEN P = gel(faF,1), E = gel(faF,2);
1068 163465 : long i, t, l = lg(P);
1069 357282 : for (i = t = 1; i < l; i++)
1070 : {
1071 193817 : long p = P[i], e = E[i], s = kross(D,p);
1072 193817 : if (e == 1) { t *= 1 + p - s; continue; }
1073 51718 : if (s == 1) { t *= upowuu(p,e); continue; }
1074 9426 : t *= 1 + usumpow(p, e) * (p - s);
1075 : }
1076 163465 : return t;
1077 : }
1078 : /* Hurwitz(D F^2)/ Hurwitz(D)
1079 : * = \sum_{f|F} f \prod_{p|f} (1-kro(D/p)/p)
1080 : * = \prod_{p^e || F} (1 + (p^e-1) / (p-1) * (p-kro(D/p))) */
1081 : GEN
1082 60721 : hclassnoF_fact(GEN P, GEN E, GEN D)
1083 : {
1084 : GEN H;
1085 60721 : long i, l = lg(P);
1086 60721 : if (l == 1) return gen_1;
1087 56164 : for (i = 1, H = NULL; i < l; i++)
1088 : {
1089 30240 : GEN t, p = gel(P,i);
1090 30240 : long e = E[i], s = kronecker(D,p);
1091 30242 : if (e == 1) t = addiu(p, 1-s);
1092 1555 : else if (s == 1) t = powiu(p, e);
1093 973 : else t = addui(1, mulii(subis(p, s), geomsum(p, e-1)));
1094 30241 : H = H? mulii(H,t): t;
1095 : }
1096 25924 : return H;
1097 : }
1098 : static GEN
1099 60718 : hclassno6_large(GEN x)
1100 : {
1101 60718 : GEN H = NULL, P, E, D = coredisc2_fact(absZ_factor(x), -1, &P, &E);
1102 60722 : long l = lg(P);
1103 :
1104 60722 : if (l > 1 && lgefint(x) == 3)
1105 : { /* F != 1, second chance */
1106 25927 : ulong h = hclassno6u_no_cache(x[2]);
1107 25927 : if (h) H = utoipos(h);
1108 : }
1109 60722 : if (!H)
1110 : {
1111 60722 : H = quadclassno(D);
1112 60726 : switch(itou_or_0(D))
1113 : {
1114 49 : case 3: H = shifti(H,1);break;
1115 7 : case 4: H = muliu(H,3); break;
1116 60670 : default:H = muliu(H,6); break;
1117 : }
1118 : }
1119 60722 : return mulii(H, hclassnoF_fact(P, E, D));
1120 : }
1121 :
1122 : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
1123 : GEN
1124 132473 : hclassno6(GEN x)
1125 : {
1126 132473 : ulong d = itou_or_0(x);
1127 132473 : if (d)
1128 : { /* create cache if d small */
1129 132472 : ulong h = d < 500000 ? hclassno6u(d): hclassno6u_no_cache(d);
1130 132472 : if (h) return utoipos(h);
1131 : }
1132 60718 : return hclassno6_large(x);
1133 : }
1134 :
1135 : GEN
1136 49 : hclassno(GEN x)
1137 : {
1138 : long a, s;
1139 49 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
1140 49 : s = signe(x);
1141 49 : if (s < 0) return gen_0;
1142 49 : if (!s) return gdivgs(gen_1, -12);
1143 49 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
1144 49 : return gdivgu(hclassno6(x), 6);
1145 : }
1146 :
1147 : /* return [N',v]; v contains all x mod N' s.t. x^2 + B x + C = 0 modulo N */
1148 : GEN
1149 2551766 : Zn_quad_roots(GEN N, GEN B, GEN C)
1150 : {
1151 2551766 : pari_sp av = avma;
1152 2551766 : GEN fa = NULL, D, w, v, P, E, F0, Q0, F, mF, A, Q, T, R, Np, N4;
1153 : long l, i, j, ct;
1154 :
1155 2551766 : if ((fa = check_arith_non0(N,"Zn_quad_roots")))
1156 : {
1157 7665 : N = typ(N) == t_VEC? gel(N,1): factorback(N);
1158 7665 : fa = clean_Z_factor(fa);
1159 : }
1160 2551766 : N = absi_shallow(N);
1161 2551766 : N4 = shifti(N,2);
1162 2551766 : D = modii(subii(sqri(B), shifti(C,2)), N4);
1163 2551766 : if (!signe(D))
1164 : { /* (x + B/2)^2 = 0 (mod N), D = B^2-4C = 0 (4N) => B even */
1165 693 : if (!fa) fa = Z_factor(N);
1166 693 : P = gel(fa,1);
1167 693 : E = ZV_to_zv(gel(fa,2));
1168 693 : l = lg(P);
1169 1519 : for (i = 1; i < l; i++) E[i] = (E[i]+1) >> 1;
1170 693 : Np = factorback2(P, E); /* x = -B mod N' */
1171 693 : B = shifti(B,-1);
1172 693 : return gerepilecopy(av, mkvec2(Np, mkvec(Fp_neg(B,Np))));
1173 : }
1174 2551073 : if (!fa)
1175 2543576 : fa = Z_factor(N4);
1176 : else /* convert to factorization of N4 = 4*N */
1177 7497 : fa = famat_reduce(famat_mulpows_shallow(fa, gen_2, 2));
1178 2551073 : P = gel(fa,1); l = lg(P);
1179 2551073 : E = ZV_to_zv(gel(fa,2));
1180 2551073 : F = cgetg(l, t_VEC);
1181 2551073 : mF= cgetg(l, t_VEC); F0 = gen_0;
1182 2551073 : Q = cgetg(l, t_VEC); Q0 = gen_1;
1183 6107395 : for (i = j = 1, ct = 0; i < l; i++)
1184 : {
1185 5596122 : GEN p = gel(P,i), q, f, mf, D0;
1186 5596122 : long t2, s = E[i], t = Z_pvalrem(D, p, &D0), d = s - t;
1187 5596124 : if (d <= 0)
1188 : {
1189 1252293 : q = powiu(p, (s+1)>>1);
1190 2162306 : Q0 = mulii(Q0, q); continue;
1191 : }
1192 : /* d > 0 */
1193 6198131 : if (odd(t)) return NULL;
1194 4158332 : t2 = t >> 1;
1195 4158332 : if (i > 1)
1196 : { /* p > 2 */
1197 2603156 : if (kronecker(D0, p) == -1) return NULL;
1198 1250879 : q = powiu(p, s - t2);
1199 1250879 : f = Zp_sqrt(D0, p, d);
1200 1250879 : if (!f) return NULL; /* p was not actually prime... */
1201 1250879 : if (t2) f = mulii(powiu(p,t2), f);
1202 1250879 : mf = Fp_neg(f, q);
1203 : }
1204 : else
1205 : { /* p = 2 */
1206 1555176 : if (d <= 3)
1207 : {
1208 1185247 : if (d == 3 && Mod8(D0) != 1) return NULL;
1209 960918 : if (d == 2 && Mod4(D0) != 1) return NULL;
1210 910014 : Q0 = int2n(1+t2); F0 = NULL; continue;
1211 : }
1212 369929 : if (Mod8(D0) != 1) return NULL;
1213 143143 : q = int2n(d - 1 + t2);
1214 143143 : f = Z2_sqrt(D0, d);
1215 143143 : if (t2) f = shifti(f, t2);
1216 143143 : mf = Fp_neg(f, q);
1217 : }
1218 1394020 : gel(Q,j) = q;
1219 1394020 : gel(F,j) = f;
1220 1394020 : gel(mF,j)= mf; j++;
1221 : }
1222 511273 : setlg(Q,j);
1223 511273 : setlg(F,j);
1224 511273 : setlg(mF,j);
1225 511273 : if (is_pm1(Q0)) A = leafcopy(F);
1226 : else
1227 : { /* append the fixed congruence (F0 mod Q0) */
1228 476350 : if (!F0) F0 = shifti(Q0,-1);
1229 476350 : A = shallowconcat(F, F0);
1230 476350 : Q = shallowconcat(Q, Q0);
1231 : }
1232 511273 : ct = 1 << (j-1);
1233 511273 : T = ZV_producttree(Q);
1234 511273 : R = ZV_chinesetree(Q,T);
1235 511272 : Np = gmael(T, lg(T)-1, 1);
1236 511272 : B = modii(B, Np);
1237 511272 : if (!signe(B)) B = NULL;
1238 511272 : Np = shifti(Np, -1); /* N' = (\prod_i Q[i]) / 2 */
1239 511272 : w = cgetg(3, t_VEC);
1240 511271 : gel(w,1) = icopy(Np);
1241 511271 : gel(w,2) = v = cgetg(ct+1, t_VEC);
1242 511273 : l = lg(F);
1243 2359177 : for (j = 1; j <= ct; j++)
1244 : {
1245 1847904 : pari_sp av2 = avma;
1246 1847904 : long m = j - 1;
1247 : GEN u;
1248 5602261 : for (i = 1; i < l; i++)
1249 : {
1250 3754357 : gel(A,i) = (m&1L)? gel(mF,i): gel(F,i);
1251 3754357 : m >>= 1;
1252 : }
1253 1847904 : u = ZV_chinese_tree(A,Q,T,R); /* u mod N' st u^2 = B^2-4C modulo 4N */
1254 1847901 : if (B) u = subii(u,B);
1255 1847901 : gel(v,j) = gerepileuptoint(av2, modii(shifti(u,-1), Np));
1256 : }
1257 511273 : return gerepileupto(av, w);
1258 : }
|
