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 22939 : isfundamental(GEN x)
57 : {
58 22939 : 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 21532 : return Z_isfundamental(x);
65 : }
66 :
67 : /* x fundamental ? */
68 : long
69 16981 : uposisfundamental(ulong x)
70 : {
71 16981 : ulong r = x & 15; /* x mod 16 */
72 16981 : if (!r) return 0;
73 16169 : switch(r & 3)
74 : { /* x mod 4 */
75 3501 : case 0: return (r == 4)? 0: uissquarefree(x >> 2);
76 6217 : case 1: return uissquarefree(x);
77 6451 : default: return 0;
78 : }
79 : }
80 : /* -x fundamental ? */
81 : long
82 36841 : unegisfundamental(ulong x)
83 : {
84 36841 : ulong r = x & 15; /* x mod 16 */
85 36841 : if (!r) return 0;
86 34279 : switch(r & 3)
87 : { /* x mod 4 */
88 10669 : case 0: return (r == 12)? 0: uissquarefree(x >> 2);
89 13281 : 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 22099 : Z_isfundamental(GEN x)
99 : {
100 : long r;
101 22099 : switch(lgefint(x))
102 : {
103 7 : case 2: return 0;
104 11417 : case 3: return signe(x) < 0? unegisfundamental(x[2])
105 31499 : : 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 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 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 = gc_upto(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 20923 : ui_pow(void *E, GEN x, GEN n)
231 20923 : { struct uimod *S = (struct uimod*)E; return FpXQ_pow(x, n, S->T, S->N); }
232 : static int
233 43953 : 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 : (void)gc_all(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 14616 : quadunit_q(GEN D, GEN d, long *pN)
321 : {
322 14616 : pari_sp av = avma;
323 : GEN p, q, q1;
324 14616 : long first = 1;
325 14616 : p = (Mod2(d) == Mod2(D))? d: subiu(d, 1);
326 14616 : q = gen_2;
327 14616 : q1 = shifti(subii(D, sqri(p)), -1);
328 : for(;;)
329 153398 : {
330 168014 : GEN r, A = dvmdii(addii(p, d), q, &r), p1 = p, t;
331 168014 : p = subii(d, r);
332 177443 : if (!first && equalii(p1, p)) { *pN = 1; return q; } /* even period */
333 162827 : first = 0;
334 162827 : t = q; q = submulii(q1, A, subii(p, p1)); q1 = t;
335 162827 : if (equalii(q, t)) { *pN = -1; return q; } /* odd period */
336 153398 : 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 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 = gc_INT(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 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 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 : (void)gc_all(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 gc_GEN(av, mkquad(quadpoly_i(D0), u, v));
518 : }
519 : long
520 69034 : quadunitnorm(GEN D)
521 : {
522 69034 : pari_sp av = avma;
523 : long s, r;
524 69034 : check_quaddisc(D, &s, &r, "quadunitnorm");
525 69027 : if (s < 0) return 1;
526 69027 : if (mod4(D) == 1 && BPSW_psp(D)) return -1;
527 61152 : if (Z_has_prime_3mod4(D)) return 1;
528 14518 : (void)quadunit_q(D, sqrti(D), &s); return gc_long(av, s);
529 : }
530 :
531 : GEN
532 49 : quadregulator(GEN x, long prec)
533 : {
534 49 : pari_sp av = avma, av2;
535 : GEN R, rsqd, u, v, sqd;
536 : long r, e;
537 :
538 49 : check_quaddisc_real(x, &r, "quadregulator");
539 49 : sqd = sqrti(x);
540 49 : rsqd = gsqrt(x,prec); av2 = avma;
541 49 : e = 0; R = real2n(1, prec); u = utoi(r); v = gen_2;
542 : for(;;)
543 140 : {
544 189 : GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
545 189 : GEN v1 = divii(subii(x,sqri(u1)),v);
546 189 : if (equalii(v,v1)) { R = mulrr(sqrr(R), divri(addir(u1,rsqd),v)); break; }
547 154 : if (equalii(u,u1)) { R = sqrr(R); break; }
548 140 : R = mulrr(R, divri(addir(u1,rsqd),v));
549 140 : e += expo(R); setexpo(R,0);
550 140 : u = u1; v = v1;
551 140 : if (e & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
552 140 : 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 49 : R = divri(R, v); e = 2*e - 1;
559 : /* avoid loss of accuracy */
560 49 : if (!((e + expo(R)) & ~EXPOBITS)) { setexpo(R, e + expo(R)); e = 0; }
561 49 : R = logr_abs(R);
562 49 : if (e) R = addrr(R, mulsr(e, mplog2(prec)));
563 49 : return gc_leaf(av, R);
564 : }
565 :
566 : /*************************************************************************/
567 : /** **/
568 : /** CLASS NUMBER **/
569 : /** **/
570 : /*************************************************************************/
571 :
572 : int
573 12551986 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
574 :
575 : static GEN
576 10887055 : qfi_pow(void *E, GEN f, GEN n)
577 10887055 : { return E? nupow(f,n,(GEN)E): qfbpow_i(f,n); }
578 : static GEN
579 7551261 : qfi_comp(void *E, GEN f, GEN g)
580 7551261 : { return E? nucomp(f,g,(GEN)E): qfbcomp_i(f,g); }
581 : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
582 : gidentical,qfb_equal1,NULL};
583 :
584 : GEN
585 3525282 : qfi_order(GEN q, GEN o)
586 3525282 : { return gen_order(q, o, NULL, &qfi_group); }
587 :
588 : GEN
589 0 : qfi_log(GEN a, GEN g, GEN o)
590 0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
591 :
592 : GEN
593 761953 : qfi_Shanks(GEN a, GEN g, long n)
594 : {
595 761953 : pari_sp av = avma;
596 : GEN T, X;
597 : long rt_n, c;
598 :
599 761953 : a = qfi_red(a);
600 761953 : g = qfi_red(g);
601 :
602 761953 : rt_n = sqrt((double)n);
603 761953 : c = n / rt_n;
604 761953 : c = (c * rt_n < n + 1) ? c + 1 : c;
605 :
606 761953 : T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
607 761953 : X = gen_Shanks(T, a, c, NULL, &qfi_group);
608 761953 : return X? gc_INT(av, X): gc_NULL(av);
609 : }
610 :
611 : GEN
612 462 : qfbclassno0(GEN x,long flag)
613 : {
614 462 : switch(flag)
615 : {
616 448 : case 0: return map_proto_G(classno,x);
617 14 : case 1: return map_proto_G(classno2,x);
618 0 : default: pari_err_FLAG("qfbclassno");
619 : }
620 : return NULL; /* LCOV_EXCL_LINE */
621 : }
622 :
623 : /* f^h = 1, return order(f). Set *pfao to its factorization */
624 : static GEN
625 5418 : find_order(void *E, GEN f, GEN h, GEN *pfao)
626 : {
627 5418 : GEN v = gen_factored_order(f, h, E, &qfi_group);
628 5418 : *pfao = gel(v,2); return gel(v,1);
629 : }
630 :
631 : static int
632 301 : ok_q(GEN q, GEN h, GEN d2, long r2)
633 : {
634 301 : if (d2)
635 : {
636 245 : if (r2 <= 2 && !mpodd(q)) return 0;
637 245 : return is_pm1(Z_ppo(q,d2));
638 : }
639 : else
640 : {
641 56 : if (r2 <= 1 && !mpodd(q)) return 0;
642 56 : return is_pm1(Z_ppo(q,h));
643 : }
644 : }
645 :
646 : /* a,b given by their factorizations. Return factorization of lcm(a,b).
647 : * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
648 : static GEN
649 182 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
650 : {
651 182 : GEN P = ZC_union_shallow(gel(Fa,1), gel(Fb,1));
652 182 : GEN A = gen_1, B = gen_1;
653 182 : long i, l = lg(P);
654 182 : GEN E = cgetg(l, t_COL);
655 707 : for (i=1; i<l; i++)
656 : {
657 525 : GEN p = gel(P,i);
658 525 : long va = Z_pval(a,p);
659 525 : long vb = Z_pval(b,p);
660 525 : if (va < vb)
661 : {
662 203 : B = mulii(B,powiu(p,vb));
663 203 : gel(E,i) = utoi(vb);
664 : }
665 : else
666 : {
667 322 : A = mulii(A,powiu(p,va));
668 322 : gel(E,i) = utoi(va);
669 : }
670 : }
671 182 : *pA = A;
672 182 : *pB = B; return mkmat2(P,E);
673 : }
674 :
675 : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
676 : static void
677 182 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
678 : {
679 182 : GEN A,B, g1 = *pg1, d1 = *pd1;
680 182 : *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
681 182 : *pg1 = gmul(qfbpow_i(g1, diviiexact(d1,A)), qfbpow_i(f, diviiexact(o,B)));
682 182 : *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
683 182 : }
684 :
685 : /* Let s = 1 or -1; D = s * d; assume Df^2 fits in an ulong
686 : * Return f / [O_{Df^2}^*:O_D^*] * \prod_{p|f} [ 1 - (D/p) p^-1 ]
687 : * The Euler product is \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
688 : ulong
689 7238 : uquadclassnoF_fact(ulong d, long s, GEN P, GEN E)
690 : {
691 7238 : long i, l = lg(P);
692 7238 : ulong H = 1;
693 9267 : for (i = 1; i < l; i++)
694 : {
695 2029 : ulong p = P[i], e = E[i];
696 2029 : long D = (long)(p == 2? d & 7: d % p), a;
697 2029 : if (s < 0) D = -D;
698 2029 : a = kross(D,p);
699 2029 : if (!a)
700 532 : H *= upowuu(p, e);
701 : else
702 : {
703 1497 : H *= p - a;
704 1497 : if (e >= 2) H *= upowuu(p, e-1);
705 : }
706 : }
707 7238 : if (l == 1) return H;
708 1532 : if (s < 0)
709 : {
710 1512 : switch(d)
711 : { /* divide by [ O_K^* : O^* ] */
712 105 : case 4: H >>= 1; break;
713 287 : case 3: H /= 3; break;
714 : }
715 : }
716 : else
717 : {
718 20 : GEN fa = mkmat2(zc_to_ZC(P), zc_to_ZC(E));
719 20 : H /= itou(quadunitindex_i(utoipos(d), factorback(fa), fa));
720 : }
721 1532 : return H;
722 : }
723 : GEN
724 48 : quadclassnoF_fact(GEN D, GEN P, GEN E)
725 : {
726 48 : GEN H = quadclassnoEuler_fact(D, P, E);
727 48 : if (lg(P) == 1) return H;
728 15 : if (signe(D) < 0)
729 : {
730 7 : switch(itou_or_0(D))
731 : { /* divide by [ O_K^* : O^* ] */
732 0 : case 4: H = shifti(H,-1); break;
733 0 : case 3: H = diviuexact(H,3); break;
734 : }
735 : }
736 : else
737 : {
738 8 : GEN fa = mkvec2(P, E);
739 8 : H = diviiexact(H, quadunitindex_i(D, factorback2(P, E), fa));
740 : }
741 15 : return H;
742 : }
743 :
744 : static ulong
745 673 : quadclassnoF_u(ulong x, long s, ulong *pD)
746 : {
747 673 : pari_sp av = avma;
748 : GEN P, E;
749 673 : ulong D = coredisc2u_fact(factoru(x), s, &P, &E);
750 673 : long H = uquadclassnoF_fact(D, s, P, E);
751 673 : *pD = D; return gc_ulong(av, H);
752 : }
753 : ulong
754 0 : unegquadclassnoF(ulong x, ulong *pD) { return quadclassnoF_u(x, -1, pD); }
755 : ulong
756 0 : uposquadclassnoF(ulong x, ulong *pD) { return quadclassnoF_u(x, 1, pD); }
757 :
758 : /* *pD = coredisc(x) */
759 : GEN
760 721 : quadclassnoF(GEN x, GEN *pD)
761 : {
762 : GEN D, P, E;
763 721 : if (lgefint(x) == 3)
764 : {
765 673 : long s = signe(x);
766 673 : ulong d, h = quadclassnoF_u(x[2], s, &d);
767 673 : if (pD) *pD = s > 0? utoipos(d): utoineg(d);
768 673 : return utoipos(h);
769 : }
770 48 : D = coredisc2_fact(absZ_factor(x), signe(x), &P, &E);
771 48 : if (pD) *pD = D;
772 48 : return quadclassnoF_fact(D, P, E);
773 : }
774 :
775 : static long
776 672 : two_rank(GEN x)
777 : {
778 672 : GEN p = gel(absZ_factor(x),1);
779 672 : long l = lg(p)-1;
780 : #if 0 /* positive disc not needed */
781 : if (signe(x) > 0)
782 : {
783 : long i;
784 : for (i=1; i<=l; i++)
785 : if (mod4(gel(p,i)) == 3) { l--; break; }
786 : }
787 : #endif
788 672 : return l-1;
789 : }
790 :
791 : static GEN
792 12733 : sqr_primeform(GEN x, ulong p) { return qfbsqr_i(primeform_u(x, p)); }
793 : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
794 : static GEN
795 672 : get_forms(GEN D, GEN *pL)
796 : {
797 672 : const long MAXFORM = 20;
798 672 : GEN L, sqrtD = gsqrt(absi_shallow(D),DEFAULTPREC);
799 672 : GEN forms = vectrunc_init(MAXFORM+1);
800 672 : long s, nforms = 0;
801 : ulong p;
802 : forprime_t S;
803 672 : L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
804 672 : s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
805 672 : if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
806 672 : if (s < 10) s = 200;
807 441 : else if (s < 20) s = 1000;
808 413 : else if (s < 5000) s = 5000;
809 672 : u_forprime_init(&S, 2, s);
810 13474321 : while ( (p = u_forprime_next(&S)) )
811 : {
812 13473649 : long d, k = kroiu(D,p);
813 : pari_sp av2;
814 13473649 : if (!k) continue;
815 13472473 : if (k > 0)
816 : {
817 6737493 : if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
818 6737493 : d = p-1;
819 : }
820 : else
821 6734980 : d = p+1;
822 13472473 : av2 = avma; affrr(divru(mulur(p,L),d), L); set_avma(av2);
823 : }
824 672 : *pL = L; return forms;
825 : }
826 :
827 : /* h ~ #G, return o = order of f, set fao = its factorization */
828 : static GEN
829 770 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
830 : {
831 770 : long s = minss(itos(sqrti(h)), 10000);
832 770 : GEN T = gen_Shanks_init(f, s, E, &qfi_group);
833 770 : GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
834 770 : return find_order(E, f, addiu(v,1), pfao);
835 : }
836 :
837 : /* if g = 1 in G/<f> ? */
838 : static int
839 5684 : equal1(void *E, GEN T, ulong N, GEN g)
840 5684 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
841 :
842 : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
843 : * FIXME: should be gen_order, but equal1 has the wrong prototype */
844 : static GEN
845 350 : relative_order(void *E, GEN a, GEN o, ulong N, GEN T)
846 : {
847 350 : pari_sp av = avma;
848 : long i, l;
849 : GEN m;
850 :
851 350 : m = get_arith_ZZM(o);
852 350 : if (!m) pari_err_TYPE("gen_order [missing order]",a);
853 350 : o = gel(m,1);
854 350 : m = gel(m,2); l = lgcols(m);
855 1148 : for (i = l-1; i; i--)
856 : {
857 798 : GEN t, y, p = gcoeff(m,i,1);
858 798 : long j, e = itos(gcoeff(m,i,2));
859 798 : if (l == 2) {
860 49 : t = gen_1;
861 49 : y = a;
862 : } else {
863 749 : t = diviiexact(o, powiu(p,e));
864 749 : y = powgi(a, t);
865 : }
866 798 : if (equal1(E, T,N,y)) o = t;
867 : else {
868 364 : for (j = 1; j < e; j++)
869 : {
870 84 : y = powgi(y, p);
871 84 : if (equal1(E, T,N,y)) break;
872 : }
873 357 : if (j < e) {
874 77 : if (j > 1) p = powiu(p, j);
875 77 : o = mulii(t, p);
876 : }
877 : }
878 : }
879 350 : return gc_GEN(av, o);
880 : }
881 :
882 : /* h(x) for x<0 using Baby Step/Giant Step.
883 : * Assumes G is not too far from being cyclic.
884 : *
885 : * Compute G^2 instead of G so as to kill most of the noncyclicity */
886 : GEN
887 819 : classno(GEN x)
888 : {
889 819 : pari_sp av = avma;
890 : long r2, k, s, i, l;
891 : GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
892 : void *E;
893 :
894 819 : if (signe(x) >= 0) return classno2(x);
895 :
896 784 : check_quaddisc(x, &s, &k, "classno");
897 784 : if (abscmpiu(x,12) <= 0) return gen_1;
898 :
899 672 : Hf = quadclassnoF(x, &D);
900 672 : if (abscmpiu(D,12) <= 0) return gc_GEN(av, Hf);
901 672 : forms = get_forms(D, &L);
902 672 : r2 = two_rank(D);
903 672 : hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
904 :
905 672 : l = lg(forms);
906 672 : order_bound = const_vec(l-1, NULL);
907 672 : E = expi(D) > 60? (void*)sqrtnint(shifti(absi_shallow(D),-2),4): NULL;
908 672 : g1 = gel(forms,1);
909 672 : gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
910 672 : q = diviiround(hin, d1); /* approximate order of G/<g1> */
911 672 : d2 = NULL; /* not computed yet */
912 672 : if (is_pm1(q)) goto END;
913 7315 : for (i=2; i < l; i++)
914 : {
915 6951 : GEN o, fao, a, F, fd, f = gel(forms,i);
916 6951 : fd = qfbpow_i(f, d1); if (is_pm1(gel(fd,1))) continue;
917 182 : F = qfbpow_i(fd, q);
918 182 : a = gel(F,1);
919 182 : o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
920 : /* f^(d1 q) = 1 */
921 182 : fao = ZM_merge_factor(fad1,fao);
922 182 : o = find_order(E, f, fao, &fao);
923 182 : gel(order_bound,i) = o;
924 : /* o = order of f, fao = factor(o) */
925 182 : update_g1(&g1,&d1,&fad1, f,o,fao);
926 182 : q = diviiround(hin, d1);
927 182 : if (is_pm1(q)) goto END;
928 : }
929 : /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
930 364 : if (expi(q) > 3)
931 : { /* q large: compute d2, 2nd elt divisor */
932 308 : ulong N, n = 2*itou(sqrti(d1));
933 308 : GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
934 308 : d2 = gen_1;
935 308 : N = itou( gceil(gdivgu(d1,n)) ); /* order(g1) <= n*N */
936 5047 : for (i = 1; i < l; i++)
937 : {
938 4802 : GEN d, f = gel(forms,i), B = gel(order_bound,i);
939 4802 : if (!B) B = find_order(E, f, fad1, /*junk*/&d);
940 4802 : f = qfbpow_i(f,d2);
941 4802 : if (equal1(E, T,N,f)) continue;
942 350 : B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
943 : /* f^B = 1 */
944 350 : d = relative_order(E, f, B, N,T);
945 350 : d2= mulii(d,d2);
946 350 : D = mulii(d1,d2);
947 350 : q = diviiround(hin,D);
948 350 : if (is_pm1(q)) { d1 = D; goto END; }
949 : }
950 : /* very probably, d2 is the 2nd elementary divisor */
951 245 : d1 = D; /* product of first two elt divisors */
952 : }
953 : /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
954 : * 2-rank */
955 301 : if (!ok_q(q,d1,d2,r2))
956 : {
957 0 : GEN q0 = q;
958 : long d;
959 0 : if (cmpii(mulii(q,d1), hin) < 0)
960 : { /* try q = q0+1,-1,+2,-2 */
961 0 : d = 1;
962 0 : do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
963 : }
964 : else
965 : { /* q0-1,+1,-2,+2 */
966 0 : d = -1;
967 0 : do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
968 : }
969 : }
970 301 : d1 = mulii(d1,q);
971 :
972 672 : END:
973 672 : return gc_INT(av, shifti(mulii(d1,Hf), r2));
974 : }
975 :
976 : /* use Euler products */
977 : GEN
978 49 : classno2(GEN x)
979 : {
980 49 : pari_sp av = avma;
981 49 : const long prec = DEFAULTPREC;
982 : long n, i, s;
983 49 : GEN p1, p2, S, p4, p5, p7, Hf, Pi, logd, sqrtd, D, half, reg = NULL;
984 :
985 49 : check_quaddisc(x, &s, NULL, "classno2");
986 49 : if (s < 0 && abscmpiu(x,12) <= 0) return gen_1;
987 :
988 49 : Hf = quadclassnoF(x, &D);
989 49 : if (s < 0 && abscmpiu(D,12) <= 0) return gc_GEN(av, Hf); /* |D| < 12*/
990 :
991 49 : Pi = mppi(prec);
992 49 : sqrtd = sqrtr_abs(itor(D, prec));
993 49 : logd = logr_abs(sqrtd); shiftr_inplace(logd,1);
994 49 : p1 = sqrtr_abs(divrr(mulir(D,logd), gmul2n(Pi,1)));
995 49 : if (s > 0)
996 : {
997 42 : GEN invlogd = invr(logd);
998 42 : reg = quadregulator(D, prec);
999 42 : p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
1000 42 : if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
1001 : }
1002 49 : n = itos_or_0( mptrunc(p1) );
1003 49 : if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
1004 :
1005 49 : p4 = divri(Pi, D); setsigne(p4, 1);
1006 49 : p7 = invr(sqrtr_abs(Pi));
1007 49 : half = real2n(-1, prec);
1008 49 : if (s > 0)
1009 : { /* i = 1, shortcut */
1010 42 : p1 = sqrtd;
1011 42 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
1012 42 : S = addrr(mulrr(p1,p5), eint1(p4,prec));
1013 1358 : for (i=2; i<=n; i++)
1014 : {
1015 1316 : long k = kroiu(D,i); if (!k) continue;
1016 1218 : p2 = mulir(sqru(i), p4);
1017 1218 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
1018 1218 : p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
1019 1218 : S = (k>0)? addrr(S,p5): subrr(S,p5);
1020 : }
1021 42 : S = shiftr(divrr(S,reg),-1);
1022 : }
1023 : else
1024 : { /* i = 1, shortcut */
1025 7 : p1 = gdiv(sqrtd, Pi);
1026 7 : p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
1027 7 : S = addrr(p5, divrr(p1, mpexp(p4)));
1028 952 : for (i=2; i<=n; i++)
1029 : {
1030 945 : long k = kroiu(D,i); if (!k) continue;
1031 945 : p2 = mulir(sqru(i), p4);
1032 945 : p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
1033 945 : p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
1034 945 : S = (k>0)? addrr(S,p5): subrr(S,p5);
1035 : }
1036 : }
1037 49 : return gc_INT(av, mulii(Hf, roundr(S)));
1038 : }
1039 :
1040 : /* 1 + q + ... + q^v, v > 0 */
1041 : static GEN
1042 1793 : geomsumu(ulong q, long v)
1043 : {
1044 1793 : GEN u = utoipos(1+q);
1045 2171 : for (; v > 1; v--) u = addui(1, mului(q, u));
1046 1793 : return u;
1047 : }
1048 : static GEN
1049 1793 : geomsum(GEN q, long v)
1050 : {
1051 : GEN u;
1052 1793 : if (lgefint(q) == 3) return geomsumu(q[2], v);
1053 0 : u = addiu(q,1);
1054 0 : for (; v > 1; v--) u = addui(1, mulii(q, u));
1055 0 : return u;
1056 : }
1057 :
1058 : /* 1+p+...+p^(e-1), e >= 1; assuming result fits in an ulong */
1059 : static ulong
1060 10387 : usumpow(ulong p, long e)
1061 : {
1062 : ulong q;
1063 : long i;
1064 10387 : if (p == 2) return (1UL << e) - 1; /* also OK if e = BITS_IN_LONG */
1065 7400 : e--; for (i = 1, q = p + 1; i < e; i++) q = p*q + 1;
1066 7345 : return q;
1067 : }
1068 : long
1069 172854 : uhclassnoF_fact(GEN faF, long D)
1070 : {
1071 172854 : GEN P = gel(faF,1), E = gel(faF,2);
1072 172854 : long i, t, l = lg(P);
1073 377788 : for (i = t = 1; i < l; i++)
1074 : {
1075 204934 : long p = P[i], e = E[i], s = kross(D,p);
1076 204934 : if (e == 1) { t *= 1 + p - s; continue; }
1077 54701 : if (s == 1) { t *= upowuu(p,e); continue; }
1078 10387 : t *= 1 + usumpow(p, e) * (p - s);
1079 : }
1080 172854 : return t;
1081 : }
1082 : /* Hurwitz(D F^2)/ Hurwitz(D)
1083 : * = \sum_{f|F} f \prod_{p|f} (1-kro(D/p)/p)
1084 : * = \prod_{p^e || F} (1 + (p^e-1) / (p-1) * (p-kro(D/p))) */
1085 : GEN
1086 112514 : hclassnoF_fact(GEN P, GEN E, GEN D)
1087 : {
1088 : GEN H;
1089 112514 : long i, l = lg(P);
1090 112514 : if (l == 1) return gen_1;
1091 104215 : for (i = 1, H = NULL; i < l; i++)
1092 : {
1093 56116 : GEN t, p = gel(P,i);
1094 56116 : long e = E[i], s = kronecker(D,p);
1095 56116 : if (e == 1) t = addiu(p, 1-s);
1096 2879 : else if (s == 1) t = powiu(p, e);
1097 1793 : else t = addui(1, mulii(subis(p, s), geomsum(p, e-1)));
1098 56116 : H = H? mulii(H,t): t;
1099 : }
1100 48099 : return H;
1101 : }
1102 : static GEN
1103 112514 : hclassno6_large(GEN x)
1104 : {
1105 112514 : GEN H = NULL, P, E, D = coredisc2_fact(absZ_factor(x), -1, &P, &E);
1106 112514 : long l = lg(P);
1107 :
1108 112514 : if (l > 1 && lgefint(x) == 3)
1109 : { /* F != 1, second chance */
1110 48098 : ulong h = hclassno6u_no_cache(x[2]);
1111 48098 : if (h) H = utoipos(h);
1112 : }
1113 112514 : if (!H)
1114 : {
1115 112514 : H = quadclassno(D);
1116 112514 : switch(itou_or_0(D))
1117 : {
1118 91 : case 3: H = shifti(H,1);break;
1119 7 : case 4: H = muliu(H,3); break;
1120 112416 : default:H = muliu(H,6); break;
1121 : }
1122 : }
1123 112514 : return mulii(H, hclassnoF_fact(P, E, D));
1124 : }
1125 :
1126 : /* x > 0, x = 0,3 (mod 4). Return 6*hclassno(x), an integer */
1127 : GEN
1128 185842 : hclassno6(GEN x)
1129 : {
1130 185842 : ulong d = itou_or_0(x);
1131 185842 : if (d)
1132 : { /* create cache if d small */
1133 185841 : ulong h = d < 500000 ? hclassno6u(d): hclassno6u_no_cache(d);
1134 185837 : if (h) return utoipos(h);
1135 : }
1136 112514 : return hclassno6_large(x);
1137 : }
1138 :
1139 : GEN
1140 49 : hclassno(GEN x)
1141 : {
1142 : long a, s;
1143 49 : if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
1144 49 : s = signe(x);
1145 49 : if (s < 0) return gen_0;
1146 49 : if (!s) return gdivgs(gen_1, -12);
1147 49 : a = mod4(x); if (a == 1 || a == 2) return gen_0;
1148 49 : return gdivgu(hclassno6(x), 6);
1149 : }
1150 :
1151 : /* return [N',v]; v contains all x mod N' s.t. x^2 + B x + C = 0 modulo N */
1152 : GEN
1153 2664249 : Zn_quad_roots(GEN N, GEN B, GEN C)
1154 : {
1155 2664249 : pari_sp av = avma;
1156 : GEN fa, D, w, v, P, E, F0, Q0, F, mF, A, Q, T, R, Np, N4;
1157 : long l, i, j, ct;
1158 :
1159 2664249 : if ((fa = check_arith_non0(N,"Zn_quad_roots")))
1160 : {
1161 7665 : N = typ(N) == t_VEC? gel(N,1): factorback(N);
1162 7665 : fa = clean_Z_factor(fa);
1163 : }
1164 2664249 : N = absi_shallow(N);
1165 2664249 : N4 = shifti(N,2);
1166 2664249 : D = modii(subii(sqri(B), shifti(C,2)), N4);
1167 2664249 : if (!signe(D))
1168 : { /* (x + B/2)^2 = 0 (mod N), D = B^2-4C = 0 (4N) => B even */
1169 763 : if (!fa) fa = Z_factor(N);
1170 763 : P = gel(fa,1);
1171 763 : E = ZV_to_zv(gel(fa,2));
1172 763 : l = lg(P);
1173 1652 : for (i = 1; i < l; i++) E[i] = (E[i]+1) >> 1;
1174 763 : Np = factorback2(P, E); /* x = -B mod N' */
1175 763 : B = shifti(B,-1);
1176 763 : return gc_GEN(av, mkvec2(Np, mkvec(Fp_neg(B,Np))));
1177 : }
1178 2663486 : if (!fa)
1179 2655989 : fa = Z_factor(N4);
1180 : else /* convert to factorization of N4 = 4*N */
1181 7497 : fa = famat_reduce(famat_mulpows_shallow(fa, gen_2, 2));
1182 2663486 : P = gel(fa,1); l = lg(P);
1183 2663486 : E = ZV_to_zv(gel(fa,2));
1184 2663486 : F = cgetg(l, t_VEC);
1185 2663486 : mF= cgetg(l, t_VEC); F0 = gen_0;
1186 2663486 : Q = cgetg(l, t_VEC); Q0 = gen_1;
1187 6391560 : for (i = j = 1, ct = 0; i < l; i++)
1188 : {
1189 5847954 : GEN p = gel(P,i), q, f, mf, D0;
1190 5847954 : long t2, s = E[i], t = Z_pvalrem(D, p, &D0), d = s - t;
1191 5847954 : if (d <= 0)
1192 : {
1193 1325947 : q = powiu(p, (s+1)>>1);
1194 2253391 : Q0 = mulii(Q0, q); continue;
1195 : }
1196 : /* d > 0 */
1197 6456240 : if (odd(t)) return NULL;
1198 4336360 : t2 = t >> 1;
1199 4336360 : if (i > 1)
1200 : { /* p > 2 */
1201 2741053 : if (kronecker(D0, p) == -1) return NULL;
1202 1325562 : q = powiu(p, s - t2);
1203 1325562 : f = Zp_sqrt(D0, p, d);
1204 1325562 : if (!f) return NULL; /* p was not actually prime... */
1205 1325562 : if (t2) f = mulii(powiu(p,t2), f);
1206 1325562 : mf = Fp_neg(f, q);
1207 : }
1208 : else
1209 : { /* p = 2 */
1210 1595307 : if (d <= 3)
1211 : {
1212 1215242 : if (d == 3 && Mod8(D0) != 1) return NULL;
1213 986706 : if (d == 2 && Mod4(D0) != 1) return NULL;
1214 927444 : Q0 = int2n(1+t2); F0 = NULL; continue;
1215 : }
1216 380065 : if (Mod8(D0) != 1) return NULL;
1217 149121 : q = int2n(d - 1 + t2);
1218 149121 : f = Z2_sqrt(D0, d);
1219 149121 : if (t2) f = shifti(f, t2);
1220 149121 : mf = Fp_neg(f, q);
1221 : }
1222 1474683 : gel(Q,j) = q;
1223 1474683 : gel(F,j) = f;
1224 1474683 : gel(mF,j)= mf; j++;
1225 : }
1226 543606 : setlg(Q,j);
1227 543606 : setlg(F,j);
1228 543606 : setlg(mF,j);
1229 543606 : if (is_pm1(Q0)) A = leafcopy(F);
1230 : else
1231 : { /* append the fixed congruence (F0 mod Q0) */
1232 504819 : if (!F0) F0 = shifti(Q0,-1);
1233 504819 : A = shallowconcat(F, F0);
1234 504819 : Q = shallowconcat(Q, Q0);
1235 : }
1236 543606 : if (j-1 >= LGnumBITS) pari_err_OVERFLOW("Zn_quad_roots");
1237 543606 : ct = 1L << (j-1);
1238 543606 : T = ZV_producttree(Q);
1239 543606 : R = ZV_chinesetree(Q,T);
1240 543606 : Np = gmael(T, lg(T)-1, 1);
1241 543606 : B = modii(B, Np);
1242 543606 : if (!signe(B)) B = NULL;
1243 543606 : Np = shifti(Np, -1); /* N' = (\prod_i Q[i]) / 2 */
1244 543606 : w = cgetg(3, t_VEC);
1245 543606 : gel(w,1) = icopy(Np);
1246 543606 : gel(w,2) = v = cgetg(ct+1, t_VEC);
1247 543606 : l = lg(F);
1248 2547699 : for (j = 1; j <= ct; j++)
1249 : {
1250 2004093 : pari_sp av2 = avma;
1251 2004093 : long m = j - 1;
1252 : GEN u;
1253 6191045 : for (i = 1; i < l; i++)
1254 : {
1255 4186952 : gel(A,i) = (m&1L)? gel(mF,i): gel(F,i);
1256 4186952 : m >>= 1;
1257 : }
1258 2004093 : u = ZV_chinese_tree(A,Q,T,R); /* u mod N' st u^2 = B^2-4C modulo 4N */
1259 2004093 : if (B) u = subii(u,B);
1260 2004093 : gel(v,j) = gc_INT(av2, modii(shifti(u,-1), Np));
1261 : }
1262 543606 : return gc_upto(av, w);
1263 : }
1264 :
1265 : GEN
1266 0 : Zn_sqrtall(GEN D, GEN fa)
1267 : {
1268 0 : pari_sp av = avma;
1269 : long i, j, k, lB, aN;
1270 : GEN a, L, V, B, N;
1271 0 : a = typ(fa) == t_INT ? fa: typ(fa) == t_VEC? gel(fa,1): factorback(fa);
1272 0 : V = Zn_quad_roots(fa, gen_0, negi(D));
1273 0 : if (!V) return NULL;
1274 0 : N = gel(V,1); B = gel(V,2); lB = lg(B);
1275 0 : aN = itou(diviiexact(a, N)); /* |a|/N */
1276 0 : L = cgetg((lB-1)*aN+1, t_VEC);
1277 0 : for (k = 1, i = 1; i < lB; i++)
1278 : {
1279 0 : GEN b = gel(B,i);
1280 0 : for (j = 0;; b = addii(b, N))
1281 : {
1282 0 : gel(L, k++) = b;
1283 0 : if (++j == aN) break;
1284 : }
1285 : }
1286 0 : return gc_upto(av, ZV_sort(L));
1287 : }
|
