Annotation of OpenXM_contrib/pari-2.2/src/modules/stark.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: stark.c,v 1.39 2001/10/01 14:53:18 xavier Exp $
2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /*******************************************************************/
17: /* */
18: /* COMPUTATION OF STARK UNITS */
19: /* OF TOTALLY REAL FIELDS */
20: /* */
21: /*******************************************************************/
22: #include "pari.h"
23: #include "parinf.h"
24:
25: #define EXTRA_PREC (DEFAULTPREC-1)
26: #define ADD_PREC (DEFAULTPREC-2)*3
27:
28: extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus);
29:
30: /********************************************************************/
31: /* Miscellaneous functions */
32: /********************************************************************/
33:
34: /* Compute the image of logelt by chi as a complex number if flag = 0,
35: otherwise as a polmod, see InitChar in part 3 */
36: static GEN
37: ComputeImagebyChar(GEN chi, GEN logelt, long flag)
38: {
39: GEN gn = gmul((GEN)chi[1], logelt), x = (GEN)chi[flag? 4: 2];
40: long d = itos((GEN)chi[3]), n = smodis(gn, d);
41: /* x^d = 1 and, if d even, x^(d/2) = -1 */
42: if ((d & 1) == 0)
43: {
44: d /= 2;
45: if (n >= d) return gneg(gpowgs(x, n-d));
46: }
47: return gpowgs(x, n);
48: }
49:
50: /* Compute the conjugate character */
51: static GEN
52: ConjChar(GEN chi, GEN cyc)
53: {
54: long i, l = lg(chi);
55: GEN p1 = cgetg(l, t_VEC);
56:
57: for (i = 1; i < l; i++)
58: if (!signe((GEN)chi[i]))
59: p1[i] = zero;
60: else
61: p1[i] = lsubii((GEN)cyc[i], (GEN)chi[i]);
62:
63: return p1;
64: }
65:
66: /* compute the next element for FindEltofGroup */
67: static GEN
68: NextEltofGroup(GEN cyc, long l, long adec)
69: {
70: GEN p1;
71: long dj, j;
72:
73: p1 = cgetg(l + 1, t_COL);
74:
75: for (j = 1; j <= l; j++)
76: {
77: dj = itos((GEN)cyc[j]);
78: p1[j] = lstoi(adec%dj);
79: adec /= dj;
80: }
81:
82: return p1;
83: }
84:
85: /* Compute all the elements of a group given by its SNF */
86: static GEN
87: FindEltofGroup(long order, GEN cyc)
88: {
89: long l, i;
90: GEN rep;
91:
92: l = lg(cyc)-1;
93:
94: rep = cgetg(order + 1, t_VEC);
95:
96: for (i = 1; i <= order; i++)
97: rep[i] = (long)NextEltofGroup(cyc, l, i);
98:
99: return rep;
100: }
101:
102: /* Let dataC as given by InitQuotient0, compute a system of
103: representatives of the quotient */
104: static GEN
105: ComputeLift(GEN dataC)
106: {
107: long order, i, av = avma;
108: GEN cyc, surj, eltq, elt;
109:
110: order = itos((GEN)dataC[1]);
111: cyc = (GEN)dataC[2];
112: surj = (GEN)dataC[3];
113:
114: eltq = FindEltofGroup(order, cyc);
115:
116: elt = cgetg(order + 1, t_VEC);
117:
118: for (i = 1; i <= order; i++)
119: elt[i] = (long)inverseimage(surj, (GEN)eltq[i]);
120:
121: return gerepileupto(av, elt);
122: }
123:
124: /* Let bnr1, bnr2 be such that mod(bnr2) | mod(bnr1), compute the
125: matrix of the surjective map Cl(bnr1) ->> Cl(bnr2) */
126: static GEN
127: GetSurjMat(GEN bnr1, GEN bnr2)
128: {
129: long nbg, i;
130: GEN gen, M;
131:
132: gen = gmael(bnr1, 5, 3);
133: nbg = lg(gen) - 1;
134:
135: M = cgetg(nbg + 1, t_MAT);
136: for (i = 1; i <= nbg; i++)
137: M[i] = (long)isprincipalray(bnr2, (GEN)gen[i]);
138:
139: return M;
140: }
141:
142: /* A character is given by a vector [(c_i), z, d, pm] such that
143: chi(id) = z ^ sum(c_i * a_i) where
144: a_i= log(id) on the generators of bnr
145: z = exp(2i * Pi / d)
146: pm = z as a polmod */
147: static GEN
148: get_Char(GEN chi, long prec)
149: {
150: GEN p2, d, _2ipi = gmul(gi, shiftr(mppi(prec), 1));
151: p2 = cgetg(5, t_VEC); d = denom(chi);
152: p2[1] = lmul(d, chi);
153: if (egalii(d, gdeux))
154: p2[2] = lstoi(-1);
155: else
156: p2[2] = lexp(gdiv(_2ipi, d), prec);
157: p2[3] = (long)d;
158: p2[4] = lmodulcp(polx[0], cyclo(itos(d), 0));
159: return p2;
160: }
161:
162: /* Let chi a character defined over bnr and primitif over bnrc,
163: compute the corresponding primitive character and the vectors of
164: prime ideals dividing bnr but not bnr. Returns NULL if bnr = bnrc */
165: static GEN
166: GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
167: {
168: long nbg, i, j, l, av = avma, nd;
169: GEN gen, cyc, U, chic, M, s, p1, cond, condc, p2, nf;
170: GEN prdiff, Mrc;
171:
172: cond = gmael(bnr, 2, 1);
173: condc = gmael(bnrc, 2, 1);
174: if (gegal(cond, condc)) return NULL;
175:
176: gen = gmael(bnr, 5, 3);
177: nbg = lg(gen) - 1;
178: cyc = gmael(bnr, 5, 2);
179: Mrc = diagonal(gmael(bnrc, 5, 2));
180: nf = gmael(bnr, 1, 7);
181:
182: cond = (GEN)cond[1];
183: condc = (GEN)condc[1];
184:
185: M = GetSurjMat(bnr, bnrc);
186: l = lg((GEN)M[1]);
187: p1 = hnfall(concatsp(M, Mrc));
188: U = (GEN)p1[2];
189:
190: chic = cgetg(l, t_VEC);
191: for (i = 1; i < l; i++)
192: {
193: s = gzero; p1 = (GEN)U[i + nbg];
194: for (j = 1; j <= nbg; j++)
195: {
196: p2 = gdiv((GEN)p1[j], (GEN)cyc[j]);
197: s = gadd(s,gmul(p2,(GEN)chi[j]));
198: }
199: chic[i] = (long)s;
200: }
201:
202: p2 = (GEN)idealfactor(nf, cond)[1];
203: l = lg(p2);
204:
205: prdiff = cgetg(l, t_COL);
206: for (nd=1, i=1; i < l; i++)
207: if (!idealval(nf, condc, (GEN)p2[i])) prdiff[nd++] = p2[i];
208: setlg(prdiff, nd);
209:
210: p1 = cgetg(3, t_VEC);
211: p1[1] = (long)get_Char(chic,prec);
212: p1[2] = lcopy(prdiff);
213:
214: return gerepileupto(av,p1);
215: }
216:
217: /* Let dataCR be a list of characters, compute the image of logelt */
218: static GEN
219: chiideal(GEN dataCR, GEN logelt, long flag)
220: {
221: long j, l = lg(dataCR);
222: GEN rep = cgetg(l, t_VEC);
223:
224: for (j = 1; j < l; j++)
225: rep[j] = (long)ComputeImagebyChar(gmael(dataCR, j, 5), logelt, flag);
226:
227: return rep;
228: }
229:
230: static GEN
231: GetDeg(GEN dataCR)
232: {
233: long i, l = lg(dataCR);
234: GEN degs = cgetg(l, t_VECSMALL);
235:
236: for (i = 1; i < l; i++)
237: degs[i] = degpol(gmael4(dataCR, i, 5, 4, 1));
238: return degs;
239: }
240:
241: /********************************************************************/
242: /* 1rst part: find the field K */
243: /********************************************************************/
244:
245: static GEN AllStark(GEN data, GEN nf, long flag, long prec);
246: static GEN InitChar0(GEN dataD, long prec);
247:
248: /* Let A be a finite abelian group given by its relation and let C
249: define a subgroup of A, compute the order of A / C, its structure and
250: the matrix expressing the generators of A on those of A / C */
251: static GEN
252: InitQuotient0(GEN DA, GEN C)
253: {
254: long i, l;
255: GEN MQ, MrC, H, snf, snf2, rep, p1;
256:
257: l = lg(DA)-1;
258: H = gcmp0(C)? DA: C;
259: MrC = gauss(H, DA);
260: snf = smith2(hnf(MrC));
261: MQ = concatsp(gmul(H, (GEN)snf[1]), DA);
262: snf2 = smith2(hnf(MQ));
263:
264: rep = cgetg(5, t_VEC);
265:
266: p1 = cgetg(l + 1, t_VEC);
267: for (i = 1; i <= l; i++)
268: p1[i] = lcopy(gcoeff((GEN)snf2[3], i, i));
269:
270: rep[1] = (long)dethnf((GEN)snf2[3]);
271: rep[2] = (long)p1;
272: rep[3] = lcopy((GEN)snf2[1]);
273: rep[4] = lcopy(C);
274:
275: return rep;
276: }
277:
278: /* Let m be a modulus et C a subgroup of Clk(m), compute all the data
279: * needed to work with the quotient Clk(m) / C namely
280: * 1) bnr(m)
281: * 2.1) its order
282: * 2.2) its structure
283: * 2.3) the matrix Clk(m) ->> Clk(m)/C
284: * 2.4) the group C */
285: static GEN
286: InitQuotient(GEN bnr, GEN C)
287: {
288: GEN Mrm, dataquo = cgetg(3, t_VEC);
289: long av;
290:
291: dataquo[1] = lcopy(bnr);
292: av = avma; Mrm = diagonal(gmael(bnr, 5, 2));
293: dataquo[2] = lpileupto(av, InitQuotient0(Mrm, C));
294:
295: return dataquo;
296: }
297:
298: /* Let s: A -> B given by P, and let DA, DB be resp. the matrix of the
299: relations of A and B, let nbA, nbB be resp. the rank of A and B,
300: compute the kernel of s. If DA = 0 then A is free */
301: static GEN
302: ComputeKernel0(GEN P, GEN DA, GEN DB, long nbA, long nbB)
303: {
304: long rk, av = avma;
305: GEN herm, mask1, mask2, U;
306:
307: herm = hnfall(concatsp(P, DB));
308: rk = nbA + nbB + 1;
309: rk -= lg((GEN)herm[1]); /* two steps: bug in pgcc 1.1.3 inlining (IS) */
310:
311: mask1 = subis(shifti(gun, nbA), 1);
312: mask2 = subis(shifti(gun, rk), 1);
313:
314: U = matextract((GEN)herm[2], mask1, mask2);
315:
316: if (!gcmp0(DA)) U = concatsp(U, DA);
317: return gerepileupto(av, hnf(U));
318: }
319:
320: /* Let m and n be two moduli such that n|m and let C be a congruence
321: group modulo n, compute the corresponding congruence group modulo m
322: ie then kernel of the map Clk(m) ->> Clk(n)/C */
323: static GEN
324: ComputeKernel(GEN bnrm, GEN dataC)
325: {
326: long av = avma, i, nbm, nbq;
327: GEN bnrn, Mrm, genm, Mrq, mgq, P;
328:
329: bnrn = (GEN)dataC[1];
330: Mrm = diagonal(gmael(bnrm, 5, 2));
331: genm = gmael(bnrm, 5, 3);
332: nbm = lg(genm) - 1;
333: Mrq = diagonal(gmael(dataC, 2, 2));
334: mgq = gmael(dataC, 2, 3);
335: nbq = lg(mgq) - 1;
336:
337: P = cgetg(nbm + 1, t_MAT);
338: for (i = 1; i <= nbm; i++)
339: P[i] = lmul(mgq, isprincipalray(bnrn, (GEN)genm[i]));
340:
341: return gerepileupto(av, ComputeKernel0(P, Mrm, Mrq, nbm, nbq));
342: }
343:
344: /* Let C a congruence group in bnr, compute its subgroups of index 2 as
345: subgroups of Clk(bnr) */
346: static GEN
347: ComputeIndex2Subgroup(GEN bnr, GEN C)
348: {
349: long nb, i, l, av = avma;
350: GEN snf, Mr, U, CU, subgrp, rep, p1;
351:
352: disable_dbg(0);
353:
354: Mr = diagonal(gmael(bnr, 5, 2));
355: snf = smith2(gmul(ginv(C), Mr));
356:
357: U = ginv((GEN)snf[1]);
358: l = lg((GEN)snf[3]);
359:
360: p1 = cgetg(l, t_VEC);
361: for (i = 1; i < l; i++)
362: p1[i] = mael3(snf, 3, i, i);
363:
364: subgrp = subgrouplist(p1, 2);
365: nb = lg(subgrp) - 1; CU = gmul(C,U);
366:
367: rep = cgetg(nb, t_VEC);
368: for (i = 1; i < nb; i++) /* skip Id which comes last */
369: rep[i] = (long)hnf(concatsp(gmul(CU, (GEN)subgrp[i]), Mr));
370:
371: disable_dbg(-1);
372: return gerepilecopy(av, rep);
373: }
374:
375: /* Let pr be a prime (pr may divide mod(bnr)), compute the indexes
376: e,f of the splitting of pr in the class field nf(bnr/subgroup) */
377: static GEN
378: GetIndex(GEN pr, GEN bnr, GEN subgroup)
379: {
380: long av = avma, v, lg, i;
381: GEN bnf, mod, mod0, mpr0, mpr, bnrpr, subpr, M, e, f, dtQ, p1;
382: GEN rep, cycpr, cycQ;
383:
384: bnf = (GEN)bnr[1];
385: mod = gmael(bnr, 2, 1);
386: mod0 = (GEN)mod[1];
387:
388: /* Compute the part of mod coprime to pr */
389: v = idealval(bnf, mod0, pr);
390: mpr0 = idealdivexact(bnf, mod0, idealpow(bnf, pr, stoi(v)));
391:
392: mpr = cgetg(3, t_VEC);
393: mpr[1] = (long)mpr0;
394: mpr[2] = mod[2];
395: if (gegal(mpr0, mod0))
396: {
397: bnrpr = bnr;
398: subpr = subgroup;
399: }
400: else
401: {
402: bnrpr = buchrayinitgen(bnf, mpr);
403: cycpr = gmael(bnrpr, 5, 2);
404: M = GetSurjMat(bnr, bnrpr);
405: M = gmul(M, subgroup);
406: subpr = hnf(concatsp(M, diagonal(cycpr)));
407: }
408:
409: /* e = #(bnr/subgroup) / #(bnrpr/subpr) */
410: e = gdiv(det(subgroup), det(subpr));
411:
412: /* f = order of [pr] in bnrpr/subpr */
413: dtQ = InitQuotient(bnrpr, subpr);
414: p1 = isprincipalray(bnrpr, pr);
415: p1 = gmul(gmael(dtQ, 2, 3), p1);
416: cycQ = gmael(dtQ, 2, 2);
417: lg = lg(cycQ) - 1;
418: f = gun;
419: for (i = 1; i <= lg; i++)
420: f = glcm(f, gdiv((GEN)cycQ[i], ggcd((GEN)cycQ[i], (GEN)p1[i])));
421:
422: rep = cgetg(3, t_VEC);
423: rep[1] = lcopy(e);
424: rep[2] = lcopy(f);
425:
426: return gerepileupto(av, rep);
427: }
428:
429:
430: /* Given a conductor and a subgroups, return the corresponding
431: complexity and precision required using quickpol */
432: static GEN
433: CplxModulus(GEN data, long *newprec, long prec)
434: {
435: long av = avma, pr, dprec;
436: GEN nf, cpl, pol, p1;
437:
438: nf = gmael3(data, 1, 1, 7);
439:
440: p1 = cgetg(6, t_VEC);
441:
442: p1[1] = data[1];
443: p1[2] = data[2];
444: p1[3] = data[3];
445: p1[4] = data[4];
446:
447: if (DEBUGLEVEL >= 2)
448: fprintferr("\nTrying modulus = %Z and subgroup = %Z\n",
449: mael3(p1, 1, 2, 1), (GEN)p1[2]);
450:
451: dprec = DEFAULTPREC;
452:
453: for (;;)
454: {
455: p1[5] = (long)InitChar0((GEN)data[3], dprec);
456: pol = AllStark(p1, nf, -1, dprec);
457: if (!gcmp0(leading_term(pol)))
458: {
459: cpl = mpabs(poldisc0(pol, 0));
460: if (!gcmp0(cpl)) break;
461: }
462: pr = gexpo(pol)>>(TWOPOTBITS_IN_LONG+1);
463: if (pr < 0) pr = 0;
464: dprec = max(dprec, pr) + EXTRA_PREC;
465:
466: if (DEBUGLEVEL >= 2) err(warnprec, "CplxModulus", dprec);
467: }
468:
469: if (DEBUGLEVEL >= 2) fprintferr("cpl = %Z\n", cpl);
470:
471: pr = gexpo(pol)>>TWOPOTBITS_IN_LONG;
472: if (pr < 0) pr = 0;
473: *newprec = max(prec, pr + EXTRA_PREC);
474:
475: return gerepileupto(av, cpl);
476: }
477:
478: /* Let f be a conductor without infinite part and let C be a
479: congruence group modulo f, compute (m,D) such that D is a
480: congruence group of conductor m where m is a multiple of f
481: divisible by all the infinite places but one, D is a subgroup of
482: index 2 of Im(C) in Clk(m), no prime dividing f splits in the
483: corresponding quadratic extension and m is of minimal norm. Return
484: bnr(m), D, quotient Ck(m) / D and Clk(m) / C */
485: /* If fl != 0, try bnd extra moduli */
486: static GEN
487: FindModulus(GEN dataC, long fl, long *newprec, long prec, long bnd)
488: {
489: long n, i, narch, nbp, maxnorm, minnorm, N, nbidnn, s, c, j, nbcand;
490: long limnorm, first = 1, pr;
491: ulong av = avma, av1, av0;
492: GEN bnr, rep, bnf, nf, f, arch, m, listid, idnormn, bnrm, ImC;
493: GEN candD, D, bpr, indpr, sgp, p1, p2, rb;
494:
495: bnr = (GEN)dataC[1];
496: sgp = gmael(dataC, 2, 4);
497: bnf = (GEN)bnr[1];
498: nf = (GEN)bnf[7];
499: N = degpol(nf[1]);
500: f = gmael3(bnr, 2, 1, 1);
501:
502: rep = cgetg(6, t_VEC);
503: for (i = 1; i <= 5; i++) rep[i] = zero;
504:
505: /* if cpl < rb, it is not necessary to try another modulus */
506: rb = powgi(gmul(gmael(bnf, 7, 3), det(f)), gmul2n(gmael(bnr, 5, 1), 3));
507:
508: bpr = (GEN)idealfactor(nf, f)[1];
509: nbp = lg(bpr) - 1;
510:
511: indpr = cgetg(nbp + 1,t_VEC);
512: for (i = 1; i <= nbp; i++)
513: {
514: p1 = GetIndex((GEN)bpr[i], bnr, sgp);
515: indpr[i] = lmulii((GEN)p1[1], (GEN)p1[2]);
516: }
517:
518: /* Initialization of the possible infinite part */
519: arch = cgetg(N+1, t_VEC);
520: for (i = 1; i <= N; i++) arch[i] = un;
521:
522: /* narch = (N == 2)? 1: N; -- if N=2, only one case is necessary */
523: narch = N;
524:
525: m = cgetg(3, t_VEC);
526: m[2] = (long)arch;
527:
528: /* we go from minnorm up to maxnorm, if necessary we increase these values.
529: If we cannot find a suitable conductor of norm < limnorm, we stop */
530: maxnorm = 50;
531: minnorm = 1;
532: limnorm = 200;
533:
534: if (DEBUGLEVEL >= 2)
535: fprintferr("Looking for a modulus of norm: ");
536:
537: av0 = avma;
538: for(;;)
539: {
540: /* compute all ideals of norm <= maxnorm */
541: disable_dbg(0);
542: listid = ideallist(nf, maxnorm);
543: disable_dbg(-1);
544: av1 = avma;
545:
546: for (n = minnorm; n <= maxnorm; n++)
547: {
548: if (DEBUGLEVEL >= 2) fprintferr(" %ld", n);
549:
550: idnormn = (GEN)listid[n];
551: nbidnn = lg(idnormn) - 1;
552: for (i = 1; i <= nbidnn; i++)
553: {
554: rep = gerepilecopy(av1, rep);
555:
556: /* finite part of the conductor */
557: m[1] = (long)idealmul(nf, f, (GEN)idnormn[i]);
558:
559: for (s = 1; s <= narch; s++)
560: {
561: /* infinite part */
562: arch[N+1-s] = zero;
563:
564: /* compute Clk(m), check if m is a conductor */
565: disable_dbg(0);
566: bnrm = buchrayinitgen(bnf, m);
567: p1 = conductor(bnrm, gzero, -1);
568: disable_dbg(-1);
569: if (signe(p1))
570: {
571: /* compute Im(C) in Clk(m)... */
572: ImC = ComputeKernel(bnrm, dataC);
573:
574: /* ... and its subgroups of index 2 */
575: candD = ComputeIndex2Subgroup(bnrm, ImC);
576: nbcand = lg(candD) - 1;
577: for (c = 1; c <= nbcand; c++)
578: {
579: /* check if m is the conductor */
580: D = (GEN)candD[c];
581: disable_dbg(0);
582: p1 = conductor(bnrm, D, -1);
583: disable_dbg(-1);
584: if (signe(p1))
585: {
586: /* check the splitting of primes */
587: for (j = 1; j <= nbp; j++)
588: {
589: p1 = GetIndex((GEN)bpr[j], bnrm, D);
590: p1 = mulii((GEN)p1[1], (GEN)p1[2]);
591: if (egalii(p1, (GEN)indpr[j])) break; /* no good */
592: }
593: if (j > nbp)
594: {
595: p2 = cgetg(6, t_VEC);
596:
597: p2[1] = lcopy(bnrm);
598: p2[2] = lcopy(D);
599: p2[3] = (long)InitQuotient((GEN)p2[1], (GEN)p2[2]);
600: p2[4] = (long)InitQuotient((GEN)p2[1], ImC);
601:
602: p1 = CplxModulus(p2, &pr, prec);
603:
604: if (first == 1 || gcmp(p1, (GEN)rep[5]) < 0)
605: {
606: *newprec = pr;
607: for (j = 1; j <= 4; j++) rep[j] = p2[j];
608: rep[5] = (long)p1;
609: }
610:
611: if (!fl || (gcmp(p1, rb) < 0))
612: {
613: rep[5] = (long)InitChar0((GEN)rep[3], *newprec);
614: return gerepilecopy(av, rep);
615: }
616: if (DEBUGLEVEL >= 2)
617: fprintferr("Trying to find another modulus...");
618: first--;
619: }
620: }
621: }
622: }
623: arch[N+1-s] = un;
624: }
625: if (first <= bnd)
626: {
627: if (DEBUGLEVEL >= 2)
628: fprintferr("No, we're done!\nModulus = %Z and subgroup = %Z\n",
629: gmael3(rep, 1, 2, 1), rep[2]);
630: rep[5] = (long)InitChar0((GEN)rep[3], *newprec);
631: return gerepilecopy(av, rep);
632: }
633: }
634: }
635: /* if necessary compute more ideals */
636: rep = gerepilecopy(av0, rep);
637:
638: minnorm = maxnorm;
639: maxnorm <<= 1;
640: if (maxnorm > limnorm)
641: err(talker, "Cannot find a suitable modulus in FindModulus");
642: }
643: }
644:
645: /********************************************************************/
646: /* 2nd part: compute W(X) */
647: /********************************************************************/
648:
649: /* compute W(chi) such that Ld(s,chi) = W(chi) Ld(1 - s, chi*),
650: if flag != 0 do not check the result */
651: static GEN
652: ComputeArtinNumber(GEN datachi, long flag, long prec)
653: {
654: long av = avma, av2, G, ms, j, i, nz, zcard, q, l, N, lim;
655: GEN chi, nc, dc, p1, cond0, cond1, elts, Msign, umod2, lambda, nf;
656: GEN sg, p2, chib, diff, vt, z, idg, mu, idh, zid, zstruc, zgen, zchi;
657: GEN classe, bnr, beta, s, tr, p3, den, muslambda, Pi, lp1, beta2;
658:
659: chi = (GEN)datachi[8];
660: /* trivial case */
661: if (cmpsi(2, (GEN)chi[3]) >= 0) return gun;
662:
663: bnr = (GEN)datachi[3];
664: nf = gmael(bnr, 1, 7);
665: diff = gmael(nf, 5, 5);
666: cond0 = gmael3(bnr, 2, 1, 1);
667: cond1 = gmael3(bnr, 2, 1, 2);
668: umod2 = gmodulcp(gun, gdeux);
669: N = degpol(nf[1]);
670: Pi = mppi(prec);
671:
672: G = - bit_accuracy(prec) >> 1;
673: nc = idealnorm(nf, cond0);
674: dc = idealmul(nf, diff, cond0);
675: den = idealnorm(nf, dc);
676: z = gexp(gdiv(gmul2n(gmul(gi, Pi), 1), den), prec);
677:
678: q = 0;
679: for (i = 1; i < lg(cond1); i++)
680: if (gcmp1((GEN)cond1[i])) q++;
681:
682: /* compute a system of elements congru to 1 mod cond0 and giving all
683: possible signatures for cond1 */
684: p1 = zarchstar(nf, cond0, cond1, q);
685: elts = (GEN)p1[2];
686: Msign = gmul((GEN)p1[3], umod2);
687: ms = lg(elts) - 1;
688:
689: /* find lambda in diff.cond such that gcd(lambda.(diff.cond)^-1,cond0) = 1
690: and lambda >(cond1)> 0 */
691: lambda = idealappr(nf, dc);
692: sg = zsigne(nf, lambda, cond1);
693: p2 = lift(gmul(Msign, sg));
694:
695: for (j = 1; j <= ms; j++)
696: if (gcmp1((GEN)p2[j])) lambda = element_mul(nf, lambda, (GEN)elts[j]);
697:
698: idg = idealdivexact(nf, lambda, dc);
699:
700: /* find mu in idg such that idh=(mu) / idg is coprime with cond0 and
701: mu >(cond1)> 0 */
702: if (!gcmp1(gcoeff(idg, 1, 1)))
703: {
704: p1 = idealfactor(nf, idg);
705: p2 = idealfactor(nf, cond0);
706:
707: l = lg((GEN)p2[1]);
708: for (i = 1; i < l; i++) coeff(p2, i, 2) = zero;
709:
710: p1 = gtrans(concatsp(gtrans(p1), gtrans(p2)));
711: mu = idealapprfact(nf, p1);
712: sg = zsigne(nf, mu, cond1);
713: p2 = lift(gmul(Msign, sg));
714:
715: for (j = 1; j <= ms; j++)
716: if (gcmp1((GEN)p2[j])) mu = element_mul(nf, mu, (GEN)elts[j]);
717:
718: idh = idealdivexact(nf, mu, idg);
719: }
720: else
721: {
722: mu = gun;
723: idh = gcopy(idg);
724: }
725:
726: muslambda = element_div(nf, mu, lambda);
727:
728: /* compute a system of generators of (Ok/cond)^* cond1-positive */
729: zid = zidealstarinitgen(nf, cond0);
730:
731: zcard = itos(gmael(zid, 2, 1));
732: zstruc = gmael(zid, 2, 2);
733: zgen = gmael(zid, 2, 3);
734: nz = lg(zgen) - 1;
735:
736: zchi = cgetg(nz + 1, t_VEC);
737: for (i = 1; i <= nz; i++)
738: {
739: p1 = (GEN)zgen[i];
740: sg = zsigne(nf, p1, cond1);
741: p2 = lift(gmul(Msign, sg));
742:
743: for (j = 1; j <= ms; j++)
744: if (gcmp1((GEN)p2[j])) p1 = element_mul(nf, p1, (GEN)elts[j]);
745:
746: classe = isprincipalray(bnr, p1);
747: zchi[i] = (long)ComputeImagebyChar(chi, classe, 0);
748: zgen[i] = (long)p1;
749: }
750:
751: /* Sum chi(beta) * exp(2i * Pi * Tr(beta * mu / lambda) where beta
752: runs through the classes of (Ok/cond0)^* and beta cond1-positive */
753:
754: p3 = cgetg(N + 1, t_COL);
755: for (i = 1; i <= N; i++) p3[i] = zero;
756:
757: vt = cgetg(N + 1, t_VEC);
758: for (i = 1; i <= N; i++)
759: {
760: p3[i] = un;
761: vt[i] = ltrace(basistoalg(nf, p3));
762: p3[i] = zero;
763: }
764:
765: lp1 = NULL;
766: s = gzero;
767:
768: av2 = avma; lim = stack_lim(av2, 1);
769:
770: for (i = 0; i < zcard; i++)
771: {
772: p1 = NextEltofGroup(zstruc, nz, i);
773:
774: /* we test if we can use the previous value
775: of beta / chib to compute the next one */
776: /* FIXME: there should be a better way of doing that! */
777: if (!lp1 || !gcmp1(gnorml2(gsub(p1, lp1))))
778: {
779: beta = gun;
780: chib = gun;
781:
782: for (j = 1; j <= nz; j++)
783: {
784: if (!gcmp0((GEN)p1[j]))
785: {
786: p2 = element_powmodideal(nf, (GEN)zgen[j], (GEN)p1[j], cond0);
787: beta = element_mulmodideal(nf, beta, p2, cond0);
788: chib = gmul(chib, powgi((GEN)zchi[j], (GEN)p1[j]));
789: }
790: }
791: }
792: else
793: {
794: /* we use the fact that NextEltofGroup is, in general,
795: obtained by adding 1 to a component of p1 */
796: for (j = 1; j <= nz; j++)
797: if (!gegal((GEN)p1[j], (GEN)lp1[j]))
798: {
799: beta = element_mulmodideal(nf, beta, (GEN)zgen[j], cond0);
800: chib = gmul(chib, (GEN)zchi[j]);
801: }
802: }
803:
804: lp1 = p1;
805: sg = zsigne(nf, beta, cond1);
806: p2 = lift(gmul(Msign, sg));
807:
808: for (j = 1; j <= ms; j++)
809: if (gcmp1((GEN)p2[j]))
810: beta = element_mul(nf, beta, (GEN)elts[j]);
811:
812: beta2 = element_mul(nf, beta, muslambda);
813: tr = gmul(vt, beta2);
814: tr = gmod(gmul(tr, den), den);
815:
816: s = gadd(s, gmul(chib, powgi(z,tr)));
817:
818: if (low_stack(lim, stack_lim(av2, 1)))
819: {
820: GEN *gptr[5];
821: gptr[0] = &s; gptr[1] = &lp1; gptr[2] = β gptr[3] = &chib;
822: if (DEBUGMEM > 1) err(warnmem,"ComputeArtinNumber");
823: gerepilemany(av2, gptr, 4);
824: }
825: }
826:
827: classe = isprincipalray(bnr, idh);
828: s = gmul(s, ComputeImagebyChar(chi, classe, 0));
829: s = gdiv(s, gsqrt(nc, prec));
830:
831: p1 = gsubgs(gabs(s, prec), 1);
832:
833: i = expo(p1);
834: if ((i > G) && !flag)
835: err(bugparier, "ComputeArtinNumber");
836:
837: return gerepileupto(av, gmul(s, gpowgs(gneg_i(gi),q)));
838: }
839:
840: /* compute the constant W of the functional equation of
841: Lambda(chi). If flag = 1 then chi is assumed to be primitive */
842: GEN
843: bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
844: {
845: long av = avma, l, i;
846: GEN cond, condc, bnrc, chic, cyc, d, p1, p2, dtcr, Pi;
847:
848: if ((flag < 0) || (flag > 1))
849: err(flagerr,"bnrrootnumber");
850:
851: checkbnr(bnr);
852:
853: cond = gmael(bnr, 2, 1);
854: l = lg(gmael(bnr, 5, 2));
855: Pi = mppi(prec);
856:
857: if ((typ(chi) != t_VEC) || (lg(chi) != l))
858: err(talker, "incorrect character in bnrrootnumber");
859:
860: if (!flag)
861: {
862: condc = bnrconductorofchar(bnr, chi);
863: if (gegal(cond, condc)) flag = 1;
864: }
865: else condc = cond;
866:
867: if (flag)
868: bnrc = bnr;
869: else
870: bnrc = buchrayinitgen((GEN)bnr[1], condc);
871:
872: chic = cgetg(l, t_VEC);
873: cyc = gmael(bnr, 5, 2);
874: for (i = 1; i < l; i++)
875: chic[i] = ldiv((GEN)chi[i], (GEN)cyc[i]);
876: d = denom(chic);
877:
878: p2 = cgetg(4, t_VEC);
879: p2[1] = lmul(d, chic);
880: if (egalii(d, gdeux))
881: p2[2] = lstoi(-1);
882: else
883: p2[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), d), prec);
884: p2[3] = (long)d;
885:
886: dtcr = cgetg(9, t_VEC);
887: dtcr[1] = (long)chi;
888: dtcr[2] = zero;
889: dtcr[3] = (long)bnrc;
890: dtcr[4] = (long)bnr;
891: dtcr[5] = (long)p2;
892: dtcr[6] = zero;
893: dtcr[7] = (long)condc;
894:
895: p1 = GetPrimChar(chi, bnr, bnrc, prec);
896:
897: if (!p1)
898: dtcr[8] = (long)p2;
899: else
900: dtcr[8] = p1[1];
901:
902: return gerepileupto(av, ComputeArtinNumber(dtcr, 0, prec));
903: }
904:
905: /********************************************************************/
906: /* 3rd part: initialize the characters */
907: /********************************************************************/
908:
909: static GEN
910: LiftChar(GEN cyc, GEN Mat, GEN chi)
911: {
912: long lm, l, i, j, av;
913: GEN lchi, s;
914:
915: lm = lg(cyc) - 1;
916: l = lg(chi) - 1;
917:
918: lchi = cgetg(lm + 1, t_VEC);
919: for (i = 1; i <= lm; i++)
920: {
921: av = avma;
922: s = gzero;
923:
924: for (j = 1; j <= l; j++)
925: s = gadd(s, gmul((GEN)chi[j], gcoeff(Mat, j, i)));
926:
927: lchi[i] = (long)gerepileupto(av, gmod(gmul(s, (GEN)cyc[i]), (GEN)cyc[i]));
928: }
929:
930: return lchi;
931: }
932:
933: /* Let chi be a character, A(chi) corresponding to the primes dividing diff
934: at s = flag. If s = 0, returns [r, A] where r is the order of vanishing
935: at s = 0 corresponding to diff. No Garbage collector */
936: static GEN
937: ComputeAChi(GEN dtcr, long flag, long prec)
938: {
939: long l, i;
940: GEN p1, ray, r, A, rep, diff, chi, bnrc;
941:
942: diff = (GEN)dtcr[6];
943: bnrc = (GEN)dtcr[3];
944: chi = (GEN)dtcr[8];
945: l = lg(diff) - 1;
946:
947: A = gun;
948: r = gzero;
949:
950: for (i = 1; i <= l; i++)
951: {
952: ray = isprincipalray(bnrc, (GEN)diff[i]);
953: p1 = ComputeImagebyChar(chi, ray, 0);
954:
955: if (flag)
956: A = gmul(A, gsub(gun, gdiv(p1, idealnorm((GEN)bnrc[1], (GEN)diff[i]))));
957: else
958: {
959: if (gcmp1(p1))
960: {
961: r = addis(r, 1);
962: A = gmul(A, glog(idealnorm((GEN)bnrc[1], (GEN)diff[i]), prec));
963: }
964: else
965: A = gmul(A, gsub(gun, p1));
966: }
967: }
968:
969: if (flag) return A;
970:
971: rep = cgetg(3, t_VEC);
972: rep[1] = (long)r;
973: rep[2] = (long)A;
974:
975: return rep;
976: }
977:
978: /* Given a list [chi, cond(chi)] of characters over Cl(bnr), compute a
979: vector dataCR containing for each character:
980: 1: chi
981: 2: the constant C(chi)
982: 3: bnr(cond(chi))
983: 4: bnr(m)
984: 5: [(c_i), z, d, pm] in bnr(m)
985: 6: diff(chi) primes dividing m but not cond(chi)
986: 7: finite part of cond(chi)
987: 8: [(c_i), z, d, pm] in bnr(cond(chi))
988: 9: [q, r1 - q, r2, rc] where
989: q = number of real places in cond(chi)
990: rc = max{r1 + r2 - q + 1, r2 + q} */
991: static GEN
992: InitChar(GEN bnr, GEN listCR, long prec)
993: {
994: GEN bnf = checkbnf(bnr), nf = checknf(bnf);
995: GEN modul, dk, C, dataCR, chi, cond, Mr, chic, gr2, p1, p2, q;
996: long N, r1, r2, prec2, h, i, j, nbg, av = avma;
997:
998: modul = gmael(bnr, 2, 1);
999: Mr = gmael(bnr, 5, 2);
1000: nbg = lg(Mr) - 1;
1001: dk = (GEN)nf[3];
1002: N = degpol(nf[1]);
1003: r1 = nf_get_r1(nf);
1004: r2 = (N - r1) >> 1; gr2 = stoi(r2);
1005: prec2 = ((prec - 2)<<1) + EXTRA_PREC;
1006: C = gmul2n(gsqrt(gdiv(absi(dk), gpowgs(mppi(prec2),N)), prec2), -r2);
1007:
1008: disable_dbg(0);
1009:
1010: h = lg(listCR) - 1;
1011: dataCR = cgetg(h + 1, t_VEC);
1012: for (i = 1; i <= h ;i++)
1013: dataCR[i] = lgetg(10, t_VEC);
1014:
1015: q = gnorml2((GEN)modul[2]);
1016: p1 = cgetg(5, t_VEC);
1017: p1[1] = (long)q;
1018: p1[2] = lsubsi(r1, q);
1019: p1[3] = (long)gr2;
1020: p1[4] = lmax(addis((GEN)p1[2], r2+1), addsi(r2, q));
1021:
1022: for (i = 1; i <= h; i++)
1023: {
1024: GEN olddata, data = (GEN)dataCR[i];
1025:
1026: chi = gmael(listCR, i, 1);
1027: cond = gmael(listCR, i, 2);
1028:
1029: /* do we already know about the invariants of chi? */
1030: olddata = NULL;
1031: for (j = 1; j < i; j++)
1032: if (gegal(cond, gmael(listCR, j, 2)))
1033: { olddata = (GEN)dataCR[j]; break; }
1034:
1035: /* if cond(chi) = cond(bnr) */
1036: if (!olddata && gegal(cond, modul))
1037: {
1038: data[2] = lmul(C, gsqrt(det((GEN)cond[1]), prec2));
1039: data[3] = (long)bnr;
1040: data[6] = lgetg(1, t_VEC);
1041: data[7] = modul[1];
1042: data[9] = (long)p1;
1043:
1044: olddata = data;
1045: }
1046:
1047: data[1] = (long)chi; /* the character */
1048: if (!olddata)
1049: {
1050: data[2] = lmul(C, gsqrt(det((GEN)cond[1]), prec2));
1051: data[3] = (long)buchrayinitgen(bnf, cond);
1052: }
1053: else
1054: {
1055: data[2] = olddata[2]; /* constant C(chi) */
1056: data[3] = olddata[3]; /* bnr(cond(chi)) */
1057: }
1058: data[4] = (long)bnr; /* bnr(m) */
1059:
1060: chic = cgetg(nbg + 1, t_VEC);
1061: for (j = 1; j <= nbg; j++)
1062: chic[j] = ldiv((GEN)chi[j], (GEN)Mr[j]);
1063: data[5] = (long)get_Char(chic,prec); /* char associated to bnr(m) */
1064:
1065: /* compute diff(chi) and the corresponding primitive character */
1066: data[7] = cond[1];
1067: p2 = GetPrimChar(chi, bnr, (GEN)data[3], prec2);
1068: if (p2)
1069: {
1070: data[6] = p2[2];
1071: data[8] = p2[1];
1072: }
1073: else
1074: {
1075: data[6] = lgetg(1, t_VEC);
1076: data[8] = data[5];
1077: }
1078:
1079: /* compute q and store [q, r1 - q, r2] */
1080: if (!olddata)
1081: {
1082: q = gnorml2((GEN)cond[2]);
1083: p2 = cgetg(5, t_VEC);
1084: p2[1] = (long)q;
1085: p2[2] = lsubsi(r1, q);
1086: p2[3] = (long)gr2;
1087: p2[4] = lmax(addis((GEN)p2[2], r2+1), addsi(r2, q));
1088: data[9] = (long)p2;
1089: }
1090: else
1091: data[9] = olddata[9];
1092: }
1093:
1094: disable_dbg(-1);
1095: return gerepilecopy(av, dataCR);
1096: }
1097:
1098: /* compute the list of characters to consider for AllStark and
1099: initialize the data to compute with them */
1100: static GEN
1101: InitChar0(GEN dataD, long prec)
1102: {
1103: GEN MrD, listCR, p1, chi, lchi, Surj, cond, bnr, p2, Mr, d, allCR;
1104: long hD, h, nc, i, j, lD, nbg, tnc, av = avma;
1105:
1106: Surj = gmael(dataD, 2, 3);
1107: MrD = gmael(dataD, 2, 2);
1108: bnr = (GEN)dataD[1];
1109: Mr = gmael(bnr, 5, 2);
1110: hD = itos(gmael(dataD, 2, 1));
1111: h = hD >> 1;
1112: lD = lg(MrD)-1;
1113: nbg = lg(Mr) - 1;
1114:
1115: disable_dbg(0);
1116:
1117: listCR = cgetg(h + 1, t_VEC); /* non-conjugate characters */
1118: nc = 1;
1119: allCR = cgetg(h + 1, t_VEC); /* all characters, including conjugates */
1120: tnc = 1;
1121:
1122: p1 = FindEltofGroup(hD, MrD);
1123:
1124: for (i = 1; tnc <= h; i++)
1125: {
1126: /* lift a character of D in Clk(m) */
1127: chi = (GEN)p1[i];
1128: for (j = 1; j <= lD; j++) chi[j] = ldiv((GEN)chi[j], (GEN)MrD[j]);
1129: lchi = LiftChar(Mr, Surj, chi);
1130:
1131: for (j = 1; j < tnc; j++)
1132: if (gegal(lchi, (GEN)allCR[j])) break;
1133: if (j != tnc) continue;
1134:
1135: cond = bnrconductorofchar(bnr, lchi);
1136: if (gcmp0((GEN)cond[2])) continue;
1137:
1138: /* the infinite part of chi is non trivial */
1139: p2 = cgetg(3, t_VEC);
1140: p2[1] = (long)lchi;
1141: p2[2] = (long)cond;
1142: listCR[nc++] = (long)p2;
1143: allCR[tnc++] = (long)lchi;
1144:
1145: /* compute the order of chi */
1146: p2 = cgetg(nbg + 1, t_VEC);
1147: for (j = 1; j <= nbg; j++)
1148: p2[j] = ldiv((GEN)lchi[j], (GEN)Mr[j]);
1149: d = denom(p2);
1150:
1151: /* if chi is not real, add its conjugate character to allCR */
1152: if (!gegal(d, gdeux))
1153: allCR[tnc++] = (long)ConjChar(lchi, Mr);
1154: }
1155:
1156: setlg(listCR, nc);
1157: disable_dbg(-1);
1158: return gerepileupto(av, InitChar(bnr, listCR, prec));
1159: }
1160:
1161: /* recompute dataCR with the new precision */
1162: static GEN
1163: CharNewPrec(GEN dataCR, GEN nf, long prec)
1164: {
1165: GEN dk, C, p1, Pi;
1166: long av = avma, N, l, j, prec2;
1167:
1168: dk = (GEN)nf[3];
1169: N = degpol((GEN)nf[1]);
1170: l = lg(dataCR) - 1;
1171: prec2 = ((prec - 2)<<1) + EXTRA_PREC;
1172:
1173: Pi = mppi(prec2);
1174:
1175: C = gsqrt(gdiv(dk, gpowgs(Pi, N)), prec2);
1176:
1177: for (j = 1; j <= l; j++)
1178: {
1179: mael(dataCR, j, 2) = lmul(C, gsqrt(det(gmael(dataCR, j, 7)), prec2));
1180:
1181: mael4(dataCR, j, 3, 1, 7) = lcopy(nf);
1182:
1183: p1 = gmael(dataCR, j, 5);
1184: p1[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), (GEN)p1[3]), prec);
1185:
1186: p1 = gmael(dataCR, j, 8);
1187: p1[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), (GEN)p1[3]), prec);
1188: }
1189:
1190: return gerepilecopy(av, dataCR);
1191: }
1192:
1193: /********************************************************************/
1194: /* 4th part: compute the coefficients an(chi) */
1195: /* */
1196: /* matan entries are arrays of ints containing the coefficients of */
1197: /* an(chi) as a polmod modulo cyclo(order(chi)) */
1198: /********************************************************************/
1199:
1200: static void
1201: _0toCoeff(int *rep, long dg)
1202: {
1203: long i;
1204: for (i=0; i<dg; i++) rep[i] = 0;
1205: }
1206:
1207: /* transform a polmod into coeff */
1208: static void
1209: Polmod2Coeff(int *rep, GEN polmod, long dg)
1210: {
1211: GEN pol = (GEN)polmod[2];
1212: long i,d = degpol(pol);
1213:
1214: pol += 2;
1215: for (i=0; i<=d; i++) rep[i] = itos((GEN)pol[i]);
1216: for ( ; i<dg; i++) rep[i] = 0;
1217: }
1218:
1219: /* initialize a cl x nmax x [degs[1], ..., degs[cl] matrix of ints */
1220: /* modified to allocate ints and pointers separately since this used to
1221: break on 64bit platforms --GN1999Sep01 */
1222: static int***
1223: InitMatAn(long cl, long nmax, GEN degs, long flag)
1224: {
1225: long si,dg,i,j,k, n = nmax+1;
1226: int *c, **pj, ***A;
1227:
1228: for (si=0, i=1; i <= cl; i++)
1229: {
1230: dg = degs[i];
1231: si += dg;
1232: }
1233: A = (int***)gpmalloc((cl+1)*sizeof(int**) + cl*n*sizeof(int*));
1234: c = (int*)gpmalloc(si*n*sizeof(int));
1235: A[0] = (int**)c; /* keep it around for FreeMat() */
1236:
1237: pj = (int**)(A + (cl+1));
1238: for (i = 1; i <= cl; i++, pj+=n)
1239: {
1240: A[i] = pj; dg = degs[i];
1241: for (j = 0; j < n; j++,c+=dg)
1242: {
1243: pj[j] = c;
1244: c[0] = (j == 1 || flag);
1245: for (k = 1; k < dg; k++) c[k] = 0;
1246: }
1247: }
1248: return A;
1249: }
1250:
1251: static void
1252: FreeMat(int ***m)
1253: {
1254: free((void *)(m[0]));
1255: free((void *)m);
1256: }
1257:
1258: /* initialize coeff reduction */
1259: /* similar change --GN1999Sep01 */
1260: static int***
1261: InitReduction(GEN dataCR, GEN degs)
1262: {
1263: long av = avma,si,sp,dg,i,j, cl = lg(dataCR)-1;
1264: int *c, **pj, ***A;
1265: GEN d,polmod,pol, x = polx[0];
1266:
1267: for (si=sp=0, i=1; i <= cl; i++)
1268: {
1269: dg = degs[i];
1270: sp += dg;
1271: si += dg*dg;
1272: }
1273: A = (int***)gpmalloc((cl+1)*sizeof(int**) + sp*sizeof(int*));
1274: c = (int*)gpmalloc(si*sizeof(int));
1275: A[0] = (int**)c; /* keep it around for FreeMat() */
1276:
1277: pj = (int**)(A + (cl+1));
1278: for (i = 1; i <= cl; i++)
1279: {
1280: A[i] = pj;
1281: d = gmael3(dataCR, i, 5, 3);
1282: pol = cyclo(itos(d), 0);
1283: dg = degs[i]; /* degree(pol) */
1284: for (j = 0; j < dg; j++,c+=dg)
1285: {
1286: pj[j] = c;
1287: polmod = gmodulcp(gpowgs(x, dg + j), pol);
1288: Polmod2Coeff(c, polmod, dg);
1289: }
1290: pj += dg;
1291: }
1292: avma = av; return A;
1293: }
1294:
1295: #if 0
1296: static void
1297: pan(int ***an,long cl, long nmax, GEN degs)
1298: {
1299: long i,j,k;
1300: for (i = 1; i <= cl; i++)
1301: {
1302: long dg = degs[i];
1303: for (j = 0; j <= nmax; j++)
1304: for (k = 0; k < dg; k++) fprintferr("%d ",an[i][j][k]);
1305: }
1306: }
1307: #endif
1308:
1309: /* multiply (with reduction) a polmod with a coeff. put result in c1 */
1310: static void
1311: MulPolmodCoeff(GEN polmod, int* c1, int** reduc, long dg)
1312: {
1313: long av,i,j;
1314: int c, *c2, *c3;
1315:
1316: if (gcmp1(polmod)) return;
1317: for (i = 0; i < dg; i++)
1318: if (c1[i]) break;
1319: if (i == dg) return;
1320: av = avma;
1321: c3 = (int*)new_chunk(2*dg);
1322: c2 = (int*)new_chunk(dg);
1323: Polmod2Coeff(c2,polmod, dg);
1324:
1325: for (i = 0; i < 2*dg; i++)
1326: {
1327: c = 0;
1328: for (j = 0; j <= i; j++)
1329: if (j < dg && j > i - dg) c += c1[j] * c2[i-j];
1330: c3[i] = c;
1331: }
1332: c2 = c1;
1333: for (i = 0; i < dg; i++)
1334: {
1335: c = c3[i];
1336: for (j = 0; j < dg; j++) c += reduc[j][i] * c3[dg+j];
1337: c2[i] = c;
1338: }
1339: /* cast necessary to work around a gcc-2.96 bug on alpha-linux (IS) */
1340: for ( ; i < (short)dg; i++) c2[i] = 0;
1341: avma = av;
1342: }
1343:
1344: /* c0 <- c0 + c2 * c1 */
1345: static void
1346: AddMulCoeff(int *c0, int *c2, int* c1, int** reduc, long dg)
1347: {
1348: long av,i,j;
1349: int c, *c3;
1350:
1351: if (!c2) /* c2 == 1 */
1352: {
1353: for (i = 0; i < dg; i++) c0[i] += c1[i];
1354: return;
1355: }
1356: for (i = 0; i <= dg; i++)
1357: if (c1[i]) break;
1358: if (i > dg) return;
1359: av = avma;
1360: c3 = (int*)new_chunk(2*dg);
1361: for (i = 0; i < 2*dg; i++)
1362: {
1363: c = 0;
1364: for (j = 0; j <= i; j++)
1365: if (j < dg && j > i - dg) c += c1[j] * c2[i-j];
1366: c3[i] = c;
1367: }
1368: for (i = 0; i < dg; i++)
1369: {
1370: c = c0[i] + c3[i];
1371: for (j = 0; j < dg; j++) c += reduc[j][i] * c3[dg+j];
1372: c0[i] = c;
1373: }
1374:
1375: avma = av;
1376: }
1377:
1378: /* returns 0 if c is zero, 1 otherwise. */
1379: static long
1380: IsZero(int* c, long dg)
1381: {
1382: long i;
1383:
1384: for (i = 0; i < dg; i++)
1385: if (c[i]) return 0;
1386: return 1;
1387: }
1388:
1389: /* evaluate the coeff. No Garbage collector */
1390: static GEN
1391: EvalCoeff(GEN z, int* c, long dg)
1392: {
1393: long i,j;
1394: GEN e, r;
1395:
1396: #if 0
1397: /* standard Horner */
1398: e = stoi(c[dg - 1]);
1399: for (i = dg - 2; i >= 0; i--)
1400: e = gadd(stoi(c[i]), gmul(z, e));
1401: #else
1402: /* specific attention to sparse polynomials */
1403: e = NULL;
1404: for (i = dg-1; i >=0; i=j-1)
1405: {
1406: for (j=i; c[j] == 0; j--)
1407: if (j==0)
1408: {
1409: if (!e) return NULL;
1410: if (i!=j) z = gpuigs(z,i-j+1);
1411: return gmul(e,z);
1412: }
1413: if (e)
1414: {
1415: r = (i==j)? z: gpuigs(z,i-j+1);
1416: e = gadd(gmul(e,r), stoi(c[j]));
1417: }
1418: else
1419: e = stoi(c[j]);
1420: }
1421: #endif
1422: return e;
1423: }
1424:
1425: /* copy the n * m array matan */
1426: static void
1427: CopyCoeff(int*** a, int*** a2, long n, long m, GEN degs)
1428: {
1429: long i,j,k;
1430:
1431: for (i = 1; i <= n; i++)
1432: {
1433: long dg = degs[i];
1434: int **b = a[i], **b2 = a2[i];
1435: for (j = 0; j <= m; j++)
1436: {
1437: int *c = b[j], *c2 = b2[j];
1438: for (k = 0; k < dg; k++) c2[k] = c[k];
1439: }
1440: }
1441: return;
1442: }
1443:
1444: /* initialize the data for GetRay */
1445: static GEN
1446: InitGetRay(GEN bnr, long nmax)
1447: {
1448: long bd, i, j, l;
1449: GEN listid, listcl, id, rep, bnf, cond;
1450:
1451: bnf = (GEN)bnr[1];
1452: cond = gmael3(bnr, 2, 1, 1);
1453:
1454: if (nmax < 1000) return NULL;
1455:
1456: rep = cgetg(4, t_VEC);
1457:
1458: disable_dbg(0);
1459: bd = min(1000, nmax / 50);
1460: listid = ideallist(bnf, bd);
1461: disable_dbg(-1);
1462:
1463: listcl = cgetg(bd + 1, t_VEC);
1464: for (i = 1; i <= bd; i++)
1465: {
1466: l = lg((GEN)listid[i]) - 1;
1467: listcl[i] = lgetg(l + 1, t_VEC);
1468:
1469: for (j = 1; j <= l; j++)
1470: {
1471: id = gmael(listid, i, j);
1472: if (gcmp1(gcoeff(idealadd(bnf, id, cond), 1, 1)))
1473: mael(listcl, i, j) = (long)isprincipalray(bnr, id);
1474: }
1475: }
1476:
1477: if (DEBUGLEVEL) msgtimer("InitGetRay");
1478:
1479: rep[1] = (long)listid;
1480: rep[2] = (long)listcl;
1481: rep[3] = nf_get_r2((GEN)bnf[7])? 0: un; /* != 0 iff nf is totally real */
1482: return rep;
1483: }
1484:
1485: /* compute the class of the prime ideal pr in cl(bnr) using dataray */
1486: static GEN
1487: GetRay(GEN bnr, GEN dataray, GEN pr, long prec)
1488: {
1489: long av = avma, N, n, bd, c;
1490: GEN id, tid, t2, u, alpha, p1, cl, listid, listcl, nf, cond;
1491:
1492: if (!dataray)
1493: return isprincipalray(bnr, pr);
1494:
1495: listid = (GEN)dataray[1];
1496: listcl = (GEN)dataray[2];
1497: cond = gmael3(bnr, 2, 1, 1);
1498: bd = lg(listid) - 1;
1499: nf = gmael(bnr, 1, 7);
1500: N = degpol(nf[1]);
1501:
1502: if (dataray[3])
1503: t2 = gmael(nf, 5, 4);
1504: else
1505: t2 = gmael(nf, 5, 3);
1506:
1507: id = prime_to_ideal(nf, pr);
1508: tid = qf_base_change(t2, id, 1);
1509:
1510: if (dataray[3])
1511: u = lllgramint(tid);
1512: else
1513: u = lllgramintern(tid,100,1,prec);
1514:
1515: if (!u) return gerepileupto(av, isprincipalray(bnr, id));
1516:
1517: c = 1; alpha = NULL;
1518: for (c=1; c<=N; c++)
1519: {
1520: p1 = gmul(id, (GEN)u[c]);
1521: if (gcmp1(gcoeff(idealadd(nf, p1, cond), 1, 1))) { alpha = p1; break; }
1522: }
1523: if (!alpha)
1524: return gerepileupto(av, isprincipalray(bnr, pr));
1525:
1526: id = idealdivexact(nf, alpha, id);
1527:
1528: n = itos(det(id));
1529: if (n > bd)
1530: cl = isprincipalray(bnr, id);
1531: else
1532: {
1533: cl = NULL;
1534: c = 1;
1535: p1 = (GEN)listid[n];
1536: while (!cl)
1537: {
1538: if (gegal((GEN)p1[c], id))
1539: cl = gmael(listcl, n, c);
1540: c++;
1541: }
1542: }
1543:
1544: return gerepileupto(av, gsub(isprincipalray(bnr, alpha), cl));
1545: }
1546:
1547: /* correct the coefficients an(chi) according with diff(chi) in place */
1548: static void
1549: CorrectCoeff(GEN dtcr, int** an, int** reduc, long nmax, long dg)
1550: {
1551: long lg, av1, j, p, q, limk, k, l, av = avma;
1552: int ***an2, **an1, *c, *c2;
1553: GEN chi, bnrc, diff, ray, ki, ki2, pr, degs;
1554:
1555: chi = (GEN)dtcr[8];
1556: bnrc = (GEN)dtcr[3];
1557: diff = (GEN)dtcr[6];
1558: lg = lg(diff) - 1;
1559: if (!lg) return;
1560:
1561: if (DEBUGLEVEL > 2) fprintferr("diff(chi) = %Z", diff);
1562:
1563: degs = cgetg(2, t_VECSMALL); degs[1] = dg;
1564: an2 = InitMatAn(1, nmax, degs, 0); an1 = an2[1];
1565: c = (int*)new_chunk(dg);
1566: av1 = avma;
1567:
1568: for (j = 1; j <= lg; j++)
1569: {
1570: for (k = 0; k <= nmax; k++)
1571: for (l = 0; l < dg; l++) an1[k][l] = an[k][l];
1572:
1573: pr = (GEN)diff[j];
1574: ray = isprincipalray(bnrc, pr);
1575: ki = ComputeImagebyChar(chi, ray, 1);
1576: ki2 = gcopy(ki);
1577:
1578: q = p = itos(powgi((GEN)pr[1], (GEN)pr[4]));
1579: limk = nmax / q;
1580:
1581: while (q <= nmax)
1582: {
1583: if (gcmp1(ki2)) c2 = NULL; else { Polmod2Coeff(c,ki2, dg); c2 = c; }
1584: for(k = 1; k <= limk; k++)
1585: AddMulCoeff(an[k*q], c2, an1[k], reduc, dg);
1586:
1587: q *= p; limk /= p;
1588: ki2 = gmul(ki2, ki);
1589: }
1590: avma = av1;
1591: }
1592: FreeMat(an2); avma = av;
1593: }
1594:
1595: /* compute the coefficients an in the general case */
1596: static int***
1597: ComputeCoeff(GEN dataCR, long nmax, long prec)
1598: {
1599: long cl, i, j, av = avma, av2, np, q, q1, limk, k, id, cpt = 10, dg, Bq;
1600: int ***matan, ***reduc, ***matan2, *c2;
1601: GEN c, degs, tabprem, bnf, pr, cond, ray, ki, ki2, prime, npg, bnr, dataray;
1602: byteptr dp = diffptr;
1603:
1604: bnr = gmael(dataCR, 1, 4);
1605: bnf = (GEN)bnr[1];
1606: cond = gmael3(bnr, 2, 1, 1);
1607: cl = lg(dataCR) - 1;
1608:
1609: dataray = InitGetRay(bnr, nmax);
1610:
1611: degs = GetDeg(dataCR);
1612: matan = InitMatAn(cl, nmax, degs, 0);
1613: matan2 = InitMatAn(cl, nmax, degs, 0);
1614: reduc = InitReduction(dataCR, degs);
1615: c = cgetg(cl + 1, t_VEC);
1616: for (i = 1; i <= cl; i++)
1617: c[i] = (long)new_chunk(degs[i]);
1618:
1619: if (DEBUGLEVEL > 1) fprintferr("p = ");
1620:
1621: prime = stoi(2); dp++;
1622: av2 = avma;
1623: while (*dp && (prime[2] <= nmax))
1624: {
1625: tabprem = primedec(bnf, prime);
1626: for (j = 1; j < lg(tabprem); j++)
1627: {
1628: pr = (GEN)tabprem[j];
1629: npg = powgi((GEN)pr[1], (GEN)pr[4]);
1630: if (is_bigint(npg) || (np=npg[2]) > nmax
1631: || idealval(bnf, cond, pr)) continue;
1632:
1633: CopyCoeff(matan, matan2, cl, nmax, degs);
1634: ray = GetRay(bnr, dataray, pr, prec);
1635: ki = chiideal(dataCR, ray, 1);
1636: ki2 = dummycopy(ki);
1637:
1638: Bq = nmax/np;
1639: for (q1 = 1; q1 <= Bq; q1 *= np)
1640: {
1641: q = q1*np;
1642: limk = nmax / q;
1643: for (id = 1; id <= cl; id++)
1644: {
1645: dg = degs[id];
1646: if (gcmp1((GEN)ki2[id]))
1647: c2 = NULL;
1648: else
1649: {
1650: c2 = (int*)c[id];
1651: Polmod2Coeff(c2, (GEN)ki2[id], dg);
1652: }
1653: for (k = 1; k <= limk; k++)
1654: AddMulCoeff(matan[id][k*q], c2, matan2[id][k], reduc[id], dg);
1655: ki2[id] = lmul((GEN)ki2[id], (GEN)ki[id]);
1656: }
1657: }
1658: }
1659: avma = av2;
1660: prime[2] += (*dp++);
1661: if (!*dp) err(primer1);
1662:
1663: if (DEBUGLEVEL > 1 && prime[2] > cpt)
1664: { fprintferr("%ld ", prime[2]); cpt += 10; }
1665: }
1666: if (DEBUGLEVEL > 1) fprintferr("\n");
1667:
1668: for (i = 1; i <= cl; i++)
1669: CorrectCoeff((GEN)dataCR[i], matan[i], reduc[i], nmax, degs[i]);
1670:
1671: FreeMat(matan2); FreeMat(reduc);
1672: avma = av; return matan;
1673: }
1674:
1675: /********************************************************************/
1676: /* 5th part: compute L-functions at s=1 */
1677: /********************************************************************/
1678:
1679: /* if flag != 0, prec means decimal digits */
1680: static GEN
1681: get_limx(long r1, long r2, long prec, GEN *pteps, long flag)
1682: {
1683: GEN eps, a, r, c0, A0, limx, Pi = mppi(prec), N, p1;
1684:
1685: N = addss(r1, 2*r2);
1686: a = gmul(gpow(gdeux, gsubgs(gdiv(stoi(r1), N), 1), DEFAULTPREC), N);
1687: r = addss(r1, r2);
1688:
1689: if (flag)
1690: *pteps = eps = gmul2n(gpowgs(dbltor(10.), -prec), -1);
1691: else
1692: *pteps = eps = gmul2n(gpowgs(dbltor(10.), (long)(-(prec-2) / pariK1)), -1);
1693:
1694: c0 = gpow(gmul2n(Pi, 1), gdiv(subis(r, 1), gdeux), DEFAULTPREC);
1695: c0 = gmul(c0, gdiv(gdeux, N));
1696: c0 = gmul(c0, gpow(gdeux, gmul(gdiv(stoi(r1), gdeux),
1697: gsubsg(1, gdiv(addis(r, 1), N))),
1698: DEFAULTPREC));
1699:
1700: A0 = glog(gdiv(gmul2n(c0, 1), eps), DEFAULTPREC);
1701:
1702: limx = gpow(gdiv(a, A0), gdiv(N, gdeux), DEFAULTPREC);
1703: p1 = gsub(glog(A0, DEFAULTPREC), glog(a, DEFAULTPREC));
1704: p1 = gmul(gmul(p1, N), addis(r, 1));
1705: p1 = gdiv(p1, gmul2n(gadd(gmul2n(A0, 1), addis(r, 1)), 1));
1706: limx = gmul(limx, gaddgs(p1, 1));
1707:
1708: return limx;
1709: }
1710:
1711: static long
1712: GetBoundN0(GEN C, long r1, long r2, long prec, long flag)
1713: {
1714: long av = avma;
1715: GEN eps, limx = get_limx(r1, r2, prec, &eps, flag);
1716:
1717: limx = gfloor(gdiv(C, limx));
1718: if (is_bigint(limx))
1719: err(talker, "Too many coefficients (%Z) needed in GetST: computation impossible", limx);
1720:
1721: avma = av; return itos(limx);
1722: }
1723:
1724: static long
1725: GetBoundi0(long r1, long r2, long prec)
1726: {
1727: long av = avma, imin, i0, itest;
1728: GEN ftest, borneps, eps, limx = get_limx(r1, r2, prec, &eps, 0);
1729: GEN Pi = mppi(DEFAULTPREC);
1730:
1731: borneps = gmul(gmul2n(gun, r2), gpow(Pi, gdiv(subss(r2, 3), gdeux),
1732: DEFAULTPREC));
1733: borneps = gdiv(gmul(borneps, gpowgs(stoi(5), r1)), eps);
1734: borneps = gdiv(borneps, gsqrt(limx, DEFAULTPREC));
1735:
1736: imin = 1;
1737: i0 = 1400;
1738: while(i0 - imin >= 4)
1739: {
1740: itest = (i0 + imin) >> 1;
1741:
1742: ftest = gpowgs(limx, itest);
1743: ftest = gmul(ftest, gpowgs(mpfactr(itest / 2, DEFAULTPREC), r1));
1744: ftest = gmul(ftest, gpowgs(mpfactr(itest, DEFAULTPREC), r2));
1745:
1746: if(gcmp(ftest, borneps) >= 0)
1747: i0 = itest;
1748: else
1749: imin = itest;
1750: }
1751: avma = av;
1752:
1753: return (i0 / 2) * 2;
1754: }
1755:
1756: /* compute the principal part at the integers s = 0, -1, -2, ..., -i0
1757: of Gamma((s+1)/2)^a Gamma(s/2)^b Gamma(s)^c / (s - z) with z = 0 and 1 */
1758: /* NOTE: this is surely not the best way to do this, but it's fast enough! */
1759: static GEN
1760: ppgamma(long a, long b, long c, long i0, long prec)
1761: {
1762: GEN cst, gamun, gamdm, an, bn, cn_evn, cn_odd, x, x2, aij, p1, cf, p2;
1763: long i, j, r, av = avma;
1764:
1765: r = max(b + c + 1, a + c);
1766:
1767: aij = cgetg(i0 + 1, t_VEC);
1768: for (i = 1; i <= i0; i++)
1769: {
1770: aij[i] = lgetg(3, t_VEC);
1771: mael(aij, i, 1) = lgetg(r + 1, t_VEC);
1772: mael(aij, i, 2) = lgetg(r + 1, t_VEC);
1773: }
1774:
1775: x = polx[0];
1776: x2 = gmul2n(x, -1);
1777:
1778: /* Euler gamma constant, values of Riemann zeta functions at
1779: positive integers */
1780: cst = cgetg(r + 2, t_VEC);
1781: cst[1] = (long)mpeuler(prec);
1782: for (i = 2; i <= r + 1; i++)
1783: cst[i] = (long)gzeta(stoi(i), prec);
1784:
1785: /* the expansion of log(Gamma(s)) at s = 1 */
1786: gamun = cgetg(r + 2, t_SER);
1787: gamun[1] = evalsigne(1) | evalvalp(0) | evalvarn(0);
1788: gamun[2] = zero;
1789: for (i = 1; i <= r; i++)
1790: {
1791: gamun[i + 2] = ldivgs((GEN)cst[i], i);
1792: if (i%2) gamun[i + 2] = lneg((GEN)gamun[i + 2]);
1793: }
1794:
1795: /* the expansion of log(Gamma(s)) at s = 1/2 */
1796: gamdm = cgetg(r + 2, t_SER);
1797: gamdm[1] = evalsigne(1) | evalvalp(0) | evalvarn(0);
1798: gamdm[2] = (long)mplog(gsqrt(mppi(prec), prec));
1799: gamdm[3] = lneg(gadd(gmul2n(glog(gdeux, prec), 1), (GEN)cst[1]));
1800: for (i = 2; i <= r; i++)
1801: gamdm[i + 2] = lmul((GEN)gamun[i + 2], subis(shifti(gun, i), 1));
1802:
1803: gamun = gexp(gamun, prec);
1804: gamdm = gexp(gamdm, prec);
1805:
1806: /* We simplify to get one of the following two expressions */
1807:
1808: /* Case 1 (b > a): sqrt{Pi}^a 2^{a - as} Gamma(s/2)^{b-a} Gamma(s)^{a + c} */
1809: if (b > a)
1810: {
1811: cf = gpui(mppi(prec), gmul2n(stoi(a), -1), prec);
1812:
1813: /* an is the expansion of Gamma(x)^{a+c} */
1814: an = gpowgs(gdiv(gamun, x), a + c);
1815:
1816: /* bn is the expansion of 2^{a - ax} */
1817: bn = gpowgs(gpow(gdeux, gsubsg(1, x), prec), a);
1818:
1819: /* cn_evn is the expansion of Gamma(x/2)^{b-a} */
1820: cn_evn = gpowgs(gdiv(gsubst(gamun, 0, x2), x2), b - a);
1821:
1822: /* cn_odd is the expansion of Gamma((x-1)/2)^{b-a} */
1823: cn_odd = gpowgs(gdiv(gsubst(gamdm, 0, x2), gsub(x2, ghalf)), b - a);
1824:
1825: for (i = 0; i < i0/2; i++)
1826: {
1827: p1 = gmul(cf, gmul(an, gmul(bn, cn_evn)));
1828:
1829: p2 = gdiv(p1, gsubgs(x, 2*i));
1830: for (j = 1; j <= r; j++)
1831: mael3(aij, 2*i + 1, 1, j) = (long)polcoeff0(p2, -j, 0);
1832:
1833: p2 = gdiv(p1, gsubgs(x, 2*i + 1));
1834: for (j = 1; j <= r; j++)
1835: mael3(aij, 2*i + 1, 2, j) = (long)polcoeff0(p2, -j, 0);
1836:
1837: /* an(x-s-1) = an(x-s) / (x-s-1)^{a+c} */
1838: an = gdiv(an, gpowgs(gsubgs(x, 2*i + 1), a + c));
1839:
1840: /* bn(x-s-1) = 2^a bn(x-s) */
1841: bn = gmul2n(bn, a);
1842:
1843: /* cn_evn(x-s-2) = cn_evn(x-s) / (x/2 - (s+2)/2)^{b-a} */
1844: cn_evn = gdiv(cn_evn, gpowgs(gsubgs(x2, i + 1), b - a));
1845:
1846: p1 = gmul(cf, gmul(an, gmul(bn, cn_odd)));
1847:
1848: p2 = gdiv(p1, gsubgs(x, 2*i + 1));
1849: for (j = 1; j <= r; j++)
1850: mael3(aij, 2*i + 2, 1, j) = (long)polcoeff0(p2, -j, 0);
1851:
1852: p2 = gdiv(p1, gsubgs(x, 2*i + 2));
1853: for (j = 1; j <= r; j++)
1854: mael3(aij, 2*i + 2, 2, j) = (long)polcoeff0(p2, -j, 0);
1855:
1856: an = gdiv(an, gpowgs(gsubgs(x, 2*i + 2), a + c));
1857: bn = gmul2n(bn, a);
1858:
1859: /* cn_odd(x-s-2) = cn_odd(x-s) / (x/2 - (s+2)/2)^{b-a} */
1860: cn_odd = gdiv(cn_odd, gpowgs(gsub(x2, gaddgs(ghalf, i + 1)), b - a));
1861: }
1862: }
1863: else
1864: /* Case 2 (b <= a): sqrt{Pi}^b 2^{b - bs} Gamma((s+1)/2)^{a-b}
1865: Gamma(s)^{b + c) */
1866: {
1867: cf = gpui(mppi(prec), gmul2n(stoi(b), -1), prec);
1868:
1869: /* an is the expansion of Gamma(x)^{b+c} */
1870: an = gpowgs(gdiv(gamun, x), b + c);
1871:
1872: /* bn is the expansion of 2^{b - bx} */
1873: bn = gpowgs(gpow(gdeux, gsubsg(1, x), prec), b);
1874:
1875: /* cn_evn is the expansion of Gamma((x+1)/2)^{a-b} */
1876: cn_evn = gpowgs(gsubst(gamdm, 0, x2), a - b);
1877:
1878: /* cn_odd is the expansion of Gamma(x/2)^{a-b} */
1879: cn_odd = gpowgs(gdiv(gsubst(gamun, 0, x2), x2), a - b);
1880:
1881: for (i = 0; i < i0/2; i++)
1882: {
1883: p1 = gmul(cf, gmul(an, gmul(bn, cn_evn)));
1884:
1885: p2 = gdiv(p1, gsubgs(x, 2*i));
1886: for (j = 1; j <= r; j++)
1887: mael3(aij, 2*i + 1, 1, j) = (long)polcoeff0(p2, -j, 0);
1888:
1889: p2 = gdiv(p1, gsubgs(x, 2*i + 1));
1890: for (j = 1; j <= r; j++)
1891: mael3(aij, 2*i + 1, 2, j) = (long)polcoeff0(p2, -j, 0);
1892:
1893: /* an(x-s-1) = an(x-s) / (x-s-1)^{b+c} */
1894: an = gdiv(an, gpowgs(gsubgs(x, 2*i + 1), b + c));
1895:
1896: /* bn(x-s-1) = 2^b bn(x-s) */
1897: bn = gmul2n(bn, b);
1898:
1899: /* cn_evn(x-s-2) = cn_evn(x-s) / (x/2 - (s+1)/2)^{a-b} */
1900: cn_evn = gdiv(cn_evn, gpowgs(gsub(x2, gaddgs(ghalf, i)), a - b));
1901:
1902: p1 = gmul(cf, gmul(an, gmul(bn, cn_odd)));
1903:
1904: p2 = gdiv(p1, gsubgs(x, 2*i + 1));
1905: for (j = 1; j <= r; j++)
1906: mael3(aij, 2*i + 2, 1, j) = (long)polcoeff0(p2, -j, 0);
1907:
1908: p2 = gdiv(p1, gsubgs(x, 2*i + 2));
1909: for (j = 1; j <= r; j++)
1910: mael3(aij, 2*i + 2, 2, j) = (long)polcoeff0(p2, -j, 0);
1911:
1912: an = gdiv(an, gpowgs(gsubgs(x, 2*i + 2), b + c));
1913: bn = gmul2n(bn, b);
1914:
1915: /* cn_odd(x-s-2) = cn_odd(x-s) / (x/2 - (s+1)/2)^{a-b} */
1916: cn_odd = gdiv(cn_odd, gpowgs(gsubgs(x2, i + 1), a - b));
1917: }
1918: }
1919:
1920: return gerepilecopy(av, aij);
1921: }
1922:
1923: static GEN
1924: GetST(GEN dataCR, long prec)
1925: {
1926: GEN N0, CC, bnr, bnf, nf, Pi, racpi, C, cond, aij, B, S, T, csurn, lncsurn;
1927: GEN degs, p1, p2, nsurc, an, rep, powlncn, powracpi;
1928: long i, j, k, n, av = avma, av1, av2, hk, fj, id, prec2, i0, nmax;
1929: long a, b, c, rc1, rc2, r, r1, r2;
1930: int ***matan;
1931:
1932: if (DEBUGLEVEL) timer2();
1933: bnr = gmael(dataCR, 1, 4);
1934: bnf = checkbnf(bnr);
1935: nf = checknf(bnf);
1936: r1 = nf_get_r1(nf);
1937: r2 = nf_get_r2(nf);
1938: hk = lg(dataCR) - 1;
1939: prec2 = ((prec - 2)<<1) + EXTRA_PREC;
1940:
1941: Pi = mppi(prec2);
1942: racpi = gsqrt(Pi, prec2);
1943:
1944: C = cgetg(hk + 1, t_VEC);
1945: cond = cgetg(hk + 1, t_VEC);
1946: N0 = new_chunk(hk+1);
1947: CC = new_chunk(hk+1);
1948: nmax = 0;
1949: for (i = 1; i <= hk; i++)
1950: {
1951: C[i] = mael(dataCR, i, 2);
1952:
1953: p1 = cgetg(3, t_VEC);
1954: p1[1] = mael(dataCR, i, 7);
1955: p1[2] = mael(dataCR, i, 9);
1956: cond[i] = (long)p1;
1957:
1958: N0[i] = GetBoundN0((GEN)C[i], r1, r2, prec, 0);
1959: if (nmax < N0[i]) nmax = N0[i];
1960: }
1961:
1962: if ((ulong)nmax > maxprime())
1963: err(talker, "Not enough precomputed primes (need all primes up to %ld)", nmax);
1964:
1965: i0 = GetBoundi0(r1, r2, prec);
1966:
1967: if(DEBUGLEVEL > 1) fprintferr("nmax = %ld and i0 = %ld\n", nmax, i0);
1968:
1969: matan = ComputeCoeff(dataCR, nmax, prec);
1970: degs = GetDeg(dataCR);
1971: if (DEBUGLEVEL) msgtimer("Compute an");
1972:
1973: p1 = cgetg(3, t_COMPLEX);
1974: p1[1] = lgetr(prec2);
1975: p1[2] = lgetr(prec2);
1976: gaffect(gzero, p1);
1977:
1978: S = cgetg(hk + 1, t_VEC);
1979: T = cgetg(hk + 1, t_VEC);
1980: for (id = 1; id <= hk; id++)
1981: {
1982: S[id] = lcopy(p1);
1983: T[id] = lcopy(p1);
1984: for (k = 1; k < id; k++)
1985: if (gegal((GEN)cond[id], (GEN)cond[k])) break;
1986: CC[id] = k;
1987: }
1988:
1989: powracpi = cgetg(hk + 1, t_VEC);
1990: for (j = 1; j <= hk; j++)
1991: powracpi[j] = (long)gpow(racpi, gmael3(dataCR, j, 9, 2), prec2);
1992:
1993: av1 = avma;
1994: if (DEBUGLEVEL > 1) fprintferr("n = ");
1995:
1996: for (id = 1; id <= hk; id++)
1997: {
1998: if (CC[id] != id) continue;
1999: p2 = gmael(dataCR, id, 9);
2000: a = itos((GEN)p2[1]);
2001: b = itos((GEN)p2[2]);
2002: c = itos((GEN)p2[3]);
2003: aij = ppgamma(a, b, c, i0, prec2);
2004: rc1 = a + c;
2005: rc2 = b + c; r = max(rc2 + 1, rc1);
2006: av2 = avma;
2007:
2008: for (n = 1; n <= N0[id]; n++)
2009: {
2010: if (DEBUGLEVEL > 1 && n%100 == 0) fprintferr("%ld ", n);
2011:
2012: for (k = 1; k <= hk; k++)
2013: if (CC[k] == id && !IsZero(matan[k][n], degs[k])) break;
2014: if (k > hk) continue;
2015:
2016: csurn = gdivgs((GEN)C[id], n);
2017: nsurc = ginv(csurn);
2018:
2019: B = cgetg(r + 1, t_VEC);
2020: lncsurn = glog(csurn, prec2);
2021: powlncn = gun;
2022: fj = 1;
2023:
2024: p1 = gzero;
2025: p2 = gzero;
2026: for (j = 1; j <= r; j++)
2027: {
2028: if (j > 2) fj = fj * (j - 1);
2029:
2030: B[j] = ldivgs(powlncn, fj);
2031: p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, i0, 2, j)));
2032: p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, i0, 1, j)));
2033:
2034: powlncn = gmul(powlncn, lncsurn);
2035: }
2036: for (i = i0 - 1; i > 1; i--)
2037: {
2038: p1 = gmul(p1, nsurc);
2039: p2 = gmul(p2, nsurc);
2040: for (j = i%2? rc2: rc1; j; j--)
2041: {
2042: p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, i, 2, j)));
2043: p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, i, 1, j)));
2044: }
2045: }
2046: p1 = gmul(p1, nsurc);
2047: p2 = gmul(p2, nsurc);
2048: for (j = 1; j <= r; j++)
2049: {
2050: p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, 1, 2, j)));
2051: p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, 1, 1, j)));
2052: }
2053:
2054: p1 = gadd(p1, gmul(csurn, (GEN)powracpi[id]));
2055:
2056: for (j = 1; j <= hk; j++)
2057: if (CC[j] == id &&
2058: (an = EvalCoeff(gmael3(dataCR, j, 5, 2), matan[j][n], degs[j])))
2059: {
2060: gaffect(gadd((GEN)S[j], gmul(p1, an)), (GEN)S[j]);
2061: gaffect(gadd((GEN)T[j], gmul(p2, gconj(an))), (GEN)T[j]);
2062: }
2063: avma = av2;
2064: }
2065: avma = av1;
2066: }
2067: FreeMat(matan);
2068:
2069: if (DEBUGLEVEL > 1) fprintferr("\n");
2070: if (DEBUGLEVEL) msgtimer("Compute S&T");
2071:
2072: rep = cgetg(3, t_VEC);
2073: rep[1] = (long)S;
2074: rep[2] = (long)T;
2075: return gerepilecopy(av, rep);
2076: }
2077:
2078: /* Given datachi, S(chi) and T(chi), return L(1, chi) if fl = 1,
2079: or [r(chi), c(chi)] where r(chi) is the rank of chi and c(chi)
2080: is given by L(s, chi) = c(chi).s^r(chi) at s = 0 if fl = 0.
2081: if fl2 = 1, adjust the value to get L_S(s, chi). */
2082: static GEN
2083: GetValue(GEN datachi, GEN S, GEN T, long fl, long fl2, long prec)
2084: {
2085: GEN W, A, q, b, c, d, rchi, cf, VL, rep, racpi, nS, nT;
2086: long av = avma;
2087:
2088: racpi = gsqrt(mppi(prec), prec);
2089: W = ComputeArtinNumber(datachi, 0, prec);
2090: A = ComputeAChi(datachi, fl, prec);
2091:
2092: d = gmael(datachi, 8, 3);
2093:
2094: q = gmael(datachi, 9, 1);
2095: b = gmael(datachi, 9, 2);
2096: c = gmael(datachi, 9, 3);
2097:
2098: rchi = addii(b, c);
2099:
2100: if (!fl)
2101: {
2102: cf = gmul2n(gpow(racpi, q, 0), itos(b));
2103:
2104: nS = gdiv(gconj(S), cf);
2105: nT = gdiv(gconj(T), cf);
2106:
2107: /* VL = W(chi).S(conj(chi)) + T(chi)) / (sqrt(Pi)^q 2^{r1 - q}) */
2108: VL = gadd(gmul(W, nS), nT);
2109: if (cmpis(d, 3) < 0) VL = greal(VL);
2110:
2111: if (fl2)
2112: {
2113: VL = gmul((GEN)A[2], VL);
2114: rchi = gadd(rchi, (GEN)A[1]);
2115: }
2116:
2117: rep = cgetg(3, t_VEC);
2118: rep[1] = (long)rchi;
2119: rep[2] = (long)VL;
2120: }
2121: else
2122: {
2123: cf = gmul((GEN)datachi[2], gpow(racpi, b, 0));
2124:
2125: /* VL = S(chi) + W(chi).T(chi)) / (C(chi) sqrt(Pi)^{r1 - q}) */
2126: rep = gdiv(gadd(S, gmul(W, T)), cf);
2127: if (cmpis(d, 3) < 0) rep = greal(rep);
2128:
2129: if (fl2) rep = gmul(A, rep);
2130: }
2131:
2132: return gerepilecopy(av, rep);
2133: }
2134:
2135: /* return the order and the first non-zero term of L(s, chi0)
2136: at s = 0. If flag = 1, adjust the value to get L_S(s, chi0). */
2137: static GEN
2138: GetValue1(GEN bnr, long flag, long prec)
2139: {
2140: GEN bnf = checkbnf(bnr), nf = checknf(bnf);
2141: GEN hk, Rk, wk, c, rep, mod0, diff;
2142: long i, l, r, r1, r2, av = avma;
2143:
2144: r1 = nf_get_r1(nf);
2145: r2 = nf_get_r2(nf);
2146: hk = gmael3(bnf, 8, 1, 1);
2147: Rk = gmael(bnf, 8, 2);
2148: wk = gmael3(bnf, 8, 4, 1);
2149:
2150: c = gneg_i(gdiv(gmul(hk, Rk), wk));
2151: r = r1 + r2 - 1;
2152:
2153: if (flag)
2154: {
2155: mod0 = gmael3(bnr, 2, 1, 1);
2156: diff = (GEN)idealfactor(nf, mod0)[1];
2157:
2158: l = lg(diff) - 1; r += l;
2159: for (i = 1; i <= l; i++)
2160: c = gmul(c, glog(idealnorm(nf, (GEN)diff[i]), prec));
2161: }
2162:
2163: rep = cgetg(3, t_VEC);
2164: rep[1] = lstoi(r);
2165: rep[2] = (long)c;
2166:
2167: return gerepilecopy(av, rep);
2168: }
2169:
2170: /********************************************************************/
2171: /* 6th part: recover the coefficients */
2172: /********************************************************************/
2173:
2174: static long
2175: TestOne(GEN plg, GEN beta, GEN B, long v, long G, long N)
2176: {
2177: long j;
2178: GEN p1;
2179:
2180: p1 = gsub(beta, (GEN)plg[v]);
2181: if (expo(p1) >= G) return 0;
2182:
2183: for (j = 1; j <= N; j++)
2184: if (j != v)
2185: {
2186: p1 = gabs((GEN)plg[j], DEFAULTPREC);
2187: if (gcmp(p1, B) > 0) return 0;
2188: }
2189: return 1;
2190: }
2191:
2192: /* Using linear dependance relations */
2193: static GEN
2194: RecCoeff2(GEN nf, GEN beta, GEN B, long v, long prec)
2195: {
2196: long N, G, i, bacmin, bacmax, av = avma, av2;
2197: GEN vec, velt, p1, cand, M, plg, pol, cand2;
2198:
2199: M = gmael(nf, 5, 1);
2200: pol = (GEN)nf[1];
2201: N = degpol(pol);
2202: vec = gtrans((GEN)gtrans(M)[v]);
2203: velt = (GEN)nf[7];
2204:
2205: G = min( - 20, - bit_accuracy(prec) >> 4);
2206:
2207: p1 = cgetg(2, t_VEC);
2208:
2209: p1[1] = lneg(beta);
2210: vec = concat(p1, vec);
2211:
2212: p1[1] = zero;
2213: velt = concat(p1, velt);
2214:
2215: bacmin = (long)(.225 * bit_accuracy(prec));
2216: bacmax = (long)(.315 * bit_accuracy(prec));
2217:
2218: av2 = avma;
2219:
2220: for (i = bacmax; i >= bacmin; i--)
2221: {
2222: p1 = lindep2(vec, i);
2223:
2224: if (signe((GEN)p1[1]))
2225: {
2226: p1 = ground(gdiv(p1, (GEN)p1[1]));
2227: cand = gmodulcp(gmul(velt, gtrans(p1)), pol);
2228: cand2 = algtobasis(nf, cand);
2229: plg = gmul(M, cand2);
2230:
2231: if (TestOne(plg, beta, B, v, G, N))
2232: return gerepilecopy(av, cand);
2233: }
2234: avma = av2;
2235: }
2236: return NULL;
2237: }
2238:
2239: GEN
2240: chk_reccoeff_init(FP_chk_fun *chk, GEN nf, GEN gram, GEN mat, long *ptprec)
2241: {
2242: GEN data = chk->data;
2243: data[6] = (long)mat;
2244: chk->data = data;
2245: return (GEN)data[7];
2246: }
2247:
2248: GEN
2249: chk_reccoeff(GEN data, GEN x)
2250: {
2251: GEN M = (GEN)data[0], beta = (GEN)data[1], B = (GEN)data[2];
2252: long v = data[3], G = data[4], N = data[5], j;
2253: GEN U = (GEN)data[6], p1 = gmul(U, x), sol, plg;
2254:
2255: if (!gcmp1((GEN)p1[1])) return NULL;
2256:
2257: sol = cgetg(N + 1, t_COL);
2258: for (j = 1; j <= N; j++)
2259: sol[j] = lmulii((GEN)p1[1], (GEN)p1[j + 1]);
2260: plg = gmul(M, sol);
2261:
2262: if (TestOne(plg, beta, B, v, G, N)) return sol;
2263: return NULL;
2264: }
2265:
2266: GEN
2267: chk_reccoeff_post(GEN data, GEN res)
2268: {
2269: return res;
2270: }
2271:
2272: /* Using Cohen's method */
2273: static GEN
2274: RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
2275: {
2276: GEN A, M, nB, cand, p1, B2, C2, data, tB, beta2, eps, nf2, Bd;
2277: long N, G, i, j, k, l, ct = 0, av = avma, prec2;
2278: FP_chk_fun *chk;
2279:
2280: N = degpol(nf[1]);
2281: G = min(-10, -bit_accuracy(prec) >> 4);
2282: eps = gpowgs(stoi(10), min(-8, (G >> 1)));
2283: tB = gpow(gmul2n(eps, N), gdivgs(gun, 1-N), DEFAULTPREC);
2284:
2285: Bd = gmin(B, tB);
2286: p1 = gceil(gdiv(glog(Bd, DEFAULTPREC), dbltor(2.3026)));
2287: prec2 = max((prec << 1) - 2, (long)(itos(p1) * pariK1 + BIGDEFAULTPREC));
2288: nf2 = nfnewprec(nf, prec2);
2289: beta2 = gprec_w(beta, prec2);
2290:
2291: LABrcf: ct++;
2292: B2 = sqri(Bd);
2293: C2 = gdiv(B2, gsqr(eps));
2294:
2295: M = gmael(nf2, 5, 1);
2296:
2297: A = cgetg(N+2, t_MAT);
2298: for (i = 1; i <= N+1; i++)
2299: A[i] = lgetg(N+2, t_COL);
2300:
2301: coeff(A, 1, 1) = ladd(gmul(C2, gsqr(beta2)), B2);
2302: for (j = 2; j <= N+1; j++)
2303: {
2304: p1 = gmul(C2, gmul(gneg_i(beta2), gcoeff(M, v, j-1)));
2305: coeff(A, 1, j) = coeff(A, j, 1) = (long)p1;
2306: }
2307: for (i = 2; i <= N+1; i++)
2308: for (j = 2; j <= N+1; j++)
2309: {
2310: p1 = gzero;
2311: for (k = 1; k <= N; k++)
2312: {
2313: GEN p2 = gmul(gcoeff(M, k, j-1), gcoeff(M, k, i-1));
2314: if (k == v) p2 = gmul(C2, p2);
2315: p1 = gadd(p1,p2);
2316: }
2317: coeff(A, i, j) = coeff(A, j, i) = (long)p1;
2318: }
2319:
2320: nB = mulsi(N+1, B2);
2321:
2322: data = new_chunk(8);
2323: data[0] = (long)M;
2324: data[1] = (long)beta;
2325: data[2] = (long)B;
2326: data[3] = v;
2327: data[4] = G;
2328: data[5] = N;
2329: data[6] = (long)NULL;
2330: data[7] = (long)nB;
2331:
2332: chk = (FP_chk_fun*)new_chunk(sizeof(FP_chk_fun));
2333: chk->f = &chk_reccoeff;
2334: chk->f_init = &chk_reccoeff_init;
2335: chk->f_post = &chk_reccoeff_post;
2336: chk->data = data;
2337: chk->skipfirst = 0;
2338:
2339: cand = fincke_pohst(A, nB, 20000, 3, prec2, chk);
2340:
2341: if (!cand)
2342: {
2343: if (ct > 3) { avma = av; return NULL; }
2344:
2345: prec2 = (prec2 << 1) - 2;
2346: if (DEBUGLEVEL >= 2) err(warnprec,"RecCoeff", prec2);
2347: nf2 = nfnewprec(nf2, prec2);
2348: beta2 = gprec_w(beta2, prec2);
2349: goto LABrcf;
2350: }
2351:
2352: cand = (GEN)cand[1];
2353: l = lg(cand) - 1;
2354:
2355: if (l == 1) return gerepileupto(av, basistoalg(nf, (GEN)cand[1]));
2356:
2357: if (DEBUGLEVEL >= 2)
2358: fprintferr("RecCoeff3: no solution found!\n");
2359:
2360: avma = av; return NULL;
2361: }
2362:
2363: /* Attempts to find a polynomial with coefficients in nf such that
2364: its coefficients are close to those of pol at the place v and
2365: less than B at all the other places */
2366: GEN
2367: RecCoeff(GEN nf, GEN pol, long v, long prec)
2368: {
2369: long av = avma, j, md, G, cl = degpol(pol);
2370: GEN p1, beta;
2371:
2372: /* if precision(pol) is too low, abort */
2373: for (j = 2; j <= cl+1; j++)
2374: {
2375: p1 = (GEN)pol[j];
2376: G = bit_accuracy(gprecision(p1)) - gexpo(p1);
2377: if (G < 34) { avma = av; return NULL; }
2378: }
2379:
2380: md = cl/2;
2381: pol = dummycopy(pol);
2382:
2383: for (j = 1; j <= cl; j++)
2384: {
2385: /* start with the coefficients in the middle,
2386: since they are the harder to recognize! */
2387: long cf = md + (j%2? j/2: -j/2);
2388: GEN bound = binome(stoi(cl), cf);
2389:
2390: bound = shifti(bound, cl - cf);
2391:
2392: if (DEBUGLEVEL > 1) fprintferr("In RecCoeff with cf = %ld and B = %Z\n",
2393: cf, bound);
2394:
2395: beta = greal((GEN)pol[cf+2]);
2396: p1 = RecCoeff2(nf, beta, bound, v, prec);
2397: if (!p1)
2398: {
2399: p1 = RecCoeff3(nf, beta, bound, v, prec);
2400: if (!p1) return NULL;
2401: }
2402: pol[cf+2] = (long)p1;
2403: }
2404: pol[cl+2] = un;
2405: return gerepilecopy(av, pol);
2406: }
2407:
2408: /*******************************************************************/
2409: /*******************************************************************/
2410: /* */
2411: /* Computation of class fields of */
2412: /* real quadratic fields using Stark units */
2413: /* */
2414: /*******************************************************************/
2415: /*******************************************************************/
2416:
2417: /* compute the coefficients an for the quadratic case */
2418: static int***
2419: computean(GEN dtcr, long nmax, long prec)
2420: {
2421: long i, j, cl, q, cp, al, v1, v2, v, fldiv, av, av1;
2422: int ***matan, ***reduc;
2423: GEN bnf, ideal, dk, degs, idno, p1, prime, chi, qg, chi1, chi2;
2424: GEN chi11, chi12, bnr, pr, pr1, pr2, xray, xray1, xray2, dataray;
2425: byteptr dp = diffptr;
2426:
2427: av = avma;
2428:
2429: cl = lg(dtcr) - 1;
2430: degs = GetDeg(dtcr);
2431:
2432: matan = InitMatAn(cl, nmax, degs, 1);
2433: reduc = InitReduction(dtcr, degs);
2434:
2435: bnr = gmael(dtcr, 1, 4); bnf = (GEN)bnr[1];
2436: dataray = InitGetRay(bnr, nmax);
2437:
2438: ideal = gmael3(bnr, 2, 1, 1);
2439: idno = idealnorm(bnf, ideal);
2440: dk = gmael(bnf, 7, 3);
2441:
2442: prime = stoi(2);
2443: dp++;
2444:
2445: av1 = avma;
2446:
2447: chi = chi1 = chi2 = NULL; /* gcc -Wall */
2448: while (*dp && prime[2] <= nmax)
2449: {
2450: qg = prime;
2451: al = 1;
2452:
2453: switch (krogs(dk, prime[2]))
2454: {
2455: /* prime is inert */
2456: case -1:
2457: fldiv = divise(idno, prime);
2458:
2459: if (!fldiv)
2460: {
2461: xray = GetRay(bnr, dataray, prime, prec);
2462: chi = chiideal(dtcr, xray, 1);
2463: chi1 = dummycopy(chi);
2464: }
2465:
2466: while(cmpis(qg, nmax) <= 0)
2467: {
2468: q = qg[2];
2469:
2470: for (cp = 1, i = q; i <= nmax; i += q, cp++)
2471: if(cp % prime[2])
2472: {
2473: if (fldiv || al%2)
2474: for (j = 1; j <= cl; j++)
2475: _0toCoeff(matan[j][i], degs[j]);
2476: else
2477: for (j = 1; j <= cl; j++)
2478: MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]);
2479: }
2480:
2481: qg = mulsi(q, prime);
2482: al++;
2483:
2484: if (al%2 && !fldiv)
2485: for (j = 1; j <= cl; j++)
2486: chi[j] = lmul((GEN)chi[j], (GEN)chi1[j]);
2487: }
2488: break;
2489:
2490: /* prime is ramified */
2491: case 0:
2492: fldiv = divise(idno, prime);
2493:
2494: if (!fldiv)
2495: {
2496: pr = (GEN)primedec(bnf, prime)[1];
2497: xray = GetRay(bnr, dataray, pr, prec);
2498: chi = chiideal(dtcr, xray, 1);
2499: chi2 = dummycopy(chi);
2500: }
2501:
2502: while(cmpis(qg, nmax) <= 0)
2503: {
2504: q = qg[2];
2505:
2506: for (cp = 1, i = q; i <= nmax; i += q, cp++)
2507: if(cp % prime[2])
2508: {
2509: if (fldiv)
2510: for(j = 1; j <= cl; j++)
2511: _0toCoeff(matan[j][i], degs[j]);
2512: else
2513: {
2514: for(j = 1; j <= cl; j++)
2515: MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]);
2516: }
2517: }
2518:
2519: qg = mulsi(q, prime);
2520: al++;
2521:
2522: if (cmpis(qg, nmax) <= 0 && !fldiv)
2523: for (j = 1; j <= cl; j++)
2524: chi[j] = lmul((GEN)chi[j], (GEN)chi2[j]);
2525: }
2526: break;
2527:
2528: /* prime is split */
2529: default: /* case 1: */
2530: p1 = primedec(bnf, prime);
2531: pr1 = (GEN)p1[1];
2532: pr2 = (GEN)p1[2];
2533: v1 = idealval(bnf, ideal, pr1);
2534: v2 = idealval(bnf, ideal, pr2);
2535:
2536: if (v1 + v2 == 0)
2537: {
2538: xray1 = GetRay(bnr, dataray, pr1, prec);
2539: xray2 = GetRay(bnr, dataray, pr2, prec);
2540: chi11 = chiideal(dtcr, xray1, 1);
2541: chi12 = chiideal(dtcr, xray2, 1);
2542:
2543: chi1 = gadd(chi11, chi12);
2544: chi2 = dummycopy(chi12);
2545:
2546: while(cmpis(qg, nmax) <= 0)
2547: {
2548: q = qg[2];
2549:
2550: for (cp = 1, i = q; i <= nmax; i += q, cp++)
2551: if(cp % prime[2])
2552: for(j = 1; j <= cl; j++)
2553: MulPolmodCoeff((GEN)chi1[j], matan[j][i], reduc[j], degs[j]);
2554:
2555: qg = mulsi(q, prime);
2556: al++;
2557:
2558: if(cmpis(qg, nmax) <= 0)
2559: for (j = 1; j <= cl; j++)
2560: {
2561: chi2[j] = lmul((GEN)chi2[j], (GEN)chi12[j]);
2562: chi1[j] = ladd((GEN)chi2[j], gmul((GEN)chi1[j], (GEN)chi11[j]));
2563: }
2564: }
2565: }
2566: else
2567: {
2568: if (v1) { v = v2; pr = pr2; } else { v = v1; pr = pr1; }
2569:
2570: if (v == 0)
2571: {
2572: xray = GetRay(bnr, dataray, pr, prec);
2573: chi1 = chiideal(dtcr, xray, 1);
2574: chi = gcopy(chi1);
2575: }
2576:
2577: while(cmpis(qg, nmax) <= 0)
2578: {
2579: q = qg[2];
2580: for (cp = 1, i = q; i <= nmax; i += q, cp++)
2581: if(cp % prime[2])
2582: {
2583: if (v)
2584: for (j = 1; j <= cl; j++)
2585: _0toCoeff(matan[j][i], degs[j]);
2586: else
2587: for (j = 1; j <= cl; j++)
2588: MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]);
2589: }
2590:
2591: qg = mulii(qg, prime);
2592: al++;
2593:
2594: if (!v && (cmpis(qg, nmax) <= 0))
2595: for (j = 1; j <= cl; j++)
2596: chi[j] = lmul((GEN)chi[j], (GEN)chi1[j]);
2597: }
2598: }
2599: break;
2600: }
2601:
2602: prime[2] += (*dp++);
2603:
2604: avma = av1;
2605: }
2606:
2607: for (i = 1; i <= cl; i++)
2608: CorrectCoeff((GEN)dtcr[i], matan[i], reduc[i], nmax, degs[i]);
2609:
2610: FreeMat(reduc);
2611: avma = av; return matan;
2612: }
2613:
2614: /* compute S and T for the quadratic case */
2615: static GEN
2616: QuadGetST(GEN data, long prec)
2617: {
2618: long av = avma, n, j, nmax, cl, av1, av2, k;
2619: int ***matan;
2620: GEN nn, C, p1, p2, c2, cexp, cn, v, veclprime2, veclprime1;
2621: GEN dtcr, cond, rep, an, cf, degs, veint1;
2622:
2623: dtcr = (GEN)data[5];
2624: cl = lg(dtcr) - 1;
2625: degs = GetDeg(dtcr);
2626:
2627: cf = gmul2n(mpsqrt(mppi(prec)), 1);
2628: C = cgetg(cl+1, t_VEC);
2629: cond = cgetg(cl+1, t_VEC);
2630: c2 = cgetg(cl + 1, t_VEC);
2631: nn = new_chunk(cl+1);
2632: nmax = 0;
2633: for (j = 1; j <= cl; j++)
2634: {
2635: C[j] = mael(dtcr, j, 2);
2636: c2[j] = ldivsg(2, (GEN)C[j]);
2637: cond[j] = mael(dtcr, j, 7);
2638: nn[j] = (long)(bit_accuracy(prec) * gtodouble((GEN)C[j]) * 0.35);
2639:
2640: nmax = max(nmax, nn[j]);
2641: }
2642:
2643: if (nmax >= VERYBIGINT)
2644: err(talker, "Too many coefficients (%ld) in QuadGetST: computation impossible", nmax);
2645:
2646: if (DEBUGLEVEL >= 2)
2647: fprintferr("nmax = %ld\n", nmax);
2648:
2649: /* compute the coefficients */
2650: matan = computean(dtcr, nmax, prec);
2651: if (DEBUGLEVEL) msgtimer("Compute an");
2652:
2653: /* allocate memory for the answer */
2654: rep = cgetg(3, t_VEC);
2655:
2656: /* allocate memory for veclprime1 */
2657: veclprime1 = cgetg(cl + 1, t_VEC);
2658: for (j = 1; j <= cl; j++)
2659: {
2660: v = cgetg(3, t_COMPLEX);
2661: v[1] = lgetr(prec);
2662: v[2] = lgetr(prec); gaffect(gzero, v);
2663: veclprime1[j] = (long)v;
2664: }
2665:
2666: av1 = avma;
2667: cn = cgetr(prec);
2668: p1 = gmul2n(cf, -1);
2669:
2670: /* compute veclprime1 */
2671: for (j = 1; j <= cl; j++)
2672: {
2673: long n0 = 0;
2674: p2 = gmael3(dtcr, j, 5, 2);
2675: cexp = gexp(gneg_i((GEN)c2[j]), prec);
2676: av2 = avma; affsr(1, cn); v = (GEN)veclprime1[j];
2677: for (n = 1; n <= nn[j]; n++)
2678: if ( (an = EvalCoeff(p2, matan[j][n], degs[j])) )
2679: {
2680: affrr(gmul(cn, gpowgs(cexp, n - n0)), cn);
2681: n0 = n;
2682: gaffect(gadd(v, gmul(divrs(cn,n), an)), v);
2683: avma = av2;
2684: }
2685: gaffect(gmul(p1, gmul(v, (GEN)C[j])), v);
2686: avma = av2;
2687: }
2688: avma = av1;
2689: rep[1] = (long)veclprime1;
2690: if (DEBUGLEVEL) msgtimer("Compute V1");
2691:
2692: /* allocate memory for veclprime2 */
2693: veclprime2 = cgetg(cl + 1, t_VEC);
2694: for (j = 1; j <= cl; j++)
2695: {
2696: v = cgetg(3, t_COMPLEX);
2697: v[1] = lgetr(prec);
2698: v[2] = lgetr(prec); gaffect(gzero, v);
2699: veclprime2[j] = (long)v;
2700: }
2701:
2702: /* compute f1(C/n) */
2703: av1 = avma;
2704:
2705: veint1 = cgetg(cl + 1, t_VEC);
2706: for (j = 1; j <= cl; j++)
2707: {
2708: p1 = NULL;
2709: for (k = 1; k < j; k++)
2710: if (gegal((GEN)cond[j], (GEN)cond[k])) { p1 = (GEN)veint1[k]; break; }
2711: if (p1 == NULL)
2712: {
2713: p1 = veceint1((GEN)c2[j], stoi(nn[j]), prec);
2714: veint1[j] = (long)p1;
2715: }
2716: av2 = avma; p2 = gmael3(dtcr, j, 5, 2);
2717: v = (GEN)veclprime2[j];
2718: for (n = 1; n <= nn[j]; n++)
2719: if ( (an = EvalCoeff(p2, matan[j][n], degs[j])) )
2720: {
2721: gaffect(gadd(v, gmul((GEN)p1[n], an)), v);
2722: avma = av2;
2723: }
2724: gaffect(gmul(cf, gconj(v)), v);
2725: avma = av2;
2726: }
2727: avma = av1;
2728: rep[2] = (long)veclprime2;
2729: if (DEBUGLEVEL) msgtimer("Compute V2");
2730: FreeMat(matan); return gerepileupto(av, rep);
2731: }
2732:
2733: #if 0
2734: /* recover a quadratic integer by an exhaustive search */
2735: static GEN
2736: recbeta2(GEN nf, GEN beta, GEN bound, long prec)
2737: {
2738: long av = avma, av2, tetpil, i, range, G, e, m;
2739: GEN om, om1, om2, dom, p1, a, b, rom, bom2, *gptr[2];
2740:
2741: G = min( - 20, - bit_accuracy(prec) >> 4);
2742:
2743: if (DEBUGLEVEL > 3)
2744: fprintferr("\n Precision needed: %ld", G);
2745:
2746: om = gmael(nf, 7, 2);
2747: rom = (GEN)nf[6];
2748: om1 = poleval(om, (GEN)rom[2]);
2749: om2 = poleval(om, (GEN)rom[1]);
2750: dom = subrr(om1, om2);
2751:
2752: /* b will run from b to b + range */
2753: p1 = gaddgs(gmul2n(gceil(absr(divir(bound, dom))), 1), 2);
2754: range = VERYBIGINT;
2755: if (cmpis(p1, VERYBIGINT) < 0)
2756: range = itos(p1);
2757:
2758: av2 = avma;
2759:
2760: b = gdiv(gsub(bound, beta), dom);
2761: if (gsigne(b) < 0)
2762: b = subis(negi(gcvtoi(gneg_i(b), &e)), 1);
2763: else
2764: b=gcvtoi(b, &e);
2765:
2766: if (e > 0) /* precision is lost in truncation */
2767: {
2768: avma = av;
2769: return NULL;
2770: }
2771:
2772: bom2 = mulir(b, om2);
2773: m = 0;
2774:
2775: for (i = 0; i <= range; i++)
2776: {
2777: /* for every b, we construct a and test it */
2778: a = grndtoi(gsub(beta, bom2), &e);
2779:
2780: if (e > 0) /* precision is lost in truncation */
2781: {
2782: avma = av;
2783: return NULL;
2784: }
2785:
2786: p1 = gsub(mpadd(a, bom2), beta);
2787:
2788: if ((DEBUGLEVEL > 3) && (expo(p1)<m))
2789: {
2790: m = expo(p1);
2791: fprintferr("\n Precision found: %ld", expo(p1));
2792: }
2793:
2794: if (gcmp0(p1) || (expo(p1) < G)) /* result found */
2795: {
2796: p1 = gadd(a, gmul(b, om));
2797: return gerepileupto(av, gmodulcp(p1, (GEN)nf[1]));
2798: }
2799:
2800: tetpil = avma;
2801:
2802: b = gaddgs(b, 1);
2803: bom2 = gadd(bom2, om2);
2804:
2805: gptr[0] = &b;
2806: gptr[1] = &bom2;
2807: gerepilemanysp(av2, tetpil, gptr, 2);
2808: }
2809:
2810: /* if it fails... */
2811: return NULL;
2812: }
2813: #endif
2814:
2815: /* return 1 if the absolute polynomial pol (over Q) defines the
2816: Hilbert class field of the real quadratic field bnf */
2817: int
2818: define_hilbert(GEN bnf, GEN pol)
2819: {
2820: long cl;
2821: GEN dk;
2822:
2823: cl = itos(gmael3(bnf, 8, 1, 1));
2824: dk = gmael(bnf, 7, 3);
2825:
2826: if (degpol(pol) == cl)
2827: if ((cl%2) || !egalii(discf(pol), gpowgs(dk,cl>>1))) return 1;
2828:
2829: return 0;
2830: }
2831:
2832: /* let polrel define Hk/k, find L/Q such that Hk=Lk and L and k are
2833: disjoint */
2834: static GEN
2835: makescind(GEN bnf, GEN polabs, long cl, long prec)
2836: {
2837: long av = avma, i, l;
2838: GEN pol, p1, nf2, dabs, dk, bas;
2839:
2840: /* check the result (a little): signature and discriminant */
2841: bas = allbase4(polabs,0,&dabs,NULL);
2842: dk = gmael(bnf,7,3);
2843: if (!egalii(dabs, gpowgs(dk,cl)) || sturm(polabs) != 2*cl)
2844: err(bugparier, "quadhilbert");
2845:
2846: /* attempt to find the subfields using polred */
2847: p1 = cgetg(3,t_VEC); p1[1]=(long)polabs; p1[2]=(long)bas;
2848: pol = polredfirstpol(p1, (prec<<1) - 2, &define_hilbert, bnf);
2849: if (DEBUGLEVEL) msgtimer("polred");
2850:
2851: if (!pol)
2852: {
2853: nf2 = nfinit0(polabs, 1, prec);
2854: p1 = subfields(nf2, stoi(cl));
2855: l = lg(p1);
2856: if (DEBUGLEVEL) msgtimer("subfields");
2857:
2858: for (i = 1; i < l; i++)
2859: {
2860: pol = gmael(p1, i, 1);
2861: if ((cl%2) || !gegal(sqri(discf(pol)), (GEN)nf2[3])) break;
2862: }
2863: if (i == l)
2864: for (i = 1; i < l; i++)
2865: {
2866: pol = gmael(p1, i, 1);
2867: if (degpol(gcoeff(nffactor(bnf, pol), 1, 1)) == cl) break;
2868: }
2869: if (i == l)
2870: err(bugparier, "makescind (no polynomial found)");
2871: }
2872: pol = polredabs(pol, prec);
2873: return gerepileupto(av, pol);
2874: }
2875:
2876: /* compute the Hilbert class field using genus class field theory when
2877: the exponent of the class group is 2 */
2878: static GEN
2879: GenusField(GEN bnf, long prec)
2880: {
2881: long hk, c, l, av = avma;
2882: GEN disc, quat, x2, pol, div, d;
2883:
2884: hk = itos(gmael3(bnf, 8, 1, 1));
2885: disc = gmael(bnf, 7, 3);
2886: quat = stoi(4);
2887: x2 = gsqr(polx[0]);
2888:
2889: if (gcmp0(modii(disc, quat))) disc = divii(disc, quat);
2890:
2891: div = divisors(disc);
2892: c = 1;
2893: l = 0;
2894: pol = NULL; /* gcc -Wall */
2895:
2896: while(l < hk)
2897: {
2898: c++;
2899: d = (GEN)div[c];
2900:
2901: if (gcmp1(modii(d, quat)))
2902: {
2903: if (!l)
2904: pol = gsub(x2, d);
2905: else
2906: pol=(GEN)compositum(pol, gsub(x2, d))[1];
2907:
2908: l = degpol(pol);
2909: }
2910: }
2911:
2912: return gerepileupto(av, polredabs(pol, prec));
2913: }
2914:
2915: /* if flag = 0 returns the reduced polynomial, flag = 1 returns the
2916: non-reduced polynomial, flag = 2 returns an absolute reduced
2917: polynomial, flag = 3 returns the polynomial of the Stark's unit,
2918: flag = -1 computes a fast and crude approximation of the result */
2919: static GEN
2920: AllStark(GEN data, GEN nf, long flag, long newprec)
2921: {
2922: long cl, i, j, cpt = 0, av, av2, N, h, v, n, bnd = 300, sq = 1, r1, r2;
2923: int ***matan;
2924: GEN p0, p1, p2, S, T, polrelnum, polrel, Lp, W, A, veczeta, sig, valchi;
2925: GEN degs, ro, C, Cmax, dataCR, cond1, L1, *gptr[2], an, Pi;
2926:
2927: N = degpol(nf[1]);
2928: r1 = nf_get_r1(nf);
2929: r2 = (N - r1)>>1;
2930: cond1 = gmael4(data, 1, 2, 1, 2);
2931: Pi = mppi(newprec);
2932:
2933: v = 1;
2934: while(gcmp1((GEN)cond1[v])) v++;
2935:
2936: LABDOUB:
2937:
2938: av = avma;
2939:
2940: dataCR = (GEN)data[5];
2941: cl = lg(dataCR)-1;
2942: degs = GetDeg(dataCR);
2943: h = itos(gmul2n(det((GEN)data[2]), -1));
2944:
2945: if (flag >= 0)
2946: {
2947: /* compute S,T differently if nf is quadratic */
2948: if (N == 2)
2949: p1 = QuadGetST(data, newprec);
2950: else
2951: p1 = GetST(dataCR, newprec);
2952:
2953: S = (GEN)p1[1];
2954: T = (GEN)p1[2];
2955:
2956: Lp = cgetg(cl + 1, t_VEC);
2957: for (i = 1; i <= cl; i++)
2958: Lp[i] = GetValue((GEN)dataCR[i], (GEN)S[i], (GEN)T[i], 0, 1, newprec)[2];
2959:
2960: if (DEBUGLEVEL) msgtimer("Compute W");
2961: }
2962: else
2963: {
2964: /* compute a crude approximation of the result */
2965: C = cgetg(cl + 1, t_VEC);
2966: for (i = 1; i <= cl; i++) C[i] = mael(dataCR, i, 2);
2967: Cmax = vecmax(C);
2968:
2969: n = GetBoundN0(Cmax, r1, r2, newprec, 0);
2970: if (n > bnd) n = bnd;
2971: if (DEBUGLEVEL) fprintferr("nmax in QuickPol: %ld \n", n);
2972:
2973: matan = ComputeCoeff(dataCR, n, newprec);
2974:
2975: p0 = cgetg(3, t_COMPLEX);
2976: p0[1] = lgetr(newprec); affsr(0, (GEN)p0[1]);
2977: p0[2] = lgetr(newprec); affsr(0, (GEN)p0[2]);
2978:
2979: L1 = cgetg(cl+1, t_VEC);
2980: /* we use the formulae L(1) = sum (an / n) */
2981: for (i = 1; i <= cl; i++)
2982: {
2983: av2 = avma;
2984: p1 = p0; p2 = gmael3(dataCR, i, 5, 2);
2985: for (j = 1; j <= n; j++)
2986: if ( (an = EvalCoeff(p2, matan[i][j], degs[i])) )
2987: p1 = gadd(p1, gdivgs(an, j));
2988: L1[i] = lpileupto(av2, p1);
2989: }
2990: FreeMat(matan);
2991:
2992: p1 = gmul2n(gpowgs(mpsqrt(Pi), N - 2), 1);
2993:
2994: Lp = cgetg(cl+1, t_VEC);
2995: for (i = 1; i <= cl; i++)
2996: {
2997: W = ComputeArtinNumber((GEN)dataCR[i], 1, newprec);
2998: A = (GEN)ComputeAChi((GEN)dataCR[i], 0, newprec)[2];
2999: W = gmul((GEN)C[i], gmul(A, W));
3000:
3001: Lp[i] = ldiv(gmul(W, gconj((GEN)L1[i])), p1);
3002: }
3003: }
3004:
3005: p1 = ComputeLift(gmael(data, 4, 2));
3006:
3007: veczeta = cgetg(h + 1, t_VEC);
3008: for (i = 1; i <= h; i++)
3009: {
3010: GEN z = gzero;
3011:
3012: sig = (GEN)p1[i];
3013: valchi = chiideal(dataCR, sig, 0);
3014:
3015: for (j = 1; j <= cl; j++)
3016: {
3017: GEN p2 = greal(gmul((GEN)Lp[j], (GEN)valchi[j]));
3018: if (!gegal(gdeux, gmael3(dataCR, j, 5, 3)))
3019: p2 = gmul2n(p2, 1); /* character not real */
3020: z = gadd(z,p2);
3021: }
3022: veczeta[i] = ldivgs(z, 2 * h);
3023: }
3024: if (DEBUGLEVEL >= 2) fprintferr("zetavalues = %Z\n", veczeta);
3025:
3026: if ((flag >=0) && (flag <= 3)) sq = 0;
3027:
3028: ro = cgetg(h+1, t_VEC); /* roots */
3029:
3030: for (;;)
3031: {
3032: if (!sq && (DEBUGLEVEL > 1))
3033: fprintferr("Checking the square-root of the Stark unit...\n");
3034:
3035: for (j = 1; j <= h; j++)
3036: {
3037: p1 = gexp(gmul2n((GEN)veczeta[j], sq), newprec);
3038: ro[j] = ladd(p1, ginv(p1));
3039: }
3040: polrelnum = roots_to_pol_intern(realun(newprec),ro, 0,0);
3041: if (DEBUGLEVEL)
3042: {
3043: if (DEBUGLEVEL >= 2) fprintferr("polrelnum = %Z\n", polrelnum);
3044: msgtimer("Compute %s", (flag < 0)? "quickpol": "polrelnum");
3045: }
3046:
3047: if (flag < 0)
3048: return gerepilecopy(av, polrelnum);
3049:
3050: /* we try to recognize this polynomial */
3051: polrel = RecCoeff(nf, polrelnum, v, newprec);
3052:
3053: if (polrel || (sq++ == 1)) break;
3054: }
3055:
3056: if (!polrel) /* if it fails... */
3057: {
3058: long pr;
3059: if (++cpt >= 3) err(precer, "stark (computation impossible)");
3060:
3061: /* we compute the precision that we need */
3062: pr = 1 + (gexpo(polrelnum)>>TWOPOTBITS_IN_LONG);
3063: if (pr < 0) pr = 0;
3064: newprec = ADD_PREC + max(newprec,pr);
3065:
3066: if (DEBUGLEVEL) err(warnprec, "AllStark", newprec);
3067:
3068: nf = nfnewprec(nf, newprec);
3069: data[5] = (long)CharNewPrec((GEN)data[5], nf, newprec);
3070:
3071: gptr[0] = &data;
3072: gptr[1] = &nf;
3073: gerepilemany(av, gptr, 2);
3074:
3075: goto LABDOUB;
3076: }
3077:
3078: /* and we compute the polynomial of eps if flag = 3 */
3079: if (flag == 3)
3080: {
3081: n = fetch_var();
3082: p1 = gsub(polx[0], gadd(polx[n], ginv(polx[n])));
3083: polrel = polresultant0(polrel, p1, 0, 0);
3084: polrel = gmul(polrel, gpowgs(polx[n], h));
3085: polrel = gsubst(polrel, n, polx[0]);
3086: polrel = gmul(polrel, leading_term(polrel));
3087: delete_var();
3088: }
3089:
3090: if (DEBUGLEVEL >= 2) fprintferr("polrel = %Z\n", polrel);
3091: if (DEBUGLEVEL) msgtimer("Recpolnum");
3092:
3093: /* we want a reduced relative polynomial */
3094: if (!flag) return gerepileupto(av, rnfpolredabs(nf, polrel, 0, newprec));
3095:
3096: /* we just want the polynomial computed */
3097: if (flag!=2) return gerepilecopy(av, polrel);
3098:
3099: /* we want a reduced absolute polynomial */
3100: return gerepileupto(av, rnfpolredabs(nf, polrel, 2, newprec));
3101: }
3102:
3103: /********************************************************************/
3104: /* Main functions */
3105: /********************************************************************/
3106:
3107: /* compute the polynomial over Q of the Hilbert class field of
3108: Q(sqrt(D)) where D is a positive fundamental discriminant */
3109: GEN
3110: quadhilbertreal(GEN D, GEN flag, long prec)
3111: {
3112: VOLATILE long av = avma, cl;
3113: long newprec;
3114: VOLATILE GEN pol, bnf, bnr, dataC, bnrh, nf, exp;
3115: void *catcherr = NULL;
3116:
3117: if (DEBUGLEVEL) timer2();
3118:
3119: disable_dbg(0);
3120: /* quick computation of the class number */
3121:
3122: cl = itos((GEN)quadclassunit0(D, 0, NULL, prec)[1]);
3123: if (cl == 1)
3124: {
3125: disable_dbg(-1);
3126: avma = av; return polx[0];
3127: }
3128:
3129: /* initialize the polynomial defining Q(sqrt{D}) as a polynomial in y */
3130: pol = quadpoly(D);
3131: setvarn(pol, fetch_var());
3132:
3133: START:
3134: /* compute the class group */
3135: bnf = bnfinit0(pol, 1, NULL, prec);
3136: nf = (GEN)bnf[7];
3137: disable_dbg(-1);
3138:
3139: if (DEBUGLEVEL) msgtimer("Compute Cl(k)");
3140:
3141: /* if the exponent of the class group is 2, use rather Genus Field Theory */
3142: exp = gmael4(bnf, 8, 1, 2, 1);
3143: if (gegal(exp, gdeux)) { delete_var(); return GenusField(bnf, prec); }
3144:
3145: { /* catch precision problems (precision too small) */
3146: jmp_buf env;
3147: if (setjmp(env))
3148: {
3149: prec += EXTRA_PREC;
3150: err (warnprec, "quadhilbertreal", prec);
3151: goto START;
3152: }
3153: catcherr = err_catch(precer, env, NULL);
3154: }
3155:
3156: /* find the modulus defining N */
3157: bnr = buchrayinitgen(bnf, gun);
3158: dataC = InitQuotient(bnr, gzero);
3159: bnrh = FindModulus(dataC, 1, &newprec, prec, gcmp0(flag)? 0: -10);
3160:
3161: if (DEBUGLEVEL) msgtimer("FindModulus");
3162:
3163: if (newprec > prec)
3164: {
3165: if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec);
3166: nf = nfnewprec(nf, newprec);
3167: }
3168:
3169: /* use the generic function AllStark */
3170: pol = AllStark(bnrh, nf, 2, newprec);
3171: delete_var();
3172: return gerepileupto(av, makescind(bnf, pol, cl, prec));
3173: }
3174:
3175: GEN
3176: bnrstark(GEN bnr, GEN subgroup, long flag, long prec)
3177: {
3178: long cl, N, newprec, av = avma, bnd = 0;
3179: GEN bnf, dataS, p1, Mcyc, nf, data;
3180:
3181: if (flag >= 4)
3182: {
3183: bnd = -10;
3184: flag -= 4;
3185: }
3186:
3187: if (flag < 0 || flag > 3) err(flagerr,"bnrstark");
3188:
3189: /* check the bnr */
3190: checkbnrgen(bnr);
3191:
3192: bnf = (GEN)bnr[1];
3193: nf = (GEN)bnf[7];
3194: Mcyc = diagonal(gmael(bnr, 5, 2));
3195: N = degpol(nf[1]);
3196: if (N == 1)
3197: err(talker, "the ground field must be distinct from Q");
3198:
3199: /* check the bnf */
3200: if (!varn(gmael(bnf, 7, 1)))
3201: err(talker, "main variable in bnrstark must not be x");
3202:
3203: if (cmpis(gmael3(bnf, 7, 2, 1), N))
3204: err(talker, "not a totally real ground base field in bnrstark");
3205:
3206: /* check the subgroup */
3207: if (gcmp0(subgroup))
3208: subgroup = Mcyc;
3209: else
3210: {
3211: p1 = gauss(subgroup, Mcyc);
3212: if (!gcmp1(denom(p1)))
3213: err(talker, "incorrect subgroup in bnrstark");
3214: }
3215:
3216: /* compute bnr(conductor) */
3217: p1 = conductor(bnr, subgroup, 2);
3218: bnr = (GEN)p1[2];
3219: subgroup = (GEN)p1[3];
3220:
3221: /* check the class field */
3222: if (!gcmp0(gmael3(bnr, 2, 1, 2)))
3223: err(talker, "not a totally real class field in bnrstark");
3224:
3225: cl = itos(det(subgroup));
3226: if (cl == 1) return polx[0];
3227:
3228: timer2();
3229:
3230: /* find a suitable extension N */
3231: dataS = InitQuotient(bnr, subgroup);
3232: data = FindModulus(dataS, 1, &newprec, prec, bnd);
3233:
3234: if (newprec > prec)
3235: {
3236: if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec);
3237: nf = nfnewprec(nf, newprec);
3238: }
3239:
3240: return gerepileupto(av, AllStark(data, nf, flag, newprec));
3241: }
3242:
3243: /* For each character of Cl(bnr)/sbgrp, compute L(1, chi) (or equivalently
3244: the first non-zero term c(chi) of the expansion at s = 0). The binary
3245: digits of flag mean 1: if 0 then compute the term c(chi) and return
3246: [r(chi), c(chi)] where r(chi) is the order of L(s, chi) at s = 0,
3247: or if 1 then compute the value at s = 1 (and in this case, only for
3248: non-trivial characters), 2: if 0 then compute the value of the
3249: primitive L-function associated to chi, if 1 then compute the value
3250: of the L-function L_S(s, chi) where S is the set of places dividing
3251: the modulus of bnr (and the infinite places), 3: return also the
3252: character */
3253: GEN
3254: bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
3255: {
3256: GEN bnf, nf, cyc, Mcyc, p1, L1, chi, lchi, clchi, allCR, listCR, dataCR;
3257: GEN S, T, rep, indCR, invCR, Qt;
3258: long N, cl, i, j, k, nc, lq, a, av = avma, ncc;
3259:
3260: checkbnr(bnr);
3261: bnf = (GEN)bnr[1];
3262: nf = (GEN)bnf[7];
3263: cyc = gmael(bnr, 5, 2);
3264: Mcyc = diagonal(cyc);
3265: ncc = lg(cyc) - 1;
3266: N = degpol(nf[1]);
3267:
3268: if (N == 1)
3269: err(talker, "the ground field must be distinct from Q");
3270:
3271: if ((flag < 0) || (flag > 8))
3272: err(flagerr,"bnrL1");
3273:
3274: /* check the bnr */
3275: checkbnrgen(bnr);
3276:
3277: /* compute bnr(conductor) */
3278: if (!(flag & 2))
3279: {
3280: p1 = conductor(bnr, gzero, 2);
3281: bnr = (GEN)p1[2];
3282: cyc = gmael(bnr, 5, 2);
3283: Mcyc = diagonal(cyc);
3284: }
3285:
3286: /* check the subgroup */
3287: if (gcmp0(sbgrp))
3288: sbgrp = Mcyc;
3289: else
3290: {
3291: if (lg(sbgrp) != ncc+1)
3292: err(talker, "incorrect subgroup in bnrL1");
3293: p1 = gauss(sbgrp, Mcyc);
3294: if (!gcmp1(denom(p1)))
3295: err(talker, "incorrect subgroup in bnrL1");
3296: }
3297:
3298: cl = labs(itos(det(sbgrp)));
3299: Qt = InitQuotient0(Mcyc, sbgrp);
3300: lq = lg((GEN)Qt[2]) - 1;
3301:
3302: /* compute all the characters */
3303: allCR = FindEltofGroup(cl, (GEN)Qt[2]);
3304:
3305:
3306: /* make a list of all non-trivial characters modulo conjugation */
3307: listCR = cgetg(cl, t_VEC);
3308: indCR = new_chunk(cl);
3309: invCR = new_chunk(cl);
3310:
3311: nc = 0;
3312:
3313: for (i = 1; i < cl; i++)
3314: {
3315: chi = (GEN)allCR[i];
3316:
3317: /* lift the character to a character on Cl(bnr) */
3318: lchi = cgetg(ncc + 1, t_VEC);
3319: for (j = 1; j <= ncc; j++)
3320: {
3321: p1 = gzero;
3322: for (k = 1; k <= lq; k++)
3323: p1 = gadd(p1, gdiv(mulii(gmael3(Qt, 3, j, k), (GEN)chi[k]),
3324: gmael(Qt, 2, k)));
3325: lchi[j] = lmodii(gmul(p1, (GEN)cyc[j]), (GEN)cyc[j]);
3326: }
3327: clchi = ConjChar(lchi, cyc);
3328:
3329: a = i;
3330: for (j = 1; j <= nc; j++)
3331: if (gegal(gmael(listCR, j, 1), clchi)) { a = -j; break; }
3332:
3333: if (a > 0)
3334: {
3335: nc++;
3336: listCR[nc] = lgetg(3, t_VEC);
3337: mael(listCR, nc, 1) = (long)lchi;
3338: mael(listCR, nc, 2) = (long)bnrconductorofchar(bnr, lchi);
3339:
3340: indCR[i] = nc;
3341: invCR[nc] = i;
3342: }
3343: else
3344: indCR[i] = -invCR[-a];
3345:
3346: allCR[i] = lcopy(lchi);
3347: }
3348:
3349: /* the trivial character has to be a row vector too! */
3350: allCR[cl] = ltrans((GEN)allCR[cl]);
3351:
3352: setlg(listCR, nc + 1);
3353: if (nc == 0) err(talker, "no non-trivial character in bnrL1");
3354:
3355: /* compute the data for these characters */
3356: dataCR = InitChar(bnr, listCR, prec);
3357:
3358: p1 = GetST(dataCR, prec);
3359:
3360: S = (GEN)p1[1];
3361: T = (GEN)p1[2];
3362:
3363: if (flag & 1)
3364: L1 = cgetg(cl, t_VEC);
3365: else
3366: L1 = cgetg(cl + 1, t_VEC);
3367:
3368: for (i = 1; i < cl; i++)
3369: {
3370: a = indCR[i];
3371: if (a > 0)
3372: L1[i] = (long)GetValue((GEN)dataCR[a], (GEN)S[a], (GEN)T[a], flag & 1,
3373: flag & 2, prec);
3374: }
3375:
3376: for (i = 1; i < cl; i++)
3377: {
3378: a = indCR[i];
3379: if (a < 0)
3380: L1[i] = lconj((GEN)L1[-a]);
3381: }
3382:
3383: if (!(flag & 1))
3384: L1[cl] = (long)GetValue1(bnr, flag & 2, prec);
3385: else
3386: cl--;
3387:
3388: if (flag & 4)
3389: {
3390: rep = cgetg(cl + 1, t_VEC);
3391: for (i = 1; i <= cl; i++)
3392: {
3393: p1 = cgetg(3, t_VEC);
3394: p1[1] = allCR[i];
3395: p1[2] = L1[i];
3396:
3397: rep[i] = (long)p1;
3398: }
3399: }
3400: else
3401: rep = L1;
3402:
3403: return gerepilecopy(av, rep);
3404: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to stark.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / modules