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