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 /*******************************************************************/
16 /* */
17 /* COMPUTATION OF STARK UNITS OF TOTALLY REAL FIELDS */
18 /* */
19 /*******************************************************************/
20 #include "pari.h"
21 #include "paripriv.h"
22
23 /* ComputeCoeff */
24 typedef struct {
25 GEN L0, L1, L11, L2; /* VECSMALL of p */
26 GEN L1ray, L11ray; /* precomputed isprincipalray(pr), pr | p */
27 GEN rayZ; /* precomputed isprincipalray(i), i < condZ */
28 long condZ; /* generates cond(bnr) \cap Z, assumed small */
29 } LISTray;
30
31 /* Char evaluation */
32 typedef struct {
33 long ord;
34 GEN *val, chi;
35 } CHI_t;
36
37 /* RecCoeff */
38 typedef struct {
39 GEN M, beta, B, U, nB;
40 long v, G, N;
41 } RC_data;
42
43 /********************************************************************/
44 /* Miscellaneous functions */
45 /********************************************************************/
46 static GEN
chi_get_c(GEN chi)47 chi_get_c(GEN chi) { return gmael(chi,1,2); }
48 static long
chi_get_deg(GEN chi)49 chi_get_deg(GEN chi) { return itou(gmael(chi,1,1)); }
50
51 /* Compute the image of logelt by character chi, zeta_ord(chi)^n; return n */
52 static ulong
CharEval_n(GEN chi,GEN logelt)53 CharEval_n(GEN chi, GEN logelt)
54 {
55 GEN gn = ZV_dotproduct(chi_get_c(chi), logelt);
56 return umodiu(gn, chi_get_deg(chi));
57 }
58 /* Compute the image of logelt by character chi, as a complex number */
59 static GEN
CharEval(GEN chi,GEN logelt)60 CharEval(GEN chi, GEN logelt)
61 {
62 ulong n = CharEval_n(chi, logelt), d = chi_get_deg(chi);
63 long nn = Fl_center(n,d,d>>1);
64 GEN x = gel(chi,2);
65 x = gpowgs(x, labs(nn));
66 if (nn < 0) x = conj_i(x);
67 return x;
68 }
69
70 /* return n such that C(elt) = z^n */
71 static ulong
CHI_eval_n(CHI_t * C,GEN logelt)72 CHI_eval_n(CHI_t *C, GEN logelt)
73 {
74 GEN n = ZV_dotproduct(C->chi, logelt);
75 return umodiu(n, C->ord);
76 }
77 /* return C(elt) */
78 static GEN
CHI_eval(CHI_t * C,GEN logelt)79 CHI_eval(CHI_t *C, GEN logelt)
80 {
81 return C->val[CHI_eval_n(C, logelt)];
82 }
83
84 static void
init_CHI(CHI_t * c,GEN CHI,GEN z)85 init_CHI(CHI_t *c, GEN CHI, GEN z)
86 {
87 long i, d = chi_get_deg(CHI);
88 GEN *v = (GEN*)new_chunk(d);
89 v[0] = gen_1;
90 if (d != 1)
91 {
92 v[1] = z;
93 for (i=2; i<d; i++) v[i] = gmul(v[i-1], z);
94 }
95 c->chi = chi_get_c(CHI);
96 c->ord = d;
97 c->val = v;
98 }
99 /* as t_POLMOD */
100 static void
init_CHI_alg(CHI_t * c,GEN CHI)101 init_CHI_alg(CHI_t *c, GEN CHI) {
102 long d = chi_get_deg(CHI);
103 GEN z;
104 switch(d)
105 {
106 case 1: z = gen_1; break;
107 case 2: z = gen_m1; break;
108 default: z = mkpolmod(pol_x(0), polcyclo(d,0));
109 }
110 init_CHI(c,CHI, z);
111 }
112 /* as t_COMPLEX */
113 static void
init_CHI_C(CHI_t * c,GEN CHI)114 init_CHI_C(CHI_t *c, GEN CHI) {
115 init_CHI(c,CHI, gel(CHI,2));
116 }
117
118 typedef struct {
119 long r; /* rank = lg(gen) */
120 GEN j; /* current elt is gen[1]^j[1] ... gen[r]^j[r] */
121 GEN cyc; /* t_VECSMALL of elementary divisors */
122 } GROUP_t;
123
124 static int
NextElt(GROUP_t * G)125 NextElt(GROUP_t *G)
126 {
127 long i = 1;
128 if (G->r == 0) return 0; /* no more elt */
129 while (++G->j[i] == G->cyc[i]) /* from 0 to cyc[i]-1 */
130 {
131 G->j[i] = 0;
132 if (++i > G->r) return 0; /* no more elt */
133 }
134 return i; /* we have multiplied by gen[i] */
135 }
136
137 /* enumerate all group elements; trivial elt comes last */
138 GEN
cyc2elts(GEN cyc)139 cyc2elts(GEN cyc)
140 {
141 long i, n;
142 GEN z;
143 GROUP_t G;
144
145 G.cyc = typ(cyc)==t_VECSMALL? cyc: gtovecsmall(cyc);
146 n = zv_prod(G.cyc);
147 G.r = lg(cyc)-1;
148 G.j = zero_zv(G.r);
149
150 z = cgetg(n+1, t_VEC);
151 gel(z,n) = leafcopy(G.j); /* trivial elt comes last */
152 for (i = 1; i < n; i++)
153 {
154 (void)NextElt(&G);
155 gel(z,i) = leafcopy(G.j);
156 }
157 return z;
158 }
159
160 /* nchi: a character given by a vector [d, (c_i)], e.g. from char_normalize
161 * such that chi(x) = e((c . log(x)) / d) where log(x) on bnr.gen */
162 static GEN
get_Char(GEN nchi,long prec)163 get_Char(GEN nchi, long prec)
164 { return mkvec2(nchi, rootsof1_cx(gel(nchi,1), prec)); }
165
166 /* prime divisors of conductor */
167 static GEN
divcond(GEN bnr)168 divcond(GEN bnr) {GEN bid = bnr_get_bid(bnr); return gel(bid_get_fact(bid),1);}
169
170 /* vector of prime ideals dividing bnr but not bnrc */
171 static GEN
get_prdiff(GEN D,GEN Dc)172 get_prdiff(GEN D, GEN Dc)
173 {
174 long n, i, l = lg(D);
175 GEN diff = cgetg(l, t_COL);
176 for (n = i = 1; i < l; i++)
177 if (!tablesearch(Dc, gel(D,i), &cmp_prime_ideal)) gel(diff,n++) = gel(D,i);
178 setlg(diff, n); return diff;
179 }
180
181 #define ch_prec(x) realprec(gel(x,1))
182 #define ch_C(x) gel(x,1)
183 #define ch_bnr(x) gel(x,2)
184 #define ch_3(x) gel(x,3)
185 #define ch_q(x) gel(x,3)[1]
186 #define ch_CHI(x) gel(x,4)
187 #define ch_diff(x) gel(x,5)
188 #define ch_CHI0(x) gel(x,6)
189 #define ch_small(x) gel(x,7)
190 #define ch_comp(x) gel(x,7)[1]
191 #define ch_phideg(x) gel(x,7)[2]
192 static long
ch_deg(GEN dtcr)193 ch_deg(GEN dtcr) { return chi_get_deg(ch_CHI(dtcr)); }
194
195 /********************************************************************/
196 /* 1rst part: find the field K */
197 /********************************************************************/
198 static GEN AllStark(GEN data, long flag, long prec);
199
200 /* Columns of C [HNF] give the generators of a subgroup of the finite abelian
201 * group A [ in terms of implicit generators ], compute data to work in A/C:
202 * 1) order
203 * 2) structure
204 * 3) the matrix A ->> A/C
205 * 4) the subgroup C */
206 static GEN
InitQuotient(GEN C)207 InitQuotient(GEN C)
208 {
209 GEN U, D = ZM_snfall_i(C, &U, NULL, 1), h = ZV_prod(D);
210 return mkvec5(h, D, U, C, cyc_normalize(D));
211 }
212
213 /* lift chi character on A/C [Qt from InitQuotient] to character on A [cyc]*/
214 static GEN
LiftChar(GEN Qt,GEN cyc,GEN chi)215 LiftChar(GEN Qt, GEN cyc, GEN chi)
216 {
217 GEN ncyc = gel(Qt,5), U = gel(Qt,3), nchi = char_normalize(chi, ncyc);
218 GEN c = ZV_ZM_mul(gel(nchi,2), U), d = gel(nchi,1);
219 return char_denormalize(cyc, d, c);
220 }
221
222 /* Let C be a subgroup, system of representatives of the quotient */
223 static GEN
ag_subgroup_classes(GEN C)224 ag_subgroup_classes(GEN C)
225 {
226 GEN U, D = ZM_snfall_i(C, &U, NULL, 1), e = cyc2elts(D);
227 long i, l = lg(e);
228
229 if (ZM_isidentity(U))
230 for (i = 1; i < l; i++) (void)vecsmall_to_vec_inplace(gel(e,i));
231 else
232 {
233 GEN Ui = ZM_inv(U,NULL);
234 for (i = 1; i < l; i++) gel(e,i) = ZM_zc_mul(Ui, gel(e,i));
235 }
236 return e;
237 }
238
239 /* Let s: A -> B given by [P,cycA,cycB] A and B, compute the kernel of s. */
240 GEN
ag_kernel(GEN S)241 ag_kernel(GEN S)
242 {
243 GEN U, P = gel(S,1), cycA = gel(S,2), DB = diagonal_shallow(gel(S,3));
244 long nA = lg(cycA)-1, rk;
245
246 rk = nA + lg(DB) - lg(ZM_hnfall_i(shallowconcat(P, DB), &U, 1));
247 return ZM_hnfmodid(matslice(U, 1,nA, 1,rk), cycA);
248 }
249 /* let H be a subgroup of A; return s(H) */
250 GEN
ag_subgroup_image(GEN S,GEN H)251 ag_subgroup_image(GEN S, GEN H)
252 { return ZM_hnfmodid(ZM_mul(gel(S,1), H), gel(S,3)); }
253
254 /* Let m and n be two moduli such that n|m and let C be a congruence
255 group modulo n, compute the corresponding congruence group modulo m
256 ie the kernel of the map Clk(m) ->> Clk(n)/C */
257 static GEN
ComputeKernel(GEN bnrm,GEN bnrn,GEN dtQ)258 ComputeKernel(GEN bnrm, GEN bnrn, GEN dtQ)
259 {
260 pari_sp av = avma;
261 GEN S = bnrsurjection(bnrm, bnrn);
262 GEN P = ZM_mul(gel(dtQ,3), gel(S,1));
263 return gerepileupto(av, ag_kernel(mkvec3(P, gel(S,2), gel(dtQ,2))));
264 }
265
266 static long
cyc_is_cyclic(GEN cyc)267 cyc_is_cyclic(GEN cyc) { return lg(cyc) <= 2 || equali1(gel(cyc,2)); }
268
269 /* Let H be a subgroup of cl(bnr)/sugbroup, return 1 if
270 cl(bnr)/subgoup/H is cyclic and the signature of the
271 corresponding field is equal to sig and no finite prime
272 dividing cond(bnr) is totally split in this field. Return 0
273 otherwise. */
274 static long
IsGoodSubgroup(GEN H,GEN bnr,GEN map)275 IsGoodSubgroup(GEN H, GEN bnr, GEN map)
276 {
277 pari_sp av = avma;
278 GEN S, mod, modH, p1, U, P, PH, bnrH, iH, qH;
279 long j;
280
281 p1 = InitQuotient(H);
282 /* quotient is non cyclic */
283 if (!cyc_is_cyclic(gel(p1,2))) return gc_long(av,0);
284
285 (void)ZM_hnfall_i(shallowconcat(map,H), &U, 0);
286 setlg(U, lg(H));
287 for (j = 1; j < lg(U); j++) setlg(gel(U,j), lg(H));
288 p1 = ZM_hnfmodid(U, bnr_get_cyc(bnr)); /* H as a subgroup of bnr */
289 modH = bnrconductor_raw(bnr, p1);
290 mod = bnr_get_mod(bnr);
291
292 /* is the signature correct? */
293 if (!gequal(gel(modH,2), gel(mod,2))) return gc_long(av, 0);
294
295 /* finite part are the same: OK */
296 if (gequal(gel(modH,1), gel(mod,1))) return gc_long(av, 1);
297
298 /* need to check the splitting of primes dividing mod but not modH */
299 bnrH = Buchray(bnr, modH, nf_INIT);
300 P = divcond(bnr);
301 PH = divcond(bnrH);
302 S = bnrsurjection(bnr, bnrH);
303 /* H as a subgroup of bnrH */
304 iH = ag_subgroup_image(S, p1);
305 qH = InitQuotient(iH);
306 for (j = 1; j < lg(P); j++)
307 {
308 GEN pr = gel(P, j), e;
309 /* if pr divides modH, it is ramified, so it's good */
310 if (tablesearch(PH, pr, cmp_prime_ideal)) continue;
311 /* inertia degree of pr in bnr(modH)/H is charorder(e, cycH) */
312 e = ZM_ZC_mul(gel(qH,3), isprincipalray(bnrH, pr));
313 e = vecmodii(e, gel(qH,2));
314 if (ZV_equal0(e)) return gc_long(av,0); /* f = 1 */
315 }
316 return gc_long(av,1);
317 }
318
319 /* compute list of nontrivial characters trivial on H, modulo complex
320 * conjugation. If flag is set, impose a nontrivial conductor at infinity */
321 static GEN
AllChars(GEN bnr,GEN dtQ,long flag)322 AllChars(GEN bnr, GEN dtQ, long flag)
323 {
324 GEN v, vchi, cyc = bnr_get_cyc(bnr);
325 long n, i, hD = itos(gel(dtQ,1));
326 hashtable *S;
327
328 v = cgetg(hD+1, t_VEC); /* nonconjugate chars */
329 vchi = cyc2elts(gel(dtQ,2));
330 S = hash_create(hD, (ulong(*)(void*))&hash_GEN,
331 (int(*)(void*,void*))&ZV_equal, 1);
332 for (i = n = 1; i < hD; i++) /* remove i = hD: trivial char */
333 { /* lift a character of D in Clk(m) */
334 GEN F, lchi = LiftChar(dtQ, cyc, zv_to_ZV(gel(vchi,i))), cchi = NULL;
335
336 if (hash_search(S, lchi)) continue;
337 F = bnrconductor_raw(bnr, lchi);
338 if (flag && gequal0(gel(F,2))) continue; /* f_oo(chi) trivial ? */
339
340 if (abscmpiu(charorder(cyc,lchi), 2) > 0)
341 { /* nonreal chi: add its conjugate character to S */
342 cchi = charconj(cyc, lchi);
343 hash_insert(S, cchi, (void*)1);
344 }
345 gel(v, n++) = cchi? mkvec3(lchi, F, cchi): mkvec2(lchi, F);
346 }
347 setlg(v, n); return v;
348 }
349
350 static GEN InitChar(GEN bnr, GEN CR, long flag, long prec);
351 static void CharNewPrec(GEN data, long prec);
352
353 /* Given a conductor and a subgroups, return the corresponding complexity and
354 * precision required using quickpol. Fill data[5] with dataCR */
355 static long
CplxModulus(GEN data,long * newprec)356 CplxModulus(GEN data, long *newprec)
357 {
358 long dprec = DEFAULTPREC;
359 pari_sp av = avma;
360 for (;;)
361 {
362 GEN cpl, pol = AllStark(data, -1, dprec);
363 cpl = RgX_fpnorml2(pol, LOWDEFAULTPREC);
364 dprec = maxss(dprec, nbits2extraprec(gexpo(pol))) + EXTRAPREC64;
365 if (!gequal0(cpl)) { *newprec = dprec; return gexpo(cpl); }
366 set_avma(av);
367 if (DEBUGLEVEL>1) pari_warn(warnprec, "CplxModulus", dprec);
368 CharNewPrec(data, dprec);
369 }
370 }
371
372 /* return A \cap B in abelian group defined by cyc. NULL = whole group */
373 static GEN
subgp_intersect(GEN cyc,GEN A,GEN B)374 subgp_intersect(GEN cyc, GEN A, GEN B)
375 {
376 GEN H, U;
377 long k, lH;
378 if (!A) return B;
379 if (!B) return A;
380 H = ZM_hnfall_i(shallowconcat(A,B), &U, 1);
381 setlg(U, lg(A)); lH = lg(H);
382 for (k = 1; k < lg(U); k++) setlg(gel(U,k), lH);
383 return ZM_hnfmodid(ZM_mul(A,U), cyc);
384 }
385
386 static void CharNewPrec(GEN dataCR, long prec);
387 /* Let (f,C) be a conductor without infinite part and a congruence group mod f.
388 * Compute (m,D) such that D is a congruence group of conductor m, f | m,
389 * divisible by all the infinite places but one, D is a subgroup of index 2 of
390 * Cm = Ker: Clk(m) -> Clk(f)/C. Consider further the subgroups H of Clk(m)/D
391 * with cyclic quotient Clk(m)/H such that no place dividing m is totally split
392 * in the extension KH corresponding to H: we want their intersection to be
393 * trivial. These H correspond to (the kernels of Galois orbits of) characters
394 * chi of Clk(m)/D such that chi(log_gen_arch(m_oo)) != 1 and for all pr | m
395 * we either have
396 * - chi(log_gen_pr(pr,1)) != 1 [pr | cond(chi) => ramified in KH]
397 * - or [pr \nmid cond(chi)] chi lifted to Clk(m/pr^oo) is not trivial at pr.
398 * We want the map from Clk(m)/D given by the vector of such caracters to have
399 * trivial kernel. Return bnr(m), D, Ck(m)/D and Clk(m)/Cm */
400 static GEN
FindModulus(GEN bnr,GEN dtQ,long * newprec)401 FindModulus(GEN bnr, GEN dtQ, long *newprec)
402 {
403 const long LIMNORM = 400;
404 long n, i, maxnorm, minnorm, N, pr, rb, iscyc, olde = LONG_MAX;
405 pari_sp av = avma;
406 GEN bnf, nf, f, varch, m, rep = NULL;
407
408 bnf = bnr_get_bnf(bnr);
409 nf = bnf_get_nf(bnf);
410 N = nf_get_degree(nf);
411 f = gel(bnr_get_mod(bnr), 1);
412
413 /* if cpl < rb, it is not necessary to try another modulus */
414 rb = expi( powii(mulii(nf_get_disc(nf), ZM_det_triangular(f)),
415 gmul2n(bnr_get_no(bnr), 3)) );
416
417 /* Initialization of the possible infinite part */
418 varch = cgetg(N+1,t_VEC);
419 for (i = 1; i <= N; i++)
420 {
421 GEN a = const_vec(N,gen_1);
422 gel(a, N+1-i) = gen_0;
423 gel(varch, i) = a;
424 }
425 m = cgetg(3, t_VEC);
426
427 /* Go from minnorm up to maxnorm; if necessary, increase these values.
428 * If the extension is cyclic then a suitable conductor exists and we go on
429 * until we find it. Else, stop at norm LIMNORM. */
430 minnorm = 1;
431 maxnorm = 50;
432 iscyc = cyc_is_cyclic(gel(dtQ,2));
433
434 if (DEBUGLEVEL>1)
435 err_printf("Looking for a modulus of norm: ");
436
437 for(;;)
438 {
439 GEN listid = ideallist0(nf, maxnorm, 4+8); /* ideals of norm <= maxnorm */
440 pari_sp av1 = avma;
441 for (n = minnorm; n <= maxnorm; n++, set_avma(av1))
442 {
443 GEN idnormn = gel(listid,n);
444 long nbidnn = lg(idnormn) - 1;
445 if (DEBUGLEVEL>1) err_printf(" %ld", n);
446 for (i = 1; i <= nbidnn; i++)
447 { /* finite part of the conductor */
448 long s;
449
450 gel(m,1) = idealmul(nf, f, gel(idnormn,i));
451 for (s = 1; s <= N; s++)
452 { /* infinite part */
453 GEN candD, Cm, bnrm;
454 long lD, c;
455
456 gel(m,2) = gel(varch,s);
457 /* compute Clk(m), check if m is a conductor */
458 bnrm = Buchray(bnf, m, nf_INIT);
459 if (!bnrisconductor(bnrm, NULL)) continue;
460
461 /* compute Im(C) in Clk(m)... */
462 Cm = ComputeKernel(bnrm, bnr, dtQ);
463 /* ... and its subgroups of index 2 with conductor m */
464 candD = subgrouplist_cond_sub(bnrm, Cm, mkvec(gen_2));
465 lD = lg(candD);
466 for (c = 1; c < lD; c++)
467 {
468 GEN data, CR, D = gel(candD,c), QD = InitQuotient(D);
469 GEN ord = gel(QD,1), cyc = gel(QD,2), map = gel(QD,3);
470 long e;
471
472 if (!cyc_is_cyclic(cyc)) /* cyclic => suitable, else test */
473 {
474 GEN lH = subgrouplist(cyc, NULL), IK = NULL;
475 long j, ok = 0;
476 for (j = 1; j < lg(lH); j++)
477 {
478 GEN H = gel(lH, j), IH = subgp_intersect(cyc, IK, H);
479 /* if H > IK, no need to test H */
480 if (IK && gidentical(IH, IK)) continue;
481 if (IsGoodSubgroup(H, bnrm, map))
482 {
483 IK = IH; /* intersection of all good subgroups */
484 if (equalii(ord, ZM_det_triangular(IK))) { ok = 1; break; }
485 }
486 }
487 if (!ok) continue;
488 }
489 CR = InitChar(bnrm, AllChars(bnrm, QD, 1), 0, DEFAULTPREC);
490 data = mkvec4(bnrm, D, ag_subgroup_classes(Cm), CR);
491 if (DEBUGLEVEL>1)
492 err_printf("\nTrying modulus = %Ps and subgroup = %Ps\n",
493 bnr_get_mod(bnrm), D);
494 e = CplxModulus(data, &pr);
495 if (DEBUGLEVEL>1) err_printf("cpl = 2^%ld\n", e);
496 if (e < olde)
497 {
498 guncloneNULL(rep); rep = gclone(data);
499 *newprec = pr; olde = e;
500 }
501 if (olde < rb) goto END; /* OK */
502 if (DEBUGLEVEL>1) err_printf("Trying to find another modulus...");
503 }
504 }
505 if (rep) goto END; /* OK */
506 }
507 }
508 /* if necessary compute more ideals */
509 minnorm = maxnorm;
510 maxnorm <<= 1;
511 if (!iscyc && maxnorm > LIMNORM) return NULL;
512 }
513 END:
514 if (DEBUGLEVEL>1)
515 err_printf("No, we're done!\nModulus = %Ps and subgroup = %Ps\n",
516 bnr_get_mod(gel(rep,1)), gel(rep,2));
517 CharNewPrec(rep, *newprec); return gerepilecopy(av, rep);
518 }
519
520 /********************************************************************/
521 /* 2nd part: compute W(X) */
522 /********************************************************************/
523
524 /* find ilambda s.t. Diff*f*ilambda integral and coprime to f
525 and ilambda >> 0 at foo, fa = factorization of f */
526 static GEN
get_ilambda(GEN nf,GEN fa,GEN foo)527 get_ilambda(GEN nf, GEN fa, GEN foo)
528 {
529 GEN x, w, E2, P = gel(fa,1), E = gel(fa,2), D = nf_get_diff(nf);
530 long i, l = lg(P);
531 if (l == 1) return gen_1;
532 w = cgetg(l, t_VEC);
533 E2 = cgetg(l, t_COL);
534 for (i = 1; i < l; i++)
535 {
536 GEN pr = gel(P,i), t = pr_get_tau(pr);
537 long e = itou(gel(E,i)), v = idealval(nf, D, pr);
538 if (v) { D = idealdivpowprime(nf, D, pr, utoipos(v)); e += v; }
539 gel(E2,i) = stoi(e+1);
540 if (typ(t) == t_MAT) t = gel(t,1);
541 gel(w,i) = gdiv(nfpow(nf, t, stoi(e)), powiu(pr_get_p(pr),e));
542 }
543 x = mkmat2(P, E2);
544 return idealchinese(nf, mkvec2(x, foo), w);
545 }
546 /* compute the list of W(chi) such that Ld(s,chi) = W(chi) Ld(1 - s, chi*),
547 * for all chi in LCHI. All chi have the same conductor (= cond(bnr)). */
548 static GEN
ArtinNumber(GEN bnr,GEN LCHI,long prec)549 ArtinNumber(GEN bnr, GEN LCHI, long prec)
550 {
551 long ic, i, j, nz, nChar = lg(LCHI)-1;
552 pari_sp av2;
553 GEN sqrtnc, cond, condZ, cond0, cond1, nf, T, cyc, vN, vB, diff, vt, idh;
554 GEN zid, gen, z, nchi, indW, W, classe, s0, s, den, ilambda, sarch;
555 CHI_t **lC;
556 GROUP_t G;
557
558 lC = (CHI_t**)new_chunk(nChar + 1);
559 indW = cgetg(nChar + 1, t_VECSMALL);
560 W = cgetg(nChar + 1, t_VEC);
561 for (ic = 0, i = 1; i <= nChar; i++)
562 {
563 GEN CHI = gel(LCHI,i);
564 if (chi_get_deg(CHI) <= 2) { gel(W,i) = gen_1; continue; }
565 ic++; indW[ic] = i;
566 lC[ic] = (CHI_t*)new_chunk(sizeof(CHI_t));
567 init_CHI_C(lC[ic], CHI);
568 }
569 if (!ic) return W;
570 nChar = ic;
571
572 nf = bnr_get_nf(bnr);
573 diff = nf_get_diff(nf);
574 T = nf_get_Tr(nf);
575 cond = bnr_get_mod(bnr);
576 cond0 = gel(cond,1); condZ = gcoeff(cond0,1,1);
577 cond1 = gel(cond,2);
578
579 sqrtnc = gsqrt(idealnorm(nf, cond0), prec);
580 ilambda = get_ilambda(nf, bid_get_fact(bnr_get_bid(bnr)), cond1);
581 idh = idealmul(nf, ilambda, idealmul(nf, diff, cond0)); /* integral */
582 ilambda = Q_remove_denom(ilambda, &den);
583 z = den? rootsof1_cx(den, prec): NULL;
584
585 /* compute a system of generators of (Ok/cond)^*, we'll make them
586 * cond1-positive in the main loop */
587 zid = Idealstar(nf, cond0, nf_GEN);
588 cyc = abgrp_get_cyc(zid);
589 gen = abgrp_get_gen(zid);
590 nz = lg(gen) - 1;
591 sarch = nfarchstar(nf, cond0, vec01_to_indices(cond1));
592
593 nchi = cgetg(nChar+1, t_VEC);
594 for (ic = 1; ic <= nChar; ic++) gel(nchi,ic) = cgetg(nz + 1, t_VECSMALL);
595 for (i = 1; i <= nz; i++)
596 {
597 if (is_bigint(gel(cyc,i)))
598 pari_err_OVERFLOW("ArtinNumber [conductor too large]");
599 gel(gen,i) = set_sign_mod_divisor(nf, NULL, gel(gen,i), sarch);
600 classe = isprincipalray(bnr, gel(gen,i));
601 for (ic = 1; ic <= nChar; ic++) {
602 GEN n = gel(nchi,ic);
603 n[i] = CHI_eval_n(lC[ic], classe);
604 }
605 }
606
607 /* Sum chi(beta) * exp(2i * Pi * Tr(beta * ilambda) where beta
608 runs through the classes of (Ok/cond0)^* and beta cond1-positive */
609 vt = gel(T,1); /* ( Tr(w_i) )_i */
610 if (typ(ilambda) == t_COL)
611 vt = ZV_ZM_mul(vt, zk_multable(nf, ilambda));
612 else
613 vt = ZC_Z_mul(vt, ilambda);
614 /*vt = den . (Tr(w_i * ilambda))_i */
615 G.cyc = gtovecsmall(cyc);
616 G.r = nz;
617 G.j = zero_zv(nz);
618 vN = zero_Flm_copy(nz, nChar);
619
620 av2 = avma;
621 vB = const_vec(nz, gen_1);
622 s0 = z? powgi(z, modii(gel(vt,1), den)): gen_1; /* for beta = 1 */
623 s = const_vec(nChar, s0);
624
625 while ( (i = NextElt(&G)) )
626 {
627 GEN b = gel(vB,i);
628 b = nfmuli(nf, b, gel(gen,i));
629 b = typ(b) == t_COL? FpC_red(b, condZ): modii(b, condZ);
630 for (j=1; j<=i; j++) gel(vB,j) = b;
631
632 for (ic = 1; ic <= nChar; ic++)
633 {
634 GEN v = gel(vN,ic), n = gel(nchi,ic);
635 v[i] = Fl_add(v[i], n[i], lC[ic]->ord);
636 for (j=1; j<i; j++) v[j] = v[i];
637 }
638
639 gel(vB,i) = b = set_sign_mod_divisor(nf, NULL, b, sarch);
640 if (!z)
641 s0 = gen_1;
642 else
643 {
644 b = typ(b) == t_COL? ZV_dotproduct(vt, b): mulii(gel(vt,1),b);
645 s0 = powgi(z, modii(b,den));
646 }
647 for (ic = 1; ic <= nChar; ic++)
648 {
649 GEN v = gel(vN,ic), val = lC[ic]->val[ v[i] ];
650 gel(s,ic) = gadd(gel(s,ic), gmul(val, s0));
651 }
652
653 if (gc_needed(av2, 1))
654 {
655 if (DEBUGMEM > 1) pari_warn(warnmem,"ArtinNumber");
656 gerepileall(av2, 2, &s, &vB);
657 }
658 }
659
660 classe = isprincipalray(bnr, idh);
661 z = powIs(- (lg(gel(sarch,1))-1));
662
663 for (ic = 1; ic <= nChar; ic++)
664 {
665 s0 = gmul(gel(s,ic), CHI_eval(lC[ic], classe));
666 gel(W, indW[ic]) = gmul(gdiv(s0, sqrtnc), z);
667 }
668 return W;
669 }
670
671 static GEN
AllArtinNumbers(GEN CR,long prec)672 AllArtinNumbers(GEN CR, long prec)
673 {
674 pari_sp av = avma;
675 GEN vChar = gel(CR,1), dataCR = gel(CR,2);
676 long j, k, cl = lg(dataCR) - 1, J = lg(vChar)-1;
677 GEN W = cgetg(cl+1,t_VEC), WbyCond, LCHI;
678
679 for (j = 1; j <= J; j++)
680 {
681 GEN LChar = gel(vChar,j), ldata = vecpermute(dataCR, LChar);
682 GEN dtcr = gel(ldata,1), bnr = ch_bnr(dtcr);
683 long l = lg(LChar);
684
685 if (DEBUGLEVEL>1)
686 err_printf("* Root Number: cond. no %ld/%ld (%ld chars)\n", j, J, l-1);
687 LCHI = cgetg(l, t_VEC);
688 for (k = 1; k < l; k++) gel(LCHI,k) = ch_CHI0(gel(ldata,k));
689 WbyCond = ArtinNumber(bnr, LCHI, prec);
690 for (k = 1; k < l; k++) gel(W,LChar[k]) = gel(WbyCond,k);
691 }
692 return gerepilecopy(av, W);
693 }
694
695 /* compute the constant W of the functional equation of
696 Lambda(chi). If flag = 1 then chi is assumed to be primitive */
697 GEN
bnrrootnumber(GEN bnr,GEN chi,long flag,long prec)698 bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
699 {
700 pari_sp av = avma;
701 GEN cyc, W;
702
703 if (flag < 0 || flag > 1) pari_err_FLAG("bnrrootnumber");
704 checkbnr(bnr);
705 if (flag)
706 {
707 cyc = bnr_get_cyc(bnr);
708 if (!char_check(cyc,chi)) pari_err_TYPE("bnrrootnumber [character]", chi);
709 }
710 else
711 {
712 bnr_char_sanitize(&bnr, &chi);
713 cyc = bnr_get_cyc(bnr);
714 }
715 chi = char_normalize(chi, cyc_normalize(cyc));
716 chi = get_Char(chi, prec);
717 W = ArtinNumber(bnr, mkvec(chi), prec);
718 return gerepilecopy(av, gel(W,1));
719 }
720
721 /********************************************************************/
722 /* 3rd part: initialize the characters */
723 /********************************************************************/
724
725 /* Let chi be a character, A(chi) corresponding to the primes dividing diff
726 * at s = flag. If s = 0, returns [r, A] where r is the order of vanishing
727 * at s = 0 corresponding to diff. */
728 static GEN
AChi(GEN dtcr,long * r,long flag,long prec)729 AChi(GEN dtcr, long *r, long flag, long prec)
730 {
731 GEN A, diff = ch_diff(dtcr), bnrc = ch_bnr(dtcr), chi = ch_CHI0(dtcr);
732 long i, l = lg(diff);
733
734 A = gen_1; *r = 0;
735 for (i = 1; i < l; i++)
736 {
737 GEN B, pr = gel(diff,i), z = CharEval(chi, isprincipalray(bnrc, pr));
738 if (flag)
739 B = gsubsg(1, gdiv(z, pr_norm(pr)));
740 else if (gequal1(z))
741 {
742 B = glog(pr_norm(pr), prec);
743 (*r)++;
744 }
745 else
746 B = gsubsg(1, z);
747 A = gmul(A, B);
748 }
749 return A;
750 }
751 /* simplified version of Achi: return 1 if L(0,chi) = 0 */
752 static int
L_vanishes_at_0(GEN D)753 L_vanishes_at_0(GEN D)
754 {
755 GEN diff = ch_diff(D), bnrc = ch_bnr(D), chi = ch_CHI0(D);
756 long i, l = lg(diff);
757 for (i = 1; i < l; i++)
758 {
759 GEN pr = gel(diff,i);
760 if (!CharEval_n(chi, isprincipalray(bnrc, pr))) return 1;
761 }
762 return 0;
763 }
764
765 static GEN
_data3(GEN arch,long r2)766 _data3(GEN arch, long r2)
767 {
768 GEN z = cgetg(4, t_VECSMALL);
769 long i, r1 = lg(arch) - 1, q = 0;
770 for (i = 1; i <= r1; i++) if (signe(gel(arch,i))) q++;
771 z[1] = q;
772 z[2] = r1 - q;
773 z[3] = r2; return z;
774 }
775 static void
ch_get3(GEN dtcr,long * a,long * b,long * c)776 ch_get3(GEN dtcr, long *a, long *b, long *c)
777 { GEN v = ch_3(dtcr); *a = v[1]; *b = v[2]; *c = v[3]; }
778 static GEN
get_C(GEN nf,long prec)779 get_C(GEN nf, long prec)
780 {
781 long r2 = nf_get_r2(nf), N = nf_get_degree(nf);
782 return gmul2n(sqrtr_abs(divir(nf_get_disc(nf), powru(mppi(prec),N))), -r2);
783 }
784 /* sort chars according to conductor */
785 static GEN
sortChars(GEN ch)786 sortChars(GEN ch)
787 {
788 long j, l = lg(ch);
789 GEN F = cgetg(l, t_VEC);
790 for (j = 1; j < l; j++) gel(F, j) = gmael(ch,j,2);
791 return vec_equiv(F);
792 }
793
794 /* Given a list [chi, F = cond(chi)] of characters over Cl(bnr), return
795 * [vChar, dataCR], where vChar containes the equivalence classes of
796 * characters with the same conductor, and dataCR contains for each character:
797 * - bnr(F)
798 * - the constant C(F) [t_REAL]
799 * - [q, r1 - q, r2, rc] where
800 * q = number of real places in F
801 * rc = max{r1 + r2 - q + 1, r2 + q}
802 * - diff(chi) primes dividing m but not F
803 * - chi in bnr(m)
804 * - chi in bnr(F).
805 * If all is unset, only compute characters s.t. L(chi,0) != 0 */
806 static GEN
InitChar(GEN bnr,GEN ch,long all,long prec)807 InitChar(GEN bnr, GEN ch, long all, long prec)
808 {
809 GEN bnf = checkbnf(bnr), nf = bnf_get_nf(bnf), mod = bnr_get_mod(bnr);
810 GEN C, dataCR, ncyc, vChar = sortChars(ch);
811 long n, l, r2 = nf_get_r2(nf), prec2 = precdbl(prec) + EXTRAPREC64;
812 long lv = lg(vChar);
813
814 C = get_C(nf, prec2);
815 ncyc = cyc_normalize(bnr_get_cyc(bnr));
816
817 dataCR = cgetg_copy(ch, &l);
818 for (n = 1; n < lv; n++)
819 {
820 GEN D, bnrc, v = gel(vChar, n); /* indices of chars of given conductor */
821 long a, i = v[1], lc = lg(v);
822 GEN F = gmael(ch,i,2);
823
824 gel(dataCR, i) = D = cgetg(8, t_VEC);
825 ch_C(D) = mulrr(C, gsqrt(ZM_det_triangular(gel(F,1)), prec2));
826 ch_3(D) = _data3(gel(F,2), r2);
827 if (gequal(F, mod))
828 {
829 ch_bnr(D) = bnrc = bnr;
830 ch_diff(D) = cgetg(1, t_VEC);
831 }
832 else
833 {
834 ch_bnr(D) = bnrc = Buchray(bnf, F, nf_INIT);
835 ch_diff(D) = get_prdiff(divcond(bnr), divcond(bnrc));
836 }
837 for (a = 1; a < lc; a++)
838 {
839 long i = v[a];
840 GEN chi = gmael(ch,i,1);
841
842 if (a > 1) gel(dataCR, i) = D = leafcopy(D);
843 chi = char_normalize(chi,ncyc);
844 ch_CHI(D) = get_Char(chi, prec2);
845 if (bnrc == bnr)
846 ch_CHI0(D) = ch_CHI(D);
847 else
848 {
849 chi = bnrchar_primitive(bnr, chi, bnrc);
850 ch_CHI0(D) = get_Char(chi, prec2);
851 }
852 /* set last */
853 ch_small(D) = mkvecsmall2(all || !L_vanishes_at_0(D),
854 eulerphiu(itou(gel(chi,1))));
855 }
856 }
857 return mkvec2(vChar, dataCR);
858 }
859
860 /* recompute dataCR with the new precision, modify bnr components in place */
861 static void
CharNewPrec(GEN data,long prec)862 CharNewPrec(GEN data, long prec)
863 {
864 long j, l, prec2 = precdbl(prec) + EXTRAPREC64;
865 GEN C, nf, dataCR = gmael(data,4,2), D = gel(dataCR,1);
866
867 if (ch_prec(D) >= prec2) return;
868 nf = bnr_get_nf(ch_bnr(D));
869 if (nf_get_prec(nf) < prec) nf = nfnewprec_shallow(nf, prec); /* not prec2 */
870 C = get_C(nf, prec2); l = lg(dataCR);
871 for (j = 1; j < l; j++)
872 {
873 GEN f0;
874 D = gel(dataCR,j); f0 = gel(bnr_get_mod(ch_bnr(D)), 1);
875 ch_C(D) = mulrr(C, gsqrt(ZM_det_triangular(f0), prec2));
876 gmael(ch_bnr(D), 1, 7) = nf;
877 ch_CHI(D) = get_Char(gel(ch_CHI(D),1), prec2);
878 ch_CHI0(D)= get_Char(gel(ch_CHI0(D),1), prec2);
879 }
880 }
881
882 /********************************************************************/
883 /* 4th part: compute the coefficients an(chi) */
884 /* */
885 /* matan entries are arrays of ints containing the coefficients of */
886 /* an(chi) as a polmod modulo polcyclo(order(chi)) */
887 /********************************************************************/
888
889 static void
_0toCoeff(int * rep,long deg)890 _0toCoeff(int *rep, long deg)
891 {
892 long i;
893 for (i=0; i<deg; i++) rep[i] = 0;
894 }
895
896 /* transform a polmod into Coeff */
897 static void
Polmod2Coeff(int * rep,GEN polmod,long deg)898 Polmod2Coeff(int *rep, GEN polmod, long deg)
899 {
900 long i;
901 if (typ(polmod) == t_POLMOD)
902 {
903 GEN pol = gel(polmod,2);
904 long d = degpol(pol);
905
906 pol += 2;
907 for (i=0; i<=d; i++) rep[i] = itos(gel(pol,i));
908 for ( ; i<deg; i++) rep[i] = 0;
909 }
910 else
911 {
912 rep[0] = itos(polmod);
913 for (i=1; i<deg; i++) rep[i] = 0;
914 }
915 }
916
917 /* initialize a deg * n matrix of ints */
918 static int**
InitMatAn(long n,long deg,long flag)919 InitMatAn(long n, long deg, long flag)
920 {
921 long i, j;
922 int *a, **A = (int**)pari_malloc((n+1)*sizeof(int*));
923 A[0] = NULL;
924 for (i = 1; i <= n; i++)
925 {
926 a = (int*)pari_malloc(deg*sizeof(int));
927 A[i] = a; a[0] = (i == 1 || flag);
928 for (j = 1; j < deg; j++) a[j] = 0;
929 }
930 return A;
931 }
932
933 static void
FreeMat(int ** A,long n)934 FreeMat(int **A, long n)
935 {
936 long i;
937 for (i = 0; i <= n; i++)
938 if (A[i]) pari_free((void*)A[i]);
939 pari_free((void*)A);
940 }
941
942 /* initialize Coeff reduction */
943 static int**
InitReduction(long d,long deg)944 InitReduction(long d, long deg)
945 {
946 long j;
947 pari_sp av = avma;
948 int **A;
949 GEN polmod, pol;
950
951 A = (int**)pari_malloc(deg*sizeof(int*));
952 pol = polcyclo(d, 0);
953 for (j = 0; j < deg; j++)
954 {
955 A[j] = (int*)pari_malloc(deg*sizeof(int));
956 polmod = gmodulo(pol_xn(deg+j, 0), pol);
957 Polmod2Coeff(A[j], polmod, deg);
958 }
959
960 set_avma(av); return A;
961 }
962
963 #if 0
964 void
965 pan(int **an, long n, long deg)
966 {
967 long i,j;
968 for (i = 1; i <= n; i++)
969 {
970 err_printf("n = %ld: ",i);
971 for (j = 0; j < deg; j++) err_printf("%d ",an[i][j]);
972 err_printf("\n");
973 }
974 }
975 #endif
976
977 /* returns 0 if c is zero, 1 otherwise. */
978 static int
IsZero(int * c,long deg)979 IsZero(int* c, long deg)
980 {
981 long i;
982 for (i = 0; i < deg; i++)
983 if (c[i]) return 0;
984 return 1;
985 }
986
987 /* set c0 <-- c0 * c1 */
988 static void
MulCoeff(int * c0,int * c1,int ** reduc,long deg)989 MulCoeff(int *c0, int* c1, int** reduc, long deg)
990 {
991 long i,j;
992 int c, *T;
993
994 if (IsZero(c0,deg)) return;
995
996 T = (int*)new_chunk(2*deg);
997 for (i = 0; i < 2*deg; i++)
998 {
999 c = 0;
1000 for (j = 0; j <= i; j++)
1001 if (j < deg && j > i - deg) c += c0[j] * c1[i-j];
1002 T[i] = c;
1003 }
1004 for (i = 0; i < deg; i++)
1005 {
1006 c = T[i];
1007 for (j = 0; j < deg; j++) c += reduc[j][i] * T[deg+j];
1008 c0[i] = c;
1009 }
1010 }
1011
1012 /* c0 <- c0 + c1 * c2 */
1013 static void
AddMulCoeff(int * c0,int * c1,int * c2,int ** reduc,long deg)1014 AddMulCoeff(int *c0, int *c1, int* c2, int** reduc, long deg)
1015 {
1016 long i, j;
1017 pari_sp av;
1018 int c, *t;
1019
1020 if (IsZero(c2,deg)) return;
1021 if (!c1) /* c1 == 1 */
1022 {
1023 for (i = 0; i < deg; i++) c0[i] += c2[i];
1024 return;
1025 }
1026 av = avma;
1027 t = (int*)new_chunk(2*deg); /* = c1 * c2, not reduced */
1028 for (i = 0; i < 2*deg; i++)
1029 {
1030 c = 0;
1031 for (j = 0; j <= i; j++)
1032 if (j < deg && j > i - deg) c += c1[j] * c2[i-j];
1033 t[i] = c;
1034 }
1035 for (i = 0; i < deg; i++)
1036 {
1037 c = t[i];
1038 for (j = 0; j < deg; j++) c += reduc[j][i] * t[deg+j];
1039 c0[i] += c;
1040 }
1041 set_avma(av);
1042 }
1043
1044 /* evaluate the Coeff. No Garbage collector */
1045 static GEN
EvalCoeff(GEN z,int * c,long deg)1046 EvalCoeff(GEN z, int* c, long deg)
1047 {
1048 long i,j;
1049 GEN e, r;
1050
1051 if (!c) return gen_0;
1052 #if 0
1053 /* standard Horner */
1054 e = stoi(c[deg - 1]);
1055 for (i = deg - 2; i >= 0; i--)
1056 e = gadd(stoi(c[i]), gmul(z, e));
1057 #else
1058 /* specific attention to sparse polynomials */
1059 e = NULL;
1060 for (i = deg-1; i >=0; i=j-1)
1061 {
1062 for (j=i; c[j] == 0; j--)
1063 if (j==0)
1064 {
1065 if (!e) return NULL;
1066 if (i!=j) z = gpowgs(z,i-j+1);
1067 return gmul(e,z);
1068 }
1069 if (e)
1070 {
1071 r = (i==j)? z: gpowgs(z,i-j+1);
1072 e = gadd(gmul(e,r), stoi(c[j]));
1073 }
1074 else
1075 e = stoi(c[j]);
1076 }
1077 #endif
1078 return e;
1079 }
1080
1081 /* copy the n * (m+1) array matan */
1082 static void
CopyCoeff(int ** a,int ** a2,long n,long m)1083 CopyCoeff(int** a, int** a2, long n, long m)
1084 {
1085 long i,j;
1086
1087 for (i = 1; i <= n; i++)
1088 {
1089 int *b = a[i], *b2 = a2[i];
1090 for (j = 0; j < m; j++) b2[j] = b[j];
1091 }
1092 }
1093
1094 static void
an_AddMul(int ** an,int ** an2,long np,long n,long deg,GEN chi,int ** reduc)1095 an_AddMul(int **an,int **an2, long np, long n, long deg, GEN chi, int **reduc)
1096 {
1097 GEN chi2 = chi;
1098 long q, qk, k;
1099 int *c, *c2 = (int*)new_chunk(deg);
1100
1101 CopyCoeff(an, an2, n/np, deg);
1102 for (q=np;;)
1103 {
1104 if (gequal1(chi2)) c = NULL; else { Polmod2Coeff(c2, chi2, deg); c = c2; }
1105 for(k = 1, qk = q; qk <= n; k++, qk += q)
1106 AddMulCoeff(an[qk], c, an2[k], reduc, deg);
1107 if (! (q = umuluu_le(q,np, n)) ) break;
1108
1109 chi2 = gmul(chi2, chi);
1110 }
1111 }
1112
1113 /* correct the coefficients an(chi) according with diff(chi) in place */
1114 static void
CorrectCoeff(GEN dtcr,int ** an,int ** reduc,long n,long deg)1115 CorrectCoeff(GEN dtcr, int** an, int** reduc, long n, long deg)
1116 {
1117 pari_sp av = avma;
1118 long lg, j;
1119 pari_sp av1;
1120 int **an2;
1121 GEN bnrc, diff;
1122 CHI_t C;
1123
1124 diff = ch_diff(dtcr); lg = lg(diff) - 1;
1125 if (!lg) return;
1126
1127 if (DEBUGLEVEL>2) err_printf("diff(CHI) = %Ps", diff);
1128 bnrc = ch_bnr(dtcr);
1129 init_CHI_alg(&C, ch_CHI0(dtcr));
1130
1131 an2 = InitMatAn(n, deg, 0);
1132 av1 = avma;
1133 for (j = 1; j <= lg; j++)
1134 {
1135 GEN pr = gel(diff,j);
1136 long Np = upr_norm(pr);
1137 GEN chi = CHI_eval(&C, isprincipalray(bnrc, pr));
1138 an_AddMul(an,an2,Np,n,deg,chi,reduc);
1139 set_avma(av1);
1140 }
1141 FreeMat(an2, n); set_avma(av);
1142 }
1143
1144 /* compute the coefficients an in the general case */
1145 static int**
ComputeCoeff(GEN dtcr,LISTray * R,long n,long deg)1146 ComputeCoeff(GEN dtcr, LISTray *R, long n, long deg)
1147 {
1148 pari_sp av = avma, av2;
1149 long i, l;
1150 int **an, **reduc, **an2;
1151 GEN L;
1152 CHI_t C;
1153
1154 init_CHI_alg(&C, ch_CHI(dtcr));
1155 an = InitMatAn(n, deg, 0);
1156 an2 = InitMatAn(n, deg, 0);
1157 reduc = InitReduction(C.ord, deg);
1158 av2 = avma;
1159
1160 L = R->L1; l = lg(L);
1161 for (i=1; i<l; i++, set_avma(av2))
1162 {
1163 long np = L[i];
1164 GEN chi = CHI_eval(&C, gel(R->L1ray,i));
1165 an_AddMul(an,an2,np,n,deg,chi,reduc);
1166 }
1167 FreeMat(an2, n);
1168
1169 CorrectCoeff(dtcr, an, reduc, n, deg);
1170 FreeMat(reduc, deg-1);
1171 set_avma(av); return an;
1172 }
1173
1174 /********************************************************************/
1175 /* 5th part: compute L-functions at s=1 */
1176 /********************************************************************/
1177 static void
deg11(LISTray * R,long p,GEN bnr,GEN pr)1178 deg11(LISTray *R, long p, GEN bnr, GEN pr) {
1179 GEN z = isprincipalray(bnr, pr);
1180 vecsmalltrunc_append(R->L1, p);
1181 vectrunc_append(R->L1ray, z);
1182 }
1183 static void
deg12(LISTray * R,long p,GEN bnr,GEN Lpr)1184 deg12(LISTray *R, long p, GEN bnr, GEN Lpr) {
1185 GEN z = isprincipalray(bnr, gel(Lpr,1));
1186 vecsmalltrunc_append(R->L11, p);
1187 vectrunc_append(R->L11ray, z);
1188 }
1189 static void
deg0(LISTray * R,long p)1190 deg0(LISTray *R, long p) { vecsmalltrunc_append(R->L0, p); }
1191 static void
deg2(LISTray * R,long p)1192 deg2(LISTray *R, long p) { vecsmalltrunc_append(R->L2, p); }
1193
1194 static void
InitPrimesQuad(GEN bnr,ulong N0,LISTray * R)1195 InitPrimesQuad(GEN bnr, ulong N0, LISTray *R)
1196 {
1197 pari_sp av = avma;
1198 GEN bnf = bnr_get_bnf(bnr), cond = gel(bnr_get_mod(bnr), 1);
1199 long p,i,l, condZ = itos(gcoeff(cond,1,1)), contZ = itos(content(cond));
1200 GEN prime, Lpr, nf = bnf_get_nf(bnf), dk = nf_get_disc(nf);
1201 forprime_t T;
1202
1203 l = 1 + primepi_upper_bound(N0);
1204 R->L0 = vecsmalltrunc_init(l);
1205 R->L2 = vecsmalltrunc_init(l); R->condZ = condZ;
1206 R->L1 = vecsmalltrunc_init(l); R->L1ray = vectrunc_init(l);
1207 R->L11= vecsmalltrunc_init(l); R->L11ray= vectrunc_init(l);
1208 prime = utoipos(2);
1209 u_forprime_init(&T, 2, N0);
1210 while ( (p = u_forprime_next(&T)) )
1211 {
1212 prime[2] = p;
1213 switch (kroiu(dk, p))
1214 {
1215 case -1: /* inert */
1216 if (condZ % p == 0) deg0(R,p); else deg2(R,p);
1217 break;
1218 case 1: /* split */
1219 Lpr = idealprimedec(nf, prime);
1220 if (condZ % p != 0) deg12(R, p, bnr, Lpr);
1221 else if (contZ % p == 0) deg0(R,p);
1222 else
1223 {
1224 GEN pr = idealval(nf, cond, gel(Lpr,1))? gel(Lpr,2): gel(Lpr,1);
1225 deg11(R, p, bnr, pr);
1226 }
1227 break;
1228 default: /* ramified */
1229 if (condZ % p == 0)
1230 deg0(R,p);
1231 else
1232 deg11(R, p, bnr, idealprimedec_galois(nf,prime));
1233 break;
1234 }
1235 }
1236 /* precompute isprincipalray(x), x in Z */
1237 R->rayZ = cgetg(condZ, t_VEC);
1238 for (i=1; i<condZ; i++)
1239 gel(R->rayZ,i) = (ugcd(i,condZ) == 1)? isprincipalray(bnr, utoipos(i)): gen_0;
1240 gerepileall(av, 7, &(R->L0), &(R->L2), &(R->rayZ),
1241 &(R->L1), &(R->L1ray), &(R->L11), &(R->L11ray) );
1242 }
1243
1244 static void
InitPrimes(GEN bnr,ulong N0,LISTray * R)1245 InitPrimes(GEN bnr, ulong N0, LISTray *R)
1246 {
1247 GEN bnf = bnr_get_bnf(bnr), cond = gel(bnr_get_mod(bnr), 1);
1248 long p, l, condZ, N = lg(cond)-1;
1249 GEN DL, prime, BOUND, nf = bnf_get_nf(bnf);
1250 forprime_t T;
1251
1252 R->condZ = condZ = itos(gcoeff(cond,1,1));
1253 l = primepi_upper_bound(N0) * N;
1254 DL = cgetg(N+1, t_VEC);
1255 R->L1 = vecsmalltrunc_init(l);
1256 R->L1ray = vectrunc_init(l);
1257 u_forprime_init(&T, 2, N0);
1258 prime = utoipos(2);
1259 BOUND = utoi(N0);
1260 while ( (p = u_forprime_next(&T)) )
1261 {
1262 pari_sp av = avma;
1263 long j, k, lP;
1264 GEN P;
1265 prime[2] = p;
1266 if (DEBUGLEVEL>1 && (p & 2047) == 1) err_printf("%ld ", p);
1267 P = idealprimedec_limit_norm(nf, prime, BOUND); lP = lg(P);
1268 for (j = 1; j < lP; j++)
1269 {
1270 GEN pr = gel(P,j), dl = NULL;
1271 if (condZ % p || !idealval(nf, cond, pr))
1272 {
1273 dl = gclone( isprincipalray(bnr, pr) );
1274 vecsmalltrunc_append(R->L1, upowuu(p, pr_get_f(pr)));
1275 }
1276 gel(DL,j) = dl;
1277 }
1278 set_avma(av);
1279 for (k = 1; k < j; k++)
1280 {
1281 if (!DL[k]) continue;
1282 vectrunc_append(R->L1ray, ZC_copy(gel(DL,k)));
1283 gunclone(gel(DL,k));
1284 }
1285 }
1286 }
1287
1288 static GEN /* cf polcoef */
_sercoeff(GEN x,long n)1289 _sercoeff(GEN x, long n)
1290 {
1291 long i = n - valp(x);
1292 return (i < 0)? gen_0: gel(x,i+2);
1293 }
1294
1295 static void
affect_coeff(GEN q,long n,GEN y,long t)1296 affect_coeff(GEN q, long n, GEN y, long t)
1297 {
1298 GEN x = _sercoeff(q,-n);
1299 if (x == gen_0) gel(y,n) = NULL;
1300 else { affgr(x, gel(y,n)); shiftr_inplace(gel(y,n), t); }
1301 }
1302 /* (x-i)(x-(i+1)) */
1303 static GEN
d2(long i)1304 d2(long i) { return deg2pol_shallow(gen_1, utoineg(2*i+1), muluu(i,i+1), 0); }
1305 /* x-i */
1306 static GEN
d1(long i)1307 d1(long i) { return deg1pol_shallow(gen_1, stoi(-i), 0); }
1308
1309 typedef struct {
1310 GEN c1, aij, bij, cS, cT, powracpi;
1311 long i0, a,b,c, r, rc1, rc2;
1312 } ST_t;
1313
1314 /* compute the principal part at the integers s = 0, -1, -2, ..., -i0
1315 * of Gamma((s+1)/2)^a Gamma(s/2)^b Gamma(s)^c / (s - z) with z = 0 and 1 */
1316 static void
ppgamma(ST_t * T,long prec)1317 ppgamma(ST_t *T, long prec)
1318 {
1319 GEN G, G1, G2, A, E, O, x, sqpi, aij, bij;
1320 long c = T->c, r = T->r, i0 = T->i0, i, j, s, t, dx;
1321 pari_sp av;
1322
1323 T->aij = aij = cgetg(i0+1, t_VEC);
1324 T->bij = bij = cgetg(i0+1, t_VEC);
1325 for (i = 1; i <= i0; i++)
1326 {
1327 GEN p1, p2;
1328 gel(aij,i) = p1 = cgetg(r+1, t_VEC);
1329 gel(bij,i) = p2 = cgetg(r+1, t_VEC);
1330 for (j=1; j<=r; j++) { gel(p1,j) = cgetr(prec); gel(p2,j) = cgetr(prec); }
1331 }
1332 av = avma; x = pol_x(0);
1333 sqpi = sqrtr_abs(mppi(prec)); /* Gamma(1/2) */
1334
1335 G1 = gexp(integser(psi1series(r-1, 0, prec)), prec); /* Gamma(1 + x) */
1336 G = shallowcopy(G1); setvalp(G,-1); /* Gamma(x) */
1337
1338 /* expansion of log(Gamma(u) / Gamma(1/2)) at u = 1/2 */
1339 G2 = cgetg(r+2, t_SER);
1340 G2[1] = evalsigne(1) | _evalvalp(1) | evalvarn(0);
1341 gel(G2,2) = gneg(gadd(gmul2n(mplog2(prec), 1), mpeuler(prec)));
1342 for (i = 1; i < r; i++) gel(G2,i+2) = mulri(gel(G1,i+2), int2um1(i));
1343 G2 = gmul(sqpi, gexp(G2, prec)); /* Gamma(1/2 + x) */
1344
1345 /* We simplify to get one of the following two expressions
1346 * if (b > a) : sqrt(Pi)^a 2^{a-au} Gamma(u)^{a+c} Gamma( u/2 )^{|b-a|}
1347 * if (b <= a): sqrt(Pi)^b 2^{b-bu} Gamma(u)^{b+c} Gamma((u+1)/2)^{|b-a|} */
1348 if (T->b > T->a)
1349 {
1350 t = T->a; s = T->b; dx = 1;
1351 E = ser_unscale(G, ghalf);
1352 O = gmul2n(gdiv(ser_unscale(G2, ghalf), d1(1)), 1); /* Gamma((x-1)/2) */
1353 }
1354 else
1355 {
1356 t = T->b; s = T->a; dx = 0;
1357 E = ser_unscale(G2, ghalf);
1358 O = ser_unscale(G, ghalf);
1359 }
1360 /* (sqrt(Pi) 2^{1-x})^t Gamma(x)^{t+c} */
1361 A = gmul(gmul(powru(gmul2n(sqpi,1), t), gpowgs(G, t+c)),
1362 gpow(gen_2, RgX_to_ser(gmulgs(x,-t), r+2), prec));
1363 /* A * Gamma((x - dx + 1)/2)^{s-t} */
1364 E = gmul(A, gpowgs(E, s-t));
1365 /* A * Gamma((x - dx)/2)^{s-t} */
1366 O = gdiv(gmul(A, gpowgs(O, s-t)), gpowgs(gsubgs(x, 1), t+c));
1367 for (i = 0; i < i0/2; i++)
1368 {
1369 GEN p1, q1, A1 = gel(aij,2*i+1), B1 = gel(bij,2*i+1);
1370 GEN p2, q2, A2 = gel(aij,2*i+2), B2 = gel(bij,2*i+2);
1371 long t1 = i * (s+t), t2 = t1 + t;
1372
1373 p1 = gdiv(E, d1(2*i));
1374 q1 = gdiv(E, d1(2*i+1));
1375 p2 = gdiv(O, d1(2*i+1));
1376 q2 = gdiv(O, d1(2*i+2));
1377 for (j = 1; j <= r; j++)
1378 {
1379 affect_coeff(p1, j, A1, t1); affect_coeff(q1, j, B1, t1);
1380 affect_coeff(p2, j, A2, t2); affect_coeff(q2, j, B2, t2);
1381 }
1382 E = gdiv(E, gmul(gpowgs(d1(2*i+1+dx), s-t), gpowgs(d2(2*i+1), t+c)));
1383 O = gdiv(O, gmul(gpowgs(d1(2*i+2+dx), s-t), gpowgs(d2(2*i+2), t+c)));
1384 }
1385 set_avma(av);
1386 }
1387
1388 /* chi != 1. Return L(1, chi) if fl & 1, else [r, c] where L(s, chi) ~ c s^r
1389 * at s = 0. */
1390 static GEN
GetValue(GEN dtcr,GEN W,GEN S,GEN T,long fl,long prec)1391 GetValue(GEN dtcr, GEN W, GEN S, GEN T, long fl, long prec)
1392 {
1393 GEN cf, z;
1394 long q, b, c, r;
1395 int isreal = (ch_deg(dtcr) <= 2);
1396
1397 ch_get3(dtcr, &q, &b, &c);
1398 if (fl & 1)
1399 { /* S(chi) + W(chi).T(chi)) / (C(chi) sqrt(Pi)^{r1 - q}) */
1400 cf = gmul(ch_C(dtcr), powruhalf(mppi(prec), b));
1401
1402 z = gadd(S, gmul(W, T));
1403 if (isreal) z = real_i(z);
1404 z = gdiv(z, cf);
1405 if (fl & 2) z = gmul(z, AChi(dtcr, &r, 1, prec));
1406 }
1407 else
1408 { /* (W(chi).S(conj(chi)) + T(chi)) / (sqrt(Pi)^q 2^{r1 - q}) */
1409 cf = gmul2n(powruhalf(mppi(prec), q), b);
1410
1411 z = gadd(gmul(W, conj_i(S)), conj_i(T));
1412 if (isreal) z = real_i(z);
1413 z = gdiv(z, cf); r = 0;
1414 if (fl & 2) z = gmul(z, AChi(dtcr, &r, 0, prec));
1415 z = mkvec2(utoi(b + c + r), z);
1416 }
1417 return z;
1418 }
1419
1420 /* return the order and the first nonzero term of L(s, chi0)
1421 at s = 0. If flag != 0, adjust the value to get L_S(s, chi0). */
1422 static GEN
GetValue1(GEN bnr,long flag,long prec)1423 GetValue1(GEN bnr, long flag, long prec)
1424 {
1425 GEN bnf = bnr_get_bnf(bnr), nf = bnf_get_nf(bnf);
1426 GEN h = bnf_get_no(bnf), R = bnf_get_reg(bnf);
1427 GEN c = gdivgs(mpmul(h, R), -bnf_get_tuN(bnf));
1428 long r = lg(nf_get_roots(nf)) - 2; /* r1 + r2 - 1 */;
1429 if (flag)
1430 {
1431 GEN diff = divcond(bnr);
1432 long i, l;
1433 l = lg(diff) - 1; r += l;
1434 for (i = 1; i <= l; i++) c = gmul(c, glog(pr_norm(gel(diff,i)), prec));
1435 }
1436 return mkvec2(utoi(r), c);
1437 }
1438
1439 /********************************************************************/
1440 /* 6th part: recover the coefficients */
1441 /********************************************************************/
1442 static long
TestOne(GEN plg,RC_data * d)1443 TestOne(GEN plg, RC_data *d)
1444 {
1445 long j, v = d->v;
1446 GEN z = gsub(d->beta, gel(plg,v));
1447 if (expo(z) >= d->G) return 0;
1448 for (j = 1; j < lg(plg); j++)
1449 if (j != v && mpcmp(d->B, mpabs_shallow(gel(plg,j))) < 0) return 0;
1450 return 1;
1451 }
1452
1453 static GEN
chk_reccoeff_init(FP_chk_fun * chk,GEN r,GEN mat)1454 chk_reccoeff_init(FP_chk_fun *chk, GEN r, GEN mat)
1455 {
1456 RC_data *d = (RC_data*)chk->data;
1457 (void)r; d->U = mat; return d->nB;
1458 }
1459
1460 static GEN
chk_reccoeff(void * data,GEN x)1461 chk_reccoeff(void *data, GEN x)
1462 {
1463 RC_data *d = (RC_data*)data;
1464 GEN v = gmul(d->U, x), z = gel(v,1);
1465
1466 if (!gequal1(z)) return NULL;
1467 *++v = evaltyp(t_COL) | evallg( lg(d->M) );
1468 if (TestOne(gmul(d->M, v), d)) return v;
1469 return NULL;
1470 }
1471
1472 /* Using Cohen's method */
1473 static GEN
RecCoeff3(GEN nf,RC_data * d,long prec)1474 RecCoeff3(GEN nf, RC_data *d, long prec)
1475 {
1476 GEN A, M, nB, cand, p1, B2, C2, tB, beta2, nf2, Bd;
1477 GEN beta = d->beta, B = d->B;
1478 long N = d->N, v = d->v, e, BIG;
1479 long i, j, k, ct = 0, prec2;
1480 FP_chk_fun chk = { &chk_reccoeff, &chk_reccoeff_init, NULL, NULL, 0 };
1481 chk.data = (void*)d;
1482
1483 d->G = minss(-10, -prec2nbits(prec) >> 4);
1484 BIG = maxss(32, -2*d->G);
1485 tB = sqrtnr(real2n(BIG-N,DEFAULTPREC), N-1);
1486 Bd = grndtoi(gmin_shallow(B, tB), &e);
1487 if (e > 0) return NULL; /* failure */
1488 Bd = addiu(Bd, 1);
1489 prec2 = nbits2prec( expi(Bd) + 192 );
1490 prec2 = maxss(precdbl(prec), prec2);
1491 B2 = sqri(Bd);
1492 C2 = shifti(B2, BIG<<1);
1493
1494 LABrcf: ct++;
1495 beta2 = gprec_w(beta, prec2);
1496 nf2 = nfnewprec_shallow(nf, prec2);
1497 d->M = M = nf_get_M(nf2);
1498
1499 A = cgetg(N+2, t_MAT);
1500 for (i = 1; i <= N+1; i++) gel(A,i) = cgetg(N+2, t_COL);
1501
1502 gcoeff(A, 1, 1) = gadd(gmul(C2, gsqr(beta2)), B2);
1503 for (j = 2; j <= N+1; j++)
1504 {
1505 p1 = gmul(C2, gmul(gneg_i(beta2), gcoeff(M, v, j-1)));
1506 gcoeff(A, 1, j) = gcoeff(A, j, 1) = p1;
1507 }
1508 for (i = 2; i <= N+1; i++)
1509 for (j = i; j <= N+1; j++)
1510 {
1511 p1 = gen_0;
1512 for (k = 1; k <= N; k++)
1513 {
1514 GEN p2 = gmul(gcoeff(M, k, j-1), gcoeff(M, k, i-1));
1515 if (k == v) p2 = gmul(C2, p2);
1516 p1 = gadd(p1,p2);
1517 }
1518 gcoeff(A, i, j) = gcoeff(A, j, i) = p1;
1519 }
1520
1521 nB = mului(N+1, B2);
1522 d->nB = nB;
1523 cand = fincke_pohst(A, nB, -1, prec2, &chk);
1524
1525 if (!cand)
1526 {
1527 if (ct > 3) return NULL;
1528 prec2 = precdbl(prec2);
1529 if (DEBUGLEVEL>1) pari_warn(warnprec,"RecCoeff", prec2);
1530 goto LABrcf;
1531 }
1532
1533 cand = gel(cand,1);
1534 if (lg(cand) == 2) return gel(cand,1);
1535
1536 if (DEBUGLEVEL>1) err_printf("RecCoeff3: no solution found!\n");
1537 return NULL;
1538 }
1539
1540 /* Using linear dependance relations */
1541 static GEN
RecCoeff2(GEN nf,RC_data * d,long prec)1542 RecCoeff2(GEN nf, RC_data *d, long prec)
1543 {
1544 pari_sp av;
1545 GEN vec, M = nf_get_M(nf), beta = d->beta;
1546 long bit, min, max, lM = lg(M);
1547
1548 d->G = minss(-20, -prec2nbits(prec) >> 4);
1549
1550 vec = shallowconcat(mkvec(gneg(beta)), row(M, d->v));
1551 min = (long)prec2nbits_mul(prec, 0.75);
1552 max = (long)prec2nbits_mul(prec, 0.98);
1553 av = avma;
1554 for (bit = max; bit >= min; bit-=32, set_avma(av))
1555 {
1556 long e;
1557 GEN v = lindep_bit(vec, bit), z = gel(v,1);
1558 if (!signe(z)) continue;
1559 *++v = evaltyp(t_COL) | evallg(lM);
1560 v = grndtoi(gdiv(v, z), &e);
1561 if (e > 0) break;
1562 if (TestOne(RgM_RgC_mul(M, v), d)) return v;
1563 }
1564 /* failure */
1565 return RecCoeff3(nf,d,prec);
1566 }
1567
1568 /* Attempts to find a polynomial with coefficients in nf such that
1569 its coefficients are close to those of pol at the place v and
1570 less than B at all the other places */
1571 static GEN
RecCoeff(GEN nf,GEN pol,long v,long prec)1572 RecCoeff(GEN nf, GEN pol, long v, long prec)
1573 {
1574 long j, md, cl = degpol(pol);
1575 pari_sp av = avma;
1576 RC_data d;
1577
1578 /* if precision(pol) is too low, abort */
1579 for (j = 2; j <= cl+1; j++)
1580 {
1581 GEN t = gel(pol, j);
1582 if (prec2nbits(precision(t)) - gexpo(t) < 34) return NULL;
1583 }
1584
1585 md = cl/2;
1586 pol = leafcopy(pol);
1587
1588 d.N = nf_get_degree(nf);
1589 d.v = v;
1590
1591 for (j = 1; j <= cl; j++)
1592 { /* start with the coefficients in the middle,
1593 since they are the harder to recognize! */
1594 long cf = md + (j%2? j/2: -j/2);
1595 GEN t, bound = shifti(binomial(utoipos(cl), cf), cl-cf);
1596
1597 if (DEBUGLEVEL>1) err_printf("RecCoeff (cf = %ld, B = %Ps)\n", cf, bound);
1598 d.beta = real_i( gel(pol,cf+2) );
1599 d.B = bound;
1600 if (! (t = RecCoeff2(nf, &d, prec)) ) return NULL;
1601 gel(pol, cf+2) = coltoalg(nf,t);
1602 }
1603 gel(pol,cl+2) = gen_1;
1604 return gerepilecopy(av, pol);
1605 }
1606
1607 /* an[q * i] *= chi for all (i,p)=1 */
1608 static void
an_mul(int ** an,long p,long q,long n,long deg,GEN chi,int ** reduc)1609 an_mul(int **an, long p, long q, long n, long deg, GEN chi, int **reduc)
1610 {
1611 pari_sp av;
1612 long c,i;
1613 int *T;
1614
1615 if (gequal1(chi)) return;
1616 av = avma;
1617 T = (int*)new_chunk(deg); Polmod2Coeff(T,chi, deg);
1618 for (c = 1, i = q; i <= n; i += q, c++)
1619 if (c == p) c = 0; else MulCoeff(an[i], T, reduc, deg);
1620 set_avma(av);
1621 }
1622 /* an[q * i] = 0 for all (i,p)=1 */
1623 static void
an_set0_coprime(int ** an,long p,long q,long n,long deg)1624 an_set0_coprime(int **an, long p, long q, long n, long deg)
1625 {
1626 long c,i;
1627 for (c = 1, i = q; i <= n; i += q, c++)
1628 if (c == p) c = 0; else _0toCoeff(an[i], deg);
1629 }
1630 /* an[q * i] = 0 for all i */
1631 static void
an_set0(int ** an,long p,long n,long deg)1632 an_set0(int **an, long p, long n, long deg)
1633 {
1634 long i;
1635 for (i = p; i <= n; i += p) _0toCoeff(an[i], deg);
1636 }
1637
1638 /* compute the coefficients an for the quadratic case */
1639 static int**
computean(GEN dtcr,LISTray * R,long n,long deg)1640 computean(GEN dtcr, LISTray *R, long n, long deg)
1641 {
1642 pari_sp av = avma, av2;
1643 long i, p, q, condZ, l;
1644 int **an, **reduc;
1645 GEN L, chi, chi1;
1646 CHI_t C;
1647
1648 init_CHI_alg(&C, ch_CHI(dtcr));
1649 condZ= R->condZ;
1650
1651 an = InitMatAn(n, deg, 1);
1652 reduc = InitReduction(C.ord, deg);
1653 av2 = avma;
1654
1655 /* all pr | p divide cond */
1656 L = R->L0; l = lg(L);
1657 for (i=1; i<l; i++) an_set0(an,L[i],n,deg);
1658
1659 /* 1 prime of degree 2 */
1660 L = R->L2; l = lg(L);
1661 for (i=1; i<l; i++, set_avma(av2))
1662 {
1663 p = L[i];
1664 if (condZ == 1) chi = C.val[0]; /* 1 */
1665 else chi = CHI_eval(&C, gel(R->rayZ, p%condZ));
1666 chi1 = chi;
1667 for (q=p;;)
1668 {
1669 an_set0_coprime(an, p,q,n,deg); /* v_p(q) odd */
1670 if (! (q = umuluu_le(q,p, n)) ) break;
1671
1672 an_mul(an,p,q,n,deg,chi,reduc);
1673 if (! (q = umuluu_le(q,p, n)) ) break;
1674 chi = gmul(chi, chi1);
1675 }
1676 }
1677
1678 /* 1 prime of degree 1 */
1679 L = R->L1; l = lg(L);
1680 for (i=1; i<l; i++, set_avma(av2))
1681 {
1682 p = L[i];
1683 chi = CHI_eval(&C, gel(R->L1ray,i));
1684 chi1 = chi;
1685 for(q=p;;)
1686 {
1687 an_mul(an,p,q,n,deg,chi,reduc);
1688 if (! (q = umuluu_le(q,p, n)) ) break;
1689 chi = gmul(chi, chi1);
1690 }
1691 }
1692
1693 /* 2 primes of degree 1 */
1694 L = R->L11; l = lg(L);
1695 for (i=1; i<l; i++, set_avma(av2))
1696 {
1697 GEN ray1, ray2, chi11, chi12, chi2;
1698
1699 p = L[i]; ray1 = gel(R->L11ray,i); /* use pr1 pr2 = (p) */
1700 if (condZ == 1)
1701 ray2 = ZC_neg(ray1);
1702 else
1703 ray2 = ZC_sub(gel(R->rayZ, p%condZ), ray1);
1704 chi11 = CHI_eval(&C, ray1);
1705 chi12 = CHI_eval(&C, ray2);
1706
1707 chi1 = gadd(chi11, chi12);
1708 chi2 = chi12;
1709 for(q=p;;)
1710 {
1711 an_mul(an,p,q,n,deg,chi1,reduc);
1712 if (! (q = umuluu_le(q,p, n)) ) break;
1713 chi2 = gmul(chi2, chi12);
1714 chi1 = gadd(chi2, gmul(chi1, chi11));
1715 }
1716 }
1717
1718 CorrectCoeff(dtcr, an, reduc, n, deg);
1719 FreeMat(reduc, deg-1);
1720 set_avma(av); return an;
1721 }
1722
1723 /* return the vector of A^i/i for i = 1...n */
1724 static GEN
mpvecpowdiv(GEN A,long n)1725 mpvecpowdiv(GEN A, long n)
1726 {
1727 pari_sp av = avma;
1728 long i;
1729 GEN v = powersr(A, n);
1730 GEN w = cgetg(n+1, t_VEC);
1731 gel(w,1) = rcopy(gel(v,2));
1732 for (i=2; i<=n; i++) gel(w,i) = divru(gel(v,i+1), i);
1733 return gerepileupto(av, w);
1734 }
1735
1736 static void GetST0(GEN bnr, GEN *pS, GEN *pT, GEN CR, long prec);
1737 /* allocate memory for GetST answer */
1738 static void
ST_alloc(GEN * pS,GEN * pT,long l,long prec)1739 ST_alloc(GEN *pS, GEN *pT, long l, long prec)
1740 {
1741 long j;
1742 *pS = cgetg(l, t_VEC);
1743 *pT = cgetg(l, t_VEC);
1744 for (j = 1; j < l; j++)
1745 {
1746 gel(*pS,j) = cgetc(prec);
1747 gel(*pT,j) = cgetc(prec);
1748 }
1749 }
1750
1751 /* compute S and T for the quadratic case. The following cases are:
1752 * 1) bnr complex;
1753 * 2) bnr real and no infinite place divide cond_chi (TODO);
1754 * 3) bnr real and one infinite place divide cond_chi;
1755 * 4) bnr real and both infinite places divide cond_chi (TODO) */
1756 static void
QuadGetST(GEN bnr,GEN * pS,GEN * pT,GEN CR,long prec)1757 QuadGetST(GEN bnr, GEN *pS, GEN *pT, GEN CR, long prec)
1758 {
1759 pari_sp av, av1, av2;
1760 long ncond, n, j, k, n0;
1761 GEN vChar = gel(CR,1), dataCR = gel(CR,2), S, T, an, cs, N0, C;
1762 LISTray LIST;
1763
1764 /* initializations */
1765 ST_alloc(pS, pT, lg(dataCR), prec); T = *pT; S = *pS;
1766 av = avma;
1767 ncond = lg(vChar)-1;
1768 C = cgetg(ncond+1, t_VEC);
1769 N0 = cgetg(ncond+1, t_VECSMALL);
1770 cs = cgetg(ncond+1, t_VECSMALL);
1771 n0 = 0;
1772 for (j = 1; j <= ncond; j++)
1773 {
1774 GEN dtcr = gel(dataCR, mael(vChar,j,1)), c = ch_C(dtcr);
1775 long r1, r2;
1776
1777 gel(C,j) = c;
1778 nf_get_sign(bnr_get_nf(ch_bnr(dtcr)), &r1, &r2);
1779 if (r1 == 2) /* real quadratic */
1780 {
1781 cs[j] = 2 + ch_q(dtcr);
1782 if (cs[j] == 2 || cs[j] == 4)
1783 { /* NOT IMPLEMENTED YET */
1784 GetST0(bnr, pS, pT, CR, prec);
1785 return;
1786 }
1787 /* FIXME: is this value of N0 correct for the general case ? */
1788 N0[j] = (long)prec2nbits_mul(prec, 0.35 * gtodouble(c));
1789 }
1790 else /* complex quadratic */
1791 {
1792 cs[j] = 1;
1793 N0[j] = (long)prec2nbits_mul(prec, 0.7 * gtodouble(c));
1794 }
1795 if (n0 < N0[j]) n0 = N0[j];
1796 }
1797 if (DEBUGLEVEL>1) err_printf("N0 = %ld\n", n0);
1798 InitPrimesQuad(bnr, n0, &LIST);
1799
1800 av1 = avma;
1801 /* loop over conductors */
1802 for (j = 1; j <= ncond; j++)
1803 {
1804 GEN c0 = gel(C,j), c1 = divur(1, c0), c2 = divur(2, c0);
1805 GEN ec1 = mpexp(c1), ec2 = mpexp(c2), LChar = gel(vChar,j);
1806 GEN vf0, vf1, cf0, cf1;
1807 const long nChar = lg(LChar)-1, NN = N0[j];
1808
1809 if (DEBUGLEVEL>1)
1810 err_printf("* conductor no %ld/%ld (N = %ld)\n\tInit: ", j,ncond,NN);
1811 if (realprec(ec1) > prec) ec1 = rtor(ec1, prec);
1812 if (realprec(ec2) > prec) ec2 = rtor(ec2, prec);
1813 switch(cs[j])
1814 {
1815 case 1:
1816 cf0 = gen_1;
1817 cf1 = c0;
1818 vf0 = mpveceint1(rtor(c1, prec), ec1, NN);
1819 vf1 = mpvecpowdiv(invr(ec1), NN); break;
1820
1821 case 3:
1822 cf0 = sqrtr(mppi(prec));
1823 cf1 = gmul2n(cf0, 1);
1824 cf0 = gmul(cf0, c0);
1825 vf0 = mpvecpowdiv(invr(ec2), NN);
1826 vf1 = mpveceint1(rtor(c2, prec), ec2, NN); break;
1827
1828 default:
1829 cf0 = cf1 = NULL; /* FIXME: not implemented */
1830 vf0 = vf1 = NULL;
1831 }
1832 for (k = 1; k <= nChar; k++)
1833 {
1834 long u = LChar[k], d, c;
1835 GEN dtcr = gel(dataCR, u), z, s, t;
1836 int **matan;
1837
1838 if (!ch_comp(dtcr)) continue;
1839 if (DEBUGLEVEL>1) err_printf("\tchar no: %ld (%ld/%ld)\n", u,k,nChar);
1840 d = ch_phideg(dtcr);
1841 z = gel(ch_CHI(dtcr), 2); s = t = gen_0; av2 = avma;
1842 matan = computean(gel(dataCR,u), &LIST, NN, d);
1843 for (n = 1, c = 0; n <= NN; n++)
1844 if ((an = EvalCoeff(z, matan[n], d)))
1845 {
1846 s = gadd(s, gmul(an, gel(vf0,n)));
1847 t = gadd(t, gmul(an, gel(vf1,n)));
1848 if (++c == 256) { gerepileall(av2,2, &s,&t); c = 0; }
1849 }
1850 gaffect(gmul(cf0, s), gel(S,u));
1851 gaffect(gmul(cf1, conj_i(t)), gel(T,u));
1852 FreeMat(matan,NN); set_avma(av2);
1853 }
1854 if (DEBUGLEVEL>1) err_printf("\n");
1855 set_avma(av1);
1856 }
1857 set_avma(av);
1858 }
1859
1860 /* s += t*u. All 3 of them t_REAL, except we allow s or u = NULL (for 0) */
1861 static GEN
_addmulrr(GEN s,GEN t,GEN u)1862 _addmulrr(GEN s, GEN t, GEN u)
1863 {
1864 if (u)
1865 {
1866 GEN v = mulrr(t, u);
1867 return s? addrr(s, v): v;
1868 }
1869 return s;
1870 }
1871 /* s += t. Both real, except we allow s or t = NULL (for exact 0) */
1872 static GEN
_addrr(GEN s,GEN t)1873 _addrr(GEN s, GEN t)
1874 { return t? (s? addrr(s, t): t) : s; }
1875
1876 /* S & T for the general case. This is time-critical: optimize */
1877 static void
get_cS_cT(ST_t * T,long n)1878 get_cS_cT(ST_t *T, long n)
1879 {
1880 pari_sp av;
1881 GEN csurn, nsurc, lncsurn, A, B, s, t, Z, aij, bij;
1882 long i, j, r, i0;
1883
1884 if (gel(T->cS,n)) return;
1885
1886 av = avma;
1887 aij = T->aij; i0= T->i0;
1888 bij = T->bij; r = T->r;
1889 Z = cgetg(r+1, t_VEC);
1890 gel(Z,1) = NULL; /* unused */
1891
1892 csurn = divru(T->c1, n);
1893 nsurc = invr(csurn);
1894 lncsurn = logr_abs(csurn);
1895
1896 if (r > 1)
1897 {
1898 gel(Z,2) = lncsurn; /* r >= 2 */
1899 for (i = 3; i <= r; i++)
1900 gel(Z,i) = divru(mulrr(gel(Z,i-1), lncsurn), i-1);
1901 /* Z[i] = ln^(i-1)(c1/n) / (i-1)! */
1902 }
1903
1904 /* i = i0 */
1905 A = gel(aij,i0); t = _addrr(NULL, gel(A,1));
1906 B = gel(bij,i0); s = _addrr(NULL, gel(B,1));
1907 for (j = 2; j <= r; j++)
1908 {
1909 s = _addmulrr(s, gel(Z,j),gel(B,j));
1910 t = _addmulrr(t, gel(Z,j),gel(A,j));
1911 }
1912 for (i = i0 - 1; i > 1; i--)
1913 {
1914 A = gel(aij,i); if (t) t = mulrr(t, nsurc);
1915 B = gel(bij,i); if (s) s = mulrr(s, nsurc);
1916 for (j = odd(i)? T->rc2: T->rc1; j > 1; j--)
1917 {
1918 s = _addmulrr(s, gel(Z,j),gel(B,j));
1919 t = _addmulrr(t, gel(Z,j),gel(A,j));
1920 }
1921 s = _addrr(s, gel(B,1));
1922 t = _addrr(t, gel(A,1));
1923 }
1924 /* i = 1 */
1925 A = gel(aij,1); if (t) t = mulrr(t, nsurc);
1926 B = gel(bij,1); if (s) s = mulrr(s, nsurc);
1927 s = _addrr(s, gel(B,1));
1928 t = _addrr(t, gel(A,1));
1929 for (j = 2; j <= r; j++)
1930 {
1931 s = _addmulrr(s, gel(Z,j),gel(B,j));
1932 t = _addmulrr(t, gel(Z,j),gel(A,j));
1933 }
1934 s = _addrr(s, T->b? mulrr(csurn, gel(T->powracpi,T->b+1)): csurn);
1935 if (!s) s = gen_0;
1936 if (!t) t = gen_0;
1937 gel(T->cS,n) = gclone(s);
1938 gel(T->cT,n) = gclone(t); set_avma(av);
1939 }
1940
1941 static void
clear_cScT(ST_t * T,long N)1942 clear_cScT(ST_t *T, long N)
1943 {
1944 GEN cS = T->cS, cT = T->cT;
1945 long i;
1946 for (i=1; i<=N; i++)
1947 if (cS[i]) {
1948 gunclone(gel(cS,i));
1949 gunclone(gel(cT,i)); gel(cS,i) = gel(cT,i) = NULL;
1950 }
1951 }
1952
1953 static void
init_cScT(ST_t * T,GEN dtcr,long N,long prec)1954 init_cScT(ST_t *T, GEN dtcr, long N, long prec)
1955 {
1956 ch_get3(dtcr, &T->a, &T->b, &T->c);
1957 T->rc1 = T->a + T->c;
1958 T->rc2 = T->b + T->c;
1959 T->r = maxss(T->rc2+1, T->rc1); /* >= 2 */
1960 ppgamma(T, prec);
1961 clear_cScT(T, N);
1962 }
1963
1964 /* return a t_REAL */
1965 static GEN
zeta_get_limx(long r1,long r2,long bit)1966 zeta_get_limx(long r1, long r2, long bit)
1967 {
1968 pari_sp av = avma;
1969 GEN p1, p2, c0, c1, A0;
1970 long r = r1 + r2, N = r + r2;
1971
1972 /* c1 = N 2^(-2r2 / N) */
1973 c1 = mulrs(powrfrac(real2n(1, DEFAULTPREC), -2*r2, N), N);
1974
1975 p1 = powru(Pi2n(1, DEFAULTPREC), r - 1);
1976 p2 = mulir(powuu(N,r), p1); shiftr_inplace(p2, -r2);
1977 c0 = sqrtr( divrr(p2, powru(c1, r+1)) );
1978
1979 A0 = logr_abs( gmul2n(c0, bit) ); p2 = divrr(A0, c1);
1980 p1 = divrr(mulur(N*(r+1), logr_abs(p2)), addsr(2*(r+1), gmul2n(A0,2)));
1981 return gerepileuptoleaf(av, divrr(addrs(p1, 1), powruhalf(p2, N)));
1982 }
1983 /* N_0 = floor( C_K / limx ). Large */
1984 static long
zeta_get_N0(GEN C,GEN limx)1985 zeta_get_N0(GEN C, GEN limx)
1986 {
1987 long e;
1988 pari_sp av = avma;
1989 GEN z = gcvtoi(gdiv(C, limx), &e); /* avoid truncation error */
1990 if (e >= 0 || is_bigint(z))
1991 pari_err_OVERFLOW("zeta_get_N0 [need too many primes]");
1992 return gc_long(av, itos(z));
1993 }
1994
1995 static GEN
eval_i(long r1,long r2,GEN limx,long i)1996 eval_i(long r1, long r2, GEN limx, long i)
1997 {
1998 GEN t = powru(limx, i);
1999 if (!r1) t = mulrr(t, powru(mpfactr(i , DEFAULTPREC), r2));
2000 else if (!r2) t = mulrr(t, powru(mpfactr(i/2, DEFAULTPREC), r1));
2001 else {
2002 GEN u1 = mpfactr(i/2, DEFAULTPREC);
2003 GEN u2 = mpfactr(i, DEFAULTPREC);
2004 if (r1 == r2) t = mulrr(t, powru(mulrr(u1,u2), r1));
2005 else t = mulrr(t, mulrr(powru(u1,r1), powru(u2,r2)));
2006 }
2007 return t;
2008 }
2009
2010 /* "small" even i such that limx^i ( (i\2)! )^r1 ( i! )^r2 > B. */
2011 static long
get_i0(long r1,long r2,GEN B,GEN limx)2012 get_i0(long r1, long r2, GEN B, GEN limx)
2013 {
2014 long imin = 1, imax = 1400;
2015 while (mpcmp(eval_i(r1,r2,limx, imax), B) < 0) { imin = imax; imax *= 2; }
2016 while(imax - imin >= 4)
2017 {
2018 long m = (imax + imin) >> 1;
2019 if (mpcmp(eval_i(r1,r2,limx, m), B) >= 0) imax = m; else imin = m;
2020 }
2021 return imax & ~1; /* make it even */
2022 }
2023 /* limx = zeta_get_limx(r1, r2, bit), a t_REAL */
2024 static long
zeta_get_i0(long r1,long r2,long bit,GEN limx)2025 zeta_get_i0(long r1, long r2, long bit, GEN limx)
2026 {
2027 pari_sp av = avma;
2028 GEN B = gmul(sqrtr( divrr(powrs(mppi(DEFAULTPREC), r2-3), limx) ),
2029 gmul2n(powuu(5, r1), bit + r2));
2030 return gc_long(av, get_i0(r1, r2, B, limx));
2031 }
2032
2033 static void
GetST0(GEN bnr,GEN * pS,GEN * pT,GEN CR,long prec)2034 GetST0(GEN bnr, GEN *pS, GEN *pT, GEN CR, long prec)
2035 {
2036 pari_sp av, av1, av2;
2037 long n, j, k, jc, n0, prec2, i0, r1, r2, ncond;
2038 GEN nf = bnr_get_nf(bnr);
2039 GEN vChar = gel(CR,1), dataCR = gel(CR,2), N0, C, an, limx, S, T;
2040 LISTray LIST;
2041 ST_t cScT;
2042
2043 ST_alloc(pS, pT, lg(dataCR), prec); T = *pT; S = *pS;
2044 av = avma;
2045 nf_get_sign(nf,&r1,&r2);
2046 ncond = lg(vChar)-1;
2047 C = cgetg(ncond+1, t_VEC);
2048 N0 = cgetg(ncond+1, t_VECSMALL);
2049 n0 = 0;
2050 limx = zeta_get_limx(r1, r2, prec2nbits(prec));
2051 for (j = 1; j <= ncond; j++)
2052 {
2053 GEN dtcr = gel(dataCR, mael(vChar,j,1)), c = ch_C(dtcr);
2054 gel(C,j) = c;
2055 N0[j] = zeta_get_N0(c, limx);
2056 if (n0 < N0[j]) n0 = N0[j];
2057 }
2058 cScT.i0 = i0 = zeta_get_i0(r1, r2, prec2nbits(prec), limx);
2059 if (DEBUGLEVEL>1) err_printf("i0 = %ld, N0 = %ld\n",i0, n0);
2060 InitPrimes(bnr, n0, &LIST);
2061 prec2 = precdbl(prec) + EXTRAPREC64;
2062 cScT.powracpi = powersr(sqrtr(mppi(prec2)), r1);
2063 cScT.cS = cgetg(n0+1, t_VEC);
2064 cScT.cT = cgetg(n0+1, t_VEC);
2065 for (j=1; j<=n0; j++) gel(cScT.cS,j) = gel(cScT.cT,j) = NULL;
2066
2067 av1 = avma;
2068 for (jc = 1; jc <= ncond; jc++)
2069 {
2070 GEN LChar = gel(vChar,jc);
2071 long nChar = lg(LChar)-1, N = N0[jc];
2072
2073 /* Can discard completely a conductor if all chars satisfy L(0,chi) = 0
2074 * Not sure whether this is possible. */
2075 if (DEBUGLEVEL>1)
2076 err_printf("* conductor no %ld/%ld (N = %ld)\n\tInit: ", jc,ncond,N);
2077
2078 cScT.c1 = gel(C,jc);
2079 init_cScT(&cScT, gel(dataCR, LChar[1]), N, prec2);
2080 av2 = avma;
2081 for (k = 1; k <= nChar; k++)
2082 {
2083 long d, c, u = LChar[k];
2084 GEN dtcr = gel(dataCR, u), z, s, t;
2085 int **matan;
2086
2087 if (!ch_comp(dtcr)) continue;
2088 if (DEBUGLEVEL>1) err_printf("\tchar no: %ld (%ld/%ld)\n", u,k,nChar);
2089 z = gel(ch_CHI(dtcr), 2);
2090 d = ch_phideg(dtcr); s = t = gen_0;
2091 matan = ComputeCoeff(dtcr, &LIST, N, d);
2092 for (n = 1, c = 0; n <= N; n++)
2093 if ((an = EvalCoeff(z, matan[n], d)))
2094 {
2095 get_cS_cT(&cScT, n);
2096 s = gadd(s, gmul(an, gel(cScT.cS,n)));
2097 t = gadd(t, gmul(an, gel(cScT.cT,n)));
2098 if (++c == 256) { gerepileall(av2,2, &s,&t); c = 0; }
2099 }
2100 gaffect(s, gel(S,u));
2101 gaffect(conj_i(t), gel(T,u));
2102 FreeMat(matan, N); set_avma(av2);
2103 }
2104 if (DEBUGLEVEL>1) err_printf("\n");
2105 set_avma(av1);
2106 }
2107 clear_cScT(&cScT, n0);
2108 set_avma(av);
2109 }
2110
2111 static void
GetST(GEN bnr,GEN * pS,GEN * pT,GEN CR,long prec)2112 GetST(GEN bnr, GEN *pS, GEN *pT, GEN CR, long prec)
2113 {
2114 if (nf_get_degree(bnr_get_nf(bnr)) == 2)
2115 QuadGetST(bnr, pS, pT, CR, prec);
2116 else
2117 GetST0(bnr, pS, pT, CR, prec);
2118 }
2119
2120 /*******************************************************************/
2121 /* */
2122 /* Class fields of real quadratic fields using Stark units */
2123 /* */
2124 /*******************************************************************/
2125 /* compute the Hilbert class field using genus class field theory when
2126 the exponent of the class group is exactly 2 (trivial group not covered) */
2127 /* Cf Herz, Construction of class fields, LNM 21, Theorem 1 (VII-6) */
2128 static GEN
GenusFieldQuadReal(GEN disc)2129 GenusFieldQuadReal(GEN disc)
2130 {
2131 long i, i0 = 0, l;
2132 pari_sp av = avma;
2133 GEN T = NULL, p0 = NULL, P;
2134
2135 P = gel(Z_factor(disc), 1);
2136 l = lg(P);
2137 for (i = 1; i < l; i++)
2138 {
2139 GEN p = gel(P,i);
2140 if (mod4(p) == 3) { p0 = p; i0 = i; break; }
2141 }
2142 l--; /* remove last prime */
2143 if (i0 == l) l--; /* ... remove p0 and last prime */
2144 for (i = 1; i < l; i++)
2145 {
2146 GEN p = gel(P,i), d, t;
2147 if (i == i0) continue;
2148 if (absequaliu(p, 2))
2149 switch (mod32(disc))
2150 {
2151 case 8: d = gen_2; break;
2152 case 24: d = shifti(p0, 1); break;
2153 default: d = p0; break;
2154 }
2155 else
2156 d = (mod4(p) == 1)? p: mulii(p0, p);
2157 t = mkpoln(3, gen_1, gen_0, negi(d)); /* x^2 - d */
2158 T = T? ZX_compositum_disjoint(T, t): t;
2159 }
2160 return gerepileupto(av, polredbest(T, 0));
2161 }
2162 static GEN
GenusFieldQuadImag(GEN disc)2163 GenusFieldQuadImag(GEN disc)
2164 {
2165 long i, l;
2166 pari_sp av = avma;
2167 GEN T = NULL, P;
2168
2169 P = gel(absZ_factor(disc), 1);
2170 l = lg(P);
2171 l--; /* remove last prime */
2172 for (i = 1; i < l; i++)
2173 {
2174 GEN p = gel(P,i), d, t;
2175 if (absequaliu(p, 2))
2176 switch (mod32(disc))
2177 {
2178 case 24: d = gen_2; break; /* disc = 8 mod 32 */
2179 case 8: d = gen_m2; break; /* disc =-8 mod 32 */
2180 default: d = gen_m1; break;
2181 }
2182 else
2183 d = (mod4(p) == 1)? p: negi(p);
2184 t = mkpoln(3, gen_1, gen_0, negi(d)); /* x^2 - d */
2185 T = T? ZX_compositum_disjoint(T, t): t;
2186 }
2187 return gerepileupto(av, polredbest(T, 0));
2188 }
2189
2190 /* if flag != 0, computes a fast and crude approximation of the result */
2191 static GEN
AllStark(GEN data,long flag,long newprec)2192 AllStark(GEN data, long flag, long newprec)
2193 {
2194 const long BND = 300;
2195 long cl, i, j, cpt = 0, N, h, v, n, r1, r2, den;
2196 pari_sp av, av2;
2197 int **matan;
2198 GEN bnr = gel(data,1), nf = bnr_get_nf(bnr), p1, p2, S, T;
2199 GEN CR = gel(data,4), dataCR = gel(CR,2);
2200 GEN polrelnum, polrel, Lp, W, vzeta, C, cond1, L1, an;
2201 LISTray LIST;
2202 pari_timer ti;
2203
2204 nf_get_sign(nf, &r1,&r2);
2205 N = nf_get_degree(nf);
2206 cond1 = gel(bnr_get_mod(bnr), 2);
2207
2208 v = 1;
2209 while (gequal1(gel(cond1,v))) v++;
2210 cl = lg(dataCR)-1;
2211 h = itos(ZM_det_triangular(gel(data,2))) >> 1;
2212
2213 LABDOUB:
2214 if (DEBUGLEVEL) timer_start(&ti);
2215 av = avma;
2216 W = AllArtinNumbers(CR, newprec);
2217 if (DEBUGLEVEL) timer_printf(&ti,"Compute W");
2218 Lp = cgetg(cl + 1, t_VEC);
2219 if (!flag)
2220 {
2221 GetST(bnr, &S, &T, CR, newprec);
2222 if (DEBUGLEVEL) timer_printf(&ti, "S&T");
2223 for (i = 1; i <= cl; i++)
2224 {
2225 GEN chi = gel(dataCR, i), v = gen_0;
2226 if (ch_comp(chi))
2227 v = gel(GetValue(chi, gel(W,i), gel(S,i), gel(T,i), 2, newprec), 2);
2228 gel(Lp, i) = v;
2229 }
2230 }
2231 else
2232 { /* compute a crude approximation of the result */
2233 C = cgetg(cl + 1, t_VEC);
2234 for (i = 1; i <= cl; i++) gel(C,i) = ch_C(gel(dataCR, i));
2235 n = zeta_get_N0(vecmax(C), zeta_get_limx(r1, r2, prec2nbits(newprec)));
2236 if (n > BND) n = BND;
2237 if (DEBUGLEVEL) err_printf("N0 in QuickPol: %ld \n", n);
2238 InitPrimes(bnr, n, &LIST);
2239
2240 L1 = cgetg(cl+1, t_VEC);
2241 /* use L(1) = sum (an / n) */
2242 for (i = 1; i <= cl; i++)
2243 {
2244 GEN dtcr = gel(dataCR,i);
2245 long d = ch_phideg(dtcr);
2246 matan = ComputeCoeff(dtcr, &LIST, n, d);
2247 av2 = avma;
2248 p1 = real_0(newprec); p2 = gel(ch_CHI(dtcr), 2);
2249 for (j = 1; j <= n; j++)
2250 if ( (an = EvalCoeff(p2, matan[j], d)) )
2251 p1 = gadd(p1, gdivgs(an, j));
2252 gel(L1,i) = gerepileupto(av2, p1);
2253 FreeMat(matan, n);
2254 }
2255 p1 = gmul2n(powruhalf(mppi(newprec), N-2), 1);
2256
2257 for (i = 1; i <= cl; i++)
2258 {
2259 long r;
2260 GEN WW, A = AChi(gel(dataCR,i), &r, 0, newprec);
2261 WW = gmul(gel(C,i), gmul(A, gel(W,i)));
2262 gel(Lp,i) = gdiv(gmul(WW, conj_i(gel(L1,i))), p1);
2263 }
2264 }
2265
2266 p1 = gel(data,3);
2267 den = flag ? h: 2*h;
2268 vzeta = cgetg(h + 1, t_VEC);
2269 for (i = 1; i <= h; i++)
2270 {
2271 GEN z = gen_0, sig = gel(p1,i);
2272 for (j = 1; j <= cl; j++)
2273 {
2274 GEN dtcr = gel(dataCR,j), CHI = ch_CHI(dtcr);
2275 GEN t = mulreal(gel(Lp,j), CharEval(CHI, sig));
2276 if (chi_get_deg(CHI) != 2) t = gmul2n(t, 1); /* character not real */
2277 z = gadd(z, t);
2278 }
2279 gel(vzeta,i) = gmul2n(gcosh(gdivgs(z,den), newprec), 1);
2280 }
2281 polrelnum = roots_to_pol(vzeta, 0);
2282 if (DEBUGLEVEL)
2283 {
2284 if (DEBUGLEVEL>1) {
2285 err_printf("polrelnum = %Ps\n", polrelnum);
2286 err_printf("zetavalues = %Ps\n", vzeta);
2287 if (!flag)
2288 err_printf("Checking the square-root of the Stark unit...\n");
2289 }
2290 timer_printf(&ti, "Compute %s", flag? "quickpol": "polrelnum");
2291 }
2292 if (flag) return gerepilecopy(av, polrelnum);
2293
2294 /* try to recognize this polynomial */
2295 polrel = RecCoeff(nf, polrelnum, v, newprec);
2296 if (!polrel)
2297 {
2298 for (j = 1; j <= h; j++)
2299 gel(vzeta,j) = gsubgs(gsqr(gel(vzeta,j)), 2);
2300 polrelnum = roots_to_pol(vzeta, 0);
2301 if (DEBUGLEVEL)
2302 {
2303 if (DEBUGLEVEL>1) {
2304 err_printf("It's not a square...\n");
2305 err_printf("polrelnum = %Ps\n", polrelnum);
2306 }
2307 timer_printf(&ti, "Compute polrelnum");
2308 }
2309 polrel = RecCoeff(nf, polrelnum, v, newprec);
2310 }
2311 if (!polrel) /* FAILED */
2312 {
2313 const long EXTRA_BITS = 64;
2314 long incr_pr;
2315 if (++cpt >= 3) pari_err_PREC( "stark (computation impossible)");
2316 /* estimate needed precision */
2317 incr_pr = prec2nbits(gprecision(polrelnum))- gexpo(polrelnum);
2318 if (incr_pr < 0) incr_pr = -incr_pr + EXTRA_BITS;
2319 newprec += nbits2extraprec(maxss(3*EXTRA_BITS, cpt*incr_pr));
2320 if (DEBUGLEVEL) pari_warn(warnprec, "AllStark", newprec);
2321 CharNewPrec(data, newprec);
2322 nf = bnr_get_nf(ch_bnr(gel(dataCR,1)));
2323 goto LABDOUB;
2324 }
2325
2326 if (DEBUGLEVEL) {
2327 if (DEBUGLEVEL>1) err_printf("polrel = %Ps\n", polrel);
2328 timer_printf(&ti, "Recpolnum");
2329 }
2330 return gerepilecopy(av, polrel);
2331 }
2332
2333 /********************************************************************/
2334 /* Main functions */
2335 /********************************************************************/
2336 static GEN
bnrstark_cyclic(GEN bnr,GEN dtQ,long prec)2337 bnrstark_cyclic(GEN bnr, GEN dtQ, long prec)
2338 {
2339 GEN v, vH, cyc = gel(dtQ,2), U = gel(dtQ,3), M = ZM_inv(U, NULL);
2340 long i, j, l = lg(M);
2341
2342 /* M = indep. generators of Cl_f/subgp, restrict to cyclic components */
2343 vH = cgetg(l, t_VEC);
2344 for (i = j = 1; i < l; i++)
2345 {
2346 if (is_pm1(gel(cyc,i))) break;
2347 gel(vH, j++) = ZM_hnfmodid(vecsplice(M,i), cyc);
2348 }
2349 setlg(vH, j); v = cgetg(l, t_VEC);
2350 for (i = 1; i < j; i++) gel(v,i) = bnrstark(bnr, gel(vH,i), prec);
2351 return v;
2352 }
2353 GEN
bnrstark(GEN bnr,GEN subgrp,long prec)2354 bnrstark(GEN bnr, GEN subgrp, long prec)
2355 {
2356 long newprec;
2357 pari_sp av = avma;
2358 GEN nf, data, dtQ;
2359
2360 /* check the bnr */
2361 checkbnr(bnr); nf = bnr_get_nf(bnr);
2362 if (nf_get_degree(nf) == 1) return galoissubcyclo(bnr, subgrp, 0, 0);
2363 if (!nf_get_varn(nf))
2364 pari_err_PRIORITY("bnrstark", nf_get_pol(nf), "=", 0);
2365 if (nf_get_r2(nf)) pari_err_DOMAIN("bnrstark", "r2", "!=", gen_0, nf);
2366
2367 /* compute bnr(conductor) */
2368 bnr_subgroup_sanitize(&bnr, &subgrp);
2369 if (gequal1(ZM_det_triangular(subgrp))) { set_avma(av); return pol_x(0); }
2370
2371 /* check the class field */
2372 if (!gequal0(gel(bnr_get_mod(bnr), 2)))
2373 pari_err_DOMAIN("bnrstark", "r2(class field)", "!=", gen_0, bnr);
2374
2375 /* find a suitable extension N */
2376 dtQ = InitQuotient(subgrp);
2377 data = FindModulus(bnr, dtQ, &newprec);
2378 if (!data) return gerepileupto(av, bnrstark_cyclic(bnr, dtQ, prec));
2379 if (DEBUGLEVEL>1 && newprec > prec)
2380 err_printf("new precision: %ld\n", newprec);
2381 return gerepileupto(av, AllStark(data, 0, newprec));
2382 }
2383
2384 /* For each character of Cl(bnr)/subgp, compute L(1, chi) (or equivalently
2385 * the first nonzero term c(chi) of the expansion at s = 0).
2386 * If flag & 1: compute the value at s = 1 (for nontrivial characters),
2387 * else compute the term c(chi) and return [r(chi), c(chi)] where r(chi) is
2388 * the order of L(s, chi) at s = 0.
2389 * If flag & 2: compute the value of the L-function L_S(s, chi) where S is the
2390 * set of places dividing the modulus of bnr (and the infinite places),
2391 * else
2392 * compute the value of the primitive L-function attached to chi,
2393 * If flag & 4: return also the character */
2394 GEN
bnrL1(GEN bnr,GEN subgp,long flag,long prec)2395 bnrL1(GEN bnr, GEN subgp, long flag, long prec)
2396 {
2397 GEN L1, ch, Qt, z;
2398 long l, h;
2399 pari_sp av = avma;
2400
2401 checkbnr(bnr);
2402 if (flag < 0 || flag > 8) pari_err_FLAG("bnrL1");
2403
2404 subgp = bnr_subgroup_check(bnr, subgp, NULL);
2405 if (!subgp) subgp = diagonal_shallow(bnr_get_cyc(bnr));
2406
2407 Qt = InitQuotient(subgp);
2408 ch = AllChars(bnr, Qt, 0); l = lg(ch);
2409 h = itou(gel(Qt,1));
2410 L1 = cgetg((flag&1)? h: h+1, t_VEC);
2411 if (l > 1)
2412 {
2413 GEN W, S, T, CR = InitChar(bnr, ch, 1, prec), dataCR = gel(CR,2);
2414 long i, j;
2415
2416 GetST(bnr, &S, &T, CR, prec);
2417 W = AllArtinNumbers(CR, prec);
2418 for (i = j = 1; i < l; i++)
2419 {
2420 GEN chi = gel(ch,i);
2421 z = GetValue(gel(dataCR,i), gel(W,i), gel(S,i), gel(T,i), flag, prec);
2422 gel(L1,j++) = (flag & 4)? mkvec2(gel(chi,1), z): z;
2423 if (lg(chi) == 4)
2424 { /* nonreal */
2425 z = conj_i(z);
2426 gel(L1, j++) = (flag & 4)? mkvec2(gel(chi,3), z): z;
2427 }
2428 }
2429 }
2430 if (!(flag & 1))
2431 { /* trivial character */
2432 z = GetValue1(bnr, flag & 2, prec);
2433 if (flag & 4)
2434 {
2435 GEN chi = zerovec(lg(bnr_get_cyc(bnr))-1);
2436 z = mkvec2(chi, z);
2437 }
2438 gel(L1,h) = z;
2439 }
2440 return gerepilecopy(av, L1);
2441 }
2442
2443 /*******************************************************************/
2444 /* */
2445 /* Hilbert and Ray Class field using Stark */
2446 /* */
2447 /*******************************************************************/
2448 /* P in A[x,y], deg_y P < 2, return P0 and P1 in A[x] such that P = P0 + P1 y */
2449 static void
split_pol_quad(GEN P,GEN * gP0,GEN * gP1)2450 split_pol_quad(GEN P, GEN *gP0, GEN *gP1)
2451 {
2452 long i, l = lg(P);
2453 GEN P0 = cgetg(l, t_POL), P1 = cgetg(l, t_POL);
2454 P0[1] = P1[1] = P[1];
2455 for (i = 2; i < l; i++)
2456 {
2457 GEN c = gel(P,i), c0 = c, c1 = gen_0;
2458 if (typ(c) == t_POL) /* write c = c1 y + c0 */
2459 switch(degpol(c))
2460 {
2461 case -1: c0 = gen_0; break;
2462 default: c1 = gel(c,3); /* fall through */
2463 case 0: c0 = gel(c,2); break;
2464 }
2465 gel(P0,i) = c0; gel(P1,i) = c1;
2466 }
2467 *gP0 = normalizepol_lg(P0, l);
2468 *gP1 = normalizepol_lg(P1, l);
2469 }
2470
2471 /* k = nf quadratic field, P relative equation of H_k (Hilbert class field)
2472 * return T in Z[X], such that H_k / Q is the compositum of Q[X]/(T) and k */
2473 static GEN
makescind(GEN nf,GEN P)2474 makescind(GEN nf, GEN P)
2475 {
2476 GEN Pp, p, pol, G, L, a, roo, P0,P1, Ny,Try, nfpol = nf_get_pol(nf);
2477 long i, is_P;
2478
2479 P = lift_shallow(P);
2480 split_pol_quad(P, &P0, &P1);
2481 /* P = P0 + y P1, Norm_{k/Q}(P) = P0^2 + Tr y P0P1 + Ny P1^2, irreducible/Q */
2482 Ny = gel(nfpol, 2);
2483 Try = negi(gel(nfpol, 3));
2484 pol = RgX_add(RgX_sqr(P0), RgX_Rg_mul(RgX_sqr(P1), Ny));
2485 if (signe(Try)) pol = RgX_add(pol, RgX_Rg_mul(RgX_mul(P0,P1), Try));
2486 /* pol = rnfequation(nf, P); */
2487 G = galoisinit(pol, NULL);
2488 L = gal_get_group(G);
2489 p = gal_get_p(G);
2490 a = FpX_oneroot(nfpol, p);
2491 /* P mod a prime \wp above p (which splits) */
2492 Pp = FpXY_evalx(P, a, p);
2493 roo = gal_get_roots(G);
2494 is_P = gequal0( FpX_eval(Pp, remii(gel(roo,1),p), p) );
2495 /* each roo[i] mod p is a root of P or (exclusive) tau(P) mod \wp */
2496 /* record whether roo[1] is a root of P or tau(P) */
2497
2498 for (i = 1; i < lg(L); i++)
2499 {
2500 GEN perm = gel(L,i);
2501 long k = perm[1]; if (k == 1) continue;
2502 k = gequal0( FpX_eval(Pp, remii(gel(roo,k),p), p) );
2503 /* roo[k] is a root of the other polynomial */
2504 if (k != is_P)
2505 {
2506 ulong o = perm_orderu(perm);
2507 if (o != 2) perm = perm_powu(perm, o >> 1);
2508 /* perm has order two and doesn't belong to Gal(H_k/k) */
2509 return galoisfixedfield(G, perm, 1, varn(P));
2510 }
2511 }
2512 pari_err_BUG("makescind");
2513 return NULL; /*LCOV_EXCL_LINE*/
2514 }
2515
2516 /* pbnf = NULL if no bnf is needed, f = NULL may be passed for a trivial
2517 * conductor */
2518 static void
quadray_init(GEN * pD,GEN f,GEN * pbnf,long prec)2519 quadray_init(GEN *pD, GEN f, GEN *pbnf, long prec)
2520 {
2521 GEN D = *pD, nf, bnf = NULL;
2522 if (typ(D) == t_INT)
2523 {
2524 int isfund;
2525 if (pbnf) {
2526 long v = f? gvar(f): NO_VARIABLE;
2527 if (v == NO_VARIABLE) v = 1;
2528 bnf = Buchall(quadpoly0(D, v), nf_FORCE, prec);
2529 nf = bnf_get_nf(bnf);
2530 isfund = equalii(D, nf_get_disc(nf));
2531 }
2532 else
2533 isfund = Z_isfundamental(D);
2534 if (!isfund) pari_err_DOMAIN("quadray", "isfundamental(D)", "=",gen_0, D);
2535 }
2536 else
2537 {
2538 bnf = checkbnf(D);
2539 nf = bnf_get_nf(bnf);
2540 if (nf_get_degree(nf) != 2)
2541 pari_err_DOMAIN("quadray", "degree", "!=", gen_2, nf_get_pol(nf));
2542 D = nf_get_disc(nf);
2543 }
2544 if (pbnf) *pbnf = bnf;
2545 *pD = D;
2546 }
2547
2548 /* compute the polynomial over Q of the Hilbert class field of
2549 Q(sqrt(D)) where D is a positive fundamental discriminant */
2550 static GEN
quadhilbertreal(GEN D,long prec)2551 quadhilbertreal(GEN D, long prec)
2552 {
2553 pari_sp av = avma;
2554 GEN bnf, pol, bnr, dtQ, data, M;
2555 long newprec;
2556 pari_timer T;
2557
2558 quadray_init(&D, NULL, &bnf, prec);
2559 switch(itou_or_0(cyc_get_expo(bnf_get_cyc(bnf))))
2560 {
2561 case 1: set_avma(av); return pol_x(0);
2562 case 2: return gerepileupto(av, GenusFieldQuadReal(D));
2563 }
2564 bnr = Buchray(bnf, gen_1, nf_INIT);
2565 M = diagonal_shallow(bnr_get_cyc(bnr));
2566 dtQ = InitQuotient(M);
2567
2568 if (DEBUGLEVEL) timer_start(&T);
2569 data = FindModulus(bnr, dtQ, &newprec);
2570 if (DEBUGLEVEL) timer_printf(&T,"FindModulus");
2571 if (!data) return gerepileupto(av, bnrstark_cyclic(bnr, dtQ, prec));
2572 pol = AllStark(data, 0, newprec);
2573 pol = makescind(bnf_get_nf(bnf), pol);
2574 return gerepileupto(av, polredbest(pol, 0));
2575 }
2576
2577 /*******************************************************************/
2578 /* */
2579 /* Hilbert and Ray Class field using CM (Schertz) */
2580 /* */
2581 /*******************************************************************/
2582 /* form^2 = 1 ? */
2583 static int
hasexp2(GEN form)2584 hasexp2(GEN form)
2585 {
2586 GEN a = gel(form,1), b = gel(form,2), c = gel(form,3);
2587 return !signe(b) || absequalii(a,b) || equalii(a,c);
2588 }
2589 static int
uhasexp2(GEN form)2590 uhasexp2(GEN form)
2591 {
2592 long a = form[1], b = form[2], c = form[3];
2593 return !b || a == labs(b) || a == c;
2594 }
2595
2596 GEN
qfbforms(GEN D)2597 qfbforms(GEN D)
2598 {
2599 ulong d = itou(D), dover3 = d/3, t, b2, a, b, c, h;
2600 GEN L = cgetg((long)(sqrt((double)d) * log2(d)), t_VEC);
2601 b2 = b = (d&1); h = 0;
2602 if (!b) /* b = 0 treated separately to avoid special cases */
2603 {
2604 t = d >> 2; /* (b^2 - D) / 4*/
2605 for (a=1; a*a<=t; a++)
2606 if (c = t/a, t == c*a) gel(L,++h) = mkvecsmall3(a,0,c);
2607 b = 2; b2 = 4;
2608 }
2609 /* now b > 0, b = D (mod 2) */
2610 for ( ; b2 <= dover3; b += 2, b2 = b*b)
2611 {
2612 t = (b2 + d) >> 2; /* (b^2 - D) / 4*/
2613 /* b = a */
2614 if (c = t/b, t == c*b) gel(L,++h) = mkvecsmall3(b,b,c);
2615 /* b < a < c */
2616 for (a = b+1; a*a < t; a++)
2617 if (c = t/a, t == c*a)
2618 {
2619 gel(L,++h) = mkvecsmall3(a, b,c);
2620 gel(L,++h) = mkvecsmall3(a,-b,c);
2621 }
2622 /* a = c */
2623 if (a * a == t) gel(L,++h) = mkvecsmall3(a,b,a);
2624 }
2625 setlg(L,h+1); return L;
2626 }
2627
2628 /* gcd(n, 24) */
2629 static long
GCD24(long n)2630 GCD24(long n)
2631 {
2632 switch(n % 24)
2633 {
2634 case 0: return 24;
2635 case 1: return 1;
2636 case 2: return 2;
2637 case 3: return 3;
2638 case 4: return 4;
2639 case 5: return 1;
2640 case 6: return 6;
2641 case 7: return 1;
2642 case 8: return 8;
2643 case 9: return 3;
2644 case 10: return 2;
2645 case 11: return 1;
2646 case 12: return 12;
2647 case 13: return 1;
2648 case 14: return 2;
2649 case 15: return 3;
2650 case 16: return 8;
2651 case 17: return 1;
2652 case 18: return 6;
2653 case 19: return 1;
2654 case 20: return 4;
2655 case 21: return 3;
2656 case 22: return 2;
2657 case 23: return 1;
2658 default: return 0;
2659 }
2660 }
2661
2662 struct gpq_data {
2663 long p, q;
2664 GEN sqd; /* sqrt(D), t_REAL */
2665 GEN u, D;
2666 GEN pq, pq2; /* p*q, 2*p*q */
2667 GEN qfpq ; /* class of \P * \Q */
2668 };
2669
2670 /* find P and Q two non principal prime ideals (above p <= q) such that
2671 * cl(P) = cl(Q) if P,Q have order 2 in Cl(K).
2672 * Ensure that e = 24 / gcd(24, (p-1)(q-1)) = 1 */
2673 /* D t_INT, discriminant */
2674 static void
init_pq(GEN D,struct gpq_data * T)2675 init_pq(GEN D, struct gpq_data *T)
2676 {
2677 const long Np = 6547; /* N.B. primepi(50000) = 5133 */
2678 const ulong maxq = 50000;
2679 GEN listp = cgetg(Np + 1, t_VECSMALL); /* primes p */
2680 GEN listP = cgetg(Np + 1, t_VEC); /* primeform(p) if of order 2, else NULL */
2681 GEN gcd24 = cgetg(Np + 1, t_VECSMALL); /* gcd(p-1, 24) */
2682 forprime_t S;
2683 long l = 1;
2684 double best = 0.;
2685 ulong q;
2686
2687 u_forprime_init(&S, 2, ULONG_MAX);
2688 T->D = D;
2689 T->p = T->q = 0;
2690 for(;;)
2691 {
2692 GEN Q;
2693 long i, gcdq, mod;
2694 int order2, store;
2695 double t;
2696
2697 q = u_forprime_next(&S);
2698 if (best > 0 && q >= maxq)
2699 {
2700 if (DEBUGLEVEL)
2701 pari_warn(warner,"possibly suboptimal (p,q) for D = %Ps", D);
2702 break;
2703 }
2704 if (kroiu(D, q) < 0) continue; /* inert */
2705 Q = redimag(primeform_u(D, q));
2706 if (is_pm1(gel(Q,1))) continue; /* Q | q is principal */
2707
2708 store = 1;
2709 order2 = hasexp2(Q);
2710 gcd24[l] = gcdq = GCD24(q-1);
2711 mod = 24 / gcdq; /* mod must divide p-1 otherwise e > 1 */
2712 listp[l] = q;
2713 gel(listP,l) = order2 ? Q : NULL;
2714 t = (q+1)/(double)(q-1);
2715 for (i = 1; i < l; i++) /* try all (p, q), p < q in listp */
2716 {
2717 long p = listp[i], gcdp = gcd24[i];
2718 double b;
2719 /* P,Q order 2 => cl(Q) = cl(P) */
2720 if (order2 && gel(listP,i) && !gequal(gel(listP,i), Q)) continue;
2721 if (gcdp % gcdq == 0) store = 0; /* already a better one in the list */
2722 if ((p-1) % mod) continue;
2723
2724 b = (t*(p+1)) / (p-1); /* (p+1)(q+1) / (p-1)(q-1) */
2725 if (b > best) {
2726 store = 0; /* (p,q) always better than (q,r) for r >= q */
2727 best = b; T->q = q; T->p = p;
2728 if (DEBUGLEVEL>2) err_printf("p,q = %ld,%ld\n", p, q);
2729 }
2730 /* won't improve with this q as largest member */
2731 if (best > 0) break;
2732 }
2733 /* if !store or (q,r) won't improve on current best pair, forget that q */
2734 if (store && t*t > best)
2735 if (++l >= Np) pari_err_BUG("quadhilbert (not enough primes)");
2736 if (!best) /* (p,q) with p < q always better than (q,q) */
2737 { /* try (q,q) */
2738 if (gcdq >= 12 && umodiu(D, q)) /* e = 1 and unramified */
2739 {
2740 double b = (t*q) / (q-1); /* q(q+1) / (q-1)^2 */
2741 if (b > best) {
2742 best = b; T->q = T->p = q;
2743 if (DEBUGLEVEL>2) err_printf("p,q = %ld,%ld\n", q, q);
2744 }
2745 }
2746 }
2747 /* If (p1+1)(q+1) / (p1-1)(q-1) <= best, we can no longer improve
2748 * even with best p : stop */
2749 if ((listp[1]+1)*t <= (listp[1]-1)*best) break;
2750 }
2751 if (DEBUGLEVEL>1)
2752 err_printf("(p, q) = %ld, %ld; gain = %f\n", T->p, T->q, 12*best);
2753 }
2754
2755 static GEN
gpq(GEN form,struct gpq_data * T)2756 gpq(GEN form, struct gpq_data *T)
2757 {
2758 pari_sp av = avma;
2759 long a = form[1], b = form[2], c = form[3];
2760 long p = T->p, q = T->q;
2761 GEN form2, w, z;
2762 int fl, real = 0;
2763
2764 form2 = qficomp(T->qfpq, mkvec3s(a, -b, c));
2765 /* form2 and form yield complex conjugate roots : only compute for the
2766 * lexicographically smallest of the 2 */
2767 fl = cmpis(gel(form2,1), a);
2768 if (fl <= 0)
2769 {
2770 if (fl < 0) return NULL;
2771 fl = cmpis(gel(form2,2), b);
2772 if (fl <= 0)
2773 {
2774 if (fl < 0) return NULL;
2775 /* form == form2 : real root */
2776 real = 1;
2777 }
2778 }
2779
2780 if (p == 2) { /* (a,b,c) = (1,1,0) mod 2 ? */
2781 if (a % q == 0 && (a & b & 1) && !(c & 1))
2782 { /* apply S : make sure that (a,b,c) represents odd values */
2783 lswap(a,c); b = -b;
2784 }
2785 }
2786 if (a % p == 0 || a % q == 0)
2787 { /* apply T^k, look for c' = a k^2 + b k + c coprime to N */
2788 while (c % p == 0 || c % q == 0)
2789 {
2790 c += a + b;
2791 b += a << 1;
2792 }
2793 lswap(a, c); b = -b; /* apply S */
2794 }
2795 /* now (a,b,c) ~ form and (a,pq) = 1 */
2796
2797 /* gcd(2a, u) = 2, w = u mod 2pq, -b mod 2a */
2798 w = Z_chinese(T->u, stoi(-b), T->pq2, utoipos(a << 1));
2799 z = double_eta_quotient(utoipos(a), w, T->D, T->p, T->q, T->pq, T->sqd);
2800 if (real && typ(z) == t_COMPLEX) z = gcopy(gel(z, 1));
2801 return gerepileupto(av, z);
2802 }
2803
2804 /* returns an equation for the Hilbert class field of Q(sqrt(D)), D < 0
2805 * fundamental discriminant */
2806 static GEN
quadhilbertimag(GEN D)2807 quadhilbertimag(GEN D)
2808 {
2809 GEN L, P, Pi, Pr, qfp, u;
2810 pari_sp av = avma;
2811 long h, i, prec;
2812 struct gpq_data T;
2813 pari_timer ti;
2814
2815 if (DEBUGLEVEL>1) timer_start(&ti);
2816 if (lgefint(D) == 3)
2817 switch (D[2]) { /* = |D|; special cases where e > 1 */
2818 case 3:
2819 case 4:
2820 case 7:
2821 case 8:
2822 case 11:
2823 case 19:
2824 case 43:
2825 case 67:
2826 case 163: return pol_x(0);
2827 }
2828 L = qfbforms(D);
2829 h = lg(L)-1;
2830 if ((1L << vals(h)) == h) /* power of 2 */
2831 { /* check whether > |Cl|/2 elements have order <= 2 ==> 2-elementary */
2832 long lim = (h>>1) + 1;
2833 for (i=1; i <= lim; i++)
2834 if (!uhasexp2(gel(L,i))) break;
2835 if (i > lim) return GenusFieldQuadImag(D);
2836 }
2837 if (DEBUGLEVEL>1) timer_printf(&ti,"class number = %ld",h);
2838 init_pq(D, &T);
2839 qfp = primeform_u(D, T.p);
2840 T.pq = muluu(T.p, T.q);
2841 T.pq2 = shifti(T.pq,1);
2842 if (T.p == T.q)
2843 {
2844 GEN qfbp2 = qficompraw(qfp, qfp);
2845 u = gel(qfbp2,2);
2846 T.u = modii(u, T.pq2);
2847 T.qfpq = redimag(qfbp2);
2848 }
2849 else
2850 {
2851 GEN qfq = primeform_u(D, T.q), bp = gel(qfp,2), bq = gel(qfq,2);
2852 T.u = Z_chinese(bp, bq, utoipos(T.p << 1), utoipos(T.q << 1));
2853 /* T.u = bp (mod 2p), T.u = bq (mod 2q) */
2854 T.qfpq = qficomp(qfp, qfq);
2855 }
2856 /* u modulo 2pq */
2857 prec = LOWDEFAULTPREC;
2858 Pr = cgetg(h+1,t_VEC);
2859 Pi = cgetg(h+1,t_VEC);
2860 for(;;)
2861 {
2862 long ex, exmax = 0, r1 = 0, r2 = 0;
2863 pari_sp av0 = avma;
2864 T.sqd = sqrtr_abs(itor(D, prec));
2865 for (i=1; i<=h; i++)
2866 {
2867 GEN s = gpq(gel(L,i), &T);
2868 if (DEBUGLEVEL>3) err_printf("%ld ", i);
2869 if (!s) continue;
2870 if (typ(s) != t_COMPLEX) gel(Pr, ++r1) = s; /* real root */
2871 else gel(Pi, ++r2) = s;
2872 ex = gexpo(s); if (ex > 0) exmax += ex;
2873 }
2874 if (DEBUGLEVEL>1) timer_printf(&ti,"roots");
2875 setlg(Pr, r1+1);
2876 setlg(Pi, r2+1);
2877 P = roots_to_pol_r1(shallowconcat(Pr,Pi), 0, r1);
2878 P = grndtoi(P,&exmax);
2879 if (DEBUGLEVEL>1) timer_printf(&ti,"product, error bits = %ld",exmax);
2880 if (exmax <= -10) break;
2881 set_avma(av0); prec += nbits2extraprec(prec2nbits(DEFAULTPREC)+exmax);
2882 if (DEBUGLEVEL) pari_warn(warnprec,"quadhilbertimag",prec);
2883 }
2884 return gerepileupto(av,P);
2885 }
2886
2887 GEN
quadhilbert(GEN D,long prec)2888 quadhilbert(GEN D, long prec)
2889 {
2890 GEN d = D;
2891 quadray_init(&d, NULL, NULL, 0);
2892 return (signe(d)>0)? quadhilbertreal(D,prec)
2893 : quadhilbertimag(d);
2894 }
2895
2896 /* return a vector of all roots of 1 in bnf [not necessarily quadratic] */
2897 static GEN
getallrootsof1(GEN bnf)2898 getallrootsof1(GEN bnf)
2899 {
2900 GEN T, u, nf = bnf_get_nf(bnf), tu;
2901 long i, n = bnf_get_tuN(bnf);
2902
2903 if (n == 2) {
2904 long N = nf_get_degree(nf);
2905 return mkvec2(scalarcol_shallow(gen_m1, N),
2906 scalarcol_shallow(gen_1, N));
2907 }
2908 tu = poltobasis(nf, bnf_get_tuU(bnf));
2909 T = zk_multable(nf, tu);
2910 u = cgetg(n+1, t_VEC); gel(u,1) = tu;
2911 for (i=2; i <= n; i++) gel(u,i) = ZM_ZC_mul(T, gel(u,i-1));
2912 return u;
2913 }
2914 /* assume bnr has the right conductor */
2915 static GEN
get_lambda(GEN bnr)2916 get_lambda(GEN bnr)
2917 {
2918 GEN bnf = bnr_get_bnf(bnr), nf = bnf_get_nf(bnf), pol = nf_get_pol(nf);
2919 GEN f = gel(bnr_get_mod(bnr), 1), labas, lamodf, u;
2920 long a, b, f2, i, lu, v = varn(pol);
2921
2922 f2 = 2 * itos(gcoeff(f,1,1));
2923 u = getallrootsof1(bnf); lu = lg(u);
2924 for (i=1; i<lu; i++)
2925 gel(u,i) = ZC_hnfrem(gel(u,i), f); /* roots of 1, mod f */
2926 if (DEBUGLEVEL>1)
2927 err_printf("quadray: looking for [a,b] != unit mod 2f\n[a,b] = ");
2928 for (a=0; a<f2; a++)
2929 for (b=0; b<f2; b++)
2930 {
2931 GEN la = deg1pol_shallow(stoi(a), stoi(b), v); /* ax + b */
2932 if (umodiu(gnorm(mkpolmod(la, pol)), f2) != 1) continue;
2933 if (DEBUGLEVEL>1) err_printf("[%ld,%ld] ",a,b);
2934
2935 labas = poltobasis(nf, la);
2936 lamodf = ZC_hnfrem(labas, f);
2937 for (i=1; i<lu; i++)
2938 if (ZV_equal(lamodf, gel(u,i))) break;
2939 if (i < lu) continue; /* la = unit mod f */
2940 if (DEBUGLEVEL)
2941 {
2942 if (DEBUGLEVEL>1) err_printf("\n");
2943 err_printf("lambda = %Ps\n",la);
2944 }
2945 return labas;
2946 }
2947 pari_err_BUG("get_lambda");
2948 return NULL;/*LCOV_EXCL_LINE*/
2949 }
2950
2951 static GEN
to_approx(GEN nf,GEN a)2952 to_approx(GEN nf, GEN a)
2953 {
2954 GEN M = nf_get_M(nf);
2955 return gadd(gel(a,1), gmul(gcoeff(M,1,2),gel(a,2)));
2956 }
2957 /* Z-basis for a (over C) */
2958 static GEN
get_om(GEN nf,GEN a)2959 get_om(GEN nf, GEN a) {
2960 return mkvec2(to_approx(nf,gel(a,2)),
2961 to_approx(nf,gel(a,1)));
2962 }
2963
2964 /* Compute all elts in class group G = [|G|,c,g], c=cyclic factors, g=gens.
2965 * Set list[j + 1] = g1^e1...gk^ek where j is the integer
2966 * ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */
2967 static GEN
getallelts(GEN bnr)2968 getallelts(GEN bnr)
2969 {
2970 GEN nf, C, c, g, list, pows, gk;
2971 long lc, i, j, no;
2972
2973 nf = bnr_get_nf(bnr);
2974 no = itos( bnr_get_no(bnr) );
2975 c = bnr_get_cyc(bnr);
2976 g = bnr_get_gen_nocheck(bnr); lc = lg(c)-1;
2977 list = cgetg(no+1,t_VEC);
2978 gel(list,1) = matid(nf_get_degree(nf)); /* (1) */
2979 if (!no) return list;
2980
2981 pows = cgetg(lc+1,t_VEC);
2982 c = leafcopy(c); settyp(c, t_VECSMALL);
2983 for (i=1; i<=lc; i++)
2984 {
2985 long k = itos(gel(c,i));
2986 c[i] = k;
2987 gk = cgetg(k, t_VEC); gel(gk,1) = gel(g,i);
2988 for (j=2; j<k; j++)
2989 gel(gk,j) = idealmoddivisor(bnr, idealmul(nf, gel(gk,j-1), gel(gk,1)));
2990 gel(pows,i) = gk; /* powers of g[i] */
2991 }
2992
2993 C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc];
2994 for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1];
2995 /* C[i] = c(k-i+1) * ... * ck */
2996 /* j < C[i+1] <==> j only involves g(k-i)...gk */
2997 i = 1;
2998 for (j=1; j < C[1]; j++)
2999 gel(list, j+1) = gmael(pows,lc,j);
3000 while(j<no)
3001 {
3002 long k;
3003 GEN a;
3004 if (j == C[i+1]) i++;
3005 a = gmael(pows,lc-i,j/C[i]);
3006 k = j%C[i] + 1;
3007 if (k > 1) a = idealmoddivisor(bnr, idealmul(nf, a, gel(list,k)));
3008 gel(list, ++j) = a;
3009 }
3010 return list;
3011 }
3012
3013 /* x quadratic integer (approximate), recognize it. If error return NULL */
3014 static GEN
findbezk(GEN nf,GEN x)3015 findbezk(GEN nf, GEN x)
3016 {
3017 GEN a,b, M = nf_get_M(nf), u = gcoeff(M,1,2);
3018 long ea, eb;
3019
3020 /* u t_COMPLEX generator of nf.zk, write x ~ a + b u, a,b in Z */
3021 b = grndtoi(mpdiv(imag_i(x), gel(u,2)), &eb);
3022 if (eb > -20) return NULL;
3023 a = grndtoi(mpsub(real_i(x), mpmul(b,gel(u,1))), &ea);
3024 if (ea > -20) return NULL;
3025 return signe(b)? coltoalg(nf, mkcol2(a,b)): a;
3026 }
3027
3028 static GEN
findbezk_pol(GEN nf,GEN x)3029 findbezk_pol(GEN nf, GEN x)
3030 {
3031 long i, lx = lg(x);
3032 GEN y = cgetg(lx,t_POL);
3033 for (i=2; i<lx; i++)
3034 if (! (gel(y,i) = findbezk(nf,gel(x,i))) ) return NULL;
3035 y[1] = x[1]; return y;
3036 }
3037
3038 /* P approximation computed at initial precision prec. Compute needed prec
3039 * to know P with 1 word worth of trailing decimals */
3040 static long
get_prec(GEN P,long prec)3041 get_prec(GEN P, long prec)
3042 {
3043 long k = gprecision(P);
3044 if (k == 3) return precdbl(prec); /* approximation not trustworthy */
3045 k = prec - k; /* lost precision when computing P */
3046 if (k < 0) k = 0;
3047 k += nbits2prec(gexpo(P) + 128);
3048 if (k <= prec) k = precdbl(prec); /* dubious: old prec should have worked */
3049 return k;
3050 }
3051
3052 /* Compute data for ellphist */
3053 static GEN
ellphistinit(GEN om,long prec)3054 ellphistinit(GEN om, long prec)
3055 {
3056 GEN res,om1b,om2b, om1 = gel(om,1), om2 = gel(om,2);
3057
3058 if (gsigne(imag_i(gdiv(om1,om2))) < 0) { swap(om1,om2); om = mkvec2(om1,om2); }
3059 om1b = conj_i(om1);
3060 om2b = conj_i(om2); res = cgetg(4,t_VEC);
3061 gel(res,1) = gdivgs(elleisnum(om,2,0,prec),12);
3062 gel(res,2) = gdiv(PiI2(prec), gmul(om2, imag_i(gmul(om1b,om2))));
3063 gel(res,3) = om2b; return res;
3064 }
3065
3066 /* Computes log(phi^*(z,om)), using res computed by ellphistinit */
3067 static GEN
ellphist(GEN om,GEN res,GEN z,long prec)3068 ellphist(GEN om, GEN res, GEN z, long prec)
3069 {
3070 GEN u = imag_i(gmul(z, gel(res,3)));
3071 GEN zst = gsub(gmul(u, gel(res,2)), gmul(z,gel(res,1)));
3072 return gsub(ellsigma(om,z,1,prec),gmul2n(gmul(z,zst),-1));
3073 }
3074
3075 /* Computes phi^*(la,om)/phi^*(1,om) where (1,om) is an oriented basis of the
3076 ideal gf*gc^{-1} */
3077 static GEN
computeth2(GEN om,GEN la,long prec)3078 computeth2(GEN om, GEN la, long prec)
3079 {
3080 GEN p1,p2,res = ellphistinit(om,prec);
3081
3082 p1 = gsub(ellphist(om,res,la,prec), ellphist(om,res,gen_1,prec));
3083 p2 = imag_i(p1);
3084 if (gexpo(real_i(p1))>20 || gexpo(p2)> prec2nbits(minss(prec,realprec(p2)))-10)
3085 return NULL;
3086 return gexp(p1,prec);
3087 }
3088
3089 /* Computes P_2(X)=polynomial in Z_K[X] closest to prod_gc(X-th2(gc)) where
3090 the product is over the ray class group bnr.*/
3091 static GEN
computeP2(GEN bnr,long prec)3092 computeP2(GEN bnr, long prec)
3093 {
3094 long clrayno, i, first = 1;
3095 pari_sp av=avma, av2;
3096 GEN listray, P0, P, lanum, la = get_lambda(bnr);
3097 GEN nf = bnr_get_nf(bnr), f = gel(bnr_get_mod(bnr), 1);
3098 listray = getallelts(bnr);
3099 clrayno = lg(listray)-1; av2 = avma;
3100 PRECPB:
3101 if (!first)
3102 {
3103 if (DEBUGLEVEL) pari_warn(warnprec,"computeP2",prec);
3104 nf = gerepilecopy(av2, nfnewprec_shallow(checknf(bnr),prec));
3105 }
3106 first = 0; lanum = to_approx(nf,la);
3107 P = cgetg(clrayno+1,t_VEC);
3108 for (i=1; i<=clrayno; i++)
3109 {
3110 GEN om = get_om(nf, idealdiv(nf,f,gel(listray,i)));
3111 GEN s = computeth2(om,lanum,prec);
3112 if (!s) { prec = precdbl(prec); goto PRECPB; }
3113 gel(P,i) = s;
3114 }
3115 P0 = roots_to_pol(P, 0);
3116 P = findbezk_pol(nf, P0);
3117 if (!P) { prec = get_prec(P0, prec); goto PRECPB; }
3118 return gerepilecopy(av, P);
3119 }
3120
3121 #define nexta(a) (a>0 ? -a : 1-a)
3122 static GEN
do_compo(GEN A0,GEN B)3123 do_compo(GEN A0, GEN B)
3124 {
3125 long a, i, l = lg(B), v = fetch_var_higher();
3126 GEN A, z;
3127 /* now v > x = pol_x(0) > nf variable */
3128 B = leafcopy(B); setvarn(B, v);
3129 for (i = 2; i < l; i++) gel(B,i) = monomial(gel(B,i), l-i-1, 0);
3130 /* B := x^deg(B) B(v/x) */
3131 A = A0 = leafcopy(A0); setvarn(A0, v);
3132 for (a = 0;; a = nexta(a))
3133 {
3134 if (a) A = RgX_translate(A0, stoi(a));
3135 z = resultant(A,B); /* in variable 0 */
3136 if (issquarefree(z)) break;
3137 }
3138 (void)delete_var(); return z;
3139 }
3140 #undef nexta
3141
3142 static GEN
galoisapplypol(GEN nf,GEN s,GEN x)3143 galoisapplypol(GEN nf, GEN s, GEN x)
3144 {
3145 long i, lx = lg(x);
3146 GEN y = cgetg(lx,t_POL);
3147
3148 for (i=2; i<lx; i++) gel(y,i) = galoisapply(nf,s,gel(x,i));
3149 y[1] = x[1]; return y;
3150 }
3151 /* x quadratic, write it as ua + v, u,v rational */
3152 static GEN
findquad(GEN a,GEN x,GEN p)3153 findquad(GEN a, GEN x, GEN p)
3154 {
3155 long tu, tv;
3156 pari_sp av = avma;
3157 GEN u,v;
3158 if (typ(x) == t_POLMOD) x = gel(x,2);
3159 if (typ(a) == t_POLMOD) a = gel(a,2);
3160 u = poldivrem(x, a, &v);
3161 u = simplify_shallow(u); tu = typ(u);
3162 v = simplify_shallow(v); tv = typ(v);
3163 if (!is_scalar_t(tu)) pari_err_TYPE("findquad", u);
3164 if (!is_scalar_t(tv)) pari_err_TYPE("findquad", v);
3165 x = deg1pol(u, v, varn(a));
3166 if (typ(x) == t_POL) x = gmodulo(x,p);
3167 return gerepileupto(av, x);
3168 }
3169 static GEN
findquad_pol(GEN p,GEN a,GEN x)3170 findquad_pol(GEN p, GEN a, GEN x)
3171 {
3172 long i, lx = lg(x);
3173 GEN y = cgetg(lx,t_POL);
3174 for (i=2; i<lx; i++) gel(y,i) = findquad(a, gel(x,i), p);
3175 y[1] = x[1]; return y;
3176 }
3177 static GEN
compocyclo(GEN nf,long m,long d)3178 compocyclo(GEN nf, long m, long d)
3179 {
3180 GEN sb,a,b,s,p1,p2,p3,res,polL,polLK,nfL, D = nf_get_disc(nf);
3181 long ell,vx;
3182
3183 p1 = quadhilbertimag(D);
3184 p2 = polcyclo(m,0);
3185 if (d==1) return do_compo(p1,p2);
3186
3187 ell = m&1 ? m : (m>>2);
3188 if (absequalui(ell,D)) /* ell = |D| */
3189 {
3190 p2 = gcoeff(nffactor(nf,p2),1,1);
3191 return do_compo(p1,p2);
3192 }
3193 if (ell%4 == 3) ell = -ell;
3194 /* nf = K = Q(a), L = K(b) quadratic extension = Q(t) */
3195 polLK = quadpoly(stoi(ell)); /* relative polynomial */
3196 res = rnfequation2(nf, polLK);
3197 vx = nf_get_varn(nf);
3198 polL = gsubst(gel(res,1),0,pol_x(vx)); /* = charpoly(t) */
3199 a = gsubst(lift_shallow(gel(res,2)), 0,pol_x(vx));
3200 b = gsub(pol_x(vx), gmul(gel(res,3), a));
3201 nfL = nfinit(polL, DEFAULTPREC);
3202 p1 = gcoeff(nffactor(nfL,p1),1,1);
3203 p2 = gcoeff(nffactor(nfL,p2),1,1);
3204 p3 = do_compo(p1,p2); /* relative equation over L */
3205 /* compute non trivial s in Gal(L / K) */
3206 sb= gneg(gadd(b, RgX_coeff(polLK,1))); /* s(b) = Tr(b) - b */
3207 s = gadd(pol_x(vx), gsub(sb, b)); /* s(t) = t + s(b) - b */
3208 p3 = gmul(p3, galoisapplypol(nfL, s, p3));
3209 return findquad_pol(nf_get_pol(nf), a, p3);
3210 }
3211
3212 /* I integral ideal in HNF. (x) = I, x small in Z ? */
3213 static long
isZ(GEN I)3214 isZ(GEN I)
3215 {
3216 GEN x = gcoeff(I,1,1);
3217 if (signe(gcoeff(I,1,2)) || !equalii(x, gcoeff(I,2,2))) return 0;
3218 return is_bigint(x)? -1: itos(x);
3219 }
3220
3221 /* Treat special cases directly. return NULL if not special case */
3222 static GEN
treatspecialsigma(GEN bnr)3223 treatspecialsigma(GEN bnr)
3224 {
3225 GEN bnf = bnr_get_bnf(bnr), nf = bnf_get_nf(bnf);
3226 GEN f = gel(bnr_get_mod(bnr), 1), D = nf_get_disc(nf);
3227 GEN p1, p2;
3228 long Ds, fl, tryf, i = isZ(f);
3229
3230 if (i == 1) return quadhilbertimag(D); /* f = 1 */
3231
3232 if (absequaliu(D,3)) /* Q(j) */
3233 {
3234 if (i == 4 || i == 5 || i == 7) return polcyclo(i,0);
3235 if (!absequaliu(gcoeff(f,1,1),9) || !absequaliu(Z_content(f),3)) return NULL;
3236 /* f = P_3^3 */
3237 p1 = mkpolmod(bnf_get_tuU(bnf), nf_get_pol(nf));
3238 return gadd(pol_xn(3,0), p1); /* x^3+j */
3239 }
3240 if (absequaliu(D,4)) /* Q(i) */
3241 {
3242 if (i == 3 || i == 5) return polcyclo(i,0);
3243 if (i != 4) return NULL;
3244 p1 = mkpolmod(bnf_get_tuU(bnf), nf_get_pol(nf));
3245 return gadd(pol_xn(2,0), p1); /* x^2+i */
3246 }
3247 Ds = smodis(D,48);
3248 if (i)
3249 {
3250 if (i==2 && Ds%16== 8) return compocyclo(nf, 4,1);
3251 if (i==3 && Ds% 3== 1) return compocyclo(nf, 3,1);
3252 if (i==4 && Ds% 8== 1) return compocyclo(nf, 4,1);
3253 if (i==6 && Ds ==40) return compocyclo(nf,12,1);
3254 return NULL;
3255 }
3256
3257 p1 = gcoeff(f,1,1); /* integer > 0 */
3258 tryf = itou_or_0(p1); if (!tryf) return NULL;
3259 p2 = gcoeff(f,2,2); /* integer > 0 */
3260 if (is_pm1(p2)) fl = 0;
3261 else {
3262 if (Ds % 16 != 8 || !absequaliu(Z_content(f),2)) return NULL;
3263 fl = 1; tryf >>= 1;
3264 }
3265 if (tryf <= 3 || umodiu(D, tryf) || !uisprime(tryf)) return NULL;
3266 if (fl) tryf <<= 2;
3267 return compocyclo(nf,tryf,2);
3268 }
3269
3270 GEN
quadray(GEN D,GEN f,long prec)3271 quadray(GEN D, GEN f, long prec)
3272 {
3273 GEN bnr, y, bnf, H = NULL;
3274 pari_sp av = avma;
3275
3276 if (isint1(f)) return quadhilbert(D, prec);
3277 quadray_init(&D, f, &bnf, prec);
3278 bnr = Buchray(bnf, f, nf_INIT|nf_GEN);
3279 if (is_pm1(bnr_get_no(bnr))) { set_avma(av); return pol_x(0); }
3280 if (signe(D) > 0)
3281 y = bnrstark(bnr, H, prec);
3282 else
3283 {
3284 bnr_subgroup_sanitize(&bnr, &H);
3285 y = treatspecialsigma(bnr);
3286 if (!y) y = computeP2(bnr, prec);
3287 }
3288 return gerepileupto(av, y);
3289 }
3290