File: [local] / OpenXM_contrib / pari / src / modules / Attic / stark.c (download)
Revision 1.1.1.1 (vendor branch), Sun Jan 9 17:35:33 2000 UTC (24 years, 5 months ago) by maekawa
Branch: PARI_GP
CVS Tags: maekawa-ipv6, VERSION_2_0_17_BETA, RELEASE_20000124, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3, RELEASE_1_1_2
Changes since 1.1: +0 -0
lines
Import PARI/GP 2.0.17 beta.
|
/*******************************************************************/
/* */
/* COMPUTATION OF STARK UNITS */
/* OF TOTALLY REAL FIELDS */
/* */
/*******************************************************************/
/* $Id: stark.c,v 1.2 1999/09/23 17:50:57 karim Exp $ */
#include "pari.h"
#define EXTRA_PREC (DEFAULTPREC-1)
GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus);
static int*** computean(GEN dtcr, long nmax, long prec);
/********************************************************************/
/* 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 */
static GEN
ComputeImagebyChar(GEN chi, GEN logelt, long flag)
{
GEN gn = gmul((GEN)chi[1], logelt), x = (GEN)chi[flag? 4: 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)
{
d /= 2;
if (n >= d) return gneg(gpowgs(x, n-d));
}
return gpowgs(x, n);
}
/* Compute the conjugate character */
static GEN
ConjChar(GEN chi, GEN cyc)
{
long i, l = lg(chi);
GEN p1 = cgetg(l, t_COL);
for (i = 1; i < l; i++)
if (!signe((GEN)chi[i]))
p1[i] = zero;
else
p1[i] = lsubii((GEN)cyc[i], (GEN)chi[i]);
return p1;
}
/* Compute all the elements of a group given by its SNF */
static GEN
FindEltofGroup(long order, GEN cyc)
{
long l, i, adec, j, dj;
GEN rep, p1;
l = lg(cyc)-1;
rep = cgetg(order + 1, t_VEC);
for (i = 1; i <= order; i++)
{
p1 = cgetg(l + 1, t_COL);
rep[i] = (long)p1;
adec = i;
for (j = l; j; j--)
{
dj = itos((GEN)cyc[j]);
p1[j] = lstoi(adec%dj);
adec /= dj;
}
}
return rep;
}
/* Let dataC as given by InitQuotient0, compute a system of
representatives of the quotient */
static GEN
ComputeLift(GEN dataC)
{
long order, i, av = avma;
GEN cyc, surj, eltq, elt;
order = itos((GEN)dataC[1]);
cyc = (GEN)dataC[2];
surj = (GEN)dataC[3];
eltq = FindEltofGroup(order, cyc);
elt = cgetg(order + 1, t_VEC);
for (i = 1; i <= order; i++)
elt[i] = (long)inverseimage(surj, (GEN)eltq[i]);
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)
{
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;
}
/* 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 */
static GEN
get_Char(GEN chi, 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;
}
/* Let chi a character defined over bnr and primitif over bnrc,
compute the corresponding primitive character and the vectors of
prime ideals dividing bnr but not bnr. 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;
GEN prdiff, Mrc;
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);
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];
chic = cgetg(l, t_VEC);
for (i = 1; i < l; i++)
{
s = gzero; p1 = (GEN)U[i + nbg];
for (j = 1; j <= nbg; j++)
{
p2 = gdiv((GEN)p1[j], (GEN)cyc[j]);
s = gadd(s,gmul(p2,(GEN)chi[j]));
}
chic[i] = (long)s;
}
p2 = (GEN)idealfactor(nf, cond)[1];
l = lg(p2);
prdiff = cgetg(l, t_COL);
for (nd=1, i=1; i < l; i++)
if (!idealval(nf, condc, (GEN)p2[i])) prdiff[nd++] = p2[i];
setlg(prdiff, nd);
p1 = cgetg(3, t_VEC);
p1[1] = (long)get_Char(chic,prec);
p1[2] = lcopy(prdiff);
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] = lgef(gmael4(dataCR, i, 5, 4, 1)) - 3;
return degs;
}
/********************************************************************/
/* 1rst part: find the field K */
/********************************************************************/
static GEN AllStark(GEN data, GEN nf, long flag, long prec);
static GEN InitChar0(GEN dataD, 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
the matrix expressing the generators of A on those of A / C */
static GEN
InitQuotient0(GEN DA, GEN C)
{
long i, l;
GEN MQ, MrC, H, snf, snf2, rep, p1;
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));
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[4] = lcopy(C);
return rep;
}
/* Let m be a modulus et C a subgroup of Clk(m), compute all the data
needed to work with the quotient Clk(m) / C namely 1) bnr(m), 2.1)
its order, 2.2) its structure, 2.3) the matrix Clk(m) ->> Clk(m)/C
and 2.4) the group C */
static GEN
InitQuotient(GEN bnr, GEN C)
{
GEN Mrm, dataquo = cgetg(3, t_VEC);
long av;
dataquo[1] = lcopy(bnr);
av = avma; Mrm = diagonal(gmael(bnr, 5, 2));
dataquo[2] = lpileupto(av, InitQuotient0(Mrm, C));
return dataquo;
}
/* 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 */
static GEN
ComputeKernel0(GEN P, GEN DA, GEN DB, long nbA, long nbB)
{
long rk, av = avma;
GEN herm, mask1, mask2, 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 */
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 */
static GEN
ComputeKernel(GEN bnrm, GEN dataC)
{
long av = avma, i, nbm, nbq;
GEN bnrn, Mrm, genm, Mrq, mgq, P;
bnrn = (GEN)dataC[1];
Mrm = diagonal(gmael(bnrm, 5, 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));
}
/* Let C a congruence group in bnr, compute its subgroups of index 2 as
subgroups of Clk(bnr) */
static GEN
ComputeIndex2Subgroup(GEN bnr, GEN C)
{
long nb, i, l, av = avma;
GEN snf, Mr, U, CU, subgrp, rep, p1;
disable_dbg(0);
Mr = diagonal(gmael(bnr, 5, 2));
snf = smith2(gmul(ginv(C), Mr));
U = ginv((GEN)snf[1]);
l = lg((GEN)snf[3]);
p1 = cgetg(l, t_VEC);
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 gerepileupto(av, gcopy(rep));
}
/* Let pr be a prime (pr may divide mod(bnr)), compute the indexes
e,f of the splitting of pr in the class field nf(bnr/subgroup) */
static GEN
GetIndex(GEN pr, GEN bnr, GEN subgroup, long prec)
{
long av = avma, v, lg, i;
GEN bnf, mod, mod0, mpr0, mpr, bnrpr, subpr, M, e, f, 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))
{
bnrpr = bnr;
subpr = subgroup;
}
else
{
bnrpr = buchrayinitgen(bnf, mpr, prec);
cycpr = gmael(bnrpr, 5, 2);
M = GetSurjMat(bnr, bnrpr);
M = gmul(M, subgroup);
subpr = hnf(concatsp(M, diagonal(cycpr)));
}
/* 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])));
rep = cgetg(3, t_VEC);
rep[1] = lcopy(e);
rep[2] = lcopy(f);
return gerepileupto(av, rep);
}
/* Given a conductor and a subgroups, return the corresponding
complexity and precision required using quickpol */
static GEN
CplxModulus(GEN data, long *newprec, long prec)
{
long av = avma, pr, dprec;
GEN nf, cpl, pol, p1;
nf = gmael3(data, 1, 1, 7);
p1 = cgetg(6, t_VEC);
p1[1] = data[1];
p1[2] = data[2];
p1[3] = data[3];
p1[4] = data[4];
if (DEBUGLEVEL >= 2)
fprintferr("\nTrying modulus = %Z and subgroup = %Z\n",
mael3(p1, 1, 2, 1), (GEN)p1[2]);
dprec = DEFAULTPREC;
for (;;)
{
p1[5] = (long)InitChar0((GEN)data[3], dprec);
pol = AllStark(p1, nf, -1, dprec);
cpl = mpabs(poldisc0(pol, 0));
if (!gcmp0(cpl)) break;
if (DEBUGLEVEL >= 2) err(warnprec, "CplxModulus", dprec);
dprec++;
}
if (DEBUGLEVEL >= 2) fprintferr("cpl = %Z\n", cpl);
pr = gexpo(pol)>>TWOPOTBITS_IN_LONG;
if (pr < 0) pr = 0;
*newprec = max(prec, pr + DEFAULTPREC);
return gerepileupto(av, cpl);
}
/* Let f be a conductor without infinite part and let C be a
congruence group modulo f, compute (m,D) such that D is a
congruence group of conductor m where m is a multiple of f
divisible by all the infinite places but one, D is a subgroup of
index 2 of Im(C) in Clk(m), no prime dividing f splits in the
corresponding quadratic extension and m is of minimal norm. Return
bnr(m), D, quotient Ck(m) / D and Clk(m) / C */
/* If fl != 0, try a second modulus is the first one seems too "bad" */
static GEN
FindModulus(GEN dataC, long fl, long *newprec, long prec)
{
long n, i, narch, nbp, maxnorm, minnorm, N, nbidnn, s, c, j, nbcand;
long av = avma, av1, av0, limnorm, tetpil, first = 1, pr;
GEN bnr, rep, bnf, nf, f, arch, m, listid, idnormn, bnrm, ImC;
GEN candD, D, bpr, indpr, sgp, p1, p2, rb;
bnr = (GEN)dataC[1];
sgp = gmael(dataC, 2, 4);
bnf = (GEN)bnr[1];
nf = (GEN)bnf[7];
N = degree((GEN)nf[1]);
f = gmael3(bnr, 2, 1, 1);
rep = cgetg(6, t_VEC);
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));
bpr = (GEN)idealfactor(nf, f)[1];
nbp = lg(bpr) - 1;
indpr = cgetg(nbp + 1,t_VEC);
for (i = 1; i <= nbp; i++)
{
p1 = GetIndex((GEN)bpr[i], bnr, sgp, prec);
indpr[i] = lmulii((GEN)p1[1], (GEN)p1[2]);
}
/* Initialization of the possible infinite part */
arch = cgetg(N+1, t_VEC);
for (i = 1; i <= N; i++) arch[i] = un;
/* narch = (N == 2)? 1: N; -- if N=2, only one case is necessary */
narch = N;
m = cgetg(3, t_VEC);
m[2] = (long)arch;
/* we go from minnorm up to maxnorm, if necessary we increase these values.
If we cannot find a suitable conductor of norm < limnorm, we stop */
maxnorm = 50;
minnorm = 1;
limnorm = 200;
if (DEBUGLEVEL >= 2)
fprintferr("Looking for a modulus of norm: ");
av0 = avma;
for(;;)
{
/* compute all ideals of norm <= maxnorm */
disable_dbg(0);
listid = ideallist(nf, maxnorm);
disable_dbg(-1);
av1 = avma;
for (n = minnorm; n <= maxnorm; n++)
{
if (DEBUGLEVEL >= 2) fprintferr(" %ld", n);
idnormn = (GEN)listid[n];
nbidnn = lg(idnormn) - 1;
for (i = 1; i <= nbidnn; i++)
{
tetpil = avma;
rep = gerepile(av1, tetpil, gcopy(rep));
/* finite part of the conductor */
m[1] = (long)idealmul(nf, f, (GEN)idnormn[i]);
for (s = 1; s <= narch; s++)
{
/* infinite part */
arch[N+1-s] = zero;
/* compute Clk(m), check if m is a conductor */
disable_dbg(0);
bnrm = buchrayinitgen(bnf, m, prec);
p1 = conductor(bnrm, gzero, -1, prec);
disable_dbg(-1);
if (signe(p1))
{
/* compute Im(C) in Clk(m)... */
ImC = ComputeKernel(bnrm, dataC);
/* ... 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, prec);
disable_dbg(-1);
if (signe(p1))
{
/* check the splitting of primes */
for (j = 1; j <= nbp; j++)
{
p1 = GetIndex((GEN)bpr[j], bnrm, D, prec);
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);
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);
p1 = CplxModulus(p2, &pr, prec);
if (first == 1)
{
*newprec = pr;
for (j = 1; j <= 4; j++) rep[j] = p2[j];
rep[5] = (long)p1;
}
else
if (gcmp(p1, (GEN)rep[5]) < 0)
{
*newprec = pr;
for (j = 1; j <= 5; j++) rep[j] = p2[j];
rep[5] = (long)p1;
}
if (!fl || (gcmp(p1, rb) < 0))
{
rep[5] = (long)InitChar0((GEN)rep[3], *newprec);
return gerepileupto(av, gcopy(rep));
}
if (DEBUGLEVEL >= 2)
fprintferr("Trying to find another modulus...");
first = 0;
}
}
}
}
arch[N+1-s] = un;
}
if (!first)
{
if (DEBUGLEVEL >= 2)
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 gerepileupto(av, gcopy(rep));
}
}
}
/* if necessary compute more ideals */
tetpil = avma;
rep = gerepile(av0, tetpil, gcopy(rep));
minnorm = maxnorm;
maxnorm <<= 1;
if (maxnorm > limnorm)
err(talker, "Cannot find a suitable modulus in FindModulus");
}
}
/********************************************************************/
/* 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 */
static GEN
ComputeArtinNumber(GEN datachi, long flag, long prec)
{
long av = avma, av2, G, ms, j, i, nz, zcard, q, l, N;
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 allclass, classe, bnr, beta, s, tr, p3, den, muslambda, Pi;
chi = (GEN)datachi[8];
/* trivial case */
if (cmpsi(2, (GEN)chi[3]) >= 0) return gun;
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);
N = degree((GEN)nf[1]);
Pi = mppi(prec);
G = - bit_accuracy(prec) >> 1;
nc = idealnorm(nf, cond0);
dc = idealmul(nf, diff, cond0);
den = idealnorm(nf, dc);
z = gexp(gdiv(gmul2n(gmul(gi, Pi), 1), den), prec);
q = 0;
for (i = 1; i < lg(cond1); i++)
if (gcmp1((GEN)cond1[i])) q++;
/* 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;
/* find lambda in diff.cond such that gcd(lambda.(diff.cond)^-1,cond0) = 1
and lambda >(cond1)> 0 */
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]);
idg = idealdivexact(nf, lambda, dc);
/* find mu in idg such that idh=(mu) / idg is coprime with cond0 and
mu >(cond1)> 0 */
if (!gcmp1(gcoeff(idg, 1, 1)))
{
p1 = idealfactor(nf, idg);
p2 = idealfactor(nf, cond0);
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]);
idh = idealdivexact(nf, mu, idg);
}
else
{
mu = gun;
idh = gcopy(idg);
}
muslambda = element_div(nf, mu, lambda);
/* compute a system of generators of (Ok/cond)^* cond1-positive */
zid = zidealstarinitgen(nf, cond0);
zcard = itos(gmael(zid, 2, 1));
zstruc = gmael(zid, 2, 2);
zgen = gmael(zid, 2, 3);
nz = lg(zgen) - 1;
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;
}
/* Sum chi(beta) * exp(2i * Pi * Tr(beta * mu / lambda) where beta
runs through the classes of (Ok/cond0)^* and beta cond1-positive */
allclass = FindEltofGroup(zcard, zstruc);
p3 = cgetg(N + 1, t_COL);
for (i = 1; i <= N; i++) p3[i] = zero;
vt = cgetg(N + 1, t_VEC);
for (i = 1; i <= N; i++)
{
p3[i] = un;
vt[i] = ltrace(basistoalg(nf, p3));
p3[i] = zero;
}
s = cgetg(3, t_COMPLEX);
s[1] = lgetr(prec);
s[2] = lgetr(prec);
gaffect(gzero, s);
av2 = avma;
for (i = 1; i <= zcard; i++)
{
beta = gun;
chib = gun;
p1 = (GEN)allclass[i];
for (j = 1; j <= nz; j++)
{
p2 = element_powmodideal(nf, (GEN)zgen[j], (GEN)p1[j], cond0);
beta = element_mul(nf, beta, p2);
chib = gmul(chib, powgi((GEN)zchi[j], (GEN)p1[j]));
}
sg = zsigne(nf, beta, cond1);
p2 = lift(gmul(Msign, sg));
for (j = 1; j <= ms; j++)
if (gcmp1((GEN)p2[j]))
beta = element_mul(nf, beta, (GEN)elts[j]);
beta = element_mul(nf, beta, muslambda);
tr = gmul(vt, beta);
tr = gmod(gmul(tr, den), den);
gaffect(gadd(s, gmul(chib, powgi(z,tr))), s);
avma = av2;
}
classe = isprincipalray(bnr, idh);
s = gmul(s, ComputeImagebyChar(chi, classe, 0));
s = gdiv(s, gsqrt(nc, prec));
p1 = gsubgs(gabs(s, prec), 1);
i = expo(p1);
if ((i > G) && !flag)
err(bugparier, "ComputeArtinNumber");
return gerepileupto(av, gmul(s, gpowgs(gneg_i(gi),q)));
}
/* compute the constant W of the functional equation of
Lambda(chi). If flag = 1 then chi is assumed to be primitive */
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;
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))
err(talker, "incorrect character in bnrrootnumber");
if (!flag)
{
condc = bnrconductorofchar(bnr, chi, prec);
if (gegal(cond, condc)) flag = 1;
}
else condc = cond;
if (flag)
bnrc = bnr;
else
bnrc = buchrayinitgen((GEN)bnr[1], condc, prec);
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));
}
/********************************************************************/
/* 3rd part: initialize the characters */
/********************************************************************/
static GEN
LiftChar(GEN cyc, GEN Mat, GEN chi)
{
long lm, l, i, j, av;
GEN lchi, s;
lm = lg(cyc) - 1;
l = lg(chi) - 1;
lchi = cgetg(lm + 1, t_COL);
for (i = 1; i <= lm; i++)
{
av = avma;
s = gzero;
for (j = 1; j <= l; j++)
s = gadd(s, gmul((GEN)chi[j], gcoeff(Mat, j, i)));
lchi[i] = (long)gerepileupto(av, gmod(gmul(s, (GEN)cyc[i]), (GEN)cyc[i]));
}
return lchi;
}
/* Let chi be a character, A(chi) corresponding to the primes dividing diff
at s = flag. If s = 0, returns [r, A] where r is the order of vanishing
at s = 0 corresponding to diff. No Garbage collector */
static GEN
ComputeAChi(GEN dtcr, long flag, long prec)
{
long l, i;
GEN p1, ray, r, A, rep, diff, chi, bnrc;
diff = (GEN)dtcr[6];
bnrc = (GEN)dtcr[3];
chi = (GEN)dtcr[8];
l = lg(diff) - 1;
A = gun;
r = gzero;
for (i = 1; i <= l; i++)
{
ray = isprincipalray(bnrc, (GEN)diff[i]);
p1 = ComputeImagebyChar(chi, ray, 0);
if (flag)
A = gmul(A, gsub(gun, gdiv(p1, idealnorm((GEN)bnrc[1], (GEN)diff[i]))));
else
{
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));
}
}
if (flag) return A;
rep = cgetg(3, t_VEC);
rep[1] = (long)r;
rep[2] = (long)A;
return rep;
}
/* Given a list [chi, cond(chi)] of characters over Cl(bnr), compute a
vector dataCR containing for each character:
1: chi
2: the constant C(chi)
3: bnr(cond(chi))
4: bnr(m)
5: [(c_i), z, d, pm] in bnr(m)
6: diff(chi) primes dividing m but not cond(chi)
7: finite part of cond(chi)
8: [(c_i), z, d, pm] in bnr(cond(chi))
9: [q, r1 - q, r2, rc] where
q = number of real places in cond(chi)
rc = max{r1 + r2 - q + 1, r2 + q} */
static GEN
InitChar(GEN bnr, GEN listCR, long prec)
{
GEN modul, bnf, dk, r1, r2, C, dataCR, chi, cond, Mr, chic;
GEN p1, p2, q;
long N, prec2, h, i, j, nbg, av = avma;
modul = gmael(bnr, 2, 1);
Mr = gmael(bnr, 5, 2);
nbg = lg(Mr) - 1;
bnf = (GEN)bnr[1];
dk = gmael(bnf, 7, 3);
r1 = gmael3(bnf, 7, 2, 1);
r2 = gmael3(bnf, 7, 2, 2);
N = degree(gmael(bnf, 7, 1));
prec2 = ((prec - 2)<<1) + EXTRA_PREC;
C = gmul2n(gsqrt(gdiv(absi(dk), gpowgs(mppi(prec2),N)), prec2), -itos(r2));
disable_dbg(0);
h = lg(listCR) - 1;
dataCR = cgetg(h + 1, t_VEC);
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] = lsub(r1, q);
p1[3] = (long)r2;
p1[4] = lmax(gaddgs(gadd((GEN)p1[2], r2), 1), gadd(r2, q));
for (i = 1; i <= h; i++)
{
GEN olddata, data = (GEN)dataCR[i];
chi = gmael(listCR, i, 1);
cond = gmael(listCR, i, 2);
/* do we already know about the invariants of chi? */
olddata = NULL;
for (j = 1; j < i; j++)
if (gegal(cond, gmael(listCR, j, 2)))
{ olddata = (GEN)dataCR[j]; break; }
/* if cond(chi) = cond(bnr) */
if (!olddata && gegal(cond, modul))
{
data[2] = lmul(C, gsqrt(det((GEN)cond[1]), prec2));
data[3] = (long)bnr;
data[6] = lgetg(1, t_VEC);
data[7] = modul[1];
data[9] = (long)p1;
olddata = data;
}
data[1] = (long)chi; /* the character */
if (!olddata)
{
data[2] = lmul(C, gsqrt(det((GEN)cond[1]), prec2));
data[3] = (long)buchrayinitgen(bnf, cond, prec);
}
else
{
data[2] = olddata[2]; /* constant C(chi) */
data[3] = olddata[3]; /* bnr(cond(chi)) */
}
data[4] = (long)bnr; /* 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)
{
data[6] = p2[2];
data[8] = p2[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] = lsubii(r1, q);
p2[3] = (long)r2;
p2[4] = lmax(addis(addii((GEN)p2[2], r2), 1), addii(r2, q));
data[9] = (long)p2;
}
else
data[9] = olddata[9];
}
disable_dbg(-1);
return gerepileupto(av, gcopy(dataCR));
}
/* compute the list of characters to consider for AllStark and
initialize the data to compute with them */
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;
Surj = gmael(dataD, 2, 3);
MrD = gmael(dataD, 2, 2);
bnr = (GEN)dataD[1];
Mr = gmael(bnr, 5, 2);
hD = itos(gmael(dataD, 2, 1));
h = hD >> 1;
lD = lg(MrD)-1;
nbg = lg(Mr) - 1;
disable_dbg(0);
listCR = cgetg(h + 1, t_VEC); /* non-conjugate characters */
nc = 1;
allCR = cgetg(h + 1, t_VEC); /* all characters, including conjugates */
tnc = 1;
p1 = FindEltofGroup(hD, MrD);
for (i = 1; tnc <= h; i++)
{
/* lift a character of D in Clk(m) */
chi = (GEN)p1[i];
for (j = 1; j <= lD; j++) chi[j] = ldiv((GEN)chi[j], (GEN)MrD[j]);
lchi = LiftChar(Mr, Surj, chi);
for (j = 1; j < tnc; j++)
if (gegal(lchi, (GEN)allCR[j])) break;
if (j != tnc) continue;
cond = bnrconductorofchar(bnr, lchi, prec);
if (gcmp0((GEN)cond[2])) continue;
/* the infinite part of chi is non trivial */
p2 = cgetg(3, t_VEC);
p2[1] = (long)lchi;
p2[2] = (long)cond;
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))
allCR[tnc++] = (long)ConjChar(lchi, Mr);
}
setlg(listCR, nc);
disable_dbg(-1);
return gerepileupto(av, InitChar(bnr, listCR, prec));
}
/* recompute dataCR with the new precision */
static GEN
CharNewPrec(GEN dataCR, GEN nf, long prec)
{
GEN dk, C, p1, Pi;
long av = avma, N, l, j, prec2;
dk = (GEN)nf[3];
N = degree((GEN)nf[1]);
l = lg(dataCR) - 1;
prec2 = ((prec - 2)<<1) + EXTRA_PREC;
Pi = mppi(prec2);
C = gsqrt(gdiv(dk, gpowgs(Pi, N)), prec2);
for (j = 1; j <= l; j++)
{
mael(dataCR, j, 2) = lmul(C, gsqrt(det(gmael(dataCR, j, 7)), prec2));
mael4(dataCR, j, 3, 1, 7) = lcopy(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);
}
return gerepileupto(av, gcopy(dataCR));
}
/********************************************************************/
/* 4th part: compute the coefficients an(chi) */
/* */
/* matan entries are arrays of ints containing the coefficients of */
/* an(chi) as a polmod modulo cyclo(order(chi)) */
/********************************************************************/
static void
_0toCoeff(int *rep, long dg)
{
long i;
for (i=0; i<dg; i++) rep[i] = 0;
}
/* transform a polmod into coeff */
static void
Polmod2Coeff(int *rep, GEN polmod, long dg)
{
GEN pol = (GEN)polmod[2];
long i,d = lgef(pol)-3;
pol += 2;
for (i=0; i<=d; i++) rep[i] = itos((GEN)pol[i]);
for ( ; i<dg; i++) rep[i] = 0;
}
/* initialize a cl x nmax x [degs[1], ..., degs[cl] matrix of ints */
/* modified to allocate ints and pointers separately since this used to
break on 64bit platforms --GN1999Sep01 */
static int***
InitMatAn(long cl, long nmax, GEN degs, long flag)
{
long si,dg,i,j,k, n = nmax+1;
int *c, **pj, ***A;
for (si=0, i=1; i <= cl; i++)
{
dg = degs[i];
si += dg;
}
A = (int***)gpmalloc((cl+1)*sizeof(int**) + cl*n*sizeof(int*));
c = (int*)gpmalloc(si*n*sizeof(int));
A[0] = (int**)c; /* keep it around for FreeMat() */
pj = (int**)(A + (cl+1));
for (i = 1; i <= cl; i++, pj+=n)
{
A[i] = pj; dg = degs[i];
for (j = 0; j < n; j++,c+=dg)
{
pj[j] = c;
c[0] = (j == 1 || flag);
for (k = 1; k < dg; k++) c[k] = 0;
}
}
return A;
}
static void
FreeMat(int ***m)
{
free((void *)(m[0]));
free((void *)m);
}
/* initialize coeff reduction */
/* similar change --GN1999Sep01 */
static int***
InitReduction(GEN dataCR, GEN degs)
{
long av = avma,si,sp,dg,i,j, cl = lg(dataCR)-1;
int *c, **pj, ***A;
GEN d,polmod,pol, x = polx[0];
for (si=sp=0, i=1; i <= cl; i++)
{
dg = degs[i];
sp += dg;
si += dg*dg;
}
A = (int***)gpmalloc((cl+1)*sizeof(int**) + sp*sizeof(int*));
c = (int*)gpmalloc(si*sizeof(int));
A[0] = (int**)c; /* keep it around for FreeMat() */
pj = (int**)(A + (cl+1));
for (i = 1; i <= cl; i++)
{
A[i] = pj;
d = gmael3(dataCR, i, 5, 3);
pol = cyclo(itos(d), 0);
dg = degs[i]; /* degree(pol) */
for (j = 0; j < dg; j++,c+=dg)
{
pj[j] = c;
polmod = gmodulcp(gpowgs(x, dg + j), pol);
Polmod2Coeff(c, polmod, dg);
}
pj += dg;
}
avma = av; return A;
}
#if 0
static void
pan(int ***an,long cl, long nmax, GEN degs)
{
long i,j,k;
for (i = 1; i <= cl; i++)
{
long dg = degs[i];
for (j = 0; j <= nmax; j++)
for (k = 0; k < dg; k++) fprintferr("%d ",an[i][j][k]);
}
}
#endif
/* multiply (with reduction) a polmod with a coeff. put result in c1 */
static void
MulPolmodCoeff(GEN polmod, int* c1, int** reduc, long dg)
{
long av,i,j;
int c, *c2, *c3;
if (gcmp1(polmod)) return;
for (i = 0; i < dg; i++)
if (c1[i]) break;
if (i == dg) return;
av = avma;
c3 = (int*)new_chunk(2*dg);
c2 = (int*)new_chunk(dg);
Polmod2Coeff(c2,polmod, dg);
for (i = 0; i < 2*dg; i++)
{
c = 0;
for (j = 0; j <= i; j++)
if (j < dg && j > i - dg) c += c1[j] * c2[i-j];
c3[i] = c;
}
c2 = c1;
for (i = 0; i < dg; i++)
{
c = c3[i];
for (j = 0; j < dg; j++) c += reduc[j][i] * c3[dg+j];
c2[i] = c;
}
for ( ; i < dg; i++) c2[i] = 0;
avma = av;
}
/* c0 <- c0 + c2 * c1 */
static void
AddMulCoeff(int *c0, int *c2, int* c1, int** reduc, long dg)
{
long av,i,j;
int c, *c3;
if (!c2) /* c2 == 1 */
{
for (i = 0; i < dg; i++) c0[i] += c1[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++)
{
c = 0;
for (j = 0; j <= i; j++)
if (j < dg && j > i - dg) c += c1[j] * c2[i-j];
c3[i] = c;
}
for (i = 0; i < dg; i++)
{
c = c0[i] + c3[i];
for (j = 0; j < dg; j++) c += reduc[j][i] * c3[dg+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)
{
long i,j;
GEN e, r;
#if 0
/* standard Horner */
e = stoi(c[dg - 1]);
for (i = dg - 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 (j=i; c[j] == 0; j--)
if (j==0)
{
if (!e) return NULL;
if (i!=j) z = gpuigs(z,i-j+1);
return gmul(e,z);
}
if (e)
{
r = (i==j)? z: gpuigs(z,i-j+1);
e = gadd(gmul(e,r), stoi(c[j]));
}
else
e = stoi(c[j]);
}
#endif
return e;
}
/* copy the n * m array matan */
static void
CopyCoeff(int*** a, int*** a2, long n, long m, GEN degs)
{
long i,j,k;
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];
}
}
return;
}
/* initialize the data for GetRay */
static GEN
InitGetRay(GEN bnr, long nmax)
{
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] != NULL iff the field is totally real */
if (!cmpsi(degree(gmael(bnf, 7, 1)), gmael3(bnf, 7, 2, 1)))
rep[3] = un;
else
rep[3] = 0;
return rep;
}
/* compute the class of the prime ideal pr in cl(bnr) using dataray */
static GEN
GetRay(GEN bnr, GEN dataray, GEN pr, long prec)
{
long av = avma, N, n, bd, c;
GEN id, tid, t2, u, alpha, p1, cl, listid, listcl, nf, cond;
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 = degree((GEN)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++)
{
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));
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++;
}
}
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)
{
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;
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);
degs = cgetg(2, t_VECSMALL); degs[1] = dg;
an2 = InitMatAn(1, nmax, degs, 0); an1 = an2[1];
c = (int*)new_chunk(dg);
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];
ray = isprincipalray(bnrc, pr);
ki = ComputeImagebyChar(chi, ray, 1);
ki2 = gcopy(ki);
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);
}
avma = av1;
}
FreeMat(an2); avma = av;
}
/* compute the coefficients an in the general case */
static int***
ComputeCoeff(GEN dataCR, long nmax, long prec)
{
long cl, i, j, av = avma, av2, np, q, limk, k, id, cpt = 10, dg;
int ***matan, ***reduc, ***matan2, *c2;
GEN c, degs, tabprem, bnf, pr, cond, ray, ki, ki2, prime, npg, bnr, dataray;
byteptr dp = diffptr;
bnr = gmael(dataCR, 1, 4);
bnf = (GEN)bnr[1];
cond = gmael3(bnr, 2, 1, 1);
cl = lg(dataCR) - 1;
dataray = InitGetRay(bnr, nmax);
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]);
if (DEBUGLEVEL > 1) fprintferr("p = ");
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;
CopyCoeff(matan, matan2, cl, nmax, degs);
ray = GetRay(bnr, dataray, pr, prec);
ki = chiideal(dataCR, ray, 1);
ki2 = dummycopy(ki);
for (q = np; q <= nmax; q *= 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]);
}
}
}
avma = av2;
prime[2] += (*dp++);
if (!*dp) err(primer1);
if (DEBUGLEVEL > 1 && prime[2] > cpt)
{ fprintferr("%ld ", prime[2]); cpt += 10; }
}
if (DEBUGLEVEL > 1) fprintferr("\n");
for (i = 1; i <= cl; i++)
CorrectCoeff((GEN)dataCR[i], matan[i], reduc[i], nmax, degs[i]);
FreeMat(matan2); FreeMat(reduc);
avma = av; return matan;
}
/********************************************************************/
/* 5th part: compute L functions at s=1 */
/********************************************************************/
/* if flag != 0, prec means decimal digits */
static GEN
get_limx(long N, long prec, GEN *pteps, long flag)
{
GEN gN, mu, alpha, beta, eps, A0, c1, c0, c2, lneps, limx, Pi = mppi(prec);
gN = stoi(N);
mu = gmul2n(gN, -1);
alpha = gmul2n(stoi(N + 3), -1);
beta = gpui(gdeux, gmul2n(gN, -1), 3);
if (flag)
*pteps = eps = gmul2n(gpowgs(dbltor(10.), -prec), -1);
else
*pteps = eps = gmul2n(gpowgs(dbltor(10.), (long)(-(prec-2) / pariK1)), -1);
A0 = gmul2n(gun, N);
A0 = gmul(A0, gpowgs(mu, N + 2));
A0 = gmul(A0, gpowgs(gmul2n(Pi, 1), 1 - N));
A0 = gsqrt(A0, 3);
c1 = gmul(mu, gpui(beta, ginv(mu), 3));
c0 = gdiv(gmul(A0, gpowgs(gmul(gdeux, Pi), N - 1)), mu);
c0 = gmul(c0, gpui(c1, gsub(gun, alpha), 3));
c2 = gdiv(gsub(alpha, gun), mu);
lneps = glog(gdiv(c0, eps), 3);
limx = gdiv(gsub(glog(lneps, 3), glog(c1, 3)), gadd(c2, gdiv(lneps, mu)));
limx = gmul(gpui(gdiv(c1, lneps), mu, 3),
gadd(gun, gmul(c2, gmul(mu, limx))));
return limx;
}
static long
GetBoundN0(GEN C, long N, long prec, long flag)
{
long av = avma, N0;
GEN eps, limx = get_limx(N, prec, &eps, flag);
N0 = itos(gfloor(gdiv(C, limx)));
avma = av; return N0;
}
static long
GetBoundi0(long N, long prec)
{
long av = avma, imin, i0, itest;
GEN ftest, borneps, eps, limx = get_limx(N, prec, &eps, 0);
borneps = gsqrt(gmul(limx, gpowgs(mppi(3),3)), 3);
borneps = gdiv(gpowgs(stoi(5), N), gmul(eps, borneps));
imin = 1;
i0 = 1400;
while(i0 - imin >= 4)
{
itest = (i0 + imin) >> 1;
ftest = gpowgs(limx, itest);
ftest = gmul(ftest, gpowgs(mpfactr(itest / 2, 3), N));
if(gcmp(ftest, borneps) >= 0)
i0 = itest;
else
imin = itest;
}
avma = av;
return (i0 / 2) * 2;
}
/* 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)
{
GEN cst, gamun, gamdm, an, bn, cn_evn, cn_odd, x, x2, aij, p1, cf, p2;
long i, j, r, av = avma;
r = max(b + c + 1, a + c);
aij = 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);
}
x = polx[0];
x2 = gmul2n(x, -1);
/* 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);
gamun[1] = evalsigne(1) | evalvalp(0) | evalvarn(0);
gamun[2] = zero;
for (i = 1; i <= r; i++)
{
gamun[i + 2] = ldivgs((GEN)cst[i], i);
if (i%2) gamun[i + 2] = lneg((GEN)gamun[i + 2]);
}
/* the expansion of log(Gamma(s)) at s = 1/2 */
gamdm = cgetg(r + 2, 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]));
for (i = 2; i <= r; i++)
gamdm[i + 2] = lmul((GEN)gamun[i + 2], subis(shifti(gun, i), 1));
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} */
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));
}
}
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);
/* an is the expansion of Gamma(x)^{b+c} */
an = gpowgs(gdiv(gamun, x), b + c);
/* bn is the expansion of 2^{b - bx} */
bn = gpowgs(gpow(gdeux, gsubsg(1, x), prec), b);
/* 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++)
{
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)^{b+c} */
an = gdiv(an, gpowgs(gsubgs(x, 2*i + 1), b + c));
/* bn(x-s-1) = 2^b bn(x-s) */
bn = gmul2n(bn, b);
/* 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));
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), b + c));
bn = gmul2n(bn, b);
/* 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));
}
}
return gerepileupto(av, gcopy(aij));
}
static GEN
GetST(GEN dataCR, long prec)
{
GEN N0, CC, bnr, bnf, 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, N, hk, fj, id, prec2, i0, nmax;
long a, b, c, rc1, rc2, r;
int ***matan;
if (DEBUGLEVEL) timer2();
bnr = gmael(dataCR, 1, 4);
bnf = (GEN)bnr[1];
N = degree(gmael(bnf, 7, 1));
hk = lg(dataCR) - 1;
prec2 = ((prec - 2)<<1) + EXTRA_PREC;
Pi = mppi(prec2);
racpi = gsqrt(Pi, prec2);
C = cgetg(hk + 1, t_VEC);
cond = cgetg(hk + 1, t_VEC);
N0 = new_chunk(hk+1);
CC = new_chunk(hk+1);
nmax = 0;
for (i = 1; i <= hk; i++)
{
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], N, prec, 0);
if (nmax < N0[i]) nmax = N0[i];
}
i0 = GetBoundi0(N, prec);
if (nmax >= VERYBIGINT)
err(talker, "Too many coefficients (%ld) in GetST: computation impossible", nmax);
if(DEBUGLEVEL > 1) fprintferr("nmax = %ld and i0 = %ld\n", nmax, i0);
matan = ComputeCoeff(dataCR, nmax, prec);
degs = GetDeg(dataCR);
if (DEBUGLEVEL) msgtimer("Compute an");
p1 = cgetg(3, t_COMPLEX);
p1[1] = lgetr(prec2);
p1[2] = lgetr(prec2);
gaffect(gzero, p1);
S = cgetg(hk + 1, t_VEC);
T = cgetg(hk + 1, t_VEC);
for (id = 1; id <= hk; id++)
{
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;
}
powracpi = cgetg(hk + 1, t_VEC);
for (j = 1; j <= hk; j++)
powracpi[j] = (long)gpow(racpi, gmael3(dataCR, j, 9, 2), prec2);
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);
av2 = avma;
for (n = 1; n <= N0[id]; 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--)
{
p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, i, 2, j)));
p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, i, 1, j)));
}
}
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;
if (DEBUGLEVEL > 1 && n%100 == 0) fprintferr("%ld ", 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 gerepileupto(av, gcopy(rep));
}
/* Given datachi, 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)
{
GEN W, A, q, b, c, d, rchi, cf, VL, rep, racpi, nS, nT;
long av = avma;
racpi = gsqrt(mppi(prec), prec);
W = ComputeArtinNumber(datachi, 0, prec);
A = ComputeAChi(datachi, fl, prec);
d = gmael(datachi, 8, 3);
q = gmael(datachi, 9, 1);
b = gmael(datachi, 9, 2);
c = gmael(datachi, 9, 3);
rchi = addii(b, c);
if (!fl)
{
cf = gmul2n(gpow(racpi, q, 0), itos(b));
nS = gdiv(gconj(S), cf);
nT = gdiv(gconj(T), 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)
{
VL = gmul((GEN)A[2], VL);
rchi = gadd(rchi, (GEN)A[1]);
}
rep = cgetg(3, t_VEC);
rep[1] = (long)rchi;
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}) */
rep = gdiv(gadd(S, gmul(W, T)), cf);
if (cmpis(d, 3) < 0) rep = greal(rep);
if (fl2) rep = gmul(A, rep);
}
return gerepileupto(av, gcopy(rep));
}
/* return the order and the first non-zero term of L(s, chi0)
at s = 0. If flag = 1, adjust the value to get L_S(s, chi0). */
static GEN
GetValue1(GEN bnr, long flag, long prec)
{
GEN bnf, hk, Rk, wk, c, r, r1, r2, rep, mod0, diff;
long i, l, av = avma;
bnf = (GEN)bnr[1];
r1 = gmael3(bnf, 7, 2, 1);
r2 = gmael3(bnf, 7, 2, 2);
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 = subis(addii(r1, r2), 1);
if (flag)
{
mod0 = gmael3(bnr, 2, 1, 1);
diff = (GEN)idealfactor((GEN)bnf[7], mod0)[1];
l = lg(diff) - 1;
r = addis(r, l);
for (i = 1; i <= l; i++)
c = gmul(c, glog(idealnorm((GEN)bnf[7], (GEN)diff[i]), prec));
}
rep = cgetg(3, t_VEC);
rep[1] = (long)r;
rep[2] = (long)c;
return gerepileupto(av, gcopy(rep));
}
/********************************************************************/
/* 6th part: recover the coefficients */
/********************************************************************/
static long
TestOne(GEN plg, GEN beta, GEN B, long v, long G, long N)
{
long j;
GEN p1;
p1 = gsub(beta, (GEN)plg[v]);
if (expo(p1) >= G) return 0;
for (j = 1; j <= N; j++)
if (j != v)
{
p1 = gabs((GEN)plg[j], DEFAULTPREC);
if (gcmp(p1, B) > 0) return 0;
}
return 1;
}
/* Using linear dependance relations */
static GEN
RecCoeff2(GEN nf, GEN beta, GEN B, long v, long prec)
{
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 = degree(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 gerepileupto(av, gcopy(cand));
}
avma = av2;
}
return NULL;
}
/* Using Cohen's method */
static GEN
RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
{
GEN A, M, nB, cand, sol, p1, plg, B2, C2, max = stoi(10000);
GEN beta2, eps, nf2;
long N, G, i, j, k, l, ct = 0, av = avma, prec2;
N = degree((GEN)nf[1]);
G = min( - 20, - bit_accuracy(prec) >> 4);
eps = gpowgs(stoi(10), - max(8, (G >> 1)));
p1 = gceil(gdiv(glog(B, DEFAULTPREC), dbltor(2.3026)));
prec2 = max((prec << 1) - 2, (long)(itos(p1) * pariK1 + BIGDEFAULTPREC));
nf2 = nfnewprec(nf, prec2);
beta2 = gprec_w(beta, prec2);
LABrcf: ct++;
B2 = sqri(B);
C2 = gdiv(B2, gsqr(eps));
M = gmael(nf2, 5, 1);
A = cgetg(N+2, t_MAT);
for (i = 1; i <= N+1; i++)
A[i] = lgetg(N+2, t_COL);
coeff(A, 1, 1) = ladd(gmul(C2, gsqr(beta2)), B2);
for (j = 2; j <= N+1; j++)
{
p1 = gmul(C2, gmul(gneg_i(beta2), gcoeff(M, v, j-1)));
coeff(A, 1, j) = coeff(A, j, 1) = (long)p1;
}
for (i = 2; i <= N+1; i++)
for (j = 2; j <= N+1; j++)
{
p1 = gzero;
for (k = 1; k <= N; k++)
{
GEN p2 = gmul(gcoeff(M, k, j-1), gcoeff(M, k, i-1));
if (k == v) p2 = gmul(C2, p2);
p1 = gadd(p1,p2);
}
coeff(A, i, j) = coeff(A, j, i) = (long)p1;
}
nB = mulsi(N+1, B2);
cand = fincke_pohst(A, nB, max, 3, prec2, NULL);
if (!cand)
{
if (ct > 3) { avma = av; return NULL; }
prec2 = (prec2 << 1) - 2;
if (DEBUGLEVEL >= 2) err(warnprec,"RecCoeff", prec2);
nf2 = nfnewprec(nf2, prec2);
beta2 = gprec_w(beta2, prec2);
goto LABrcf;
}
cand = (GEN)cand[3];
l = lg(cand) - 1;
if (DEBUGLEVEL >= 2)
fprintferr("Found %ld vector(s) in RecCoeff3 \n", l);
sol = cgetg(N + 1, t_COL);
for (i = 1; i <= l; i++)
{
p1 = (GEN)cand[i];
if (gcmp1(mpabs((GEN)p1[1])))
{
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 gerepileupto(av, basistoalg(nf, sol));
}
}
avma = av; return NULL;
}
/* Attempts to find an algebraic integer close to beta at the place v
and less than B at all the others */
GEN
RecCoeff(GEN nf, GEN pol, long v, long prec)
{
long av = avma, j, G, cl = lgef(pol)-3;
GEN p1, beta, Bmax = stoi(10000);
/* if precision(pol) is too low, abort */
for (j = 2; j <= cl+1; j++)
{
p1 = (GEN)pol[j];
G = bit_accuracy(gprecision(p1)) - gexpo(p1);
if (G < 34) { avma = av; return NULL; }
}
pol = dummycopy(pol);
for (j = 2; j <= cl+1; j++)
{
GEN bound = binome(stoi(cl), j - 2);
bound = shifti(bound, cl + 2 - j);
if (DEBUGLEVEL > 1) fprintferr("In RecCoeff with B = %Z\n", bound);
beta = greal((GEN)pol[j]);
p1 = RecCoeff2(nf, beta, bound, v, prec);
if (!p1)
{
if (cmpii(bound, Bmax) > 0) bound = Bmax;
p1 = RecCoeff3(nf, beta, bound, v, prec);
if (!p1) return NULL;
}
pol[j] = (long)p1;
}
pol[j] = un;
return gerepileupto(av, gcopy(pol));
}
/*******************************************************************/
/*******************************************************************/
/* */
/* Computation of 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)
{
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;
av = avma;
cl = lg(dtcr) - 1;
degs = GetDeg(dtcr);
matan = InitMatAn(cl, nmax, degs, 1);
reduc = InitReduction(dtcr, degs);
bnr = gmael(dtcr, 1, 4); bnf = (GEN)bnr[1];
dataray = InitGetRay(bnr, nmax);
ideal = gmael3(bnr, 2, 1, 1);
idno = idealnorm(bnf, ideal);
dk = gmael(bnf, 7, 3);
prime = stoi(2);
dp++;
av1 = avma;
while (*dp && prime[2] <= nmax)
{
qg = prime;
al = 1;
switch (krogs(dk, prime[2]))
{
/* prime is inert */
case -1:
fldiv = divise(idno, prime);
if (!fldiv)
{
xray = GetRay(bnr, dataray, prime, prec);
chi = chiideal(dtcr, xray, 1);
chi1 = 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 || 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 */
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;
}
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;
}
/* compute S and T for the quadratic case */
static GEN
QuadGetST(GEN data, 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;
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 */
rep = cgetg(3, t_VEC);
/* allocate memory for veclprime1 */
veclprime1 = 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);
veclprime1[j] = (long)v;
}
av1 = avma;
cn = cgetr(prec);
p1 = gmul2n(cf, -1);
/* compute veclprime1 */
for (j = 1; j <= cl; 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;
}
avma = av1;
rep[1] = (long)veclprime1;
if (DEBUGLEVEL) msgtimer("Compute V1");
/* 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;
}
/* compute f1(C/n) */
av1 = avma;
veint1 = cgetg(cl + 1, t_VEC);
for (j = 1; j <= cl; 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);
}
#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];
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 */
{
avma = av;
return NULL;
}
p1 = gsub(mpadd(a, bom2), beta);
if ((DEBUGLEVEL > 3) && (expo(p1)<m))
{
m = expo(p1);
fprintferr("\n Precision found: %ld", expo(p1));
}
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 it fails... */
return NULL;
}
#endif
/* let polrel define Hk/k, find L/Q such that Hk=Lk and L and k are
disjoint */
static GEN
makescind(GEN nf, GEN polabs, long cl, long prec)
{
long av = avma, i, l;
GEN pol, p1, nf2, dabs, dk, bas;
/* check the result (a little): signature and discriminant */
bas = allbase4(polabs,0,&dabs,NULL);
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;
p1 = polred(p1, (prec<<1) - 2);
l = lg(p1);
for (i = 1; i < l; i++)
{
pol = (GEN)p1[i];
if (degree(pol) == cl)
if (cl % 2 || !gegal(sqri(discf(pol)), dabs)) break;
}
if (DEBUGLEVEL) msgtimer("polred");
/* ... if it fails, then use nfsubfields */
if (i == l)
{
nf2 = nfinit0(polabs, 1, prec);
p1 = subfields(nf2, stoi(cl));
l = lg(p1);
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;
}
if (i == l)
for (i = 1; i < l; i++)
{
pol = gmael(p1, i, 1);
if (degree(gcoeff(nffactor(nf, pol), 1, 1)) == cl) break;
}
if (i == l)
err(bugparier, "makescind (no polynomial found)");
}
pol = polredabs(pol, prec);
return gerepileupto(av, pol);
}
/* compute the Hilbert class field using genus class field theory when
the exponent of the class group is 2 */
static GEN
GenusField(GEN bnf, long prec)
{
long hk, c, l, av = avma;
GEN disc, quat, 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);
div = divisors(disc);
c = 1;
l = 0;
while(l < hk)
{
c++;
d = (GEN)div[c];
if (gcmp1(modii(d, quat)))
{
if (!l)
pol = gsub(x2, d);
else
pol=(GEN)compositum(pol, gsub(x2, d))[1];
l = degree(pol);
}
}
return gerepileupto(av, polredabs(pol, prec));
}
/* if flag = 0 returns the reduced polynomial, flag = 1 returns the
non-reduced polynomial, flag = 2 returns an absolute reduced
polynomial, flag = 3 returns the polynomial of the Stark's unit,
flag = -1 computes a fast and crude approximation of the result */
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;
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;
N = degree((GEN)nf[1]);
cond1 = gmael4(data, 1, 2, 1, 2);
Pi = mppi(newprec);
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));
if (flag >= 0)
{
/* compute S,T differently if nf is quadratic */
if (N == 2)
p1 = QuadGetST(data, newprec);
else
p1 = GetST(dataCR, 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];
if (DEBUGLEVEL) msgtimer("Compute W");
}
else
{
/* 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, N, newprec, 0);
if (n > bnd) n = bnd;
if (DEBUGLEVEL) fprintferr("nmax in QuickPol: %ld \n", n);
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]);
L1 = cgetg(cl+1, t_VEC);
/* we use the formulae L(1) = sum (an / n) */
for (i = 1; i <= cl; 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])) )
p1 = gadd(p1, gdivgs(an, j));
L1[i] = lpileupto(av2, p1);
}
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));
Lp[i] = ldiv(gmul(W, gconj((GEN)L1[i])), p1);
}
}
p1 = ComputeLift(gmael(data, 4, 2));
veczeta = cgetg(h + 1, t_VEC);
for (i = 1; i <= h; i++)
{
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);
}
veczeta[i] = ldivgs(z, 2 * h);
}
if (DEBUGLEVEL >= 2) fprintferr("zetavalues = %Z\n", veczeta);
ro = cgetg(h+1, t_VEC); /* roots */
for (j = 1; j <= h; j++)
{
p1 = gexp(gmul2n((GEN)veczeta[j], 1), newprec);
ro[j] = ladd(p1, ginv(p1));
}
polrelnum = roots_to_pol_intern(realun(newprec),ro, 0,0);
if (DEBUGLEVEL)
{
if (DEBUGLEVEL >= 2) fprintferr("polrelnum = %Z\n", polrelnum);
msgtimer("Compute %s", (flag < 0)? "quickpol": "polrelnum");
}
if (flag < 0)
return gerepileupto(av, gcopy(polrelnum));
/* we try to recognize this polynomial */
polrel = RecCoeff(nf, polrelnum, v, newprec);
if (!polrel) /* if it fails... */
{
long pr;
if (++cpt >= 3) err(talker,
"insufficient precision: computation impossible");
/* we compute the precision that we need */
pr = 1 + (gexpo(polrelnum)>>TWOPOTBITS_IN_LONG);
if (pr < 0) pr = 0;
newprec = DEFAULTPREC + max(newprec,pr);
if (DEBUGLEVEL) err(warnprec, "AllStark", newprec);
nf = nfnewprec(nf, newprec);
data[5] = (long)CharNewPrec((GEN)data[5], nf, newprec);
gptr[0] = &data;
gptr[1] = &nf;
gerepilemany(av, gptr, 2);
goto LABDOUB;
}
/* and we compute the polynomial of eps if flag = 3 */
if (flag == 3)
{
n = fetch_var();
p1 = gsub(polx[0], gadd(polx[n], ginv(polx[n])));
polrel = polresultant0(polrel, p1, 0, 0);
polrel = gmul(polrel, gpowgs(polx[n], h));
polrel = gsubst(polrel, n, polx[0]);
polrel = gmul(polrel, leading_term(polrel));
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));
/* we just want the polynomial computed */
if (flag!=2) return gerepileupto(av, gcopy(polrel));
/* we want a reduced absolute polynomial */
return gerepileupto(av, rnfpolredabs(nf, polrel, 2, newprec));
}
/********************************************************************/
/* Main functions */
/********************************************************************/
/* compute the polynomial over Q of the Hilbert class field of
Q(sqrt(D)) where D is a positive fundamental discriminant */
GEN
quadhilbertreal(GEN D, long prec)
{
long av = avma, cl, newprec;
GEN pol, bnf, bnr, dataC, bnrh, nf, exp;
if (DEBUGLEVEL) timer2();
disable_dbg(0);
/* quick computation of the class number */
cl = itos((GEN)quadclassunit0(D, 0, NULL, prec)[1]);
if (cl == 1)
{
disable_dbg(-1);
avma = av; return polx[0];
}
/* initialize the polynomial defining Q(sqrt{D}) as a polynomial in y */
pol = quadpoly(D);
setvarn(pol, fetch_var());
/* 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); }
/* find the modulus defining N */
bnr = buchrayinitgen(bnf, gun, prec);
dataC = InitQuotient(bnr, gzero);
bnrh = FindModulus(dataC, 1, &newprec, prec);
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(nf, pol, cl, prec));
}
GEN
bnrstark(GEN bnr, GEN subgroup, long flag, long prec)
{
long cl, N, newprec, av = avma;
GEN bnf, dataS, p1, Mcyc, nf, data;
if (flag < 0 || flag > 3) err(flagerr,"bnrstark");
/* check the bnr */
checkbnrgen(bnr);
bnf = (GEN)bnr[1];
nf = (GEN)bnf[7];
Mcyc = diagonal(gmael(bnr, 5, 2));
N = degree((GEN)nf[1]);
if (N == 1)
err(talker, "the ground field must be distinct from Q");
/* check the bnf */
if (!varn(gmael(bnf, 7, 1)))
err(talker, "main variable in bnrstark must not be x");
if (cmpis(gmael3(bnf, 7, 2, 1), N))
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");
}
/* compute bnr(conductor) */
p1 = conductor(bnr, subgroup, 2, prec);
bnr = (GEN)p1[2];
subgroup = (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));
if (cl == 1) return polx[0];
timer2();
/* find a suitable extension N */
dataS = InitQuotient(bnr, subgroup);
data = FindModulus(dataS, 1, &newprec, prec);
if (newprec > prec)
{
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 bnr, 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,
or if 1 then compute the value at s = 1 (and in this case, only for
non-trivial characters), 2: if 0 then compute the value of the
primitive L-function associated to chi, if 1 then compute the value
of the L-function L_S(s, chi) where S is the set of places dividing
the modulus of bnr (and the infinite places), 3: return also the
character */
GEN
bnrL1(GEN bnr, long flag, long prec)
{
GEN bnf, nf, cyc, Mcyc, p1, L1, chi, cchi, allCR, listCR, dataCR;
GEN S, T, rep, indCR, invCR;
long N, cl, i, j, nc, a, av = avma;
bnf = (GEN)bnr[1];
nf = (GEN)bnf[7];
cyc = gmael(bnr, 5, 2);
Mcyc = diagonal(cyc);
N = degree((GEN)nf[1]);
if (N == 1)
err(talker, "the ground field must distinct from Q");
if ((flag < 0) || (flag > 8))
err(flagerr,"bnrL1");
/* check the bnr */
checkbnrgen(bnr);
/* compute bnr(conductor) */
if (!(flag & 2))
{
p1 = conductor(bnr, gzero, 2, prec);
bnr = (GEN)p1[2];
cyc = gmael(bnr, 5, 2);
Mcyc = diagonal(cyc);
}
cl = itos(det(Mcyc));
/* compute all the characters */
allCR = FindEltofGroup(cl, cyc);
/* make a list of all non-trivial characters modulo conjugation */
listCR = cgetg(cl, t_VEC);
indCR = new_chunk(cl);
invCR = new_chunk(cl);
nc = 0;
for (i = 1; i < cl; i++)
{
chi = (GEN)allCR[i];
cchi = ConjChar(chi, cyc);
a = i;
for (j = 1; j <= nc; j++)
if (gegal(gmael(listCR, j, 1), cchi)) a = -j;
if (a > 0)
{
nc++;
listCR[nc] = lgetg(3, t_VEC);
mael(listCR, nc, 1) = (long)chi;
mael(listCR, nc, 2) = (long)bnrconductorofchar(bnr, chi, prec);
indCR[i] = nc;
invCR[nc] = i;
}
else
indCR[i] = -invCR[-a];
}
setlg(listCR, nc + 1);
if (nc == 0) err(talker, "no non-trivial character in bnrL1");
/* compute the data for these characters */
dataCR = InitChar(bnr, listCR, prec);
p1 = GetST(dataCR, prec);
S = (GEN)p1[1];
T = (GEN)p1[2];
if (flag & 1)
L1 = cgetg(cl, t_VEC);
else
L1 = cgetg(cl + 1, t_VEC);
for (i = 1; i < cl; i++)
{
a = indCR[i];
if (a > 0)
L1[i] = (long)GetValue((GEN)dataCR[a], (GEN)S[a], (GEN)T[a], flag & 1,
flag & 2, prec);
}
for (i = 1; i < cl; i++)
{
a = indCR[i];
if (a < 0)
L1[i] = lconj((GEN)L1[-a]);
}
if (!(flag & 1))
L1[cl] = (long)GetValue1(bnr, flag & 2, prec);
else
cl--;
if (flag & 4)
{
rep = cgetg(cl + 1, t_VEC);
for (i = 1; i <= cl; i++)
{
p1 = cgetg(3, t_VEC);
p1[1] = allCR[i];
p1[2] = L1[i];
rep[i] = (long)p1;
}
}
else
rep = L1;
return gerepileupto(av, gcopy(rep));
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to stark.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / modules