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