allocatemem (2^29);
\r shimura.gp
sufficiently_accurate_bnfinit (pol, flag) =
/* For this file, just a wrapper around bnfinit. */
{
return (bnfinit (pol, flag));
}
init_cmfield (A, B) =
{
return (init_cmfield_basic (A, B));
}
shimura_subgroup (cmfield) =
/* cmfield: as output by init_cmfield
If unitindex==1, returns the image of the class group of Kr under the
type norm map as a vector. A third and fourth component is added as
described in shimura_cosets.
If unitindex==2, only half of the image of the type norm is computed,
the second half has entries [a, n*eps0]. */
{
my (K0, K, Kr, m, no, S:vecsmall, wt, tngen, dlog, Clr, hr, clr,
Tn, adlog, offset);
K0 = cmfield [1];
K = cmfield [2];
Kr = cmfield [7];
m = vector (Kr.no);
/* overdimension, shorten later */
/* We need a set to keep track of already covered ideal classes. Since
the built-in set functions are slow, we use a binary vector. The
variable wt contains weights so that a scalar product maps abstract
group elements to their index in the vector. (An additional shift by 1
is needed to start with the index 1.) */
S = vectorsmall (K.no);
wt = vector (#K.clgp [2]);
wt [1] = 1;
for (i = 2, #wt,
wt [i] = wt [i-1] * K.clgp [2][i-1];
);
/* precompute the type norms of the generators of Cl (Kr) and the
decomposition over Cl (K) of their ideals */
tngen = vector (#Kr.clgp [2], i,
type_norm_reduced (cmfield, Kr.clgp [3][i]));
dlog = Mat (vector (#tngen, i,
bnfisprincipal (K, tngen [i][1], 0)));
/* a small prime not dividing the class number, determines frequency of
type norm reduction */
offset = 11;
while (K.no % offset == 0,
offset = nextprime (offset + 1);
);
/* enumerate Cl (Kr) in Clr, with hr being the current length of the vector.
Code copied from explode_abstractgroup. */
/* prepare first class group entry */
Clr = vector (Kr.no);
Clr [1] = vectorv (#Kr.clgp [2]);
hr = 1;
/* handle principal ideal with generator 1 */
Tn = vector (#Clr);
/* list of class type norms */
Tn [1] = [[1, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0; 0, 0, 0, 1], Mod (1, K0.pol)];
adlog = vectorv (#K.clgp [2]);
no = 1;
/* number of current entries in m */
m [no] = [Tn [1][1], Tn [1][2], adlog, 1];
S [1] = 1;
/* enumerate */
for (i = 1, #Kr.clgp [2],
/* treat generator i and its powers */
for (j = 1, Kr.clgp [2][i] - 1,
/* treat power j of generator i */
for (k = 1, hr,
Clr [j*hr + k] = Clr [k];
Clr [j*hr + k][i] = Clr [(j-1)*hr + k][i] + 1;
clr = Clr [j*hr + k];
/* handle ideal */
Tn [j*hr + k] = vector (2);
Tn [j*hr + k][1] = idealmul (K,
Tn [(j-1)*hr + k][1], tngen [i][1]);
Tn [j*hr + k][2] = Tn [(j-1)*hr + k][2] * tngen [i][2];
adlog = dlog * clr;
for (r = 1, #adlog,
adlog [r] %= K.clgp [2][r];
);
if (!S [wt * adlog + 1],
no++;
/* reduce type norm before storing */
Tn [j*hr + k] = reduce_type_norm (cmfield, Tn [j*hr + k]);
m [no] = [Tn [j*hr+k][1], Tn [j*hr+k][2], adlog, 1];
S [wt * adlog + 1] = 1;
,
/* reduce type norm from time to time */
if ((j*hr + k) % offset == 0,
Tn [j*hr + k] = reduce_type_norm (cmfield, Tn [j*hr + k]);
);
);
);
);
hr = hr * Kr.clgp [2][i];
);
return (vector (no, i, m [i]));
};
next_shimura_element (cmfield, S:vecsmall, wt) =
/* cmfield: as output by init_cmfield
S a binary vector whose entries indicate which ideals in the class group
have already be covered
wt the weights used to map abstract class group elements to indices of S
returns an element of the Shimura group with ideal not contained in S;
a third and fourth component are added as described in shimura_cosets.
*/
{
my (K0, K, Krel, K0narrow, idealwt, Cl, a, prin, alpha, tn);
K0 = cmfield [1];
K = cmfield [2];
Krel = cmfield [3];
K0narrow = bnrinit (K0, [1, [1, 1]], 1);
idealwt = vector (#K.clgp [2], i, log (idealnorm (K, K.clgp [3][i])));
Cl = vecsort (explode_abstractgroup (K.clgp [2]),
(v1, v2) -> cmp_weight (v1, v2, idealwt));
/* sort according to norm of unreduced representative; also helps
to avoid being stuck in one cyclic component */
for (i = 1, #Cl,
/* check if ideal already covered */
if (!S [wt * Cl [i] + 1],
/* check the Shimura condition */
a = idealfactorback (K, K.clgp [3], Cl [i], 1);
prin = bnrisprincipal (K0narrow, rnfidealnormrel (Krel, K.zk*a));
if (!prin [1],
alpha = Mod (K0.zk * prin [2], K0.pol);
if (!is_totally_positive (K0, alpha),
error ("*** Error in random_shimura_element: non totally positive generator");
);
tn = reduce_type_norm (cmfield, [a, alpha]);
/* superfluous in examples, probably since a is already small
due to the choice of a small element of Cl */
return ([tn [1], tn [2], Cl [i], 1]);
);
);
);
};
shimura_cosets (cmfield) =
/* cmfield: as output by init_cmfield
returns a vector of cosets of the image of Cl (Kr) under the type norm map
each entry is a vector of four components:
The first two components are the type norm [a, n] and represent the
Shimura group element; the other two components are technical and
used in the function triple.
The third component is the decomposition of a over the class group of K.
The fourth component is meaningful only when unitindex=2; it is 1 if
n is the "original" value computed by the type norm map (currently
an integer, but this may change in the future); it is 0 if n has
been obtained by a multiplication with eps0.
*/
{
my (res, K, Kr, unitindex, eps0, noC, m, S, wt, tn, cosrep, coset);
gettime ();
K = cmfield [2];
eps0 = cmfield [4];
unitindex = cmfield [5];
Kr = cmfield [7];
noC = cmfield [10];
gettime ();
S = vectorsmall (K.no);
wt = vector (#K.clgp [2]);
wt [1] = 1;
for (i = 2, #wt,
wt [i] = wt [i-1] * K.clgp [2][i-1];
);
/* compute the image m of Cl (Kr) */
m = shimura_subgroup (cmfield);
for (i = 1, #m,
S [wt * m [i][3] + 1] = 1;
);
res = vector (noC / (#m * unitindex));
res [1] = m;
printf ("Time for Shimura subgroup: %.1f\n", gettime () / 1000.0);
print ("Number of further cosets: ", noC / (#m * unitindex) - 1);
/* compute other cosets */
for (i = 2, noC / (#m * unitindex),
/* get new coset representative */
cosrep = next_shimura_element (cmfield, S, wt);
printf ("Time for coset representative %i: %.1f\n", i-1, gettime () / 1000.0);
/* multiply with m */
coset = vector (#m);
for (j = 1, #m,
tn = reduce_type_norm (cmfield,
[idealmul (K, cosrep [1], m [j][1]), cosrep [2] * m [j][2]]);
coset [j] = [tn [1], tn [2], cosrep [3] + m [j][3], 1];
for (k = 1, #K.clgp [2],
coset [j][3][k] %= K.clgp [2][k];
);
S [wt * coset [j][3] + 1] = 1;
);
res [i] = coset;
printf ("Time for coset %i: %.1f\n", i-1, gettime () / 1000.0);
);
if (unitindex == 2,
/* add at the end all the elements corresponding to multiplication
of the generator by eps0 */
for (i = 1, #res,
res [i] = concat (res [i],
vector (#m, j, [res [i][j][1],
(eps0 * res [i][j][2]).pol,
res [i][j][3],
0]));
);
);
return (res);
};
triple (cmfield, bp, sh, dlogconj) =
/* cmfield: as output by init_cmfield
bp: a base point for the action of the Shimura group as found by
find_basepoint
sh: an element of the Shimura group
dlogconj: a matrix containing (as columns) the precomputed action of
complex conjugation on the generators of the classgroup of K
returns a triple [a, xi, flag] = sh*bp
Watch out for an error in BrGrLa11: The action of sh on bp is by
multiplying the first components, but _dividing_ the second ones.
flag can take three values:
1: A real root of the class polynomial is represented by this triple.
Happens when (a, xi) and (abar, xi) are equivalent in the Shimura
sense, that is, (a/abar, 1) is equivalent to (OK, 1) in the
Shimura group.
If unitindex==1, this means that a and abar are equivalent in the
class group.
If unitindex==2, this means that furthermore, the unit passing from
one to the other is a power of eps0^2.
2: A pair of complex-conjugate roots of the class polynomial is
represented by this triple.
If unitindex==1, a and abar are compared in the abstract class group,
and a is retained only if it is smaller.
If unitindex==2, this criterion still holds. If a and abar are
equivalent in the class group and the unit passing from one to the
other is an odd power of eps0, then we decide according to sh[2]:
The function is called twice with the same sh[1] (leading to the
same a); we retain only one of them according to sh [4].
0: As 2, but this is "the other" complex-conjugate root and should be
omitted later.
*/
{
my (K0, K, Krel, unitindex, a, xi, abar, prin, prinbar, comp,
flag, u, uubar);
K0 = cmfield [1];
K = cmfield [2];
Krel = cmfield [3];
unitindex = cmfield [5];
a = idealmul (K, bp [1], sh [1]);
xi = (bp [2] / rnfeltup (Krel, sh [2])).pol; /* .pol to remove the modulus */
prin = bp [3] + sh [3];
for (i = 1, #K.clgp [2],
prin [i] %= K.clgp [2][i];
);
prinbar = dlogconj * prin;
for (i=1, #prinbar,
prinbar [i] %= K.clgp [2][i];
);
comp = cmp_compat (prin, prinbar);
if (comp < 0,
flag = 2,
/* else */ if (comp > 0,
flag = 0,
/* else comp == 0 */
if (unitindex == 1,
flag = 1,
/* else unitindex == 2 */
abar = complex_conjugate_ideal (K, a);
u = K.zk * bnfisprincipal (K, idealdiv (K, a, abar)) [2];
uubar = rnfeltdown (Krel, (u * complex_conjugate_element (u)) % K.pol);
if (bnfisunit (K0, uubar)[1] % 2 == 0,
flag = 1,
/* else */ if (sh [4] == 1,
flag = 2,
/* else */
flag = 0;
););
);
););
return ([a, xi, flag]);
};
symplectic_bases (cmfield) =
/* cmfield a quartic CM field
returns symplectic bases that yield period matrices for factors of the
class polynomial
Precisely, the return value is a vector with as many elements as there
are cosets of the Shimura subgroup (or factors of the class polynomial);
each entry is a list with as many entries as there are period
matrices; each of these entries is a two-element vector, the first
element being the symplectic basis, the second one being 1 or 2 for a
real or a pair of complex conjugate roots of the class polynomial */
{
my (A, B, Arred, Brred, K0r, K0, K, unitindex, Phi,
D, cosets, dlogconj, bases, m, noC, bp, trip);
K0 = cmfield [1];
A = polcoeff (K0.pol, 1);
B = polcoeff (K0.pol, 0);
Arred = cmfield [14];
Brred = cmfield [15];
K0r = cmfield [6];
K = cmfield [2];
unitindex = cmfield [5];
Phi = cmfield [9];
noC = cmfield [10];
D = K0.disc;
cosets = shimura_cosets (cmfield);
gettime ();
bp = find_basepoint (cmfield, Phi);
printf ("Time for base point: %.1f\n", gettime () / 1000.0);
dlogconj = Mat (vector (#K.clgp [2], i,
bnfisprincipal (K, complex_conjugate_ideal (K, K.clgp [3][i]), 0)));
gettime ();
bases = vector (#cosets);
for (l = 1, #cosets,
m = cosets [l];
bases [l] = List ();
for (k=1, #m,
trip = triple (cmfield, bp, m [k], dlogconj);
if (trip [3] != 0,
trip [1] = switch_ideal_basis_for_symplectic_input (K, trip [1]);
listput (bases [l],
[symplectic_basis (K, trip [1], trip [2]), trip [3]]);
);
);
);
printf ("Time for symplectic bases: %.1f\n", gettime () / 1000.0);
return (bases);
};
/*
A = 35; B = 65;
A = 352; B = 6257;
\p 100
cm = init_cmfield (A, B);
print_period_matrices (cm, 1);
*/