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