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