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