=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/modules/Attic/stark.c,v retrieving revision 1.1 retrieving revision 1.2 diff -u -p -r1.1 -r1.2 --- OpenXM_contrib/pari-2.2/src/modules/Attic/stark.c 2001/10/02 11:17:12 1.1 +++ OpenXM_contrib/pari-2.2/src/modules/Attic/stark.c 2002/09/11 07:27:06 1.2 @@ -1,4 +1,4 @@ -/* $Id: stark.c,v 1.1 2001/10/02 11:17:12 noro Exp $ +/* $Id: stark.c,v 1.2 2002/09/11 07:27:06 noro Exp $ Copyright (C) 2000 The PARI group. @@ -15,8 +15,7 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, /*******************************************************************/ /* */ -/* COMPUTATION OF STARK UNITS */ -/* OF TOTALLY REAL FIELDS */ +/* COMPUTATION OF STARK UNITS OF TOTALLY REAL FIELDS */ /* */ /*******************************************************************/ #include "pari.h" @@ -26,17 +25,48 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, #define ADD_PREC (DEFAULTPREC-2)*3 extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus); +extern GEN bnrGetSurj(GEN bnr1, GEN bnr2); +/* ComputeCoeff */ +typedef struct { + GEN L0, L1, L11, L2; /* VECSMALL of p */ + GEN *L1ray, *L11ray; /* precomputed isprincipalray(pr), pr | p */ + GEN *rayZ; /* precomputed isprincipalray(i), i < condZ */ + long condZ; /* generates cond(bnr) \cap Z, assumed small */ +} LISTray; + +/* Char evaluation */ +typedef struct { + long ord; + GEN *val, chi; +} CHI_t; + +/* RecCoeff */ +typedef struct { + GEN M, beta, B, U, nB; + long v, G, N; +} RC_data; + /********************************************************************/ /* Miscellaneous functions */ /********************************************************************/ - -/* Compute the image of logelt by chi as a complex number if flag = 0, - otherwise as a polmod, see InitChar in part 3 */ +/* exp(2iPi/den), assume den a t_INT */ +GEN +InitRU(GEN den, long prec) +{ + GEN c,s, z; + if (egalii(den, gdeux)) return stoi(-1); + gsincos(divri(gmul2n(mppi(prec),1), den), &s, &c, prec); + z = cgetg(3, t_COMPLEX); + z[1] = (long)c; + z[2] = (long)s; return z; +} +/* Compute the image of logelt by chi as a complex number + see InitChar in part 3 */ static GEN -ComputeImagebyChar(GEN chi, GEN logelt, long flag) +ComputeImagebyChar(GEN chi, GEN logelt) { - GEN gn = gmul((GEN)chi[1], logelt), x = (GEN)chi[flag? 4: 2]; + GEN gn = gmul((GEN)chi[1], logelt), x = (GEN)chi[2]; long d = itos((GEN)chi[3]), n = smodis(gn, d); /* x^d = 1 and, if d even, x^(d/2) = -1 */ if ((d & 1) == 0) @@ -47,6 +77,46 @@ ComputeImagebyChar(GEN chi, GEN logelt, long flag) return gpowgs(x, n); } +/* return n such that C(elt) = z^n */ +static long +EvalChar_n(CHI_t *C, GEN logelt) +{ + GEN n = gmul(C->chi, logelt); + return smodis(n, C->ord); +} +/* return C(elt) */ +static GEN +EvalChar(CHI_t *C, GEN logelt) +{ + return C->val[EvalChar_n(C, logelt)]; +} + +static void +init_CHI(CHI_t *c, GEN CHI, GEN z) +{ + long i, d = itos((GEN)CHI[3]); + GEN *v = (GEN*)new_chunk(d); + v[1] = z; + for (i=2; ichi = (GEN)CHI[1]; + c->ord = d; + c->val = v; +} + +/* as t_COMPLEX */ +static void +init_CHI_alg(CHI_t *c, GEN CHI) { + GEN z = gmodulcp(polx[0], cyclo(itos((GEN)CHI[3]), 0)); + init_CHI(c,CHI,z); +} +/* as t_POLMOD */ +static void +init_CHI_C(CHI_t *c, GEN CHI) { + init_CHI(c,CHI, (GEN)CHI[2]); +} + /* Compute the conjugate character */ static GEN ConjChar(GEN chi, GEN cyc) @@ -59,43 +129,49 @@ ConjChar(GEN chi, GEN cyc) p1[i] = zero; else p1[i] = lsubii((GEN)cyc[i], (GEN)chi[i]); - + return p1; } -/* compute the next element for FindEltofGroup */ -static GEN -NextEltofGroup(GEN cyc, long l, long adec) -{ - GEN p1; - long dj, j; +typedef struct { + long r; /* rank = lg(gen) */ + GEN j; /* current elt is gen[1]^j[1] ... gen[r]^j[r] */ + GEN cyc; /* t_VECSMALL of elementary divisors */ +} GROUP_t; - p1 = cgetg(l + 1, t_COL); - - for (j = 1; j <= l; j++) +static int +NextElt(GROUP_t *G) +{ + long i = 1; + if (G->r == 0) return 0; /* no more elt */ + while (++G->j[i] == G->cyc[i]) /* from 0 to cyc[i]-1 */ { - dj = itos((GEN)cyc[j]); - p1[j] = lstoi(adec%dj); - adec /= dj; + G->j[i] = 0; + if (++i > G->r) return 0; /* no more elt */ } - - return p1; + return i; /* we have multiplied by gen[i] */ } /* Compute all the elements of a group given by its SNF */ static GEN -FindEltofGroup(long order, GEN cyc) +EltsOfGroup(long order, GEN cyc) { - long l, i; + long i; GEN rep; + GROUP_t G; - l = lg(cyc)-1; + G.cyc = gtovecsmall(cyc); + G.r = lg(cyc)-1; + G.j = vecsmall_const(G.r, 0); rep = cgetg(order + 1, t_VEC); + rep[order] = (long)small_to_col(G.j); - for (i = 1; i <= order; i++) - rep[i] = (long)NextEltofGroup(cyc, l, i); - + for (i = 1; i < order; i++) + { + (void)NextElt(&G); + rep[i] = (long)small_to_col(G.j); + } return rep; } @@ -104,14 +180,15 @@ FindEltofGroup(long order, GEN cyc) static GEN ComputeLift(GEN dataC) { - long order, i, av = avma; + long order, i; + gpmem_t av = avma; GEN cyc, surj, eltq, elt; order = itos((GEN)dataC[1]); cyc = (GEN)dataC[2]; surj = (GEN)dataC[3]; - eltq = FindEltofGroup(order, cyc); + eltq = EltsOfGroup(order, cyc); elt = cgetg(order + 1, t_VEC); @@ -121,72 +198,53 @@ ComputeLift(GEN dataC) return gerepileupto(av, elt); } -/* Let bnr1, bnr2 be such that mod(bnr2) | mod(bnr1), compute the - matrix of the surjective map Cl(bnr1) ->> Cl(bnr2) */ static GEN -GetSurjMat(GEN bnr1, GEN bnr2) +get_chic(GEN chi, GEN cyc) { - long nbg, i; - GEN gen, M; - - gen = gmael(bnr1, 5, 3); - nbg = lg(gen) - 1; - - M = cgetg(nbg + 1, t_MAT); - for (i = 1; i <= nbg; i++) - M[i] = (long)isprincipalray(bnr2, (GEN)gen[i]); - - return M; + long i, l = lg(chi); + GEN chic = cgetg(l, t_VEC); + for (i = 1; i < l; i++) + chic[i] = ldiv((GEN)chi[i], (GEN)cyc[i]); + return chic; } /* A character is given by a vector [(c_i), z, d, pm] such that chi(id) = z ^ sum(c_i * a_i) where a_i= log(id) on the generators of bnr - z = exp(2i * Pi / d) - pm = z as a polmod */ + z = exp(2i * Pi / d) */ static GEN -get_Char(GEN chi, long prec) +get_Char(GEN chic, long prec) { - GEN p2, d, _2ipi = gmul(gi, shiftr(mppi(prec), 1)); - p2 = cgetg(5, t_VEC); d = denom(chi); - p2[1] = lmul(d, chi); - if (egalii(d, gdeux)) - p2[2] = lstoi(-1); - else - p2[2] = lexp(gdiv(_2ipi, d), prec); - p2[3] = (long)d; - p2[4] = lmodulcp(polx[0], cyclo(itos(d), 0)); - return p2; + GEN d, C; + C = cgetg(4, t_VEC); + d = denom(chic); + C[1] = lmul(d, chic); + C[2] = (long)InitRU(d, prec); + C[3] = (long)d; return C; } -/* Let chi a character defined over bnr and primitif over bnrc, +/* Let chi a character defined over bnr and primitive over bnrc, compute the corresponding primitive character and the vectors of - prime ideals dividing bnr but not bnr. Returns NULL if bnr = bnrc */ + prime ideals dividing bnr but not bnrc. Returns NULL if bnr = bnrc */ static GEN GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec) { - long nbg, i, j, l, av = avma, nd; - GEN gen, cyc, U, chic, M, s, p1, cond, condc, p2, nf; + long nbg, i, j, l, nd; + gpmem_t av = avma; + GEN s, chic, cyc, U, M, p1, cond, condc, p2, nf; GEN prdiff, Mrc; - cond = gmael(bnr, 2, 1); + cond = gmael(bnr, 2, 1); condc = gmael(bnrc, 2, 1); if (gegal(cond, condc)) return NULL; - gen = gmael(bnr, 5, 3); - nbg = lg(gen) - 1; - cyc = gmael(bnr, 5, 2); + cyc = gmael(bnr, 5, 2); nbg = lg(cyc)-1; Mrc = diagonal(gmael(bnrc, 5, 2)); nf = gmael(bnr, 1, 7); - cond = (GEN)cond[1]; - condc = (GEN)condc[1]; - - M = GetSurjMat(bnr, bnrc); - l = lg((GEN)M[1]); - p1 = hnfall(concatsp(M, Mrc)); - U = (GEN)p1[2]; - + M = bnrGetSurj(bnr, bnrc); + U = (GEN)hnfall(concatsp(M, Mrc))[2]; + l = lg((GEN)M[1]); chic = cgetg(l, t_VEC); for (i = 1; i < l; i++) { @@ -194,11 +252,13 @@ GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec) for (j = 1; j <= nbg; j++) { p2 = gdiv((GEN)p1[j], (GEN)cyc[j]); - s = gadd(s,gmul(p2,(GEN)chi[j])); + s = gadd(s, gmul(p2,(GEN)chi[j])); } chic[i] = (long)s; } + cond = (GEN)cond[1]; + condc = (GEN)condc[1]; p2 = (GEN)idealfactor(nf, cond)[1]; l = lg(p2); @@ -214,27 +274,14 @@ GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec) return gerepileupto(av,p1); } -/* Let dataCR be a list of characters, compute the image of logelt */ static GEN -chiideal(GEN dataCR, GEN logelt, long flag) -{ - long j, l = lg(dataCR); - GEN rep = cgetg(l, t_VEC); - - for (j = 1; j < l; j++) - rep[j] = (long)ComputeImagebyChar(gmael(dataCR, j, 5), logelt, flag); - - return rep; -} - -static GEN GetDeg(GEN dataCR) { long i, l = lg(dataCR); GEN degs = cgetg(l, t_VECSMALL); for (i = 1; i < l; i++) - degs[i] = degpol(gmael4(dataCR, i, 5, 4, 1)); + degs[i] = itos(phi(gmael3(dataCR, i, 5, 3))); return degs; } @@ -244,6 +291,7 @@ GetDeg(GEN dataCR) static GEN AllStark(GEN data, GEN nf, long flag, long prec); static GEN InitChar0(GEN dataD, long prec); +static GEN QuadGetST(GEN dataCR, GEN vChar, long prec); /* Let A be a finite abelian group given by its relation and let C define a subgroup of A, compute the order of A / C, its structure and @@ -251,25 +299,18 @@ static GEN InitChar0(GEN dataD, long prec); static GEN InitQuotient0(GEN DA, GEN C) { - long i, l; - GEN MQ, MrC, H, snf, snf2, rep, p1; + GEN D, MQ, MrC, H, U, rep; - l = lg(DA)-1; H = gcmp0(C)? DA: C; MrC = gauss(H, DA); - snf = smith2(hnf(MrC)); - MQ = concatsp(gmul(H, (GEN)snf[1]), DA); - snf2 = smith2(hnf(MQ)); + (void)smithall(hnf(MrC), &U, NULL); + MQ = concatsp(gmul(H, U), DA); + D = smithall(hnf(MQ), &U, NULL); rep = cgetg(5, t_VEC); - - p1 = cgetg(l + 1, t_VEC); - for (i = 1; i <= l; i++) - p1[i] = lcopy(gcoeff((GEN)snf2[3], i, i)); - - rep[1] = (long)dethnf((GEN)snf2[3]); - rep[2] = (long)p1; - rep[3] = lcopy((GEN)snf2[1]); + rep[1] = (long)dethnf_i(D); + rep[2] = (long)mattodiagonal(D); + rep[3] = lcopy(U); rep[4] = lcopy(C); return rep; @@ -286,7 +327,7 @@ static GEN InitQuotient(GEN bnr, GEN C) { GEN Mrm, dataquo = cgetg(3, t_VEC); - long av; + gpmem_t av; dataquo[1] = lcopy(bnr); av = avma; Mrm = diagonal(gmael(bnr, 5, 2)); @@ -296,49 +337,46 @@ InitQuotient(GEN bnr, GEN C) } /* Let s: A -> B given by P, and let DA, DB be resp. the matrix of the - relations of A and B, let nbA, nbB be resp. the rank of A and B, - compute the kernel of s. If DA = 0 then A is free */ + relations of A and B, compute the kernel of s. If DA = 0 then A is free */ static GEN -ComputeKernel0(GEN P, GEN DA, GEN DB, long nbA, long nbB) +ComputeKernel0(GEN P, GEN DA, GEN DB) { - long rk, av = avma; - GEN herm, mask1, mask2, U; + gpmem_t av = avma; + long nbA = lg(DA)-1, rk; + GEN herm, U; herm = hnfall(concatsp(P, DB)); - rk = nbA + nbB + 1; - rk -= lg((GEN)herm[1]); /* two steps: bug in pgcc 1.1.3 inlining (IS) */ + rk = nbA + lg(DB) - lg(herm[1]); + U = (GEN)herm[2]; + U = vecextract_i(U, 1,rk); + U = rowextract_i(U, 1,nbA); - mask1 = subis(shifti(gun, nbA), 1); - mask2 = subis(shifti(gun, rk), 1); - - U = matextract((GEN)herm[2], mask1, mask2); - if (!gcmp0(DA)) U = concatsp(U, DA); return gerepileupto(av, hnf(U)); } /* Let m and n be two moduli such that n|m and let C be a congruence group modulo n, compute the corresponding congruence group modulo m - ie then kernel of the map Clk(m) ->> Clk(n)/C */ + ie the kernel of the map Clk(m) ->> Clk(n)/C */ static GEN ComputeKernel(GEN bnrm, GEN dataC) { - long av = avma, i, nbm, nbq; + long i, nbm; + gpmem_t av = avma; GEN bnrn, Mrm, genm, Mrq, mgq, P; bnrn = (GEN)dataC[1]; Mrm = diagonal(gmael(bnrm, 5, 2)); + Mrq = diagonal(gmael(dataC, 2, 2)); genm = gmael(bnrm, 5, 3); nbm = lg(genm) - 1; - Mrq = diagonal(gmael(dataC, 2, 2)); mgq = gmael(dataC, 2, 3); - nbq = lg(mgq) - 1; P = cgetg(nbm + 1, t_MAT); for (i = 1; i <= nbm; i++) P[i] = lmul(mgq, isprincipalray(bnrn, (GEN)genm[i])); - return gerepileupto(av, ComputeKernel0(P, Mrm, Mrq, nbm, nbq)); + return gerepileupto(av, ComputeKernel0(P, Mrm, Mrq)); } /* Let C a congruence group in bnr, compute its subgroups of index 2 as @@ -346,30 +384,39 @@ ComputeKernel(GEN bnrm, GEN dataC) static GEN ComputeIndex2Subgroup(GEN bnr, GEN C) { - long nb, i, l, av = avma; - GEN snf, Mr, U, CU, subgrp, rep, p1; + gpmem_t av = avma; + long nb, i; + GEN D, Mr, U, T, subgrp; - disable_dbg(0); + disable_dbg(0); Mr = diagonal(gmael(bnr, 5, 2)); - snf = smith2(gmul(ginv(C), Mr)); + D = smithall(gauss(C, Mr), &U, NULL); + T = gmul(C,ginv(U)); + subgrp = subgrouplist(D, gdeux); + nb = lg(subgrp) - 1; + setlg(subgrp, nb); /* skip Id which comes last */ + for (i = 1; i < nb; i++) + subgrp[i] = (long)hnf(concatsp(gmul(T, (GEN)subgrp[i]), Mr)); - U = ginv((GEN)snf[1]); - l = lg((GEN)snf[3]); + disable_dbg(-1); + return gerepilecopy(av, subgrp); +} - p1 = cgetg(l, t_VEC); +GEN +Order(GEN cyc, GEN x) +{ + gpmem_t av = avma; + long i, l = lg(cyc); + GEN c,o,f = gun; for (i = 1; i < l; i++) - p1[i] = mael3(snf, 3, i, i); - - subgrp = subgrouplist(p1, 2); - nb = lg(subgrp) - 1; CU = gmul(C,U); - - rep = cgetg(nb, t_VEC); - for (i = 1; i < nb; i++) /* skip Id which comes last */ - rep[i] = (long)hnf(concatsp(gmul(CU, (GEN)subgrp[i]), Mr)); - - disable_dbg(-1); - return gerepilecopy(av, rep); + { + o = (GEN)cyc[i]; + c = mppgcd(o, (GEN)x[i]); + if (!is_pm1(c)) o = diviiexact(o,c); + f = mpppcm(f, o); + } + return gerepileuptoint(av, f); } /* Let pr be a prime (pr may divide mod(bnr)), compute the indexes @@ -377,51 +424,47 @@ ComputeIndex2Subgroup(GEN bnr, GEN C) static GEN GetIndex(GEN pr, GEN bnr, GEN subgroup) { - long av = avma, v, lg, i; - GEN bnf, mod, mod0, mpr0, mpr, bnrpr, subpr, M, e, f, dtQ, p1; + long v, e, f; + gpmem_t av = avma; + GEN bnf, mod, mod0, bnrpr, subpr, M, dtQ, p1; GEN rep, cycpr, cycQ; bnf = (GEN)bnr[1]; mod = gmael(bnr, 2, 1); mod0 = (GEN)mod[1]; - /* Compute the part of mod coprime to pr */ v = idealval(bnf, mod0, pr); - mpr0 = idealdivexact(bnf, mod0, idealpow(bnf, pr, stoi(v))); - - mpr = cgetg(3, t_VEC); - mpr[1] = (long)mpr0; - mpr[2] = mod[2]; - if (gegal(mpr0, mod0)) + if (v == 0) { bnrpr = bnr; subpr = subgroup; + e = 1; } else { + GEN mpr = cgetg(3, t_VEC); + GEN mpr0 = idealdivpowprime(bnf, mod0, pr, stoi(v)); + mpr[1] = (long)mpr0; /* part of mod coprime to pr */ + mpr[2] = mod[2]; bnrpr = buchrayinitgen(bnf, mpr); cycpr = gmael(bnrpr, 5, 2); - M = GetSurjMat(bnr, bnrpr); + M = bnrGetSurj(bnr, bnrpr); M = gmul(M, subgroup); subpr = hnf(concatsp(M, diagonal(cycpr))); + /* e = #(bnr/subgroup) / #(bnrpr/subpr) */ + e = itos( divii(dethnf_i(subgroup), dethnf_i(subpr)) ); } - /* e = #(bnr/subgroup) / #(bnrpr/subpr) */ - e = gdiv(det(subgroup), det(subpr)); - /* f = order of [pr] in bnrpr/subpr */ dtQ = InitQuotient(bnrpr, subpr); p1 = isprincipalray(bnrpr, pr); p1 = gmul(gmael(dtQ, 2, 3), p1); cycQ = gmael(dtQ, 2, 2); - lg = lg(cycQ) - 1; - f = gun; - for (i = 1; i <= lg; i++) - f = glcm(f, gdiv((GEN)cycQ[i], ggcd((GEN)cycQ[i], (GEN)p1[i]))); + f = itos( Order(cycQ, p1) ); - rep = cgetg(3, t_VEC); - rep[1] = lcopy(e); - rep[2] = lcopy(f); + rep = cgetg(3, t_VECSMALL); + rep[1] = (long)e; + rep[2] = (long)f; return gerepileupto(av, rep); } @@ -432,7 +475,8 @@ GetIndex(GEN pr, GEN bnr, GEN subgroup) static GEN CplxModulus(GEN data, long *newprec, long prec) { - long av = avma, pr, dprec; + long pr, dprec; + gpmem_t av = avma, av2; GEN nf, cpl, pol, p1; nf = gmael3(data, 1, 1, 7); @@ -450,13 +494,14 @@ CplxModulus(GEN data, long *newprec, long prec) dprec = DEFAULTPREC; + av2 = avma; for (;;) { p1[5] = (long)InitChar0((GEN)data[3], dprec); - pol = AllStark(p1, nf, -1, dprec); + pol = gerepileupto(av2, AllStark(p1, nf, -1, dprec)); if (!gcmp0(leading_term(pol))) { - cpl = mpabs(poldisc0(pol, 0)); + cpl = gnorml2(gtovec(pol)); if (!gcmp0(cpl)) break; } pr = gexpo(pol)>>(TWOPOTBITS_IN_LONG+1); @@ -488,7 +533,7 @@ FindModulus(GEN dataC, long fl, long *newprec, long pr { long n, i, narch, nbp, maxnorm, minnorm, N, nbidnn, s, c, j, nbcand; long limnorm, first = 1, pr; - ulong av = avma, av1, av0; + gpmem_t av = avma, av1, av0; GEN bnr, rep, bnf, nf, f, arch, m, listid, idnormn, bnrm, ImC; GEN candD, D, bpr, indpr, sgp, p1, p2, rb; @@ -503,16 +548,16 @@ FindModulus(GEN dataC, long fl, long *newprec, long pr for (i = 1; i <= 5; i++) rep[i] = zero; /* if cpl < rb, it is not necessary to try another modulus */ - rb = powgi(gmul(gmael(bnf, 7, 3), det(f)), gmul2n(gmael(bnr, 5, 1), 3)); + rb = powgi(gmul((GEN)nf[3], det(f)), gmul2n(gmael(bnr, 5, 1), 3)); bpr = (GEN)idealfactor(nf, f)[1]; nbp = lg(bpr) - 1; - indpr = cgetg(nbp + 1,t_VEC); + indpr = cgetg(nbp + 1,t_VECSMALL); for (i = 1; i <= nbp; i++) { p1 = GetIndex((GEN)bpr[i], bnr, sgp); - indpr[i] = lmulii((GEN)p1[1], (GEN)p1[2]); + indpr[i] = p1[1] * p1[2]; } /* Initialization of the possible infinite part */ @@ -529,7 +574,7 @@ FindModulus(GEN dataC, long fl, long *newprec, long pr If we cannot find a suitable conductor of norm < limnorm, we stop */ maxnorm = 50; minnorm = 1; - limnorm = 200; + limnorm = 400; if (DEBUGLEVEL >= 2) fprintferr("Looking for a modulus of norm: "); @@ -566,69 +611,59 @@ FindModulus(GEN dataC, long fl, long *newprec, long pr bnrm = buchrayinitgen(bnf, m); p1 = conductor(bnrm, gzero, -1); disable_dbg(-1); - if (signe(p1)) - { - /* compute Im(C) in Clk(m)... */ - ImC = ComputeKernel(bnrm, dataC); + arch[N+1-s] = un; + if (!signe(p1)) continue; - /* ... and its subgroups of index 2 */ - candD = ComputeIndex2Subgroup(bnrm, ImC); - nbcand = lg(candD) - 1; - for (c = 1; c <= nbcand; c++) - { - /* check if m is the conductor */ - D = (GEN)candD[c]; - disable_dbg(0); - p1 = conductor(bnrm, D, -1); - disable_dbg(-1); - if (signe(p1)) - { - /* check the splitting of primes */ - for (j = 1; j <= nbp; j++) - { - p1 = GetIndex((GEN)bpr[j], bnrm, D); - p1 = mulii((GEN)p1[1], (GEN)p1[2]); - if (egalii(p1, (GEN)indpr[j])) break; /* no good */ - } - if (j > nbp) - { - p2 = cgetg(6, t_VEC); + /* compute Im(C) in Clk(m)... */ + ImC = ComputeKernel(bnrm, dataC); - p2[1] = lcopy(bnrm); - p2[2] = lcopy(D); - p2[3] = (long)InitQuotient((GEN)p2[1], (GEN)p2[2]); - p2[4] = (long)InitQuotient((GEN)p2[1], ImC); + /* ... and its subgroups of index 2 */ + candD = ComputeIndex2Subgroup(bnrm, ImC); + nbcand = lg(candD) - 1; + for (c = 1; c <= nbcand; c++) + { + /* check if m is the conductor */ + D = (GEN)candD[c]; + disable_dbg(0); + p1 = conductor(bnrm, D, -1); + disable_dbg(-1); + if (!signe(p1)) continue; - p1 = CplxModulus(p2, &pr, prec); + /* check the splitting of primes */ + for (j = 1; j <= nbp; j++) + { + p1 = GetIndex((GEN)bpr[j], bnrm, D); + if (p1[1] * p1[2] == indpr[j]) break; /* no good */ + } + if (j <= nbp) continue; - if (first == 1 || gcmp(p1, (GEN)rep[5]) < 0) - { - *newprec = pr; - for (j = 1; j <= 4; j++) rep[j] = p2[j]; - rep[5] = (long)p1; - } + p2 = cgetg(6, t_VEC); + p2[1] = (long)bnrm; + p2[2] = (long)D; + p2[3] = (long)InitQuotient(bnrm, D); + p2[4] = (long)InitQuotient(bnrm, ImC); - if (!fl || (gcmp(p1, rb) < 0)) - { - rep[5] = (long)InitChar0((GEN)rep[3], *newprec); - return gerepilecopy(av, rep); - } - if (DEBUGLEVEL >= 2) - fprintferr("Trying to find another modulus..."); - first--; - } - } - } - } - arch[N+1-s] = un; + p1 = CplxModulus(p2, &pr, prec); + + if (first == 1 || gcmp(p1, (GEN)rep[5]) < 0) + { + *newprec = pr; + for (j = 1; j <= 4; j++) rep[j] = p2[j]; + rep[5] = (long)p1; + } + + if (!fl || gcmp(p1, rb) < 0) goto END; /* OK */ + + if (DEBUGLEVEL>1) fprintferr("Trying to find another modulus..."); + first--; + } } if (first <= bnd) { - if (DEBUGLEVEL >= 2) - fprintferr("No, we're done!\nModulus = %Z and subgroup = %Z\n", + if (DEBUGLEVEL>1) + fprintferr("No, we're done!\nModulus = %Z and subgroup = %Z\n", gmael3(rep, 1, 2, 1), rep[2]); - rep[5] = (long)InitChar0((GEN)rep[3], *newprec); - return gerepilecopy(av, rep); + goto END; /* OK */ } } } @@ -640,40 +675,55 @@ FindModulus(GEN dataC, long fl, long *newprec, long pr if (maxnorm > limnorm) err(talker, "Cannot find a suitable modulus in FindModulus"); } +END: + rep[5] = (long)InitChar0((GEN)rep[3], *newprec); + return gerepilecopy(av, rep); } /********************************************************************/ /* 2nd part: compute W(X) */ /********************************************************************/ -/* compute W(chi) such that Ld(s,chi) = W(chi) Ld(1 - s, chi*), - if flag != 0 do not check the result */ +/* compute the list of W(chi) such that Ld(s,chi) = W(chi) Ld(1 - s, chi*), + * for all chi in LCHI. All chi have the same conductor (= cond(bnr)). + * if check == 0 do not check the result */ static GEN -ComputeArtinNumber(GEN datachi, long flag, long prec) +ComputeArtinNumber(GEN bnr, GEN LCHI, long check, long prec) { - long av = avma, av2, G, ms, j, i, nz, zcard, q, l, N, lim; - GEN chi, nc, dc, p1, cond0, cond1, elts, Msign, umod2, lambda, nf; - GEN sg, p2, chib, diff, vt, z, idg, mu, idh, zid, zstruc, zgen, zchi; - GEN classe, bnr, beta, s, tr, p3, den, muslambda, Pi, lp1, beta2; + long ic, i, j, nz, q, N, nChar = lg(LCHI)-1; + gpmem_t av = avma, av2, lim; + GEN sqrtnc, dc, cond, condZ, cond0, cond1, lambda, nf, T; + GEN cyc, *vN, *vB, diff, vt, idg, mu, idh, zid, *gen, z, *nchi; + GEN indW, W, CHI, classe, s0, *s, tr, den, muslambda, beta2, sarch; + CHI_t **lC; + GROUP_t G; - chi = (GEN)datachi[8]; - /* trivial case */ - if (cmpsi(2, (GEN)chi[3]) >= 0) return gun; + lC = (CHI_t**)new_chunk(nChar + 1); + indW = cgetg(nChar + 1, t_VECSMALL); + W = cgetg(nChar + 1, t_VEC); + for (ic = 0, i = 1; i <= nChar; i++) + { + CHI = (GEN)LCHI[i]; + if (cmpsi(2, (GEN)CHI[3]) >= 0) { W[i] = un; continue; } /* trivial case */ + ic++; indW[ic] = i; + lC[ic] = (CHI_t*)new_chunk(sizeof(CHI_t)); + init_CHI_C(lC[ic], CHI); + } + if (!ic) return W; + nChar = ic; - bnr = (GEN)datachi[3]; nf = gmael(bnr, 1, 7); diff = gmael(nf, 5, 5); - cond0 = gmael3(bnr, 2, 1, 1); - cond1 = gmael3(bnr, 2, 1, 2); - umod2 = gmodulcp(gun, gdeux); + T = gmael(nf,5,4); + cond = gmael(bnr, 2, 1); + cond0 = (GEN)cond[1]; condZ = gcoeff(cond0,1,1); + cond1 = (GEN)cond[2]; N = degpol(nf[1]); - Pi = mppi(prec); - G = - bit_accuracy(prec) >> 1; - nc = idealnorm(nf, cond0); + sqrtnc = gsqrt(idealnorm(nf, cond0), prec); dc = idealmul(nf, diff, cond0); den = idealnorm(nf, dc); - z = gexp(gdiv(gmul2n(gmul(gi, Pi), 1), den), prec); + z = InitRU(den, prec); q = 0; for (i = 1; i < lg(cond1); i++) @@ -681,160 +731,146 @@ ComputeArtinNumber(GEN datachi, long flag, long prec) /* compute a system of elements congru to 1 mod cond0 and giving all possible signatures for cond1 */ - p1 = zarchstar(nf, cond0, cond1, q); - elts = (GEN)p1[2]; - Msign = gmul((GEN)p1[3], umod2); - ms = lg(elts) - 1; + sarch = zarchstar(nf, cond0, cond1, q); /* find lambda in diff.cond such that gcd(lambda.(diff.cond)^-1,cond0) = 1 - and lambda >(cond1)> 0 */ + and lambda >> 0 at cond1 */ lambda = idealappr(nf, dc); - sg = zsigne(nf, lambda, cond1); - p2 = lift(gmul(Msign, sg)); - - for (j = 1; j <= ms; j++) - if (gcmp1((GEN)p2[j])) lambda = element_mul(nf, lambda, (GEN)elts[j]); - + lambda = set_sign_mod_idele(nf, NULL, lambda, cond,sarch); idg = idealdivexact(nf, lambda, dc); /* find mu in idg such that idh=(mu) / idg is coprime with cond0 and - mu >(cond1)> 0 */ + mu >> 0 at cond1 */ if (!gcmp1(gcoeff(idg, 1, 1))) { - p1 = idealfactor(nf, idg); - p2 = idealfactor(nf, cond0); + GEN p1 = idealfactor(nf, idg); + GEN p2 = idealfactor(nf, cond0); + p2[2] = (long)zerocol(lg(p2[1])-1); + p1 = concat_factor(p1,p2); - l = lg((GEN)p2[1]); - for (i = 1; i < l; i++) coeff(p2, i, 2) = zero; - - p1 = gtrans(concatsp(gtrans(p1), gtrans(p2))); mu = idealapprfact(nf, p1); - sg = zsigne(nf, mu, cond1); - p2 = lift(gmul(Msign, sg)); - - for (j = 1; j <= ms; j++) - if (gcmp1((GEN)p2[j])) mu = element_mul(nf, mu, (GEN)elts[j]); - + mu = set_sign_mod_idele(nf, NULL,mu, cond,sarch); idh = idealdivexact(nf, mu, idg); } else { mu = gun; - idh = gcopy(idg); + idh = idg; } - muslambda = element_div(nf, mu, lambda); + muslambda = gmul(den, element_div(nf, mu, lambda)); /* compute a system of generators of (Ok/cond)^* cond1-positive */ zid = zidealstarinitgen(nf, cond0); + cyc = gmael(zid, 2, 2); + gen = (GEN*)gmael(zid, 2, 3); + nz = lg(gen) - 1; - zcard = itos(gmael(zid, 2, 1)); - zstruc = gmael(zid, 2, 2); - zgen = gmael(zid, 2, 3); - nz = lg(zgen) - 1; + nchi = (GEN*)cgetg(nChar+1, t_VEC); + for (ic = 1; ic <= nChar; ic++) nchi[ic] = cgetg(nz + 1, t_VECSMALL); - zchi = cgetg(nz + 1, t_VEC); for (i = 1; i <= nz; i++) { - p1 = (GEN)zgen[i]; - sg = zsigne(nf, p1, cond1); - p2 = lift(gmul(Msign, sg)); - - for (j = 1; j <= ms; j++) - if (gcmp1((GEN)p2[j])) p1 = element_mul(nf, p1, (GEN)elts[j]); - - classe = isprincipalray(bnr, p1); - zchi[i] = (long)ComputeImagebyChar(chi, classe, 0); - zgen[i] = (long)p1; + if (is_bigint(cyc[i])) + err(talker,"conductor too large in ComputeArtinNumber"); + gen[i] = set_sign_mod_idele(nf, NULL,gen[i], cond,sarch); + classe = isprincipalray(bnr, gen[i]); + for (ic = 1; ic <= nChar; ic++) + nchi[ic][i] = (long)EvalChar_n(lC[ic], classe); } /* Sum chi(beta) * exp(2i * Pi * Tr(beta * mu / lambda) where beta runs through the classes of (Ok/cond0)^* and beta cond1-positive */ - p3 = cgetg(N + 1, t_COL); - for (i = 1; i <= N; i++) p3[i] = zero; + vt = cgetg(N + 1, t_VEC); /* Tr(w_i) */ + for (i = 1; i <= N; i++) vt[i] = coeff(T,i,1); - vt = cgetg(N + 1, t_VEC); - for (i = 1; i <= N; i++) - { - p3[i] = un; - vt[i] = ltrace(basistoalg(nf, p3)); - p3[i] = zero; - } + G.cyc = gtovecsmall(cyc); + G.r = nz; + G.j = vecsmall_const(nz, 0); - lp1 = NULL; - s = gzero; + vN = (GEN*)cgetg(nChar+1, t_VEC); + for (ic = 1; ic <= nChar; ic++) vN[ic] = vecsmall_const(nz, 0); av2 = avma; lim = stack_lim(av2, 1); + vB = (GEN*)cgetg(nz+1, t_VEC); + for (i=1; i<=nz; i++) vB[i] = gun; - for (i = 0; i < zcard; i++) + tr = gmod(gmul(vt, muslambda), den); /* for beta = 1 */ + s0 = powgi(z, tr); + s = (GEN*)cgetg(nChar+1, t_VEC); + for (ic = 1; ic <= nChar; ic++) s[ic] = s0; + + for (;;) { - p1 = NextEltofGroup(zstruc, nz, i); + if (! (i = NextElt(&G)) ) break; - /* we test if we can use the previous value - of beta / chib to compute the next one */ - /* FIXME: there should be a better way of doing that! */ - if (!lp1 || !gcmp1(gnorml2(gsub(p1, lp1)))) - { - beta = gun; - chib = gun; + vB[i] = FpV_red(element_mul(nf, vB[i], gen[i]), condZ); + for (j=1; jord); + for (j=1; jval[ind]; + s[ic] = gadd(s[ic], gmul(val, s0)); + } - beta2 = element_mul(nf, beta, muslambda); - tr = gmul(vt, beta2); - tr = gmod(gmul(tr, den), den); - - s = gadd(s, gmul(chib, powgi(z,tr))); - if (low_stack(lim, stack_lim(av2, 1))) { - GEN *gptr[5]; - gptr[0] = &s; gptr[1] = &lp1; gptr[2] = β gptr[3] = &chib; + GEN *gptr[2]; gptr[0] = (GEN*)&s; gptr[1] = (GEN*)&vB; if (DEBUGMEM > 1) err(warnmem,"ComputeArtinNumber"); - gerepilemany(av2, gptr, 4); + gerepilemany(av2, gptr, 2); } } classe = isprincipalray(bnr, idh); - s = gmul(s, ComputeImagebyChar(chi, classe, 0)); - s = gdiv(s, gsqrt(nc, prec)); + z = gpowgs(gneg_i(gi),q); + for (ic = 1; ic <= nChar; ic++) + { + s0 = gmul(s[ic], EvalChar(lC[ic], classe)); + s0 = gdiv(s0, sqrtnc); + if (check && - expo(subrs(gnorm(s0), 1)) < bit_accuracy(prec) >> 1) + err(bugparier, "ComputeArtinNumber"); + W[indW[ic]] = lmul(s0, z); + } + return gerepilecopy(av, W); +} - p1 = gsubgs(gabs(s, prec), 1); +static GEN +ComputeAllArtinNumbers(GEN dataCR, GEN vChar, int check, long prec) +{ + const long cl = lg(dataCR) - 1; + long ncond = lg(vChar)-1; + long jc, k; + GEN W = cgetg(cl+1,t_VEC), WbyCond, LCHI; - i = expo(p1); - if ((i > G) && !flag) - err(bugparier, "ComputeArtinNumber"); + for (jc = 1; jc <= ncond; jc++) + { + const GEN LChar = (GEN)vChar[jc]; + const long nChar = lg(LChar)-1; + GEN *ldata = (GEN*)vecextract_p(dataCR, LChar); + GEN dtcr = ldata[1], bnr = (GEN)dtcr[3]; - return gerepileupto(av, gmul(s, gpowgs(gneg_i(gi),q))); + if (DEBUGLEVEL>1) + fprintferr("* Root Number: cond. no %ld/%ld (%ld chars)\n", + jc,ncond,nChar); + LCHI = cgetg(nChar + 1, t_VEC); + for (k = 1; k <= nChar; k++) LCHI[k] = ldata[k][8]; + WbyCond = ComputeArtinNumber(bnr, LCHI, check, prec); + for (k = 1; k <= nChar; k++) W[LChar[k]] = WbyCond[k]; + } + return W; } /* compute the constant W of the functional equation of @@ -842,64 +878,38 @@ ComputeArtinNumber(GEN datachi, long flag, long prec) GEN bnrrootnumber(GEN bnr, GEN chi, long flag, long prec) { - long av = avma, l, i; - GEN cond, condc, bnrc, chic, cyc, d, p1, p2, dtcr, Pi; + long l; + gpmem_t av = avma; + GEN cond, condc, bnrc, CHI; - if ((flag < 0) || (flag > 1)) - err(flagerr,"bnrrootnumber"); + if (flag < 0 || flag > 1) err(flagerr,"bnrrootnumber"); checkbnr(bnr); - cond = gmael(bnr, 2, 1); l = lg(gmael(bnr, 5, 2)); - Pi = mppi(prec); - if ((typ(chi) != t_VEC) || (lg(chi) != l)) + if (typ(chi) != t_VEC || lg(chi) != l) err(talker, "incorrect character in bnrrootnumber"); - if (!flag) + if (flag) condc = cond; + else { condc = bnrconductorofchar(bnr, chi); if (gegal(cond, condc)) flag = 1; } - else condc = cond; if (flag) + { + GEN cyc = gmael(bnr, 5, 2); bnrc = bnr; + CHI = get_Char(get_chic(chi,cyc), prec); + } else + { bnrc = buchrayinitgen((GEN)bnr[1], condc); - - chic = cgetg(l, t_VEC); - cyc = gmael(bnr, 5, 2); - for (i = 1; i < l; i++) - chic[i] = ldiv((GEN)chi[i], (GEN)cyc[i]); - d = denom(chic); - - p2 = cgetg(4, t_VEC); - p2[1] = lmul(d, chic); - if (egalii(d, gdeux)) - p2[2] = lstoi(-1); - else - p2[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), d), prec); - p2[3] = (long)d; - - dtcr = cgetg(9, t_VEC); - dtcr[1] = (long)chi; - dtcr[2] = zero; - dtcr[3] = (long)bnrc; - dtcr[4] = (long)bnr; - dtcr[5] = (long)p2; - dtcr[6] = zero; - dtcr[7] = (long)condc; - - p1 = GetPrimChar(chi, bnr, bnrc, prec); - - if (!p1) - dtcr[8] = (long)p2; - else - dtcr[8] = p1[1]; - - return gerepileupto(av, ComputeArtinNumber(dtcr, 0, prec)); + CHI = (GEN)GetPrimChar(chi, bnr, bnrc, prec)[1]; + } + return gerepilecopy(av, (GEN)ComputeArtinNumber(bnrc, _vec(CHI), 1, prec)[1]); } /********************************************************************/ @@ -909,7 +919,8 @@ bnrrootnumber(GEN bnr, GEN chi, long flag, long prec) static GEN LiftChar(GEN cyc, GEN Mat, GEN chi) { - long lm, l, i, j, av; + long lm, l, i, j; + gpmem_t av; GEN lchi, s; lm = lg(cyc) - 1; @@ -936,45 +947,57 @@ LiftChar(GEN cyc, GEN Mat, GEN chi) static GEN ComputeAChi(GEN dtcr, long flag, long prec) { - long l, i; - GEN p1, ray, r, A, rep, diff, chi, bnrc; + long l, i, r; + GEN p1, ray, A, rep, diff, chi, bnrc, nf; diff = (GEN)dtcr[6]; - bnrc = (GEN)dtcr[3]; + bnrc = (GEN)dtcr[3]; nf = checknf(bnrc); chi = (GEN)dtcr[8]; l = lg(diff) - 1; A = gun; - r = gzero; + r = 0; for (i = 1; i <= l; i++) { + GEN B; ray = isprincipalray(bnrc, (GEN)diff[i]); - p1 = ComputeImagebyChar(chi, ray, 0); + p1 = ComputeImagebyChar(chi, ray); if (flag) - A = gmul(A, gsub(gun, gdiv(p1, idealnorm((GEN)bnrc[1], (GEN)diff[i])))); - else + B = gsub(gun, gdiv(p1, idealnorm(nf, (GEN)diff[i]))); + else if (gcmp1(p1)) { - if (gcmp1(p1)) - { - r = addis(r, 1); - A = gmul(A, glog(idealnorm((GEN)bnrc[1], (GEN)diff[i]), prec)); - } - else - A = gmul(A, gsub(gun, p1)); + B = glog(idealnorm(nf, (GEN)diff[i]), prec); + r++; } + else + B = gsub(gun, p1); + A = gmul(A, B); } if (flag) return A; rep = cgetg(3, t_VEC); - rep[1] = (long)r; + rep[1] = lstoi(r); rep[2] = (long)A; - return rep; } +static GEN +_data9(GEN arch, long r1, long r2) +{ + GEN z = cgetg(5, t_VECSMALL); + long i, b, q = 0; + + for (i=1; i<=r1; i++) if (signe(arch[i])) q++; + z[1] = q; b = r1 - q; + z[2] = b; + z[3] = r2; + z[4] = max(b+r2+1, r2+q); + return z; +} + /* Given a list [chi, cond(chi)] of characters over Cl(bnr), compute a vector dataCR containing for each character: 1: chi @@ -992,17 +1015,16 @@ static GEN InitChar(GEN bnr, GEN listCR, long prec) { GEN bnf = checkbnf(bnr), nf = checknf(bnf); - GEN modul, dk, C, dataCR, chi, cond, Mr, chic, gr2, p1, p2, q; - long N, r1, r2, prec2, h, i, j, nbg, av = avma; + GEN modul, dk, C, dataCR, chi, cond, Mr, z1, p1; + long N, r1, r2, prec2, h, i, j; + gpmem_t av = avma; modul = gmael(bnr, 2, 1); Mr = gmael(bnr, 5, 2); - nbg = lg(Mr) - 1; dk = (GEN)nf[3]; N = degpol(nf[1]); - r1 = nf_get_r1(nf); - r2 = (N - r1) >> 1; gr2 = stoi(r2); - prec2 = ((prec - 2)<<1) + EXTRA_PREC; + nf_get_sign(nf, &r1,&r2); + prec2 = ((prec-2) << 1) + EXTRA_PREC; C = gmul2n(gsqrt(gdiv(absi(dk), gpowgs(mppi(prec2),N)), prec2), -r2); disable_dbg(0); @@ -1012,12 +1034,7 @@ InitChar(GEN bnr, GEN listCR, long prec) for (i = 1; i <= h ;i++) dataCR[i] = lgetg(10, t_VEC); - q = gnorml2((GEN)modul[2]); - p1 = cgetg(5, t_VEC); - p1[1] = (long)q; - p1[2] = lsubsi(r1, q); - p1[3] = (long)gr2; - p1[4] = lmax(addis((GEN)p1[2], r2+1), addsi(r2, q)); + z1 = _data9((GEN)modul[2],r1,r2); for (i = 1; i <= h; i++) { @@ -1039,7 +1056,7 @@ InitChar(GEN bnr, GEN listCR, long prec) data[3] = (long)bnr; data[6] = lgetg(1, t_VEC); data[7] = modul[1]; - data[9] = (long)p1; + data[9] = (long)z1; olddata = data; } @@ -1056,39 +1073,22 @@ InitChar(GEN bnr, GEN listCR, long prec) data[3] = olddata[3]; /* bnr(cond(chi)) */ } data[4] = (long)bnr; /* bnr(m) */ + data[5] = (long)get_Char(get_chic(chi,Mr),prec); /* associated to bnr(m) */ - chic = cgetg(nbg + 1, t_VEC); - for (j = 1; j <= nbg; j++) - chic[j] = ldiv((GEN)chi[j], (GEN)Mr[j]); - data[5] = (long)get_Char(chic,prec); /* char associated to bnr(m) */ - /* compute diff(chi) and the corresponding primitive character */ data[7] = cond[1]; - p2 = GetPrimChar(chi, bnr, (GEN)data[3], prec2); - if (p2) + p1 = GetPrimChar(chi, bnr, (GEN)data[3], prec2); + if (p1) { - data[6] = p2[2]; - data[8] = p2[1]; + data[6] = p1[2]; + data[8] = p1[1]; } else { data[6] = lgetg(1, t_VEC); data[8] = data[5]; } - - /* compute q and store [q, r1 - q, r2] */ - if (!olddata) - { - q = gnorml2((GEN)cond[2]); - p2 = cgetg(5, t_VEC); - p2[1] = (long)q; - p2[2] = lsubsi(r1, q); - p2[3] = (long)gr2; - p2[4] = lmax(addis((GEN)p2[2], r2+1), addsi(r2, q)); - data[9] = (long)p2; - } - else - data[9] = olddata[9]; + data[9] = olddata? olddata[9]: (long)_data9((GEN)cond[2],r1,r2); } disable_dbg(-1); @@ -1101,7 +1101,8 @@ static GEN InitChar0(GEN dataD, long prec) { GEN MrD, listCR, p1, chi, lchi, Surj, cond, bnr, p2, Mr, d, allCR; - long hD, h, nc, i, j, lD, nbg, tnc, av = avma; + long hD, h, nc, i, j, lD, tnc; + gpmem_t av = avma; Surj = gmael(dataD, 2, 3); MrD = gmael(dataD, 2, 2); @@ -1110,7 +1111,6 @@ InitChar0(GEN dataD, long prec) hD = itos(gmael(dataD, 2, 1)); h = hD >> 1; lD = lg(MrD)-1; - nbg = lg(Mr) - 1; disable_dbg(0); @@ -1119,7 +1119,7 @@ InitChar0(GEN dataD, long prec) allCR = cgetg(h + 1, t_VEC); /* all characters, including conjugates */ tnc = 1; - p1 = FindEltofGroup(hD, MrD); + p1 = EltsOfGroup(hD, MrD); for (i = 1; tnc <= h; i++) { @@ -1142,14 +1142,9 @@ InitChar0(GEN dataD, long prec) listCR[nc++] = (long)p2; allCR[tnc++] = (long)lchi; - /* compute the order of chi */ - p2 = cgetg(nbg + 1, t_VEC); - for (j = 1; j <= nbg; j++) - p2[j] = ldiv((GEN)lchi[j], (GEN)Mr[j]); - d = denom(p2); - /* if chi is not real, add its conjugate character to allCR */ - if (!gegal(d, gdeux)) + d = Order(Mr, lchi); + if (!egalii(d, gdeux)) allCR[tnc++] = (long)ConjChar(lchi, Mr); } @@ -1162,32 +1157,28 @@ InitChar0(GEN dataD, long prec) static GEN CharNewPrec(GEN dataCR, GEN nf, long prec) { - GEN dk, C, p1, Pi; - long av = avma, N, l, j, prec2; + GEN dk, C, p1; + long N, l, j, prec2; dk = (GEN)nf[3]; N = degpol((GEN)nf[1]); - l = lg(dataCR) - 1; prec2 = ((prec - 2)<<1) + EXTRA_PREC; - Pi = mppi(prec2); + C = mpsqrt(gdiv(dk, gpowgs(mppi(prec2), N))); - C = gsqrt(gdiv(dk, gpowgs(Pi, N)), prec2); - - for (j = 1; j <= l; j++) + l = lg(dataCR); + for (j = 1; j < l; j++) { - mael(dataCR, j, 2) = lmul(C, gsqrt(det(gmael(dataCR, j, 7)), prec2)); + GEN dtcr = (GEN)dataCR[j]; + dtcr[2] = lmul(C, gsqrt(dethnf_i((GEN)dtcr[7]), prec2)); - mael4(dataCR, j, 3, 1, 7) = lcopy(nf); + mael3(dtcr, 3, 1, 7) = (long)nf; - p1 = gmael(dataCR, j, 5); - p1[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), (GEN)p1[3]), prec); - - p1 = gmael(dataCR, j, 8); - p1[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), (GEN)p1[3]), prec); + p1 = (GEN)dtcr[5]; p1[2] = (long)InitRU((GEN)p1[3], prec2); + p1 = (GEN)dtcr[8]; p1[2] = (long)InitRU((GEN)p1[3], prec2); } - return gerepilecopy(av, dataCR); + return dataCR; } /********************************************************************/ @@ -1198,210 +1189,169 @@ CharNewPrec(GEN dataCR, GEN nf, long prec) /********************************************************************/ static void -_0toCoeff(int *rep, long dg) +_0toCoeff(int *rep, long deg) { long i; - for (i=0; i i - dg) c += c1[j] * c2[i-j]; - c3[i] = c; + if (j < deg && j > i - deg) c += c0[j] * c1[i-j]; + T[i] = c; } - c2 = c1; - for (i = 0; i < dg; i++) + for (i = 0; i < deg; i++) { - c = c3[i]; - for (j = 0; j < dg; j++) c += reduc[j][i] * c3[dg+j]; - c2[i] = c; + c = T[i]; + for (j = 0; j < deg; j++) c += reduc[j][i] * T[deg+j]; + c0[i] = c; } - /* cast necessary to work around a gcc-2.96 bug on alpha-linux (IS) */ - for ( ; i < (short)dg; i++) c2[i] = 0; - avma = av; } -/* c0 <- c0 + c2 * c1 */ +/* c0 <- c0 + c1 * c2 */ static void -AddMulCoeff(int *c0, int *c2, int* c1, int** reduc, long dg) +AddMulCoeff(int *c0, int *c1, int* c2, int** reduc, long deg) { - long av,i,j; - int c, *c3; + long i, j; + gpmem_t av; + int c, *t; - if (!c2) /* c2 == 1 */ + if (IsZero(c2,deg)) return; + if (!c1) /* c1 == 1 */ { - for (i = 0; i < dg; i++) c0[i] += c1[i]; + for (i = 0; i < deg; i++) c0[i] += c2[i]; return; } - for (i = 0; i <= dg; i++) - if (c1[i]) break; - if (i > dg) return; av = avma; - c3 = (int*)new_chunk(2*dg); - for (i = 0; i < 2*dg; i++) + t = (int*)new_chunk(2*deg); /* = c1 * c2, not reduced */ + for (i = 0; i < 2*deg; i++) { c = 0; for (j = 0; j <= i; j++) - if (j < dg && j > i - dg) c += c1[j] * c2[i-j]; - c3[i] = c; + if (j < deg && j > i - deg) c += c1[j] * c2[i-j]; + t[i] = c; } - for (i = 0; i < dg; i++) + for (i = 0; i < deg; i++) { - c = c0[i] + c3[i]; - for (j = 0; j < dg; j++) c += reduc[j][i] * c3[dg+j]; - c0[i] = c; + c = t[i]; + for (j = 0; j < deg; j++) c += reduc[j][i] * t[deg+j]; + c0[i] += c; } - avma = av; } -/* returns 0 if c is zero, 1 otherwise. */ -static long -IsZero(int* c, long dg) -{ - long i; - - for (i = 0; i < dg; i++) - if (c[i]) return 0; - return 1; -} - /* evaluate the coeff. No Garbage collector */ static GEN -EvalCoeff(GEN z, int* c, long dg) +EvalCoeff(GEN z, int* c, long deg) { long i,j; GEN e, r; + if (!c) return gzero; #if 0 /* standard Horner */ - e = stoi(c[dg - 1]); - for (i = dg - 2; i >= 0; i--) + e = stoi(c[deg - 1]); + for (i = deg - 2; i >= 0; i--) e = gadd(stoi(c[i]), gmul(z, e)); #else /* specific attention to sparse polynomials */ e = NULL; - for (i = dg-1; i >=0; i=j-1) + for (i = deg-1; i >=0; i=j-1) { for (j=i; c[j] == 0; j--) if (j==0) @@ -1422,297 +1372,265 @@ EvalCoeff(GEN z, int* c, long dg) return e; } -/* copy the n * m array matan */ +/* copy the n * (m+1) array matan */ static void -CopyCoeff(int*** a, int*** a2, long n, long m, GEN degs) +CopyCoeff(int** a, int** a2, long n, long m) { - long i,j,k; + long i,j; for (i = 1; i <= n; i++) { - long dg = degs[i]; - int **b = a[i], **b2 = a2[i]; - for (j = 0; j <= m; j++) - { - int *c = b[j], *c2 = b2[j]; - for (k = 0; k < dg; k++) c2[k] = c[k]; - } + int *b = a[i], *b2 = a2[i]; + for (j = 0; j < m; j++) b2[j] = b[j]; } - return; } -/* initialize the data for GetRay */ -static GEN -InitGetRay(GEN bnr, long nmax) +/* return q*p if <= n. Beware overflow */ +static long +next_pow(long q, long p, long n) { - long bd, i, j, l; - GEN listid, listcl, id, rep, bnf, cond; - - bnf = (GEN)bnr[1]; - cond = gmael3(bnr, 2, 1, 1); - - if (nmax < 1000) return NULL; - - rep = cgetg(4, t_VEC); - - disable_dbg(0); - bd = min(1000, nmax / 50); - listid = ideallist(bnf, bd); - disable_dbg(-1); - - listcl = cgetg(bd + 1, t_VEC); - for (i = 1; i <= bd; i++) - { - l = lg((GEN)listid[i]) - 1; - listcl[i] = lgetg(l + 1, t_VEC); - - for (j = 1; j <= l; j++) - { - id = gmael(listid, i, j); - if (gcmp1(gcoeff(idealadd(bnf, id, cond), 1, 1))) - mael(listcl, i, j) = (long)isprincipalray(bnr, id); - } - } - - if (DEBUGLEVEL) msgtimer("InitGetRay"); - - rep[1] = (long)listid; - rep[2] = (long)listcl; - rep[3] = nf_get_r2((GEN)bnf[7])? 0: un; /* != 0 iff nf is totally real */ - return rep; + const GEN x = muluu((ulong)q, (ulong)p); + const ulong qp = (ulong)x[2]; + return (lgefint(x) > 3 || qp > (ulong)n)? 0: qp; } -/* compute the class of the prime ideal pr in cl(bnr) using dataray */ -static GEN -GetRay(GEN bnr, GEN dataray, GEN pr, long prec) +static void +an_AddMul(int **an,int **an2, long np, long n, long deg, GEN chi, int **reduc) { - long av = avma, N, n, bd, c; - GEN id, tid, t2, u, alpha, p1, cl, listid, listcl, nf, cond; + GEN chi2 = chi; + long q, qk, k; + int *c, *c2 = (int*)new_chunk(deg); - if (!dataray) - return isprincipalray(bnr, pr); - - listid = (GEN)dataray[1]; - listcl = (GEN)dataray[2]; - cond = gmael3(bnr, 2, 1, 1); - bd = lg(listid) - 1; - nf = gmael(bnr, 1, 7); - N = degpol(nf[1]); - - if (dataray[3]) - t2 = gmael(nf, 5, 4); - else - t2 = gmael(nf, 5, 3); - - id = prime_to_ideal(nf, pr); - tid = qf_base_change(t2, id, 1); - - if (dataray[3]) - u = lllgramint(tid); - else - u = lllgramintern(tid,100,1,prec); - - if (!u) return gerepileupto(av, isprincipalray(bnr, id)); - - c = 1; alpha = NULL; - for (c=1; c<=N; c++) + CopyCoeff(an, an2, n/np, deg); + for (q=np;;) { - p1 = gmul(id, (GEN)u[c]); - if (gcmp1(gcoeff(idealadd(nf, p1, cond), 1, 1))) { alpha = p1; break; } - } - if (!alpha) - return gerepileupto(av, isprincipalray(bnr, pr)); + if (gcmp1(chi2)) c = NULL; else { Polmod2Coeff(c2, chi2, deg); c = c2; } + for(k = 1, qk = q; qk <= n; k++, qk += q) + AddMulCoeff(an[qk], c, an2[k], reduc, deg); + if (! (q = next_pow(q,np, n)) ) break; - id = idealdivexact(nf, alpha, id); - - n = itos(det(id)); - if (n > bd) - cl = isprincipalray(bnr, id); - else - { - cl = NULL; - c = 1; - p1 = (GEN)listid[n]; - while (!cl) - { - if (gegal((GEN)p1[c], id)) - cl = gmael(listcl, n, c); - c++; - } + chi2 = gmul(chi2, chi); } - - return gerepileupto(av, gsub(isprincipalray(bnr, alpha), cl)); } /* correct the coefficients an(chi) according with diff(chi) in place */ static void -CorrectCoeff(GEN dtcr, int** an, int** reduc, long nmax, long dg) +CorrectCoeff(GEN dtcr, int** an, int** reduc, long n, long deg) { - long lg, av1, j, p, q, limk, k, l, av = avma; - int ***an2, **an1, *c, *c2; - GEN chi, bnrc, diff, ray, ki, ki2, pr, degs; + gpmem_t av = avma; + long lg, j, np; + gpmem_t av1; + int **an2; + GEN bnrc, diff, ray, chi, CHI, pr; + CHI_t C; - chi = (GEN)dtcr[8]; - bnrc = (GEN)dtcr[3]; diff = (GEN)dtcr[6]; lg = lg(diff) - 1; if (!lg) return; - if (DEBUGLEVEL > 2) fprintferr("diff(chi) = %Z", diff); + bnrc = (GEN)dtcr[3]; + CHI = (GEN)dtcr[8]; init_CHI_alg(&C, CHI); - degs = cgetg(2, t_VECSMALL); degs[1] = dg; - an2 = InitMatAn(1, nmax, degs, 0); an1 = an2[1]; - c = (int*)new_chunk(dg); - av1 = avma; + if (DEBUGLEVEL > 2) fprintferr("diff(CHI) = %Z", diff); + an2 = InitMatAn(n, deg, 0); + av1 = avma; for (j = 1; j <= lg; j++) { - for (k = 0; k <= nmax; k++) - for (l = 0; l < dg; l++) an1[k][l] = an[k][l]; - pr = (GEN)diff[j]; + np = itos(powgi((GEN)pr[1], (GEN)pr[4])); + ray = isprincipalray(bnrc, pr); - ki = ComputeImagebyChar(chi, ray, 1); - ki2 = gcopy(ki); + chi = EvalChar(&C,ray); - q = p = itos(powgi((GEN)pr[1], (GEN)pr[4])); - limk = nmax / q; - - while (q <= nmax) - { - if (gcmp1(ki2)) c2 = NULL; else { Polmod2Coeff(c,ki2, dg); c2 = c; } - for(k = 1; k <= limk; k++) - AddMulCoeff(an[k*q], c2, an1[k], reduc, dg); - - q *= p; limk /= p; - ki2 = gmul(ki2, ki); - } + an_AddMul(an,an2,np,n,deg,chi,reduc); avma = av1; } - FreeMat(an2); avma = av; + FreeMat(an2, n); avma = av; } /* compute the coefficients an in the general case */ -static int*** -ComputeCoeff(GEN dataCR, long nmax, long prec) +static int** +ComputeCoeff(GEN dtcr, LISTray *R, long n, long deg) { - long cl, i, j, av = avma, av2, np, q, q1, limk, k, id, cpt = 10, dg, Bq; - int ***matan, ***reduc, ***matan2, *c2; - GEN c, degs, tabprem, bnf, pr, cond, ray, ki, ki2, prime, npg, bnr, dataray; - byteptr dp = diffptr; + gpmem_t av = avma, av2; + long i, l, np; + int **an, **reduc, **an2; + GEN L, CHI, ray, chi; + CHI_t C; - bnr = gmael(dataCR, 1, 4); - bnf = (GEN)bnr[1]; - cond = gmael3(bnr, 2, 1, 1); - cl = lg(dataCR) - 1; + CHI = (GEN)dtcr[5]; init_CHI_alg(&C, CHI); - dataray = InitGetRay(bnr, nmax); + an = InitMatAn(n, deg, 0); + an2 = InitMatAn(n, deg, 0); + reduc = InitReduction(CHI, deg); + av2 = avma; - degs = GetDeg(dataCR); - matan = InitMatAn(cl, nmax, degs, 0); - matan2 = InitMatAn(cl, nmax, degs, 0); - reduc = InitReduction(dataCR, degs); - c = cgetg(cl + 1, t_VEC); - for (i = 1; i <= cl; i++) - c[i] = (long)new_chunk(degs[i]); + L = R->L1; l = lg(L); + for (i=1; iL1ray[i]; + chi = EvalChar(&C, ray); + an_AddMul(an,an2,np,n,deg,chi,reduc); + } + FreeMat(an2, n); - if (DEBUGLEVEL > 1) fprintferr("p = "); + CorrectCoeff(dtcr, an, reduc, n, deg); + FreeMat(reduc, deg-1); + avma = av; return an; +} - prime = stoi(2); dp++; - av2 = avma; - while (*dp && (prime[2] <= nmax)) - { - tabprem = primedec(bnf, prime); - for (j = 1; j < lg(tabprem); j++) - { - pr = (GEN)tabprem[j]; - npg = powgi((GEN)pr[1], (GEN)pr[4]); - if (is_bigint(npg) || (np=npg[2]) > nmax - || idealval(bnf, cond, pr)) continue; +/********************************************************************/ +/* 5th part: compute L-functions at s=1 */ +/********************************************************************/ +static void +deg11(LISTray *R, long p, GEN bnr, GEN pr) { + GEN z = isprincipalray(bnr, pr); + appendL(R->L1, (GEN)p); + appendL((GEN)R->L1ray, z); +} +static void +deg12(LISTray *R, long p, GEN bnr, GEN Lpr) { + GEN z = isprincipalray(bnr, (GEN)Lpr[1]); + appendL(R->L11, (GEN)p); + appendL((GEN)R->L11ray, z); +} +static void +deg0(LISTray *R, long p) { + appendL(R->L0, (GEN)p); +} +static void +deg2(LISTray *R, long p) { + appendL(R->L2, (GEN)p); +} - CopyCoeff(matan, matan2, cl, nmax, degs); - ray = GetRay(bnr, dataray, pr, prec); - ki = chiideal(dataCR, ray, 1); - ki2 = dummycopy(ki); +/* pi(x) <= ?? */ +static long +PiBound(long x) +{ + double lx = log((double)x); + return 1 + (long) (x/lx * (1 + 3/(2*lx))); +} - Bq = nmax/np; - for (q1 = 1; q1 <= Bq; q1 *= np) +static void +InitPrimesQuad(GEN bnr, long nmax, LISTray *R) +{ + gpmem_t av = avma; + GEN bnf = (GEN)bnr[1], cond = gmael3(bnr,2,1,1); + long p,i,l, condZ = itos(gcoeff(cond,1,1)), contZ = itos(content(cond)); + GEN prime, pr, nf = checknf(bnf), dk = (GEN)nf[3]; + byteptr d = diffptr + 1; + GEN *gptr[7]; + + l = 1 + PiBound(nmax); + R->L0 = cget1(l, t_VECSMALL); + R->L2 = cget1(l, t_VECSMALL); R->condZ = condZ; + R->L1 = cget1(l, t_VECSMALL); R->L1ray = (GEN*)cget1(l, t_VEC); + R->L11 = cget1(l, t_VECSMALL); R->L11ray = (GEN*)cget1(l, t_VEC); + prime = stoi(2); + for (p = 2; p <= nmax; prime[2] = p) { + switch (krogs(dk, p)) + { + case -1: /* inert */ + if (condZ % p == 0) deg0(R,p); else deg2(R,p); + break; + case 1: /* split */ + pr = primedec(nf, prime); + if (condZ % p != 0) deg12(R, p, bnr, pr); + else if (contZ % p == 0) deg0(R,p); + else { - q = q1*np; - limk = nmax / q; - for (id = 1; id <= cl; id++) - { - dg = degs[id]; - if (gcmp1((GEN)ki2[id])) - c2 = NULL; - else - { - c2 = (int*)c[id]; - Polmod2Coeff(c2, (GEN)ki2[id], dg); - } - for (k = 1; k <= limk; k++) - AddMulCoeff(matan[id][k*q], c2, matan2[id][k], reduc[id], dg); - ki2[id] = lmul((GEN)ki2[id], (GEN)ki[id]); - } + pr = idealval(nf, cond, (GEN)pr[1])? (GEN)pr[2]: (GEN)pr[1]; + deg11(R, p, bnr, pr); } + break; + default: /* ramified */ + if (condZ % p == 0) deg0(R,p); + else + { + pr = (GEN)primedec(nf, prime)[1]; + deg11(R, p, bnr, pr); + } + break; } - avma = av2; - prime[2] += (*dp++); - if (!*dp) err(primer1); + NEXT_PRIME_VIADIFF(p,d); + } + /* precompute isprincipalray(x), x in Z */ + R->rayZ = (GEN*)cgetg(condZ, t_VEC); + for (i=1; irayZ[i] = (cgcd(i,condZ) == 1)? isprincipalray(bnr, stoi(i)): gzero; - if (DEBUGLEVEL > 1 && prime[2] > cpt) - { fprintferr("%ld ", prime[2]); cpt += 10; } - } - if (DEBUGLEVEL > 1) fprintferr("\n"); + gptr[0] = &(R->L0); + gptr[1] = &(R->L2); gptr[2] = (GEN*)&(R->rayZ); + gptr[3] = &(R->L1); gptr[5] = (GEN*)&(R->L1ray); + gptr[4] = &(R->L11); gptr[6] = (GEN*)&(R->L11ray); + gerepilemany(av,gptr,7); +} - for (i = 1; i <= cl; i++) - CorrectCoeff((GEN)dataCR[i], matan[i], reduc[i], nmax, degs[i]); +static void +InitPrimes(GEN bnr, long nmax, LISTray *R) +{ + gpmem_t av = avma; + GEN bnf = (GEN)bnr[1], cond = gmael3(bnr,2,1,1); + long np,p,j, condZ = itos(gcoeff(cond,1,1)); + GEN Npr, tabpr, prime, pr, nf = checknf(bnf); + byteptr d = diffptr + 1; + GEN *gptr[7]; - FreeMat(matan2); FreeMat(reduc); - avma = av; return matan; + R->condZ = condZ; + R->L1 = cget1(nmax, t_VECSMALL); + R->L1ray = (GEN*)cget1(nmax, t_VEC); + prime = stoi(2); + for (p = 2; p <= nmax; prime[2] = p) + { + tabpr = primedec(nf, prime); + for (j = 1; j < lg(tabpr); j++) + { + pr = (GEN)tabpr[j]; + Npr = powgi((GEN)pr[1], (GEN)pr[4]); + if (is_bigint(Npr) || (np = Npr[2]) > nmax) continue; + if (condZ % p == 0 && idealval(nf, cond, pr)) continue; + + appendL(R->L1, (GEN)np); + appendL((GEN)R->L1ray, isprincipalray(bnr, pr)); + } + NEXT_PRIME_VIADIFF(p,d); + } + gptr[0] = &(R->L1); gptr[1] = (GEN*)&(R->L1ray); + gerepilemany(av,gptr,2); } -/********************************************************************/ -/* 5th part: compute L-functions at s=1 */ -/********************************************************************/ - -/* if flag != 0, prec means decimal digits */ static GEN -get_limx(long r1, long r2, long prec, GEN *pteps, long flag) +get_limx(long r1, long r2, long prec, GEN *pteps) { - GEN eps, a, r, c0, A0, limx, Pi = mppi(prec), N, p1; + GEN eps, p1, a, c0, A0, limx, Pi2 = gmul2n(mppi(DEFAULTPREC), 1); + long r, N; - N = addss(r1, 2*r2); - a = gmul(gpow(gdeux, gsubgs(gdiv(stoi(r1), N), 1), DEFAULTPREC), N); - r = addss(r1, r2); + N = r1 + 2*r2; r = r1 + r2; + a = gmulgs(gpow(gdeux, gdiv(stoi(-2*r2), stoi(N)), DEFAULTPREC), N); - if (flag) - *pteps = eps = gmul2n(gpowgs(dbltor(10.), -prec), -1); - else - *pteps = eps = gmul2n(gpowgs(dbltor(10.), (long)(-(prec-2) / pariK1)), -1); - - c0 = gpow(gmul2n(Pi, 1), gdiv(subis(r, 1), gdeux), DEFAULTPREC); - c0 = gmul(c0, gdiv(gdeux, N)); - c0 = gmul(c0, gpow(gdeux, gmul(gdiv(stoi(r1), gdeux), - gsubsg(1, gdiv(addis(r, 1), N))), + eps = real2n(-bit_accuracy(prec), DEFAULTPREC); + c0 = gpowgs(mpsqrt(Pi2), r-1); + c0 = gdivgs(gmul2n(c0,1), N); + c0 = gmul(c0, gpow(gdeux, gdiv(stoi(r1 * (r2-1)), stoi(2*N)), DEFAULTPREC)); - A0 = glog(gdiv(gmul2n(c0, 1), eps), DEFAULTPREC); + A0 = glog(gdiv(gmul2n(c0,1), eps), DEFAULTPREC); - limx = gpow(gdiv(a, A0), gdiv(N, gdeux), DEFAULTPREC); + limx = gpow(gdiv(a, A0), gdiv(stoi(N), gdeux), DEFAULTPREC); p1 = gsub(glog(A0, DEFAULTPREC), glog(a, DEFAULTPREC)); - p1 = gmul(gmul(p1, N), addis(r, 1)); - p1 = gdiv(p1, gmul2n(gadd(gmul2n(A0, 1), addis(r, 1)), 1)); - limx = gmul(limx, gaddgs(p1, 1)); + p1 = gmulgs(p1, N*(r+1)); + p1 = gdiv(p1, gaddgs(gmul2n(A0, 2), 2*(r+1))); - return limx; + if (pteps) *pteps = eps; + return gmul(limx, gaddgs(p1, 1)); } static long -GetBoundN0(GEN C, long r1, long r2, long prec, long flag) +GetBoundN0(GEN C, long r1, long r2, long prec) { - long av = avma; - GEN eps, limx = get_limx(r1, r2, prec, &eps, flag); + gpmem_t av = avma; + GEN limx = get_limx(r1, r2, prec, NULL); limx = gfloor(gdiv(C, limx)); if (is_bigint(limx)) @@ -1724,12 +1642,12 @@ GetBoundN0(GEN C, long r1, long r2, long prec, long static long GetBoundi0(long r1, long r2, long prec) { - long av = avma, imin, i0, itest; - GEN ftest, borneps, eps, limx = get_limx(r1, r2, prec, &eps, 0); - GEN Pi = mppi(DEFAULTPREC); + long imin, i0, itest; + gpmem_t av = avma; + GEN ftest, borneps, eps, limx = get_limx(r1, r2, prec, &eps); + GEN sqrtPi = mpsqrt(mppi(DEFAULTPREC)); - borneps = gmul(gmul2n(gun, r2), gpow(Pi, gdiv(subss(r2, 3), gdeux), - DEFAULTPREC)); + borneps = gmul(gmul2n(gun, r2), gpowgs(sqrtPi, r2-3)); borneps = gdiv(gmul(borneps, gpowgs(stoi(5), r1)), eps); borneps = gdiv(borneps, gsqrt(limx, DEFAULTPREC)); @@ -1740,376 +1658,418 @@ GetBoundi0(long r1, long r2, long prec) itest = (i0 + imin) >> 1; ftest = gpowgs(limx, itest); - ftest = gmul(ftest, gpowgs(mpfactr(itest / 2, DEFAULTPREC), r1)); - ftest = gmul(ftest, gpowgs(mpfactr(itest, DEFAULTPREC), r2)); + ftest = gmul(ftest, gpowgs(mpfactr(itest/2, DEFAULTPREC), r1)); + ftest = gmul(ftest, gpowgs(mpfactr(itest , DEFAULTPREC), r2)); - if(gcmp(ftest, borneps) >= 0) + if (gcmp(ftest, borneps) >= 0) i0 = itest; else imin = itest; } avma = av; - return (i0 / 2) * 2; + return i0 &= ~1; /* make it even */ } +static GEN /* cf polcoeff */ +_sercoeff(GEN x, long n) +{ + long i = n - valp(x); + return (i < 0)? gzero: (GEN)x[i+2]; +} + +static void +affect_coeff(GEN q, long n, GEN y) +{ + GEN x = _sercoeff(q,-n); + if (x == gzero) y[n] = zero; else gaffect(x, (GEN)y[n]); +} + +typedef struct { + GEN c1, *aij, *bij, *powracpi, *cS, *cT; + long i0, a,b,c, r, rc1, rc2; +} ST_t; + /* compute the principal part at the integers s = 0, -1, -2, ..., -i0 of Gamma((s+1)/2)^a Gamma(s/2)^b Gamma(s)^c / (s - z) with z = 0 and 1 */ /* NOTE: this is surely not the best way to do this, but it's fast enough! */ -static GEN -ppgamma(long a, long b, long c, long i0, long prec) +static void +ppgamma(ST_t *T, long prec) { - GEN cst, gamun, gamdm, an, bn, cn_evn, cn_odd, x, x2, aij, p1, cf, p2; - long i, j, r, av = avma; + GEN eul, gam,gamun,gamdm, an,bn,cn_evn,cn_odd, x,x2,X,Y, cf, sqpi; + GEN p1, p2, *aij, *bij; + long a = T->a; + long b = T->b; + long c = T->c, r = T->r, i0 = T->i0; + long i,j, s,t; + gpmem_t av; - r = max(b + c + 1, a + c); - - aij = cgetg(i0 + 1, t_VEC); + aij = (GEN*)cgetg(i0+1, t_VEC); + bij = (GEN*)cgetg(i0+1, t_VEC); for (i = 1; i <= i0; i++) { - aij[i] = lgetg(3, t_VEC); - mael(aij, i, 1) = lgetg(r + 1, t_VEC); - mael(aij, i, 2) = lgetg(r + 1, t_VEC); + aij[i] = p1 = cgetg(r+1, t_VEC); + bij[i] = p2 = cgetg(r+1, t_VEC); + for (j=1; j<=r; j++) { p1[j] = lgetr(prec); p2[j] = lgetr(prec); } } + av = avma; x = polx[0]; - x2 = gmul2n(x, -1); + x2 = gmul2n(x, -1); /* x/2 */ + eul = mpeuler(prec); + sqpi= mpsqrt(mppi(prec)); /* Gamma(1/2) */ - /* Euler gamma constant, values of Riemann zeta functions at - positive integers */ - cst = cgetg(r + 2, t_VEC); - cst[1] = (long)mpeuler(prec); - for (i = 2; i <= r + 1; i++) - cst[i] = (long)gzeta(stoi(i), prec); - - /* the expansion of log(Gamma(s)) at s = 1 */ - gamun = cgetg(r + 2, t_SER); + /* expansion of log(Gamma(u)) at u = 1 */ + gamun = cgetg(r+3, t_SER); gamun[1] = evalsigne(1) | evalvalp(0) | evalvarn(0); gamun[2] = zero; - for (i = 1; i <= r; i++) + gamun[3] = lneg(eul); + for (i = 2; i <= r; i++) { - gamun[i + 2] = ldivgs((GEN)cst[i], i); - if (i%2) gamun[i + 2] = lneg((GEN)gamun[i + 2]); + p1 = gdivgs(gzeta(stoi(i),prec), i); + if (odd(i)) p1 = gneg(p1); + gamun[i+2] = (long)p1; } + gamun = gexp(gamun, prec); /* Gamma(1 + x) */ + gam = gdiv(gamun,x); /* Gamma(x) */ - /* the expansion of log(Gamma(s)) at s = 1/2 */ - gamdm = cgetg(r + 2, t_SER); + /* expansion of log(Gamma(u) / Gamma(1/2)) at u = 1/2 */ + gamdm = cgetg(r+3, t_SER); gamdm[1] = evalsigne(1) | evalvalp(0) | evalvarn(0); - gamdm[2] = (long)mplog(gsqrt(mppi(prec), prec)); - gamdm[3] = lneg(gadd(gmul2n(glog(gdeux, prec), 1), (GEN)cst[1])); + gamdm[2] = zero; + gamdm[3] = lneg(gadd(gmul2n(mplog2(prec), 1), eul)); for (i = 2; i <= r; i++) - gamdm[i + 2] = lmul((GEN)gamun[i + 2], subis(shifti(gun, i), 1)); + gamdm[i+2] = lmul((GEN)gamun[i+2], subis(shifti(gun,i), 1)); + gamdm = gmul(sqpi, gexp(gamdm, prec)); /* Gamma(1/2 + x) */ - gamun = gexp(gamun, prec); - gamdm = gexp(gamdm, prec); - - /* We simplify to get one of the following two expressions */ - - /* Case 1 (b > a): sqrt{Pi}^a 2^{a - as} Gamma(s/2)^{b-a} Gamma(s)^{a + c} */ + /* We simplify to get one of the following two expressions + * if (b > a) : sqrt{Pi}^a 2^{a-au} Gamma(u)^{a+c} Gamma( u/2 )^{|b-a|} + * if (b <= a): sqrt{Pi}^b 2^{b-bu} Gamma(u)^{b+c) Gamma((u+1)/2)^{|b-a|} */ if (b > a) { - cf = gpui(mppi(prec), gmul2n(stoi(a), -1), prec); - - /* an is the expansion of Gamma(x)^{a+c} */ - an = gpowgs(gdiv(gamun, x), a + c); - - /* bn is the expansion of 2^{a - ax} */ - bn = gpowgs(gpow(gdeux, gsubsg(1, x), prec), a); - - /* cn_evn is the expansion of Gamma(x/2)^{b-a} */ - cn_evn = gpowgs(gdiv(gsubst(gamun, 0, x2), x2), b - a); - - /* cn_odd is the expansion of Gamma((x-1)/2)^{b-a} */ - cn_odd = gpowgs(gdiv(gsubst(gamdm, 0, x2), gsub(x2, ghalf)), b - a); - - for (i = 0; i < i0/2; i++) - { - p1 = gmul(cf, gmul(an, gmul(bn, cn_evn))); - - p2 = gdiv(p1, gsubgs(x, 2*i)); - for (j = 1; j <= r; j++) - mael3(aij, 2*i + 1, 1, j) = (long)polcoeff0(p2, -j, 0); - - p2 = gdiv(p1, gsubgs(x, 2*i + 1)); - for (j = 1; j <= r; j++) - mael3(aij, 2*i + 1, 2, j) = (long)polcoeff0(p2, -j, 0); - - /* an(x-s-1) = an(x-s) / (x-s-1)^{a+c} */ - an = gdiv(an, gpowgs(gsubgs(x, 2*i + 1), a + c)); - - /* bn(x-s-1) = 2^a bn(x-s) */ - bn = gmul2n(bn, a); - - /* cn_evn(x-s-2) = cn_evn(x-s) / (x/2 - (s+2)/2)^{b-a} */ - cn_evn = gdiv(cn_evn, gpowgs(gsubgs(x2, i + 1), b - a)); - - p1 = gmul(cf, gmul(an, gmul(bn, cn_odd))); - - p2 = gdiv(p1, gsubgs(x, 2*i + 1)); - for (j = 1; j <= r; j++) - mael3(aij, 2*i + 2, 1, j) = (long)polcoeff0(p2, -j, 0); - - p2 = gdiv(p1, gsubgs(x, 2*i + 2)); - for (j = 1; j <= r; j++) - mael3(aij, 2*i + 2, 2, j) = (long)polcoeff0(p2, -j, 0); - - an = gdiv(an, gpowgs(gsubgs(x, 2*i + 2), a + c)); - bn = gmul2n(bn, a); - - /* cn_odd(x-s-2) = cn_odd(x-s) / (x/2 - (s+2)/2)^{b-a} */ - cn_odd = gdiv(cn_odd, gpowgs(gsub(x2, gaddgs(ghalf, i + 1)), b - a)); - } + t = a; s = b; X = x2; Y = gsub(x2,ghalf); + p1 = gsubst(gam,0,x2); + p2 = gdiv(gsubst(gamdm,0,x2), Y); /* Gamma((x-1)/2) */ } else - /* Case 2 (b <= a): sqrt{Pi}^b 2^{b - bs} Gamma((s+1)/2)^{a-b} - Gamma(s)^{b + c) */ { - cf = gpui(mppi(prec), gmul2n(stoi(b), -1), prec); + t = b; s = a; X = gadd(x2,ghalf); Y = x2; + p1 = gsubst(gamdm,0,x2); + p2 = gsubst(gam,0,x2); + } + cf = gpowgs(sqpi, t); + an = gpowgs(gpow(gdeux, gsubsg(1,x), prec), t); /* 2^{t-tx} */ + bn = gpowgs(gam, t+c); /* Gamma(x)^{t+c} */ + cn_evn = gpowgs(p1, s-t); /* Gamma(X)^{s-t} */ + cn_odd = gpowgs(p2, s-t); /* Gamma(Y)^{s-t} */ + for (i = 0; i < i0/2; i++) + { + GEN C1,q1, A1 = aij[2*i+1], B1 = bij[2*i+1]; + GEN C2,q2, A2 = aij[2*i+2], B2 = bij[2*i+2]; - /* an is the expansion of Gamma(x)^{b+c} */ - an = gpowgs(gdiv(gamun, x), b + c); + C1 = gmul(cf, gmul(bn, gmul(an, cn_evn))); + p1 = gdiv(C1, gsubgs(x, 2*i)); + q1 = gdiv(C1, gsubgs(x, 2*i+1)); - /* bn is the expansion of 2^{b - bx} */ - bn = gpowgs(gpow(gdeux, gsubsg(1, x), prec), b); + /* an(x-u-1) = 2^t an(x-u) */ + an = gmul2n(an, t); + /* bn(x-u-1) = bn(x-u) / (x-u-1)^{t+c} */ + bn = gdiv(bn, gpowgs(gsubgs(x, 2*i+1), t+c)); - /* cn_evn is the expansion of Gamma((x+1)/2)^{a-b} */ - cn_evn = gpowgs(gsubst(gamdm, 0, x2), a - b); - - /* cn_odd is the expansion of Gamma(x/2)^{a-b} */ - cn_odd = gpowgs(gdiv(gsubst(gamun, 0, x2), x2), a - b); - - for (i = 0; i < i0/2; i++) + C2 = gmul(cf, gmul(bn, gmul(an, cn_odd))); + p2 = gdiv(C2, gsubgs(x, 2*i+1)); + q2 = gdiv(C2, gsubgs(x, 2*i+2)); + for (j = 1; j <= r; j++) { - p1 = gmul(cf, gmul(an, gmul(bn, cn_evn))); + affect_coeff(p1, j, A1); affect_coeff(q1, j, B1); + affect_coeff(p2, j, A2); affect_coeff(q2, j, B2); + } - p2 = gdiv(p1, gsubgs(x, 2*i)); - for (j = 1; j <= r; j++) - mael3(aij, 2*i + 1, 1, j) = (long)polcoeff0(p2, -j, 0); + an = gmul2n(an, t); + bn = gdiv(bn, gpowgs(gsubgs(x, 2*i+2), t+c)); + /* cn_evn(x-2i-2) = cn_evn(x-2i) / (X - (i+1))^{s-t} */ + /* cn_odd(x-2i-3) = cn_odd(x-2i-1)/ (Y - (i+1))^{s-t} */ + cn_evn = gdiv(cn_evn, gpowgs(gsubgs(X,i+1), s-t)); + cn_odd = gdiv(cn_odd, gpowgs(gsubgs(Y,i+1), s-t)); + } + T->aij = aij; + T->bij = bij; avma = av; +} - p2 = gdiv(p1, gsubgs(x, 2*i + 1)); - for (j = 1; j <= r; j++) - mael3(aij, 2*i + 1, 2, j) = (long)polcoeff0(p2, -j, 0); +static GEN +_cond(GEN dtcr, int quad) +{ + GEN cond; + if (quad) cond = (GEN)dtcr[7]; + else + { + cond = cgetg(3, t_VEC); + cond[1] = dtcr[7]; + cond[2] = dtcr[9]; + } + return cond; +} - /* an(x-s-1) = an(x-s) / (x-s-1)^{b+c} */ - an = gdiv(an, gpowgs(gsubgs(x, 2*i + 1), b + c)); +/* sort chars according to conductor */ +static GEN +sortChars(GEN dataCR, int quad) +{ + const long cl = lg(dataCR) - 1; + GEN vCond = cgetg(cl+1, t_VEC); + GEN CC = cgetg(cl+1, t_VECSMALL); + GEN nvCond = cgetg(cl+1, t_VECSMALL); + long j,k, ncond; + GEN vChar; - /* bn(x-s-1) = 2^b bn(x-s) */ - bn = gmul2n(bn, b); + for (j = 1; j <= cl; j++) nvCond[j] = 0; - /* cn_evn(x-s-2) = cn_evn(x-s) / (x/2 - (s+1)/2)^{a-b} */ - cn_evn = gdiv(cn_evn, gpowgs(gsub(x2, gaddgs(ghalf, i)), a - b)); + ncond = 0; + for (j = 1; j <= cl; j++) + { + GEN cond = _cond((GEN)dataCR[j], quad); + for (k = 1; k <= ncond; k++) + if (gegal(cond, (GEN)vCond[k])) break; + if (k > ncond) vCond[++ncond] = (long)cond; + nvCond[k]++; CC[j] = k; /* char j has conductor number k */ + } + vChar = cgetg(ncond+1, t_VEC); + for (k = 1; k <= ncond; k++) + { + vChar[k] = lgetg(nvCond[k]+1, t_VECSMALL); + nvCond[k] = 0; + } + for (j = 1; j <= cl; j++) + { + k = CC[j]; nvCond[k]++; + mael(vChar,k,nvCond[k]) = j; + } + return vChar; +} - p1 = gmul(cf, gmul(an, gmul(bn, cn_odd))); +static void +get_cS_cT(ST_t *T, long n) +{ + gpmem_t av; + GEN csurn, nsurc, lncsurn; + GEN A,B,s,t,Z,*aij,*bij; + long i,j,r,i0; - p2 = gdiv(p1, gsubgs(x, 2*i + 1)); - for (j = 1; j <= r; j++) - mael3(aij, 2*i + 2, 1, j) = (long)polcoeff0(p2, -j, 0); + if (T->cS[n]) return; - p2 = gdiv(p1, gsubgs(x, 2*i + 2)); - for (j = 1; j <= r; j++) - mael3(aij, 2*i + 2, 2, j) = (long)polcoeff0(p2, -j, 0); + av = avma; + aij = T->aij; i0= T->i0; + bij = T->bij; r = T->r; + Z = cgetg(r+1, t_VEC); + Z[1] = un; - an = gdiv(an, gpowgs(gsubgs(x, 2*i + 2), b + c)); - bn = gmul2n(bn, b); + csurn = gdivgs(T->c1, n); + nsurc = ginv(csurn); + lncsurn = mplog(csurn); - /* cn_odd(x-s-2) = cn_odd(x-s) / (x/2 - (s+1)/2)^{a-b} */ - cn_odd = gdiv(cn_odd, gpowgs(gsubgs(x2, i + 1), a - b)); + Z[2] = (long)lncsurn; /* r >= 2 */ + for (i = 3; i <= r; i++) + { + s = gmul((GEN)Z[i-1], lncsurn); + Z[i] = ldivgs(s, i-1); /* Z[i] = ln^(i-1)(c1/n) / (i-1)! */ + } + + /* i = i0 */ + A = aij[i0]; t = (GEN)A[1]; + B = bij[i0]; s = (GEN)B[1]; + for (j = 2; j <= r; j++) + { + s = gadd(s, gmul((GEN)Z[j], (GEN)B[j])); + t = gadd(t, gmul((GEN)Z[j], (GEN)A[j])); } + for (i = i0 - 1; i > 1; i--) + { + A = aij[i]; t = gmul(t, nsurc); + B = bij[i]; s = gmul(s, nsurc); + for (j = odd(i)? T->rc2: T->rc1; j; j--) + { + s = gadd(s, gmul((GEN)Z[j], (GEN)B[j])); + t = gadd(t, gmul((GEN)Z[j], (GEN)A[j])); + } } + /* i = 1 */ + A = aij[1]; t = gmul(t, nsurc); + B = bij[1]; s = gmul(s, nsurc); + for (j = 1; j <= r; j++) + { + s = gadd(s, gmul((GEN)Z[j], (GEN)B[j])); + t = gadd(t, gmul((GEN)Z[j], (GEN)A[j])); + } + s = gadd(s, gmul(csurn, T->powracpi[T->b])); + T->cS[n] = gclone(s); + T->cT[n] = gclone(t); avma = av; +} - return gerepilecopy(av, aij); +static void +clear_cScT(ST_t *T, long N) +{ + GEN *cS = T->cS, *cT = T->cT; + long i; + for (i=1; i<=N; i++) + if (cS[i]) { gunclone(cS[i]); gunclone(cT[i]); cS[i] = cT[i] = NULL; } } +static void +init_cScT(ST_t *T, GEN CHI, long N, long prec) +{ + GEN p1 = (GEN)CHI[9]; + T->a = p1[1]; + T->b = p1[2]; + T->c = p1[3]; + T->rc1 = T->a + T->c; + T->rc2 = T->b + T->c; + T->r = max(T->rc2+1, T->rc1); /* >= 2 */ + ppgamma(T, prec); + clear_cScT(T, N); +} + static GEN -GetST(GEN dataCR, long prec) +GetST(GEN dataCR, GEN vChar, long prec) { - GEN N0, CC, bnr, bnf, nf, Pi, racpi, C, cond, aij, B, S, T, csurn, lncsurn; - GEN degs, p1, p2, nsurc, an, rep, powlncn, powracpi; - long i, j, k, n, av = avma, av1, av2, hk, fj, id, prec2, i0, nmax; - long a, b, c, rc1, rc2, r, r1, r2; - int ***matan; + const long cl = lg(dataCR) - 1; + gpmem_t av, av1, av2; + long ncond, n, j, k, jc, nmax, prec2, i0, r1, r2; + GEN bnr, nf, racpi, *powracpi; + GEN rep, N0, C, T, S, an, degs; + LISTray LIST; + ST_t cScT; - if (DEBUGLEVEL) timer2(); - bnr = gmael(dataCR, 1, 4); - bnf = checkbnf(bnr); - nf = checknf(bnf); - r1 = nf_get_r1(nf); - r2 = nf_get_r2(nf); - hk = lg(dataCR) - 1; - prec2 = ((prec - 2)<<1) + EXTRA_PREC; + bnr = gmael(dataCR,1,4); + nf = checknf(bnr); + /* compute S,T differently if nf is quadratic */ + if (degpol(nf[1]) == 2) return QuadGetST(dataCR, vChar, prec); - Pi = mppi(prec2); - racpi = gsqrt(Pi, prec2); + if (DEBUGLEVEL) (void)timer2(); + /* allocate memory for answer */ + rep = cgetg(3, t_VEC); + S = cgetg(cl+1, t_VEC); rep[1] = (long)S; + T = cgetg(cl+1, t_VEC); rep[2] = (long)T; + for (j = 1; j <= cl; j++) + { + S[j] = (long)cgetc(prec); + T[j] = (long)cgetc(prec); + } + av = avma; - C = cgetg(hk + 1, t_VEC); - cond = cgetg(hk + 1, t_VEC); - N0 = new_chunk(hk+1); - CC = new_chunk(hk+1); + /* initializations */ + degs = GetDeg(dataCR); + ncond = lg(vChar)-1; + nf_get_sign(nf,&r1,&r2); + + C = cgetg(ncond+1, t_VEC); + N0 = cgetg(ncond+1, t_VECSMALL); nmax = 0; - for (i = 1; i <= hk; i++) + for (j = 1; j <= ncond; j++) { - C[i] = mael(dataCR, i, 2); - - p1 = cgetg(3, t_VEC); - p1[1] = mael(dataCR, i, 7); - p1[2] = mael(dataCR, i, 9); - cond[i] = (long)p1; - - N0[i] = GetBoundN0((GEN)C[i], r1, r2, prec, 0); - if (nmax < N0[i]) nmax = N0[i]; + C[j] = mael(dataCR, mael(vChar,j,1), 2); + N0[j] = GetBoundN0((GEN)C[j], r1, r2, prec); + if (nmax < N0[j]) nmax = N0[j]; } - if ((ulong)nmax > maxprime()) - err(talker, "Not enough precomputed primes (need all primes up to %ld)", nmax); - + err(talker, "Not enough precomputed primes (need all p <= %ld)", nmax); i0 = GetBoundi0(r1, r2, prec); - if(DEBUGLEVEL > 1) fprintferr("nmax = %ld and i0 = %ld\n", nmax, i0); + if (DEBUGLEVEL>1) fprintferr("nmax = %ld, i0 = %ld\n", nmax, i0); + InitPrimes(gmael(dataCR,1,4), nmax, &LIST); - matan = ComputeCoeff(dataCR, nmax, prec); - degs = GetDeg(dataCR); - if (DEBUGLEVEL) msgtimer("Compute an"); + prec2 = ((prec-2) << 1) + EXTRA_PREC; + racpi = mpsqrt(mppi(prec2)); + powracpi = (GEN*)cgetg(r1+2,t_VEC); + powracpi++; powracpi[0] = gun; powracpi[1] = racpi; + for (j=2; j<=r1; j++) powracpi[j] = mulrr(powracpi[j-1], racpi); + cScT.powracpi = powracpi; - p1 = cgetg(3, t_COMPLEX); - p1[1] = lgetr(prec2); - p1[2] = lgetr(prec2); - gaffect(gzero, p1); + cScT.cS = (GEN*)cgetg(nmax+1, t_VEC); + cScT.cT = (GEN*)cgetg(nmax+1, t_VEC); + for (j=1; j<=nmax; j++) cScT.cS[j] = cScT.cT[j] = NULL; - S = cgetg(hk + 1, t_VEC); - T = cgetg(hk + 1, t_VEC); - for (id = 1; id <= hk; id++) + cScT.i0 = i0; + + av1 = avma; + for (jc = 1; jc <= ncond; jc++) { - S[id] = lcopy(p1); - T[id] = lcopy(p1); - for (k = 1; k < id; k++) - if (gegal((GEN)cond[id], (GEN)cond[k])) break; - CC[id] = k; - } + const GEN LChar = (GEN)vChar[jc]; + const long nChar = lg(LChar)-1, NN = N0[jc]; - powracpi = cgetg(hk + 1, t_VEC); - for (j = 1; j <= hk; j++) - powracpi[j] = (long)gpow(racpi, gmael3(dataCR, j, 9, 2), prec2); + if (DEBUGLEVEL>1) + fprintferr("* conductor no %ld/%ld (N = %ld)\n\tInit: ", jc,ncond,NN); - av1 = avma; - if (DEBUGLEVEL > 1) fprintferr("n = "); - - for (id = 1; id <= hk; id++) - { - if (CC[id] != id) continue; - p2 = gmael(dataCR, id, 9); - a = itos((GEN)p2[1]); - b = itos((GEN)p2[2]); - c = itos((GEN)p2[3]); - aij = ppgamma(a, b, c, i0, prec2); - rc1 = a + c; - rc2 = b + c; r = max(rc2 + 1, rc1); + cScT.c1 = (GEN)C[jc]; + init_cScT(&cScT, (GEN)dataCR[LChar[1]], NN, prec2); av2 = avma; - - for (n = 1; n <= N0[id]; n++) + for (k = 1; k <= nChar; k++) { - if (DEBUGLEVEL > 1 && n%100 == 0) fprintferr("%ld ", n); - - for (k = 1; k <= hk; k++) - if (CC[k] == id && !IsZero(matan[k][n], degs[k])) break; - if (k > hk) continue; - - csurn = gdivgs((GEN)C[id], n); - nsurc = ginv(csurn); - - B = cgetg(r + 1, t_VEC); - lncsurn = glog(csurn, prec2); - powlncn = gun; - fj = 1; - - p1 = gzero; - p2 = gzero; - for (j = 1; j <= r; j++) - { - if (j > 2) fj = fj * (j - 1); - - B[j] = ldivgs(powlncn, fj); - p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, i0, 2, j))); - p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, i0, 1, j))); - - powlncn = gmul(powlncn, lncsurn); - } - for (i = i0 - 1; i > 1; i--) - { - p1 = gmul(p1, nsurc); - p2 = gmul(p2, nsurc); - for (j = i%2? rc2: rc1; j; j--) + const long t = LChar[k], d = degs[t]; + const GEN z = gmael3(dataCR, t, 5, 2); + GEN p1 = gzero, p2 = gzero; + long c = 0; + int **matan; + + if (DEBUGLEVEL>1) + fprintferr("\tcharacter no: %ld (%ld/%ld)\n", t,k,nChar); + matan = ComputeCoeff((GEN)dataCR[t], &LIST, NN, d); + for (n = 1; n <= NN; n++) + if ((an = EvalCoeff(z, matan[n], d))) { - p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, i, 2, j))); - p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, i, 1, j))); + get_cS_cT(&cScT, n); + p1 = gadd(p1, gmul(an, cScT.cS[n])); + p2 = gadd(p2, gmul(gconj(an), cScT.cT[n])); + if (++c == 256) + { GEN *gptr[2]; gptr[0]=&p1; gptr[1]=&p2; + gerepilemany(av2,gptr,2); c = 0; + } } - } - p1 = gmul(p1, nsurc); - p2 = gmul(p2, nsurc); - for (j = 1; j <= r; j++) - { - p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, 1, 2, j))); - p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, 1, 1, j))); - } - - p1 = gadd(p1, gmul(csurn, (GEN)powracpi[id])); - - for (j = 1; j <= hk; j++) - if (CC[j] == id && - (an = EvalCoeff(gmael3(dataCR, j, 5, 2), matan[j][n], degs[j]))) - { - gaffect(gadd((GEN)S[j], gmul(p1, an)), (GEN)S[j]); - gaffect(gadd((GEN)T[j], gmul(p2, gconj(an))), (GEN)T[j]); - } - avma = av2; + gaffect(p1, (GEN)S[t]); + gaffect(p2, (GEN)T[t]); + FreeMat(matan, NN); avma = av2; } + if (DEBUGLEVEL>1) fprintferr("\n"); avma = av1; } - FreeMat(matan); - - if (DEBUGLEVEL > 1) fprintferr("\n"); - if (DEBUGLEVEL) msgtimer("Compute S&T"); - - rep = cgetg(3, t_VEC); - rep[1] = (long)S; - rep[2] = (long)T; - return gerepilecopy(av, rep); + if (DEBUGLEVEL) msgtimer("S&T"); + clear_cScT(&cScT, nmax); + avma = av; return rep; } -/* Given datachi, S(chi) and T(chi), return L(1, chi) if fl = 1, +/* Given W(chi) [possibly NULL], S(chi) and T(chi), return L(1, chi) if fl = 1, or [r(chi), c(chi)] where r(chi) is the rank of chi and c(chi) is given by L(s, chi) = c(chi).s^r(chi) at s = 0 if fl = 0. if fl2 = 1, adjust the value to get L_S(s, chi). */ static GEN -GetValue(GEN datachi, GEN S, GEN T, long fl, long fl2, long prec) +GetValue(GEN dtcr, GEN W, GEN S, GEN T, long fl, long fl2, long prec) { - GEN W, A, q, b, c, d, rchi, cf, VL, rep, racpi, nS, nT; - long av = avma; + gpmem_t av = avma; + GEN A, cf, VL, rep, racpi, p1; + long q, b, c; + int isreal; - racpi = gsqrt(mppi(prec), prec); - W = ComputeArtinNumber(datachi, 0, prec); - A = ComputeAChi(datachi, fl, prec); + racpi = mpsqrt(mppi(prec)); + if (!W) + W = (GEN)ComputeArtinNumber((GEN)dtcr[3], _vec((GEN)dtcr[8]), 1, prec)[1]; - d = gmael(datachi, 8, 3); + isreal = (itos(gmael(dtcr, 8, 3)) <= 2); - q = gmael(datachi, 9, 1); - b = gmael(datachi, 9, 2); - c = gmael(datachi, 9, 3); + p1 = (GEN)dtcr[9]; + q = p1[1]; + b = p1[2]; + c = p1[3]; - rchi = addii(b, c); - if (!fl) - { - cf = gmul2n(gpow(racpi, q, 0), itos(b)); + { /* VL = (W(chi).S(conj(chi)) + T(chi)) / (sqrt(Pi)^q 2^{r1 - q}) */ + GEN rchi = stoi(b + c); + cf = gmul2n(gpowgs(racpi, q), b); - nS = gdiv(gconj(S), cf); - nT = gdiv(gconj(T), cf); + VL = gadd(gmul(W, gconj(S)), gconj(T)); + if (isreal) VL = greal(VL); + VL = gdiv(VL, cf); - /* VL = W(chi).S(conj(chi)) + T(chi)) / (sqrt(Pi)^q 2^{r1 - q}) */ - VL = gadd(gmul(W, nS), nT); - if (cmpis(d, 3) < 0) VL = greal(VL); - if (fl2) { + A = ComputeAChi(dtcr, fl, prec); VL = gmul((GEN)A[2], VL); rchi = gadd(rchi, (GEN)A[1]); } @@ -2119,17 +2079,21 @@ GetValue(GEN datachi, GEN S, GEN T, long fl, long fl2, rep[2] = (long)VL; } else - { - cf = gmul((GEN)datachi[2], gpow(racpi, b, 0)); + { /* VL = S(chi) + W(chi).T(chi)) / (C(chi) sqrt(Pi)^{r1 - q}) */ + cf = gmul((GEN)dtcr[2], gpowgs(racpi, b)); - /* VL = S(chi) + W(chi).T(chi)) / (C(chi) sqrt(Pi)^{r1 - q}) */ - rep = gdiv(gadd(S, gmul(W, T)), cf); - if (cmpis(d, 3) < 0) rep = greal(rep); + rep = gadd(S, gmul(W, T)); + if (isreal) rep = greal(rep); + rep = gdiv(rep, cf); - if (fl2) rep = gmul(A, rep); + if (fl2) + { + A = ComputeAChi(dtcr, fl, prec); + rep = gmul(A, rep); + } } - return gerepilecopy(av, rep); + return gerepilecopy(av, rep); } /* return the order and the first non-zero term of L(s, chi0) @@ -2139,14 +2103,15 @@ GetValue1(GEN bnr, long flag, long prec) { GEN bnf = checkbnf(bnr), nf = checknf(bnf); GEN hk, Rk, wk, c, rep, mod0, diff; - long i, l, r, r1, r2, av = avma; + long i, l, r, r1, r2; + gpmem_t av = avma; r1 = nf_get_r1(nf); r2 = nf_get_r2(nf); hk = gmael3(bnf, 8, 1, 1); Rk = gmael(bnf, 8, 2); wk = gmael3(bnf, 8, 4, 1); - + c = gneg_i(gdiv(gmul(hk, Rk), wk)); r = r1 + r2 - 1; @@ -2163,126 +2128,72 @@ GetValue1(GEN bnr, long flag, long prec) rep = cgetg(3, t_VEC); rep[1] = lstoi(r); rep[2] = (long)c; - return gerepilecopy(av, rep); } /********************************************************************/ /* 6th part: recover the coefficients */ /********************************************************************/ - static long -TestOne(GEN plg, GEN beta, GEN B, long v, long G, long N) +TestOne(GEN plg, RC_data *d) { - long j; + long j, v = d->v; GEN p1; - p1 = gsub(beta, (GEN)plg[v]); - if (expo(p1) >= G) return 0; + p1 = gsub(d->beta, (GEN)plg[v]); + if (expo(p1) >= d->G) return 0; - for (j = 1; j <= N; j++) + for (j = 1; j <= d->N; j++) if (j != v) { p1 = gabs((GEN)plg[j], DEFAULTPREC); - if (gcmp(p1, B) > 0) return 0; + if (gcmp(p1, d->B) > 0) return 0; } return 1; } -/* Using linear dependance relations */ static GEN -RecCoeff2(GEN nf, GEN beta, GEN B, long v, long prec) +chk_reccoeff_init(FP_chk_fun *chk, GEN gram, GEN mat) { - long N, G, i, bacmin, bacmax, av = avma, av2; - GEN vec, velt, p1, cand, M, plg, pol, cand2; - - M = gmael(nf, 5, 1); - pol = (GEN)nf[1]; - N = degpol(pol); - vec = gtrans((GEN)gtrans(M)[v]); - velt = (GEN)nf[7]; - - G = min( - 20, - bit_accuracy(prec) >> 4); - - p1 = cgetg(2, t_VEC); - - p1[1] = lneg(beta); - vec = concat(p1, vec); - - p1[1] = zero; - velt = concat(p1, velt); - - bacmin = (long)(.225 * bit_accuracy(prec)); - bacmax = (long)(.315 * bit_accuracy(prec)); - - av2 = avma; - - for (i = bacmax; i >= bacmin; i--) - { - p1 = lindep2(vec, i); - - if (signe((GEN)p1[1])) - { - p1 = ground(gdiv(p1, (GEN)p1[1])); - cand = gmodulcp(gmul(velt, gtrans(p1)), pol); - cand2 = algtobasis(nf, cand); - plg = gmul(M, cand2); - - if (TestOne(plg, beta, B, v, G, N)) - return gerepilecopy(av, cand); - } - avma = av2; - } - return NULL; + RC_data *d = (RC_data*)chk->data; + (void)gram; d->U = mat; return d->nB; } -GEN -chk_reccoeff_init(FP_chk_fun *chk, GEN nf, GEN gram, GEN mat, long *ptprec) +static GEN +chk_reccoeff(void *data, GEN x) { - GEN data = chk->data; - data[6] = (long)mat; - chk->data = data; - return (GEN)data[7]; -} + RC_data *d = (RC_data*)data; + long N = d->N, j; + GEN p1 = gmul(d->U, x), sol, plg; -GEN -chk_reccoeff(GEN data, GEN x) -{ - GEN M = (GEN)data[0], beta = (GEN)data[1], B = (GEN)data[2]; - long v = data[3], G = data[4], N = data[5], j; - GEN U = (GEN)data[6], p1 = gmul(U, x), sol, plg; - if (!gcmp1((GEN)p1[1])) return NULL; - + sol = cgetg(N + 1, t_COL); for (j = 1; j <= N; j++) sol[j] = lmulii((GEN)p1[1], (GEN)p1[j + 1]); - plg = gmul(M, sol); - - if (TestOne(plg, beta, B, v, G, N)) return sol; + plg = gmul(d->M, sol); + + if (TestOne(plg, d)) return sol; return NULL; } -GEN -chk_reccoeff_post(GEN data, GEN res) -{ - return res; -} - /* Using Cohen's method */ static GEN -RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec) +RecCoeff3(GEN nf, RC_data *d, long prec) { - GEN A, M, nB, cand, p1, B2, C2, data, tB, beta2, eps, nf2, Bd; - long N, G, i, j, k, l, ct = 0, av = avma, prec2; - FP_chk_fun *chk; + GEN A, M, nB, cand, p1, B2, C2, tB, beta2, eps, nf2, Bd; + GEN beta = d->beta, B = d->B; + long N = d->N, v = d->v; + long i, j, k, l, ct = 0, prec2; + gpmem_t av = avma; + FP_chk_fun chk = { &chk_reccoeff, &chk_reccoeff_init, NULL, 0 }; + chk.data = (void*)d; - N = degpol(nf[1]); - G = min(-10, -bit_accuracy(prec) >> 4); - eps = gpowgs(stoi(10), min(-8, (G >> 1))); + d->G = min(-10, -bit_accuracy(prec) >> 4); + eps = gpowgs(stoi(10), min(-8, (d->G >> 1))); tB = gpow(gmul2n(eps, N), gdivgs(gun, 1-N), DEFAULTPREC); - Bd = gmin(B, tB); + Bd = gceil(gmin(B, tB)); p1 = gceil(gdiv(glog(Bd, DEFAULTPREC), dbltor(2.3026))); prec2 = max((prec << 1) - 2, (long)(itos(p1) * pariK1 + BIGDEFAULTPREC)); nf2 = nfnewprec(nf, prec2); @@ -2293,6 +2204,7 @@ RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec) C2 = gdiv(B2, gsqr(eps)); M = gmael(nf2, 5, 1); + d->M = M; A = cgetg(N+2, t_MAT); for (i = 1; i <= N+1; i++) @@ -2318,26 +2230,9 @@ RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec) } nB = mulsi(N+1, B2); + d->nB = nB; + cand = fincke_pohst(A, nB, 20000, 1, prec2, &chk); - data = new_chunk(8); - data[0] = (long)M; - data[1] = (long)beta; - data[2] = (long)B; - data[3] = v; - data[4] = G; - data[5] = N; - data[6] = (long)NULL; - data[7] = (long)nB; - - chk = (FP_chk_fun*)new_chunk(sizeof(FP_chk_fun)); - chk->f = &chk_reccoeff; - chk->f_init = &chk_reccoeff_init; - chk->f_post = &chk_reccoeff_post; - chk->data = data; - chk->skipfirst = 0; - - cand = fincke_pohst(A, nB, 20000, 3, prec2, chk); - if (!cand) { if (ct > 3) { avma = av; return NULL; } @@ -2360,14 +2255,62 @@ RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec) avma = av; return NULL; } -/* Attempts to find a polynomial with coefficients in nf such that - its coefficients are close to those of pol at the place v and +/* Using linear dependance relations */ +static GEN +RecCoeff2(GEN nf, RC_data *d, long prec) +{ + long i, bacmin, bacmax; + gpmem_t av = avma, av2; + GEN vec, velt, p1, cand, M, plg, pol, cand2; + GEN beta = d->beta; + + d->G = min(-20, -bit_accuracy(prec) >> 4); + M = gmael(nf, 5, 1); + pol = (GEN)nf[1]; + vec = gtrans((GEN)gtrans(M)[d->v]); + velt = (GEN)nf[7]; + + p1 = cgetg(2, t_VEC); + p1[1] = lneg(beta); + vec = concat(p1, vec); + + p1[1] = zero; + velt = concat(p1, velt); + + bacmin = (long)(.225 * bit_accuracy(prec)); + bacmax = (long)(.315 * bit_accuracy(prec)); + + av2 = avma; + + for (i = bacmax; i >= bacmin; i-=16) + { + p1 = lindep2(vec, i); + + if (signe((GEN)p1[1])) + { + p1 = ground(gdiv(p1, (GEN)p1[1])); + cand = gmodulcp(gmul(velt, gtrans(p1)), pol); + cand2 = algtobasis(nf, cand); + plg = gmul(M, cand2); + + if (TestOne(plg, d)) return gerepilecopy(av, cand); + } + avma = av2; + } + /* failure */ + return RecCoeff3(nf,d,prec); +} + +/* Attempts to find a polynomial with coefficients in nf such that + its coefficients are close to those of pol at the place v and less than B at all the other places */ GEN RecCoeff(GEN nf, GEN pol, long v, long prec) { - long av = avma, j, md, G, cl = degpol(pol); + long j, md, G, cl = degpol(pol); + gpmem_t av = avma; GEN p1, beta; + RC_data d; /* if precision(pol) is too low, abort */ for (j = 2; j <= cl+1; j++) @@ -2380,25 +2323,24 @@ RecCoeff(GEN nf, GEN pol, long v, long prec) md = cl/2; pol = dummycopy(pol); + d.N = degpol(nf[1]); + d.v = v; + for (j = 1; j <= cl; j++) - { - /* start with the coefficients in the middle, + { /* start with the coefficients in the middle, since they are the harder to recognize! */ long cf = md + (j%2? j/2: -j/2); GEN bound = binome(stoi(cl), cf); bound = shifti(bound, cl - cf); - if (DEBUGLEVEL > 1) fprintferr("In RecCoeff with cf = %ld and B = %Z\n", - cf, bound); + if (DEBUGLEVEL > 1) + fprintferr("In RecCoeff with cf = %ld and B = %Z\n", cf, bound); beta = greal((GEN)pol[cf+2]); - p1 = RecCoeff2(nf, beta, bound, v, prec); - if (!p1) - { - p1 = RecCoeff3(nf, beta, bound, v, prec); - if (!p1) return NULL; - } + d.beta = beta; + d.B = bound; + if (! (p1 = RecCoeff2(nf, &d, prec)) ) return NULL; pol[cf+2] = (long)p1; } pol[cl+2] = un; @@ -2406,470 +2348,291 @@ RecCoeff(GEN nf, GEN pol, long v, long prec) } /*******************************************************************/ -/*******************************************************************/ /* */ -/* Computation of class fields of */ -/* real quadratic fields using Stark units */ +/* Class fields of real quadratic fields using Stark units */ /* */ /*******************************************************************/ -/*******************************************************************/ -/* compute the coefficients an for the quadratic case */ -static int*** -computean(GEN dtcr, long nmax, long prec) +/* an[q * i] *= chi for all (i,p)=1 */ +static void +an_mul(int **an, long p, long q, long n, long deg, GEN chi, int **reduc) { - long i, j, cl, q, cp, al, v1, v2, v, fldiv, av, av1; - int ***matan, ***reduc; - GEN bnf, ideal, dk, degs, idno, p1, prime, chi, qg, chi1, chi2; - GEN chi11, chi12, bnr, pr, pr1, pr2, xray, xray1, xray2, dataray; - byteptr dp = diffptr; + gpmem_t av; + long c,i; + int *T; + if (gcmp1(chi)) return; av = avma; + T = (int*)new_chunk(deg); Polmod2Coeff(T,chi, deg); + for (c = 1, i = q; i <= n; i += q, c++) + if (c == p) c = 0; else MulCoeff(an[i], T, reduc, deg); + avma = av; +} +/* an[q * i] = 0 for all (i,p)=1 */ +static void +an_set0_coprime(int **an, long p, long q, long n, long deg) +{ + long c,i; + for (c = 1, i = q; i <= n; i += q, c++) + if (c == p) c = 0; else _0toCoeff(an[i], deg); +} +/* an[q * i] = 0 for all i */ +static void +an_set0(int **an, long p, long n, long deg) +{ + long i; + for (i = p; i <= n; i += p) _0toCoeff(an[i], deg); +} - cl = lg(dtcr) - 1; - degs = GetDeg(dtcr); +/* compute the coefficients an for the quadratic case */ +static int** +computean(GEN dtcr, LISTray *R, long n, long deg) +{ + gpmem_t av = avma, av2; + long i, p, q, condZ, l; + int **an, **reduc; + GEN L, CHI, chi, chi1, ray; + CHI_t C; - matan = InitMatAn(cl, nmax, degs, 1); - reduc = InitReduction(dtcr, degs); + CHI = (GEN)dtcr[5]; init_CHI_alg(&C, CHI); + condZ= R->condZ; - bnr = gmael(dtcr, 1, 4); bnf = (GEN)bnr[1]; - dataray = InitGetRay(bnr, nmax); + an = InitMatAn(n, deg, 1); + reduc = InitReduction(CHI, deg); + av2 = avma; - ideal = gmael3(bnr, 2, 1, 1); - idno = idealnorm(bnf, ideal); - dk = gmael(bnf, 7, 3); + /* all pr | p divide cond */ + L = R->L0; l = lg(L); + for (i=1; iL2; l = lg(L); + for (i=1; irayZ[p % condZ]; + chi = EvalChar(&C, ray); + } + chi1 = chi; + for (q=p;;) + { + an_set0_coprime(an, p,q,n,deg); /* v_p(q) odd */ + if (! (q = next_pow(q,p, n)) ) break; - switch (krogs(dk, prime[2])) + an_mul(an,p,q,n,deg,chi,reduc); + if (! (q = next_pow(q,p, n)) ) break; + chi = gmul(chi, chi1); + } + } + + /* 1 prime of degree 1 */ + L = R->L1; l = lg(L); + for (i=1; iL1ray[i]; + chi = EvalChar(&C, ray); + chi1 = chi; + for(q=p;;) { - /* prime is inert */ - case -1: - fldiv = divise(idno, prime); + an_mul(an,p,q,n,deg,chi,reduc); + if (! (q = next_pow(q,p, n)) ) break; + chi = gmul(chi, chi1); + } + } - if (!fldiv) - { - xray = GetRay(bnr, dataray, prime, prec); - chi = chiideal(dtcr, xray, 1); - chi1 = dummycopy(chi); - } + /* 2 primes of degree 1 */ + L = R->L11; l = lg(L); + for (i=1; iL11ray[i]; /* use pr1 pr2 = (p) */ + if (condZ == 1) + ray2 = gneg(ray1); + else + ray2 = gsub(R->rayZ[p % condZ], ray1); + chi11 = EvalChar(&C, ray1); + chi12 = EvalChar(&C, ray2); - for (cp = 1, i = q; i <= nmax; i += q, cp++) - if(cp % prime[2]) - { - if (fldiv || al%2) - for (j = 1; j <= cl; j++) - _0toCoeff(matan[j][i], degs[j]); - else - for (j = 1; j <= cl; j++) - MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]); - } - - qg = mulsi(q, prime); - al++; - - if (al%2 && !fldiv) - for (j = 1; j <= cl; j++) - chi[j] = lmul((GEN)chi[j], (GEN)chi1[j]); - } - break; - - /* prime is ramified */ - case 0: - fldiv = divise(idno, prime); - - if (!fldiv) - { - pr = (GEN)primedec(bnf, prime)[1]; - xray = GetRay(bnr, dataray, pr, prec); - chi = chiideal(dtcr, xray, 1); - chi2 = dummycopy(chi); - } - - while(cmpis(qg, nmax) <= 0) - { - q = qg[2]; - - for (cp = 1, i = q; i <= nmax; i += q, cp++) - if(cp % prime[2]) - { - if (fldiv) - for(j = 1; j <= cl; j++) - _0toCoeff(matan[j][i], degs[j]); - else - { - for(j = 1; j <= cl; j++) - MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]); - } - } - - qg = mulsi(q, prime); - al++; - - if (cmpis(qg, nmax) <= 0 && !fldiv) - for (j = 1; j <= cl; j++) - chi[j] = lmul((GEN)chi[j], (GEN)chi2[j]); - } - break; - - /* prime is split */ - default: /* case 1: */ - p1 = primedec(bnf, prime); - pr1 = (GEN)p1[1]; - pr2 = (GEN)p1[2]; - v1 = idealval(bnf, ideal, pr1); - v2 = idealval(bnf, ideal, pr2); - - if (v1 + v2 == 0) - { - xray1 = GetRay(bnr, dataray, pr1, prec); - xray2 = GetRay(bnr, dataray, pr2, prec); - chi11 = chiideal(dtcr, xray1, 1); - chi12 = chiideal(dtcr, xray2, 1); - - chi1 = gadd(chi11, chi12); - chi2 = dummycopy(chi12); - - while(cmpis(qg, nmax) <= 0) - { - q = qg[2]; - - for (cp = 1, i = q; i <= nmax; i += q, cp++) - if(cp % prime[2]) - for(j = 1; j <= cl; j++) - MulPolmodCoeff((GEN)chi1[j], matan[j][i], reduc[j], degs[j]); - - qg = mulsi(q, prime); - al++; - - if(cmpis(qg, nmax) <= 0) - for (j = 1; j <= cl; j++) - { - chi2[j] = lmul((GEN)chi2[j], (GEN)chi12[j]); - chi1[j] = ladd((GEN)chi2[j], gmul((GEN)chi1[j], (GEN)chi11[j])); - } - } - } - else - { - if (v1) { v = v2; pr = pr2; } else { v = v1; pr = pr1; } - - if (v == 0) - { - xray = GetRay(bnr, dataray, pr, prec); - chi1 = chiideal(dtcr, xray, 1); - chi = gcopy(chi1); - } - - while(cmpis(qg, nmax) <= 0) - { - q = qg[2]; - for (cp = 1, i = q; i <= nmax; i += q, cp++) - if(cp % prime[2]) - { - if (v) - for (j = 1; j <= cl; j++) - _0toCoeff(matan[j][i], degs[j]); - else - for (j = 1; j <= cl; j++) - MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]); - } - - qg = mulii(qg, prime); - al++; - - if (!v && (cmpis(qg, nmax) <= 0)) - for (j = 1; j <= cl; j++) - chi[j] = lmul((GEN)chi[j], (GEN)chi1[j]); - } - } - break; + chi1 = gadd(chi11, chi12); + chi2 = chi12; + for(q=p;;) + { + an_mul(an,p,q,n,deg,chi1,reduc); + if (! (q = next_pow(q,p, n)) ) break; + chi2 = gmul(chi2, chi12); + chi1 = gadd(chi2, gmul(chi1, chi11)); } - - prime[2] += (*dp++); - - avma = av1; } - for (i = 1; i <= cl; i++) - CorrectCoeff((GEN)dtcr[i], matan[i], reduc[i], nmax, degs[i]); - - FreeMat(reduc); - avma = av; return matan; + CorrectCoeff(dtcr, an, reduc, n, deg); + FreeMat(reduc, deg-1); + avma = av; return an; } /* compute S and T for the quadratic case */ static GEN -QuadGetST(GEN data, long prec) +QuadGetST(GEN dataCR, GEN vChar, long prec) { - long av = avma, n, j, nmax, cl, av1, av2, k; - int ***matan; - GEN nn, C, p1, p2, c2, cexp, cn, v, veclprime2, veclprime1; - GEN dtcr, cond, rep, an, cf, degs, veint1; + const long cl = lg(dataCR) - 1; + gpmem_t av, av1, av2; + long ncond, n, j, k, nmax; + GEN rep, N0, C, T, S, cf, cfh, an, degs; + LISTray LIST; - dtcr = (GEN)data[5]; - cl = lg(dtcr) - 1; - degs = GetDeg(dtcr); - - cf = gmul2n(mpsqrt(mppi(prec)), 1); - C = cgetg(cl+1, t_VEC); - cond = cgetg(cl+1, t_VEC); - c2 = cgetg(cl + 1, t_VEC); - nn = new_chunk(cl+1); - nmax = 0; - for (j = 1; j <= cl; j++) - { - C[j] = mael(dtcr, j, 2); - c2[j] = ldivsg(2, (GEN)C[j]); - cond[j] = mael(dtcr, j, 7); - nn[j] = (long)(bit_accuracy(prec) * gtodouble((GEN)C[j]) * 0.35); - - nmax = max(nmax, nn[j]); - } - - if (nmax >= VERYBIGINT) - err(talker, "Too many coefficients (%ld) in QuadGetST: computation impossible", nmax); - - if (DEBUGLEVEL >= 2) - fprintferr("nmax = %ld\n", nmax); - - /* compute the coefficients */ - matan = computean(dtcr, nmax, prec); - if (DEBUGLEVEL) msgtimer("Compute an"); - - /* allocate memory for the answer */ + /* allocate memory for answer */ rep = cgetg(3, t_VEC); - - /* allocate memory for veclprime1 */ - veclprime1 = cgetg(cl + 1, t_VEC); + S = cgetg(cl+1, t_VEC); rep[1] = (long)S; + T = cgetg(cl+1, t_VEC); rep[2] = (long)T; for (j = 1; j <= cl; j++) { - v = cgetg(3, t_COMPLEX); - v[1] = lgetr(prec); - v[2] = lgetr(prec); gaffect(gzero, v); - veclprime1[j] = (long)v; + S[j] = (long)cgetc(prec); + T[j] = (long)cgetc(prec); } + av = avma; - av1 = avma; - cn = cgetr(prec); - p1 = gmul2n(cf, -1); - - /* compute veclprime1 */ - for (j = 1; j <= cl; j++) + /* initializations */ + degs = GetDeg(dataCR); + ncond = lg(vChar)-1; + C = cgetg(ncond+1, t_VEC); + N0 = cgetg(ncond+1, t_VECSMALL); + nmax = 0; + for (j = 1; j <= ncond; j++) { - long n0 = 0; - p2 = gmael3(dtcr, j, 5, 2); - cexp = gexp(gneg_i((GEN)c2[j]), prec); - av2 = avma; affsr(1, cn); v = (GEN)veclprime1[j]; - for (n = 1; n <= nn[j]; n++) - if ( (an = EvalCoeff(p2, matan[j][n], degs[j])) ) - { - affrr(gmul(cn, gpowgs(cexp, n - n0)), cn); - n0 = n; - gaffect(gadd(v, gmul(divrs(cn,n), an)), v); - avma = av2; - } - gaffect(gmul(p1, gmul(v, (GEN)C[j])), v); - avma = av2; + C[j] = mael(dataCR, mael(vChar,j,1), 2); + N0[j] = (long)(bit_accuracy(prec) * gtodouble((GEN)C[j]) * 0.35); + if (nmax < N0[j]) nmax = N0[j]; } - avma = av1; - rep[1] = (long)veclprime1; - if (DEBUGLEVEL) msgtimer("Compute V1"); + if ((ulong)nmax > maxprime()) + err(talker, "Not enough precomputed primes (need all p <= %ld)", nmax); + if (DEBUGLEVEL>1) fprintferr("nmax = %ld\n", nmax); + InitPrimesQuad(gmael(dataCR,1,4), nmax, &LIST); - /* allocate memory for veclprime2 */ - veclprime2 = cgetg(cl + 1, t_VEC); - for (j = 1; j <= cl; j++) - { - v = cgetg(3, t_COMPLEX); - v[1] = lgetr(prec); - v[2] = lgetr(prec); gaffect(gzero, v); - veclprime2[j] = (long)v; - } + cf = gmul2n(mpsqrt(mppi(prec)), 1); + cfh = gmul2n(cf, -1); - /* compute f1(C/n) */ av1 = avma; - - veint1 = cgetg(cl + 1, t_VEC); - for (j = 1; j <= cl; j++) + /* loop over conductors */ + for (j = 1; j <= ncond; j++) { - p1 = NULL; - for (k = 1; k < j; k++) - if (gegal((GEN)cond[j], (GEN)cond[k])) { p1 = (GEN)veint1[k]; break; } - if (p1 == NULL) - { - p1 = veceint1((GEN)c2[j], stoi(nn[j]), prec); - veint1[j] = (long)p1; - } - av2 = avma; p2 = gmael3(dtcr, j, 5, 2); - v = (GEN)veclprime2[j]; - for (n = 1; n <= nn[j]; n++) - if ( (an = EvalCoeff(p2, matan[j][n], degs[j])) ) - { - gaffect(gadd(v, gmul((GEN)p1[n], an)), v); - avma = av2; - } - gaffect(gmul(cf, gconj(v)), v); - avma = av2; - } - avma = av1; - rep[2] = (long)veclprime2; - if (DEBUGLEVEL) msgtimer("Compute V2"); - FreeMat(matan); return gerepileupto(av, rep); -} + const GEN c1 = (GEN)C[j], c2 = divsr(2,c1), cexp = mpexp(gneg(c2)); + const GEN LChar = (GEN)vChar[j]; + const long nChar = lg(LChar)-1, NN = N0[j]; + GEN veint1, vcn = cgetg(NN+1, t_VEC); -#if 0 -/* recover a quadratic integer by an exhaustive search */ -static GEN -recbeta2(GEN nf, GEN beta, GEN bound, long prec) -{ - long av = avma, av2, tetpil, i, range, G, e, m; - GEN om, om1, om2, dom, p1, a, b, rom, bom2, *gptr[2]; + if (DEBUGLEVEL>1) + fprintferr("* conductor no %ld/%ld (N = %ld)\n\tInit: ", j,ncond,NN); + veint1 = veceint1(c2, stoi(NN), prec); + vcn[1] = (long)cexp; + for (n=2; n<=NN; n++) vcn[n] = lmulrr((GEN)vcn[n-1], cexp); + av2 = avma; + for (n=2; n<=NN; n++, avma = av2) + affrr(divrs((GEN)vcn[n],n), (GEN)vcn[n]); - G = min( - 20, - bit_accuracy(prec) >> 4); - - if (DEBUGLEVEL > 3) - fprintferr("\n Precision needed: %ld", G); - - om = gmael(nf, 7, 2); - rom = (GEN)nf[6]; - om1 = poleval(om, (GEN)rom[2]); - om2 = poleval(om, (GEN)rom[1]); - dom = subrr(om1, om2); - - /* b will run from b to b + range */ - p1 = gaddgs(gmul2n(gceil(absr(divir(bound, dom))), 1), 2); - range = VERYBIGINT; - if (cmpis(p1, VERYBIGINT) < 0) - range = itos(p1); - - av2 = avma; - - b = gdiv(gsub(bound, beta), dom); - if (gsigne(b) < 0) - b = subis(negi(gcvtoi(gneg_i(b), &e)), 1); - else - b=gcvtoi(b, &e); - - if (e > 0) /* precision is lost in truncation */ - { - avma = av; - return NULL; - } - - bom2 = mulir(b, om2); - m = 0; - - for (i = 0; i <= range; i++) - { - /* for every b, we construct a and test it */ - a = grndtoi(gsub(beta, bom2), &e); - - if (e > 0) /* precision is lost in truncation */ + for (k = 1; k <= nChar; k++) { - avma = av; - return NULL; - } + const long t = LChar[k], d = degs[t]; + const GEN z = gmael3(dataCR, t, 5, 2); + GEN p1 = gzero, p2 = gzero; + int **matan; + long c = 0; - p1 = gsub(mpadd(a, bom2), beta); - - if ((DEBUGLEVEL > 3) && (expo(p1)1) + fprintferr("\tcharacter no: %ld (%ld/%ld)\n", t,k,nChar); + matan = computean((GEN)dataCR[t], &LIST, NN, d); + for (n = 1; n <= NN; n++) + if ((an = EvalCoeff(z, matan[n], d))) + { + p1 = gadd(p1, gmul(an, (GEN)vcn[n])); + p2 = gadd(p2, gmul(an, (GEN)veint1[n])); + if (++c == 256) + { GEN *gptr[2]; gptr[0]=&p1; gptr[1]=&p2; + gerepilemany(av2,gptr,2); c = 0; + } + } + gaffect(gmul(cfh, gmul(p1,c1)), (GEN)S[t]); + gaffect(gmul(cf, gconj(p2)), (GEN)T[t]); + FreeMat(matan,NN); avma = av2; } - - if (gcmp0(p1) || (expo(p1) < G)) /* result found */ - { - p1 = gadd(a, gmul(b, om)); - return gerepileupto(av, gmodulcp(p1, (GEN)nf[1])); - } - - tetpil = avma; - - b = gaddgs(b, 1); - bom2 = gadd(bom2, om2); - - gptr[0] = &b; - gptr[1] = &bom2; - gerepilemanysp(av2, tetpil, gptr, 2); + if (DEBUGLEVEL>1) fprintferr("\n"); + avma = av1; } - - /* if it fails... */ - return NULL; + if (DEBUGLEVEL) msgtimer("S & T"); + avma = av; return rep; } -#endif +typedef struct { + long cl; + GEN dkpow; +} DH_t; + /* return 1 if the absolute polynomial pol (over Q) defines the Hilbert class field of the real quadratic field bnf */ -int -define_hilbert(GEN bnf, GEN pol) +static GEN +define_hilbert(void *S, GEN pol) { - long cl; - GEN dk; + DH_t *T = (DH_t*)S; + GEN d = modulargcd(derivpol(pol), pol); - cl = itos(gmael3(bnf, 8, 1, 1)); - dk = gmael(bnf, 7, 3); - - if (degpol(pol) == cl) - if ((cl%2) || !egalii(discf(pol), gpowgs(dk,cl>>1))) return 1; - - return 0; + if (degpol(pol) != T->cl + degpol(d)) return NULL; + pol = gdivexact(pol, d); + return (T->cl & 1 || !egalii(smalldiscf(pol), T->dkpow))? pol: NULL; } /* let polrel define Hk/k, find L/Q such that Hk=Lk and L and k are disjoint */ static GEN -makescind(GEN bnf, GEN polabs, long cl, long prec) +makescind(GEN nf, GEN polrel, long cl) { - long av = avma, i, l; - GEN pol, p1, nf2, dabs, dk, bas; + long i, l; + gpmem_t av = avma; + GEN pol, polabs, L, BAS, nf2, dabs, dk, bas; + DH_t T; + FP_chk_fun CHECK; - /* check the result (a little): signature and discriminant */ - bas = allbase4(polabs,0,&dabs,NULL); - dk = gmael(bnf,7,3); + BAS = rnfpolredabs(nf, polrel, nf_ABSOLUTE|nf_ADDZK); + polabs = (GEN)BAS[1]; + bas = (GEN)BAS[2]; + dabs = gmul(ZX_disc(polabs), gsqr(det2(vecpol_to_mat(bas, 2*cl)))); + + /* check result (a little): signature and discriminant */ + dk = (GEN)nf[3]; if (!egalii(dabs, gpowgs(dk,cl)) || sturm(polabs) != 2*cl) err(bugparier, "quadhilbert"); /* attempt to find the subfields using polred */ - p1 = cgetg(3,t_VEC); p1[1]=(long)polabs; p1[2]=(long)bas; - pol = polredfirstpol(p1, (prec<<1) - 2, &define_hilbert, bnf); + T.cl = cl; + T.dkpow = (cl & 1) ? NULL: gpowgs(dk, cl>>1); + CHECK.f = &define_hilbert; + CHECK.data = (void*)&T; + pol = polredfirstpol(BAS, 0, &CHECK); if (DEBUGLEVEL) msgtimer("polred"); if (!pol) { - nf2 = nfinit0(polabs, 1, prec); - p1 = subfields(nf2, stoi(cl)); - l = lg(p1); + nf2 = nfinit0(BAS, 0, DEFAULTPREC); + L = subfields(nf2, stoi(cl)); + l = lg(L); if (DEBUGLEVEL) msgtimer("subfields"); for (i = 1; i < l; i++) { - pol = gmael(p1, i, 1); - if ((cl%2) || !gegal(sqri(discf(pol)), (GEN)nf2[3])) break; + pol = gmael(L, i, 1); + if (cl & 1 || !egalii(smalldiscf(pol), T.dkpow)) break; } if (i == l) for (i = 1; i < l; i++) { - pol = gmael(p1, i, 1); - if (degpol(gcoeff(nffactor(bnf, pol), 1, 1)) == cl) break; + pol = gmael(L, i, 1); + if (degpol(gcoeff(nffactor(nf, pol), 1, 1)) == cl) break; } if (i == l) err(bugparier, "makescind (no polynomial found)"); } - pol = polredabs(pol, prec); + pol = polredabs0(pol, nf_PARTIALFACT); return gerepileupto(av, pol); } @@ -2878,38 +2641,35 @@ makescind(GEN bnf, GEN polabs, long cl, long prec) static GEN GenusField(GEN bnf, long prec) { - long hk, c, l, av = avma; - GEN disc, quat, x2, pol, div, d; + long hk, c, l; + gpmem_t av = avma; + GEN disc, x2, pol, div, d; hk = itos(gmael3(bnf, 8, 1, 1)); disc = gmael(bnf, 7, 3); - quat = stoi(4); x2 = gsqr(polx[0]); - if (gcmp0(modii(disc, quat))) disc = divii(disc, quat); - + if (mod4(disc) == 0) disc = divis(disc, 4); div = divisors(disc); - c = 1; l = 0; - pol = NULL; /* gcc -Wall */ + pol = NULL; - while(l < hk) + for (c = 2; l < hk; c++) /* skip c = 1 : div[1] = gun */ { - c++; d = (GEN)div[c]; - - if (gcmp1(modii(d, quat))) + if (mod4(d) == 1) { - if (!l) - pol = gsub(x2, d); + GEN t = gsub(x2, d); + if (!pol) + pol = t; else - pol=(GEN)compositum(pol, gsub(x2, d))[1]; + pol = (GEN)compositum(pol, t)[1]; l = degpol(pol); } } - return gerepileupto(av, polredabs(pol, prec)); + return gerepileupto(av, polredabs0(pol, nf_PARTIALFACT)); } /* if flag = 0 returns the reduced polynomial, flag = 1 returns the @@ -2919,86 +2679,83 @@ GenusField(GEN bnf, long prec) static GEN AllStark(GEN data, GEN nf, long flag, long newprec) { - long cl, i, j, cpt = 0, av, av2, N, h, v, n, bnd = 300, sq = 1, r1, r2; - int ***matan; - GEN p0, p1, p2, S, T, polrelnum, polrel, Lp, W, A, veczeta, sig, valchi; - GEN degs, ro, C, Cmax, dataCR, cond1, L1, *gptr[2], an, Pi; + long cl, i, j, cpt = 0, N, h, v, n, bnd = 300, sq = 1, r1, r2; + gpmem_t av, av2; + int **matan; + GEN p0, p1, p2, S, T, polrelnum, polrel, Lp, W, A, veczeta, sig; + GEN vChar, degs, ro, C, Cmax, dataCR, cond1, L1, *gptr[2], an, Pi; + LISTray LIST; N = degpol(nf[1]); r1 = nf_get_r1(nf); r2 = (N - r1)>>1; cond1 = gmael4(data, 1, 2, 1, 2); Pi = mppi(newprec); + dataCR = (GEN)data[5]; v = 1; while(gcmp1((GEN)cond1[v])) v++; -LABDOUB: - - av = avma; - - dataCR = (GEN)data[5]; cl = lg(dataCR)-1; degs = GetDeg(dataCR); h = itos(gmul2n(det((GEN)data[2]), -1)); +LABDOUB: + + av = avma; + vChar = sortChars(dataCR, N==2); + W = ComputeAllArtinNumbers(dataCR, vChar, (flag >= 0), newprec); if (flag >= 0) { - /* compute S,T differently if nf is quadratic */ - if (N == 2) - p1 = QuadGetST(data, newprec); - else - p1 = GetST(dataCR, newprec); - + p1 = GetST(dataCR, vChar, newprec); S = (GEN)p1[1]; T = (GEN)p1[2]; - Lp = cgetg(cl + 1, t_VEC); for (i = 1; i <= cl; i++) - Lp[i] = GetValue((GEN)dataCR[i], (GEN)S[i], (GEN)T[i], 0, 1, newprec)[2]; + Lp[i] = GetValue((GEN)dataCR[i], (GEN)W[i], (GEN)S[i], (GEN)T[i], + 0, 1, newprec)[2]; if (DEBUGLEVEL) msgtimer("Compute W"); } else - { - /* compute a crude approximation of the result */ + { /* compute a crude approximation of the result */ C = cgetg(cl + 1, t_VEC); for (i = 1; i <= cl; i++) C[i] = mael(dataCR, i, 2); Cmax = vecmax(C); - n = GetBoundN0(Cmax, r1, r2, newprec, 0); + n = GetBoundN0(Cmax, r1, r2, newprec); if (n > bnd) n = bnd; if (DEBUGLEVEL) fprintferr("nmax in QuickPol: %ld \n", n); + InitPrimes(gmael(dataCR,1,4), n, &LIST); - matan = ComputeCoeff(dataCR, n, newprec); - p0 = cgetg(3, t_COMPLEX); - p0[1] = lgetr(newprec); affsr(0, (GEN)p0[1]); - p0[2] = lgetr(newprec); affsr(0, (GEN)p0[2]); + p0[1] = (long)realzero(newprec); + p0[2] = (long)realzero(newprec); L1 = cgetg(cl+1, t_VEC); /* we use the formulae L(1) = sum (an / n) */ for (i = 1; i <= cl; i++) { + matan = ComputeCoeff((GEN)dataCR[i], &LIST, n, degs[i]); av2 = avma; p1 = p0; p2 = gmael3(dataCR, i, 5, 2); for (j = 1; j <= n; j++) - if ( (an = EvalCoeff(p2, matan[i][j], degs[i])) ) + if ( (an = EvalCoeff(p2, matan[j], degs[i])) ) p1 = gadd(p1, gdivgs(an, j)); L1[i] = lpileupto(av2, p1); + FreeMat(matan, n); } - FreeMat(matan); p1 = gmul2n(gpowgs(mpsqrt(Pi), N - 2), 1); Lp = cgetg(cl+1, t_VEC); for (i = 1; i <= cl; i++) { - W = ComputeArtinNumber((GEN)dataCR[i], 1, newprec); - A = (GEN)ComputeAChi((GEN)dataCR[i], 0, newprec)[2]; - W = gmul((GEN)C[i], gmul(A, W)); + GEN dtcr = (GEN)dataCR[i], WW; + A = (GEN)ComputeAChi(dtcr, 0, newprec)[2]; + WW= gmul((GEN)C[i], gmul(A, (GEN)W[i])); - Lp[i] = ldiv(gmul(W, gconj((GEN)L1[i])), p1); + Lp[i] = ldiv(gmul(WW, gconj((GEN)L1[i])), p1); } } @@ -3010,26 +2767,25 @@ LABDOUB: GEN z = gzero; sig = (GEN)p1[i]; - valchi = chiideal(dataCR, sig, 0); - for (j = 1; j <= cl; j++) { - GEN p2 = greal(gmul((GEN)Lp[j], (GEN)valchi[j])); - if (!gegal(gdeux, gmael3(dataCR, j, 5, 3))) - p2 = gmul2n(p2, 1); /* character not real */ - z = gadd(z,p2); + GEN CHI = gmael(dataCR,j,5); + GEN val = ComputeImagebyChar(CHI, sig); + GEN p2 = greal(gmul((GEN)Lp[j], val)); + if (itos((GEN)CHI[3]) != 2) p2 = gmul2n(p2, 1); /* character not real */ + z = gadd(z, p2); } veczeta[i] = ldivgs(z, 2 * h); } if (DEBUGLEVEL >= 2) fprintferr("zetavalues = %Z\n", veczeta); - if ((flag >=0) && (flag <= 3)) sq = 0; + if (flag >=0 && flag <= 3) sq = 0; ro = cgetg(h+1, t_VEC); /* roots */ - + for (;;) { - if (!sq && (DEBUGLEVEL > 1)) + if (DEBUGLEVEL > 1 && !sq) fprintferr("Checking the square-root of the Stark unit...\n"); for (j = 1; j <= h; j++) @@ -3043,10 +2799,10 @@ LABDOUB: if (DEBUGLEVEL >= 2) fprintferr("polrelnum = %Z\n", polrelnum); msgtimer("Compute %s", (flag < 0)? "quickpol": "polrelnum"); } - + if (flag < 0) return gerepilecopy(av, polrelnum); - + /* we try to recognize this polynomial */ polrel = RecCoeff(nf, polrelnum, v, newprec); @@ -3066,11 +2822,10 @@ LABDOUB: if (DEBUGLEVEL) err(warnprec, "AllStark", newprec); nf = nfnewprec(nf, newprec); - data[5] = (long)CharNewPrec((GEN)data[5], nf, newprec); + dataCR = CharNewPrec(dataCR, nf, newprec); - gptr[0] = &data; - gptr[1] = &nf; - gerepilemany(av, gptr, 2); + gptr[0] = &dataCR; + gptr[1] = &nf; gerepilemany(av, gptr, 2); goto LABDOUB; } @@ -3084,20 +2839,20 @@ LABDOUB: polrel = gmul(polrel, gpowgs(polx[n], h)); polrel = gsubst(polrel, n, polx[0]); polrel = gmul(polrel, leading_term(polrel)); - delete_var(); + (void)delete_var(); } if (DEBUGLEVEL >= 2) fprintferr("polrel = %Z\n", polrel); if (DEBUGLEVEL) msgtimer("Recpolnum"); /* we want a reduced relative polynomial */ - if (!flag) return gerepileupto(av, rnfpolredabs(nf, polrel, 0, newprec)); + if (!flag) return gerepileupto(av, rnfpolredabs(nf, polrel, 0)); /* we just want the polynomial computed */ if (flag!=2) return gerepilecopy(av, polrel); /* we want a reduced absolute polynomial */ - return gerepileupto(av, rnfpolredabs(nf, polrel, 2, newprec)); + return gerepileupto(av, rnfpolredabs(nf, polrel, nf_ABSOLUTE)); } /********************************************************************/ @@ -3109,12 +2864,12 @@ LABDOUB: GEN quadhilbertreal(GEN D, GEN flag, long prec) { - VOLATILE long av = avma, cl; - long newprec; - VOLATILE GEN pol, bnf, bnr, dataC, bnrh, nf, exp; - void *catcherr = NULL; + gpmem_t av = avma; + long cl, newprec; + GEN pol, bnf, bnr, dataC, bnrh, nf, exp; - if (DEBUGLEVEL) timer2(); + (void)≺ /* prevent longjmp clobbering it */ + if (DEBUGLEVEL) (void)timer2(); disable_dbg(0); /* quick computation of the class number */ @@ -3130,55 +2885,57 @@ quadhilbertreal(GEN D, GEN flag, long prec) pol = quadpoly(D); setvarn(pol, fetch_var()); - START: +START: /* compute the class group */ bnf = bnfinit0(pol, 1, NULL, prec); nf = (GEN)bnf[7]; disable_dbg(-1); - if (DEBUGLEVEL) msgtimer("Compute Cl(k)"); /* if the exponent of the class group is 2, use rather Genus Field Theory */ exp = gmael4(bnf, 8, 1, 2, 1); - if (gegal(exp, gdeux)) { delete_var(); return GenusField(bnf, prec); } + if (gegal(exp, gdeux)) { (void)delete_var(); return GenusField(bnf, prec); } - { /* catch precision problems (precision too small) */ - jmp_buf env; - if (setjmp(env)) + CATCH(precer) { + prec += EXTRA_PREC; pol = NULL; + err (warnprec, "quadhilbertreal", prec); + } TRY { + /* find the modulus defining N */ + bnr = buchrayinitgen(bnf, gun); + dataC = InitQuotient(bnr, gzero); + bnrh = FindModulus(dataC, 1, &newprec, prec, gcmp0(flag)? 0: -10); + + if (DEBUGLEVEL) msgtimer("FindModulus"); + + if (newprec > prec) { - prec += EXTRA_PREC; - err (warnprec, "quadhilbertreal", prec); - goto START; + if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec); + nf = nfnewprec(nf, newprec); } - catcherr = err_catch(precer, env, NULL); - } + pol = AllStark(bnrh, nf, 1, newprec); + } ENDCATCH; + if (!pol) goto START; - /* find the modulus defining N */ - bnr = buchrayinitgen(bnf, gun); - dataC = InitQuotient(bnr, gzero); - bnrh = FindModulus(dataC, 1, &newprec, prec, gcmp0(flag)? 0: -10); + pol = makescind(nf, pol, cl); + (void)delete_var(); + return gerepileupto(av, pol); +} - if (DEBUGLEVEL) msgtimer("FindModulus"); - - if (newprec > prec) - { - if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec); - nf = nfnewprec(nf, newprec); - } - - /* use the generic function AllStark */ - pol = AllStark(bnrh, nf, 2, newprec); - delete_var(); - return gerepileupto(av, makescind(bnf, pol, cl, prec)); +static GEN +get_subgroup(GEN subgp, GEN cyc) +{ + if (!subgp || gcmp0(subgp)) return cyc; + return gcmp1(denom(gauss(subgp, cyc)))? subgp: NULL; } GEN -bnrstark(GEN bnr, GEN subgroup, long flag, long prec) +bnrstark(GEN bnr, GEN subgrp, long flag, long prec) { - long cl, N, newprec, av = avma, bnd = 0; + long cl, N, newprec, bnd = 0; + gpmem_t av = avma; GEN bnf, dataS, p1, Mcyc, nf, data; - if (flag >= 4) + if (flag >= 4) { bnd = -10; flag -= 4; @@ -3188,59 +2945,52 @@ bnrstark(GEN bnr, GEN subgroup, long flag, long pre /* check the bnr */ checkbnrgen(bnr); - - bnf = (GEN)bnr[1]; - nf = (GEN)bnf[7]; - Mcyc = diagonal(gmael(bnr, 5, 2)); + bnf = checkbnf(bnr); + nf = checknf(bnf); N = degpol(nf[1]); - if (N == 1) - err(talker, "the ground field must be distinct from Q"); + if (N == 1) return galoissubcyclo(bnr, subgrp, 0, 0); + + Mcyc = diagonal(gmael(bnr, 5, 2)); /* check the bnf */ - if (!varn(gmael(bnf, 7, 1))) + if (!varn(nf[1])) err(talker, "main variable in bnrstark must not be x"); - if (cmpis(gmael3(bnf, 7, 2, 1), N)) + if (nf_get_r2(nf)) err(talker, "not a totally real ground base field in bnrstark"); - /* check the subgroup */ - if (gcmp0(subgroup)) - subgroup = Mcyc; - else - { - p1 = gauss(subgroup, Mcyc); - if (!gcmp1(denom(p1))) - err(talker, "incorrect subgroup in bnrstark"); - } + /* check the subgrp */ + if (! (subgrp = get_subgroup(subgrp,Mcyc)) ) + err(talker, "incorrect subgrp in bnrstark"); /* compute bnr(conductor) */ - p1 = conductor(bnr, subgroup, 2); + p1 = conductor(bnr, subgrp, 2); bnr = (GEN)p1[2]; - subgroup = (GEN)p1[3]; + subgrp = (GEN)p1[3]; /* check the class field */ if (!gcmp0(gmael3(bnr, 2, 1, 2))) err(talker, "not a totally real class field in bnrstark"); - cl = itos(det(subgroup)); + cl = itos(det(subgrp)); if (cl == 1) return polx[0]; - timer2(); + if (DEBUGLEVEL) (void)timer2(); /* find a suitable extension N */ - dataS = InitQuotient(bnr, subgroup); + dataS = InitQuotient(bnr, subgrp); data = FindModulus(dataS, 1, &newprec, prec, bnd); if (newprec > prec) { - if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec); + if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec); nf = nfnewprec(nf, newprec); } return gerepileupto(av, AllStark(data, nf, flag, newprec)); } -/* For each character of Cl(bnr)/sbgrp, compute L(1, chi) (or equivalently +/* For each character of Cl(bnr)/subgp, compute L(1, chi) (or equivalently the first non-zero term c(chi) of the expansion at s = 0). The binary digits of flag mean 1: if 0 then compute the term c(chi) and return [r(chi), c(chi)] where r(chi) is the order of L(s, chi) at s = 0, @@ -3251,13 +3001,14 @@ bnrstark(GEN bnr, GEN subgroup, long flag, long pre the modulus of bnr (and the infinite places), 3: return also the character */ GEN -bnrL1(GEN bnr, GEN sbgrp, long flag, long prec) +bnrL1(GEN bnr, GEN subgp, long flag, long prec) { GEN bnf, nf, cyc, Mcyc, p1, L1, chi, lchi, clchi, allCR, listCR, dataCR; - GEN S, T, rep, indCR, invCR, Qt; - long N, cl, i, j, k, nc, lq, a, av = avma, ncc; + GEN S, T, rep, indCR, invCR, Qt, vChar; + long N, cl, i, j, k, nc, lq, a, ncc; + gpmem_t av = avma; - checkbnr(bnr); + checkbnrgen(bnr); bnf = (GEN)bnr[1]; nf = (GEN)bnf[7]; cyc = gmael(bnr, 5, 2); @@ -3265,15 +3016,9 @@ bnrL1(GEN bnr, GEN sbgrp, long flag, long prec) ncc = lg(cyc) - 1; N = degpol(nf[1]); - if (N == 1) - err(talker, "the ground field must be distinct from Q"); + if (N == 1) err(talker, "the ground field must be distinct from Q"); + if (flag < 0 || flag > 8) err(flagerr,"bnrL1"); - if ((flag < 0) || (flag > 8)) - err(flagerr,"bnrL1"); - - /* check the bnr */ - checkbnrgen(bnr); - /* compute bnr(conductor) */ if (!(flag & 2)) { @@ -3284,24 +3029,15 @@ bnrL1(GEN bnr, GEN sbgrp, long flag, long prec) } /* check the subgroup */ - if (gcmp0(sbgrp)) - sbgrp = Mcyc; - else - { - if (lg(sbgrp) != ncc+1) - err(talker, "incorrect subgroup in bnrL1"); - p1 = gauss(sbgrp, Mcyc); - if (!gcmp1(denom(p1))) - err(talker, "incorrect subgroup in bnrL1"); - } + if (! (subgp = get_subgroup(subgp,Mcyc)) ) + err(talker, "incorrect subgroup in bnrL1"); - cl = labs(itos(det(sbgrp))); - Qt = InitQuotient0(Mcyc, sbgrp); + cl = labs(itos(det(subgp))); + Qt = InitQuotient0(Mcyc, subgp); lq = lg((GEN)Qt[2]) - 1; - /* compute all the characters */ - allCR = FindEltofGroup(cl, (GEN)Qt[2]); - + /* compute all characters */ + allCR = EltsOfGroup(cl, (GEN)Qt[2]); /* make a list of all non-trivial characters modulo conjugation */ listCR = cgetg(cl, t_VEC); @@ -3316,11 +3052,11 @@ bnrL1(GEN bnr, GEN sbgrp, long flag, long prec) /* lift the character to a character on Cl(bnr) */ lchi = cgetg(ncc + 1, t_VEC); - for (j = 1; j <= ncc; j++) + for (j = 1; j <= ncc; j++) { p1 = gzero; for (k = 1; k <= lq; k++) - p1 = gadd(p1, gdiv(mulii(gmael3(Qt, 3, j, k), (GEN)chi[k]), + p1 = gadd(p1, gdiv(mulii(gmael3(Qt, 3, j, k), (GEN)chi[k]), gmael(Qt, 2, k))); lchi[j] = lmodii(gmul(p1, (GEN)cyc[j]), (GEN)cyc[j]); } @@ -3355,8 +3091,8 @@ bnrL1(GEN bnr, GEN sbgrp, long flag, long prec) /* compute the data for these characters */ dataCR = InitChar(bnr, listCR, prec); - p1 = GetST(dataCR, prec); - + vChar= sortChars(dataCR,0); + p1 = GetST(dataCR, vChar, prec); S = (GEN)p1[1]; T = (GEN)p1[2]; @@ -3369,8 +3105,8 @@ bnrL1(GEN bnr, GEN sbgrp, long flag, long prec) { a = indCR[i]; if (a > 0) - L1[i] = (long)GetValue((GEN)dataCR[a], (GEN)S[a], (GEN)T[a], flag & 1, - flag & 2, prec); + L1[i] = (long)GetValue((GEN)dataCR[a], NULL, (GEN)S[a], (GEN)T[a], + flag & 1, flag & 2, prec); } for (i = 1; i < cl; i++)