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