allocatemem(2^28);
install (extendedgcd, G);
install (Q_denom, G);
'x; 'y; 'z;
/* order of variables:
z for the real subfields
y for the quartic CM fields
x for their Galois closure */
/* {{{ Functions for dealing with the Shimura group structure */
/* We manipulate two groups:
* CK: The Shimura group.
* S: The subgroup of CK which is generated by the image of the class
* group of the reflex field under the type norm map.
*
* CK enjoys the exact sequence already quoted above:
*
* 1 \to U0+/N_{K/K0}(U) \to CK \to Cl(OK) \to Cl+(OK0) \to 1
*
* The group U0+/N_{K/K0}(U) has already been documented above.
*/
/* {{{ In order to obtain the structure of CK, the procedure goes as follows:
*
* 0: Needed basic equipment on CK:
* - basic product operation.
* - identity element.
* - equality test.
* - reduction operation.
* - factorback routine.
* 1: Compute the kernel of the map Cl(OK) \to Cl+(OK0):
* - apply bnfisprincipal in the narrow class group Cl+(OK0) on the
* images of the generators of Cl(OK).
* - knowing the Smith invariants of Cl+(OK0), deduce combinations of
* these generators, which yield zero image (hnf computation). These
* combinations thus generate the kernel.
* 2: Get relations among the ideals generating this kernel.
* - use the Smith invariants of Cl(OK) to deduce relations among the
* above generators. (hnf again).
* 3: Obtain generators for CK:
* - Lift arbitrarily the above generators to pairs in CK (pick
* arbitrary totally positive generators). If the group
* U0+/N_{K/K0}(U) is not trivial, add an extra generator (OK,eps0).
* Together, the corresponding pairs are called the ``primary''
* generators of CK, and denoted Cgens.
* - The relation matrix computed in (2) is not complete, because not
* everything which maps to 0 in Cl(OK) is necessarily 0 in CK. Thus
* compute all combinations, and possibly add an extra fixing
* (OK,eps0) component. We thus form an expanded relation matrix.
* 4: Finish:
* - Use the above relation matrix to compute the Smith invariants of
* CK, as well as its generators.
* }}} */
/* {{{ An important by-product of this computation is a generalized dlog
* procedure. Obtaining the dlog of a given element a of CK goes as
* follows:
*
* 0: map to Cl(OK)
* - strip off the second component.
* - apply bnfisprincipal.
* 1: convert to combination of the generators of the kernel.
* - a simple linear system. We expect the existence of an integer
* solution.
* 2: fix to get equality in CK.
* - use the previous exponents to form a combination of the elected
* (primary) generators of CK.
* - add a possible (OK,eps0) component if equality is not reached.
* 3: convert to Smith generators.
* - This uses one of the change of basis matrices from the Smith form
* computation.
* }}} */
/* {{{ The computation of S is not considerably simpler, although it proceeds
* through classical group theory operations.
* 1: Compute generators:
* - The type norm implementation is completely sufficient.
* 2: Get relations.
* - Compute the dlogs of the images. Form a matrix.
* - Add a diagonal block with the Smith invariants.
* - Do an hnf computation to get relations.
* 3: Get invariants.
* - Compute the Smith form of that matrix to get the Smith invariants
* of S.
* - Keep track of the change of basis matrix for later use.
* 4: Get coset representatives.
* - The hnf computation from step 2 actually produces generators for
* the quotient group. These are given by the columns of the change
* of basis matrix, on the part where the dlog of the result is not
* zero.
* }}} */
/***** complex conjugation and norm *****/
complex_conjugate_element (alpha) =
/* alpha an element of K or Kr in terms of y, or a vector of such elements
returns the complex conjugate of alpha, also in terms of y */
{
if (type(alpha) == "t_POLMOD",
return (Mod (subst (alpha.pol, y, -y), alpha.mod));
);
return (subst (alpha, y, -y));
}
complex_norm (alpha) =
/* alpha an element of K or Kr as a t_POLMOD in terms of y
returns the complex norm of alpha as a polynomial in z, the generator
of K or Kr */
{
my (n);
my (complex_conjugate_element = complex_conjugate_element);
n = lift (alpha * complex_conjugate_element (alpha));
return (z * polcoeff (n, 2) + polcoeff (n, 0));
}
complex_conjugate_ideal (K, af) =
/* K a CM field
af an ideal of K
returns the image (in HNF) of af under the complex conjugation
automorphism y \mapsto -y. */
{
/* switch to polynomial basis, conjugate, switch back to integral basis,
hnf */
return (idealhnf (K, Mat (matalgtobasis(K,
complex_conjugate_element (K.zk*af)))));
};
/***** output and file handling *****/
delete_file (file) =
/* file the name of a file to be deleted, or 0 if nothing is to be done */
{
if (file != 0,
system (Strprintf ("rm -f %s", file));
);
};
output_element (file, v) =
/* v a scalar to be printed
file the name as a string of the file to print to, or 0 to print to
stdout */
{
if (file == 0,
print1 (v),
/* else */
write1 (file, v);
);
};
output (file, v) =
/* v a vector of elements to be printed
file the name as a string of the file to print to, or 0 to print to
stdout */
{
for (i = 1, #v,
output_element (file, v [i]);
);
};
print_group (G) =
/* G an abelian group given by its cycle structure
pretty prints the structure of G
*/
{
if (G == 0,
print ("C1"),
/* else */
print1 ("C", G [#G]);
for (i = -(#G-1), -1, print1 (" x C", G [-i]));
print ();
);
};
print_period_matrices (cmfield, tofile) =
/* cmfield a quartic CM field
If tofile==false, prints on screen the (reduced) period matrices in the
format defined in README.format.
Otherwise, prints them into the file named D_A_B.in, where D is the
discriminant of the real-quadratic subfield. */
{
my (K0, A, B, Arred, Brred, K0r, D, file, Omegavec, m, Omega, den, v);
K0 = cmfield [1];
A = polcoeff (K0.pol, 1);
B = polcoeff (K0.pol, 0);
Arred = cmfield [14];
Brred = cmfield [15];
K0r = cmfield [6];
D = K0.disc;
if (tofile,
file = concat ([D, "_", A, "_", B, ".in"]);
file = concat (Strexpand("$DATA_PREFIX"), file);
delete_file (file),
/* else */
file = 0;
);
Omegavec = period_matrices (cmfield);
/* starting output */
output (file, [-3, "\n"]); /* version string */
output (file, [D, " ", A, " ", B, "\n"]);
output (file, [K0r.disc, " ", Arred, " ", Brred, "\n"]);
output (file, [#Omegavec, "\n"]);
for (l = 1, #Omegavec,
m = Omegavec [l];
output (file, ["\n", sum (i = 1, #m, m [i][2]), "\n"]);
for (k=1, #m,
output (file, [m [k][2], "\n"]);
Omega = m [k][1];
den = Q_denom (Omega);
output (file, [den, "\n"]);
v = [Vecrev_compat (Omega [1,1], 4), Vecrev_compat (Omega [1,2], 4),
Vecrev_compat (Omega [2,2], 4)];
for (i = 1, 3,
for (j = 1, 3,
output (file, [den * v [i][j], " "]);
);
output (file, [den * v [i][4], "\n"]);
);
);
);
}
/***** compatibility functions for pari < 2.6 *****/
/* xxx_compat may be replaced by xxx after switching to pari >= 2.6. */
cmp_compat (f, g) =
/* f, g two objects to be compared (vectors, matrices, scalars, polynomials)
returns a comparison value */
{
my (tmp);
if ((type (f) == "t_INT" || type (f) == "t_FRAC") &&
(type (g) == "t_INT" || type (g) == "t_FRAC"),
return (lex (f, g));
);
if (type (f) != type (g),
error ("*** Different types in cmp_compat: ", type(f), " ", type(g));
);
if (type (f) == "t_POL",
return (cmp_compat (Vec (f), Vec (g)));
);
if (type (f) == "t_VEC" || type (f) == "t_COL",
if (#f < #g,
return (-1);
);
if (#f > #g,
return (1);
);
for (i = 1, #f,
tmp = cmp_compat (f [i], g [i]);
if (tmp != 0,
return (tmp);
);
);
return (0);
);
if (type (f) == "t_MAT",
if (#f < #g,
return (-1);
);
if (#f > #g,
return (1);
);
for (i = 1, #f,
tmp = cmp_compat (f [,i], g [,i]);
if (tmp != 0,
return (tmp);
);
);
return (0);
);
error (["*** Type ", type (f), " not implemented in cmp_compat"]);
}
Vecrev_compat (p, n) =
/* In this function, p is supposed to be a polynomial; a vector with n
coefficients is returned even if deg (p) < n-1.
*/
{
my (v);
v = Vecrev (p);
if (#v < n, v = vector (n, i, if (i <= #v, v [i], 0)));
return (v);
}
matinv (M) =
/* Function working around slow inversion of 2x2 matrices in pari */
{
return ((1/matdet (M)) * [M [2,2], -M [1,2]; -M [2,1], M [1,1]]);
}
/***** basic CM field handling *****/
canonical_cm_parameters (A, B) =
/* A, B positive integers such that A^2 - 4*B > 0 and
z^4 + A*z^2 + B is irreducible
returns A' > 0, B' > 0 such that z^4 + A'*z^2 + B' defines the same
field, and A' minimal.
*/
{
my (K, b, M, m, param);
K = bnfinit (y^4+A*y^2+B);
b = K.zk;
/* Write a basis for the Q-vector space of purely imaginary algebraic
numbers, that is, numbers of the form alpha-\bar alpha, expressed
by a matrix with respect to the integral basis of K. */
M = Mat (matalgtobasis (K,
vector (4, i, b [i] - complex_conjugate_element (b [i]))));
/* Intersect with OK and write as algebraic integers again. */
b = matbasistoalg (K, Vec (matrixqz (M, -2)));
/* Minimise the trace form on the Z-module. */
m= qfminim (matrix (#b, #b, i, j, -nfelttrace (K, b[i]*b[j])/2));
/* If there are several minimal integer directions, choose the optimal
one. Example: 95 532 versus 95 2128. */
param = vecsort (vector (#m [3], i,
[-1/2 * nfelttrace (K, (b * m[3][,i])^2), nfeltnorm (K, b * m[3][,i])]),
cmp_compat);
return (param [1]);
}
init_cmfield_basic (A, B) =
/* A, B positive integers such that A^2 - 4*B > 0 and
z^2 + A*z + B is irreducible
returns a CM field structure as a vector containing
1 K0: the real subfield of K
2 K: the CM field
3 Krel: K as a relative extension of K0
4 eps0: a fundamental unit of K0 that is not totally negative
5 unitindex: the index of N_{K/K0}(U) in U0+
6 K0r: the real subfield of the reflex field
7 Kr: the reflex field
8 Krrel: Kr as a relative extension of K0r
9 Phi: an element of K whose square is a totally negative number
in K0. A CM type is implicitly defined from Phi, by choosing the
embeddings which send this element to complex numbers on the
positive imaginary axis (there are only one or two possible
values for this element, modulo equivalence of CM types,
and multiplication by rational constants).
10 noC: size of the Shimura group
11 Lrel: the Galois closure L of K and Kr, as a relative extension of K
12 Lrrel: L as a relative extension of Kr
If K is already Galois, Lrel and Lrrel are set to 0.
13 sigm: If K is galois, a generator of the Galois group of Kr.
0 otherwise.
14 Arred,
15 Brred: reduced, canonical parameters for the reflex field
*/
{
my (K0, Krel, K, eps0, eps, unitindex, Ar, Br, K0r, Krrel, Kr,
Phi, noC, Lrel, Lrrel, sigm, Arred, Brred, tmp);
if (A^2 - 4*B <= 0,
error ("*** Error in main: Equation does not define a quartic CM field");
);
tmp = canonical_cm_parameters (A, B);
if ((A != tmp [1] || B != tmp [2]),
warning ([A, B], " not minimal, use ", tmp, " instead!");
);
/* create CM field and subfields */
K0 = sufficiently_accurate_bnfinit (z^2 + A*z + B, 1);
Krel = rnfinit (K0, y^2-z);
K = bnfinit (Krel.pol, 1); /* K = bnfinit (y^4 + A*y^2 + B, 1); */
/* units */
eps0 = K0.fu [1];
eps = K.fu [1];
if (norm (eps0) == -1,
unitindex = 1,
/* else */
if (bnfsignunit (K0)[1,1] == -1, eps0 = -eps0);
/* forces eps0 to be totally positive */
tmp = rnfeltup (Krel, eps0);
if (tmp == eps || tmp == -eps,
unitindex = 2,
/* else */
unitindex = 1;
);
);
/* size of Shimura group */
noC = K.no / K0.no;
if (norm (eps0) == 1,
noC = noC * 2 / unitindex;
);
/* cm types, reflex fields and Galois closures */
Phi = y;
/* The other possible CM type is defined by
-y*(y^2+(-Tr_{K0/Q} (y^2))/2)). */
if (galoisinit (K.pol),
K0r = K0;
Krrel = Krel;
Kr = K;
Lrel = 0;
Lrrel = 0;
sigm = nfgaloisconj (Kr)[3];
if (sigm == y || sigm == -y,
error ("*** Error in init_cmfield: sigma of order 1 or 2");
),
/* else */
/* Do not reduce the parameters for the reflex field, since it
complicates the Galois closure computations. */
Ar = 2*A;
Br = A^2-4*B;
tmp = canonical_cm_parameters (Ar, Br);
Arred = tmp [1];
Brred = tmp [2];
K0r = bnfinit (z^2 + Ar*z + Br);
Krrel = rnfinit (K0r, y^2-z);
Kr = bnfinit (Krrel.pol);
/* Force the definition of the closure which is consistent with
* our notational choice */
Lrel = rnfinit (K, x^2-Mod(2*y,K.pol)*x+Mod(2*y^2+A,K.pol));
Lrrel = rnfinit (Kr, x^2-Mod(y,Kr.pol)*x+Mod((y^2+A)/2,Kr.pol));
/* forces Lrel.pol == Lrrel.pol, so we can move consistently up and down */
if (Lrel.pol != Lrrel.pol,
error ("*** Error in init_cmfield: inconsistent Galois closure");
);
sigm = 0;
);
return ([K0, K, Krel, eps0, unitindex, K0r, Kr, Krrel, Phi, noC,
Lrel, Lrrel, sigm, Arred, Brred]);
};
print_cmfield (cm) =
/* cm a cm field
prints some information on class groups and so on */
{
my (K0, K, eps0, unitindex, K0r, Kr, noC, Lrel, Arred, Brred,
CKdata, Sdata);
K0 = cm [1];
K = cm [2];
eps0 = cm [4];
unitindex = cm [5];
K0r = cm [6];
Kr = cm [7];
noC = cm [10];
Lrel = cm [11];
Arred = cm [14];
Brred = cm [15];
if (#cm > 15,
CKdata = cm [16];
Sdata = cm [17];
);
print ("G = ", if (Lrel == 0, "C4", "D4"));
print ("DAB = ", K0.disc, " ", polcoeff (K0.pol, 1), " ", polcoeff (K0.pol, 0));
print ("DABr = ", K0r.disc, " ", Arred, " ", Brred);
print ("h0 = ", K0.no);
print ("h1 = ", noC);
print1 ("Cl0 = ");
print_group (K0.clgp [2]);
print1 ("Cl = ");
print_group (K.clgp [2]);
print ("h0r = ", K0r.no);
print ("h1r = ", Kr.no / K0r.no);
print1 ("Cl0r = ");
print_group (K0r.clgp [2]);
print1 ("Clr = ");
print_group (Kr.clgp [2]);
if (#cm > 15,
print1 ("CK = ");
print_group (CKdata [1]);
print1 ("S = ");
print_group (Sdata [1]);
);
print ("eps0norm = ", nfeltnorm (K0, eps0));
print ("unitindex = ", unitindex);
};
explode_cmfield (Aloc, Bloc) =
/* calls init_cmfield and creates global variables from its result;
useful for testing purposes */
{
cm = init_cmfield (Aloc, Bloc);
K0 = cm [1];
K = cm [2];
Krel = cm [3];
eps0 = cm [4];
unitindex = cm [5];
K0r = cm [6];
Kr = cm [7];
Krrel = cm [8];
Phi = cm [9];
noC = cm [10];
Lrel = cm [11];
Lrrel = cm [12];
sigm = cm [13];
Arred = cm [14];
Brred = cm [15];
if (#cm > 15,
CKdata = cm [16];
Sdata = cm [17];
);
A = Aloc;
B = Bloc;
Ar = 2*A;
Br = A^2 - 4*B;
};
/***** Symplectic forms *****/
bezout_vector (v) =
/* v a column vector of integers
returns gcd (v) and a column vector u such that u~ * v = gcd (v) */
{
my (tmp);
tmp = extendedgcd (v);
return ([tmp [1], tmp [2][,#v]]);
};
switch_ideal_basis_for_symplectic_input (K, af) =
/* K a quartic CM field
af an ideal of K (usually in HNF) with respect to the integral basis
returns af with respect to the polynomial basis {1, y, y^2, y^3},
such that it is of HNF in this basis
and the first two elements form a totally isotropic subspace
with respect to the symplectic form
*/
{
my (P, bf, den);
P = [1, 0, 0, 0; 0, 0, 1, 0; 0, 1, 0, 0; 0, 0, 0, 1];
bf = P * K.nf [8]^-1 * af;
/* af in the basis {1, y^2, y, y^3} */
den = denominator (bf);
bf = mathnf (den * bf);
/* Since the defining polynomial of K contains only even powers of y,
multiplication by y^2 of an element in K with only odd resp. even
powers of y preserves this property. The same holds for complex
conjugation, which replaces y by -y.
Now b[,1] is a constant and b[,2] a linear combination of 1 and y^2.
Since xi defining the symplectic form E contains only odd powers
of y and Tr (y) = Tr (y^3) = 0, we have
E (b[,1], b[,1]) = E (b[,1], b[,2]) = E (b[,2], b[,2]) = 0.
Moreover, small coefficients with respect to the polynomial basis
will hopefully yield smaller non-zero entries. */
bf = P * bf / den;
/* switch back to the basis {1, y, y^2, y^3} */
return (bf);
}
symplectic_transformation (A) =
/* A the matrix of a symplectic form
returns the transformation matrix towards a symplectic basis
follows Alg. 4.2 on p. 50 of Streng10 */
{
my (n, E, V, EV, k, ep, e, vp, v);
n = #A;
E = matrix (n, 0);
V = matrix (n, 0);
/* matrices of column vectors e_1,...,e_{i-1} and v_1,...,v_{i-1} */
EV = matrix (n, 0);
/* E and V concatenated in an interleaved fashion: e_1, v_1, e_2, ... */
for (i = 1, n/2,
/* look for vector eip, independent of EV, among the unit vectors */
k = 0;
until (matrank (concat (EV, ep)) == 2*i-1,
k = k + 1;
ep = vectorv (n, l, l==k);
);
/* normalise ep to e */
e = ep - sum (j = 1, i-1, E [, j]~ * A * ep * V [, j])
+ sum (j = 1, i-1, V [, j]~ * A * ep * E [, j]);
e /= content (e);
/* choose vp such that e~ * A * vp == 1 */
vp = bezout_vector (e~ * A);
if (vp [1] != 1,
error ("*** Error in symplectic_transformation",
"*** gcd not 1 for vp=", vp);
);
vp = vp [2];
/* normalise vp to v */
v = vp - sum (j = 1, i-1, E [, j]~ * A * vp * V [, j])
+ sum (j = 1, i-1, V [, j]~ * A * vp * E [, j]);
/* add vectors to matrices */
E = concat (E, e);
V = concat (V, v);
EV = concat ([EV, e, v]);
);
return (concat (E, V));
};
symplectic_form (K, af, xi) =
/* K a quartic CM field
af an ideal over the polynomial basis {1, y, y^2, y^3}
xi such that (Phi, af, xi) forms an admissible triple
for some CM type Phi
returns the matrix of the symplectic form associated with af and xi */
{
my (res, basis);
basis = [1, y, y^2, y^3] * af;
res = matrix (4, 4, i, j,
nfelttrace (K,
(xi * complex_conjugate_element (basis [i]) * basis [j]) % K.pol));
return (res);
};
symplectic_basis (K, af, xi) =
/* K a quartic CM field
af an ideal over the polynomial basis {1, y, y^2, y^3}
xi such that (Phi, af, xi) forms an admissible triple
for some CM type Phi
returns a symplectic basis of af with respect to the symplectic form
associated with af and xi (in terms of the polynomial basis of K) */
{
return (af * symplectic_transformation (symplectic_form (K, af, xi)));
};
/***** Type norms *****/
reduce_type_norm (cm, tn) =
/* cm as output by init_cmfield
tn a tuple [a, n] as output by type_norm
returns an equivalent, but smaller representative of the given class
in the Shimura group
*/
{
my (K0, K, tnred, alpha, u, pow);
K0 = cm [1];
K = cm [2];
/* First reduce the ideal and keep track of the second component. */
tnred = idealred (K, [tn [1], 1]);
alpha = nfbasistoalg (K, tnred [2]);
alpha = (alpha * complex_conjugate_element (alpha)).pol;
alpha = if (poldegree (alpha) > 0, polcoeff (alpha, 2) * z, 0)
+ polcoeff (alpha, 0);
tnred [2] = Mod (tn [2] / alpha, K0.pol);
/* Then reduce the second component modulo N_{K/K0}(U). The latter
contains ; both are equal when Norm(eps0)=+1, and
conjecturally, when Norm(eps0)=-1 and K/Q is not Galois of type C2xC2,
see the discussion at the top of this file. */
/* determine the optimal unit for reduction */
u = tnred [2]
/ nfbasistoalg (K0,
bnfisprincipal (K0, idealhnf (K0, tnred [2])) [2]);
/* In the git version of gp, the call to idealhnf is not
needed any more. */
/* determine the power of K0.fu [1] and round to a multiple of 2 */
pow = bnfisunit (K0, u) [1];
pow = sign (pow) * (abs (pow) - pow % 2);
if (pow != 0,
tnred [2] /= K0.fu [1]^pow;
);
return (tnred);
}
type_norm (cmfield, ar) =
/* cmfield: as output by init_cmfield
ar: an ideal of Kr
returns m(ar) (in the notation of BrGrLa11), that is, a pair [a, n]
with a the type norm of ar and n = N_{Kr/Q} (ar) */
{
my (Kr, Lrel, Lrrel, sigm, a, n);
Kr = cmfield [7];
Lrel = cmfield [11];
Lrrel = cmfield [12];
sigm = cmfield [13];
if (Lrel == 0, /* Galois case */
a = idealmul (Kr, ar, nfgaloisapply (Kr, sigm, ar)),
/* else non Galois case */
a = rnfidealnormrel (Lrel, rnfidealup (Lrrel, ar));
);
n = idealnorm (Kr, ar);
return ([a, n]);
};
type_norm_reduced (cm, ar) =
/* input as for type_norm
output is additionally reduced by integers
*/
{
return (reduce_type_norm (cm, type_norm (cm, ar)));
}
/***** symbolic treatment of embedding of K0r */
to_standard_basis_K0r (f, Ar) =
/* f an element of K0r given as a linear polynomial in z
Ar the linear coefficient of K0r.pol
returns [a, b] such that f = a + b*sqrt (B),
where B = (Ar^2 - 4*B) / 16 = K0r.pol.disc / 16,
and z is interpreted as (-Ar - sqrt (K0r.pol.disc)/2.
Then when the complex embedding of y is that of maximal, positive
imaginary part, the embedding of a + b*sqrt (B) is that with the
positive root of B.
Notice also that if f has integral coefficients, then the return vector
is also integral.
*/
{
my (f0, f1);
f0 = polcoeff (f, 0);
f1 = polcoeff (f, 1);
return ([f0 - (Ar/2) * f1, -2*f1]);
}
cmp_K0r (f, g, Ar, B) =
/* f, g two elements of K0r given as linear polynomials in z
Ar the linear coefficient of K0r.pol in front of z
B the constant coefficient of K0.pol
returns -1, 0, 1 depending on whether f <, =, > g after the embedding of
K0r into R that sends z to the smaller root of K0r.pol
*/
{
my (h);
if (f == g,
return (0);
);
h = to_standard_basis_K0r (f - g, Ar);
return (cmp_zero_K0r (h, B));
}
cmp_zero_K0r (h, B) =
/* h an element of K0r given as a two-element vector such that
h = h [1] + h [2] * sqrt (B)
B the constant coefficient of K0.pol
returns -1, 0, 1 depending on whether h <, =, > 0 after the embedding of
K0r into R that sends B to its positive root
*/
{
if (h [1] <= 0 && h [2] <= 0,
if (h [1] == 0 && h [2] == 0,
return (0)
,
return (-1);
);
);
if (h [1] >= 0 && h [2] >= 0,
return (1);
);
/* one of h1 and h2 is positive, the other one negative */
if (h [1]^2 > B * h [2]^2,
return (sign (h [1]));
,
return (sign (h [2]));
);
}
floor_K0r (f, Ar, B) =
/* f an element of K0r given as a linear polynomial in z
Ar the linear coefficient of K0r.pol in front of z
B the constant coefficient of K0.pol
returns the rational integer, rounded down, obtained when replacing z in f
by the smaller root of K0r.pol
*/
{
my (h, a, b, res);
h = to_standard_basis_K0r (f, Ar);
/* now f = h[1] + h[2] * sqrt (B); we use the positive root */
/* We use floor(a) + floor(b) <= floor(a+b) <= floor(a) + floor(b) + 1. */
a = h [1] \ 1;
b = sqrtint ((h [2]^2 * B) \ 1);
if (h [2] < 0,
b = -b-1;
);
res = a + b;
if (cmp_zero_K0r ([res + 1 - h [1], -h[2]], B) <= 0,
return (res + 1);
,
return (res);
);
}
/***** Siegel space *****/
/* Matrices in Sp_4 (\Z) are given as 2x2 matrices with entries
that are 2x2 matrices */
minkowski_reduction_K0r (M) =
/* M a real, symmetric, positive definite 2x2 matrix with entries in K0r,
given as t_POLMOD in z
returns U such that U*M*transpose (U) = [a, b; b, c] satisfies
0 <= 2*b <= a <= c */
{
my (K0rpol, Ar, B, S, U, R, ok, q);
K0rpol = M [1, 1].mod;
Ar = polcoeff (K0rpol, 1);
B = poldisc (K0rpol) / 16;
S = [0, 1; 1, 0];
/* Exchanges a and c. */
U = matid (2);
ok = 0;
until (ok,
ok = 1;
/* Obtain -a < 2*b <= a. */
q = floor_K0r (lift (M [1,2] / M [1, 1]) + 1/2, Ar, B);
R = [1, 0; -q, 1];
M = R * M * mattranspose (R);
U = R * U;
/* Obtain a <= c. */
if (cmp_K0r (lift (M [1, 1] - M [2,2]), 0, Ar, B) > 0,
ok = 0;
M = S * M * S;
U = S * U;
);
);
/* Force b >= 0. */
if (cmp_K0r (lift (M [1, 2]), 0, Ar, B) < 0,
R = [1, 0; 0, -1];
M = R * M * R;
U = R * U;
);
return (U);
}
real_reduction_K0r (M, Ar, B) =
/* M a symmetric 2x2 matrix with entries in K0r, given as linear polynomials
in z
Ar the linear coefficient of K0r.pol
B the constant coefficient of K0.pol
returns T such that M + T = [a, b; b, c] satisfies
-1/2 < a, b, c <= 1/2 after the embedding that sends z to the smaller
real root of K0r.pol
*/
{
my (t1, t2, t3);
t1 = -floor_K0r (M [1, 1] + 1/2, Ar, B);
t2 = -floor_K0r (M [1, 2] + 1/2, Ar, B);
t3 = -floor_K0r (M [2, 2] + 1/2, Ar, B);
return ([t1, t2; t2, t3]);
}
siegel_transform_apply (T, M) =
/* T = [A, B; C, D] a matrix in Sp_4 (\Z)
M a complex, symmetric 2x2 matrix
returns T \circ M = (A M + B)*(C M + D)^(-1) */
{
return ((T [1,1] * M + T [1,2]) * (T [2,1] * M + T [2,2])^(-1));
}
gottschling_reduction_Kr (M, Ar, B) =
/* M a period matrix with entries in Kr, given as t_POLMOD
Ar the linear coefficient of K0r.pol
B the constant coefficient of K0.pol
returns [N, minnormdet], where N = [A, B; C, D] is the one of the 19
matrices in Sp_4 (\Z) determined in Gottschling59 (and given
explicitly on page 134 of Dupont06) such that |det (C*M+D)| is minimal;
and minnormdet is the square of this minimal value as a t_POLMOD
representing an element of K0r */
{
my (a, b, c, d, normdet, minnormdet, index);
a = vector (19, i, matrix (2, 2));
a [16] = matid (2);
a [17] = matid (2);
a [18] = matid (2);
a [19] = -matid (2);
b = vector (19, i, -matid (2));
b [18] = matrix (2, 2);
b [19] = matrix (2, 2);
c = vector (19, i, matid (2));
c [16][2,2] = 0;
c [17][1,1] = 0;
c [18] = [1, -1; -1, 1];
c [19] = [1, -1; -1, 1];
d = vector (19, i, matrix (2, 2));
d [ 2][1,1] = 1;
d [ 3][1,1] = -1;
d [ 4][2,2] = 1;
d [ 5][2,2] = -1;
d [ 6] = matid (2);
d [ 7] = -matid (2);
d [ 8] = [-1, 0; 0, 1];
d [ 9] = [ 1, 0; 0, -1];
d [10] = [ 0, 1; 1, 0];
d [11] = [ 0, -1; -1, 0];
d [12] = [ 1, 1; 1, 0];
d [13] = [-1, -1; -1, 0];
d [14] = [ 0, 1; 1, 1];
d [15] = [ 0, -1; -1, -1];
d [16][2,2] = 1;
d [17][1,1] = 1;
d [18] = matid (2);
d [19] = -matid (2);
for (i = 1, 19,
normdet = complex_norm (matdet (c [i] * M + d [i]));
if (i == 1 || cmp_K0r (normdet, minnormdet, Ar, B) < 0,
minnormdet = normdet;
index = i;
);
);
return ([[a [index], b [index]; c [index], d [index]], minnormdet]);
}
siegel_reduction_Kr (M, Krpol) =
/* M a period matrix with entries in Kr
Krpol the polynomial defining Kr
returns the reduced representative of M in the fundamental domain of
Sp_4 (Z) with respect to the embedding that sends y to the largest root
of Krpol
*/
{
my (Ar, K0rpol, B, U, T, Nmin, ok, Mreal, Mimag, tmp);
Ar = polcoeff (Krpol, 2);
K0rpol = z^2 + Ar * z + polcoeff (Krpol, 0);
B = poldisc (K0rpol) / 16;
M = Mod (M, Krpol);
until (ok,
tmp = lift (M);
Mimag = matrix (2, 2, i, j,
Mod (polcoeff (tmp [i, j], 3) * z + polcoeff (tmp [i, j], 1),
K0rpol));
/* the matrix with entries in K0r such that y*Mimag = i*imag (MKr) */
U = minkowski_reduction_K0r (Mimag);
M = U * M * mattranspose (U);
tmp = lift (M);
Mreal = matrix (2, 2, i, j,
polcoeff (tmp [i, j], 2) * z + polcoeff (tmp [i, j], 0));
T = real_reduction_K0r (Mreal, Ar, B);
M = M + T;
Nmin = gottschling_reduction_Kr (M, Ar, B);
ok = (cmp_K0r (Nmin [2], 1, Ar, B) >= 0);
if (!ok,
M = siegel_transform_apply (Nmin [1], M);
);
);
return (lift (M));
}
/***** Period matrices *****/
period_matrix (cmfield, basis) =
/* cmfield a quartic CM field
basis the symplectic basis of an ideal of K
returns a period matrix with entries in Kr
FIXME: so far, only the dihedral case is implemented */
{
my (Kr, Lrrel, aK, a, asigma, E, V, Z, den);
Kr = cmfield [7];
Lrrel = cmfield [12];
aK = [1, y, y^2, y^3] * basis;
/* We take x as the generating element of Lrrel over Kr. We used to
* have x==-y, but newer text in the article now sets x==y.
*
* For expressing y^sigma, we have y^sigma = y^r - y ;
* previously this was the pari object y+x (y representing the
* generating element of Kr, hence the mathematical object y^r ; x
* representing the generating element of Lrrel over Kr, i.e. -y).
* Now that we change the meaning, it's rather y-x.
*/
/* basis as elements of K */
if(Lrrel==0,
/* Galois */
sigm = cmfield [13];
a = vector(4,i,Mod(aK[i],K.pol));
asigma = vector(4,i,Mod(subst(aK[i],y,sigm),K.pol));
,
/* non-Galois */
a = vector (4, i, Mod (subst (aK [i], y, x), Lrrel [1]));
/* basis as elements of Lrrel = L/Kr */
asigma = vector (4, i, Mod (subst (aK [i], y, y-x), Lrrel [1]));
/* image of a under the non-trivial Galois automorphism of L/Kr */
);
E = matrix (2, 2, i, j, if (i == 1, a [j], asigma [j]));
V = matrix (2, 2, i, j, if (i == 1, a [j+2], asigma [j+2]));
Z = simplify (lift (matinv (V) * E));
den = denominator (Z);
Z = (((1 / den) % Kr.pol) * (den * Z)) % Kr.pol;
/* In the cyclic case, note that we still have something wrt Kr.pol
* which is the minimal polynomial of y, not of y+y^sigma ! */
return (Z);
}
reduced_period_matrix (cmfield, basis) =
/* cmfield a quartic CM field
basis the symplectic basis of an ideal of K
returns a period matrix with entries in Kr that is reduced into
the fundamental domain, where we assume that the CM type is
given by y
FIXME: so far, only the dihedral case is implemented */
{
my (K0, A, B, Kr, Z);
K0 = cmfield [1];
A = polcoeff (K0.pol, 1);
B = polcoeff (K0.pol, 0);
Kr = cmfield [7];
if (cmfield [9] != y,
error ("*** reduced_period_matrix only works for the standard CM type");
);
Z = period_matrix (cmfield, basis);
return (siegel_reduction_Kr (Z, Kr.pol));
}
unique_period_matrix_in_conjugation_orbit (Omega) =
/* Omega a period matrix with entries in Kr, given as polynomials in y
* returns one of Omega or -\bar Omega (which yields theta function values
* that are complex conjugates of those corresponding to Omega),
* according to an arbitrary criterion: the first non-zero coefficient of
* even exponent encountered in Omega is made positive.
*/
{
my (coeff);
for (i = 1, 2,
for (j = 1, 2,
for (k = 0, 1,
coeff = polcoeff (Omega [i,j], 2*k);
if (coeff != 0,
if (coeff > 0,
return (Omega);
,
return (-subst (Omega, y, -y));
);
);
);
);
);
return (Omega);
}
period_matrices (cmfield) =
/* cmfield a quartic CM field
returns 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 (reduced) period matrix, the second one being 1 or 2
for a real or a pair of complex conjugate roots of the class polynomial
*/
{
my (bases, Omegavec);
bases = symplectic_bases (cmfield);
gettime ();
Omegavec = vector (#bases, i, vector (#bases [i], j,
[unique_period_matrix_in_conjugation_orbit
(reduced_period_matrix (cmfield, bases [i][j][1])),
bases [i][j][2]]));
printf ("Time for period matrices: %.1f\n", gettime () / 1000.0);
/* sort to obtain well-defined output */
Omegavec = vecsort (vector (#Omegavec, i,
vecsort (Omegavec [i], cmp_compat)), cmp_compat);
return (Omegavec);
}
/***** Other CM computations *****/
is_totally_positive (K0, xi) =
/* K0: a real quadratic field in the variable z
xi: an element of K0, given as an expression in z
returns 1 if xi is totally positive, 0 otherwise */
{
return (nfelttrace (K0, xi) >= 0 && nfeltnorm (K0, xi) >= 0);
};
is_totally_positive_imaginary (K0, Krel, Phi, xi) =
/* K0 the real subfield of the CM field
Krel the CM field as a relative extension of K0
Phi an element defining a CM type phi
xi an element of K
checks whether phi (xi) lies on the positive imaginary axis
for all phi \in Phi -- which is equivalent to saying that xi/Phi is a
totally positive element of K0.
*/
{
return (is_totally_positive (K0, rnfeltdown (Krel, xi/Phi)));
};
explode_classgroup (K) =
/* K a CM field
returns a vector containing all class group elements of K
each entry is itself a vector with two components: the ideal and its
decomposition (as a column vector) over the class group generators */
{
my (res, h, orders, gens, pow);
res = vector (K.clgp [1]);
res [1] = [[1,0,0,0; 0,1,0,0; 0,0,1,0; 0,0,0,1],
vectorv (# K.clgp [2])];
h = 1; /* current cardinality of res */
orders = K.clgp [2];
gens = K.clgp [3];
for (i = 1, #orders,
/* treat generator i and its powers */
pow = gens [i];
for (j = 1, orders [i] - 1,
/* treat power j of generator i */
for (k = 1, h,
res [j*h + k] = [idealred (K, idealmul (K, res [k][1], pow)),
res [k][2]];
res [j*h + k][2][i] = j;
);
pow = idealred (K, idealmul (K, pow, gens [i]));
);
h = h * orders [i];
);
return (res);
};
explode_abstractgroup (v) =
/* v a vector of the elementary divisors of an abstract group
returns a vector containing all elements of the abstract group
each entry is itself a column vector w of the same length as v,
with 0 <= w [i] < v [i] */
{
my (res, h, i);
res = vector (prod (i = 1, #v, v [i]));
res [1] = vectorv (#v);
h = 1; /* current cardinality of res */
for (i = 1, #v,
/* treat generator i and its powers */
for (j = 1, v [i] - 1,
/* treat power j of generator i */
for (k = 1, h,
res [j*h + k] = res [k];
res [j*h + k][i] = j;
);
);
h = h * v [i];
);
return (res);
};
cmp_norm (v1, v2) =
/* v1, v2 entries of a vector as returned by explode_classgroup
returns -1, 0, 1 depending on whether v1 <, =, > v2,
where the order is given first by the norm of the ideal, and,
for equal norm, lexicographically with respect to the class group
decomposition
*/
{
my (n1, n2);
n1 = prod (i = 1, 4, v1 [1][i, i]);
n2 = prod (i = 1, 4, v2 [1][i, i]);
if (n1 < n2,
return (-1);
, if (n1 > n2,
return (1);
););
for (i = 1, #v1 [2],
if (v1 [2][i] < v2 [2][i],
return (-1);
,
if (v1 [2][i] > v2 [2][i],
return (1);
););
);
return (0);
};
cmp_weight (v1, v2, w) =
/* v1, v2 entries of a vector as returned by explode_abstractgroup
returns -1, 0, 1 depending on whether v1 <, =, > v2,
where the order is given first by the product of the entries weighted
by w (a vector of the same lenght as v1 and v2), and, for equal weight,
lexicographically
*/
{
my (n1, n2);
n1 = sum (i = 1, #v1, v1 [i] * w [i]);
n2 = sum (i = 1, #v2, v2 [i] * w [i]);
if (n1 < n2,
return (-1);
, if (n1 > n2,
return (1);
););
for (i = 1, #v1,
if (v1 [i] < v2 [i],
return (-1);
,
if (v1 [i] > v2 [i],
return (1);
););
);
return (0);
};
find_basepoint (cmfield, Phi) =
/* cmfield as output by init_cmfield
Phi a CM type of K given by two complex embeddings of x
returns a vector [b, xi, bdlog] such that b is an ideal of OK,
xi as in Section 3.1, p. 45 of Streng10,
and bdlog the decomposition of b over the classgroup of OK
*/
{
my (K0, K, Krel, dlogconj, dlogdiff, Clexp,
b, bdlog, prin, xi0, units, xi);
K0 = cmfield [1];
Krel = cmfield [3];
K = cmfield [2];
Krel = cmfield [3];
dlogconj = Mat (vector (#K.clgp [2], i,
bnfisprincipal (K, complex_conjugate_ideal (K, K.clgp [3][i]), 0)));
dlogdiff = bnfisprincipal (K, K.diff, 0);
Clexp = vecsort (explode_classgroup (K), cmp_norm);
for (i = 1, #Clexp,
bdlog = Clexp [i][2];
prin = bdlog + dlogconj * bdlog + dlogdiff;
for (j = 1, #prin,
prin [j] %= K.clgp [2][j];
);
if (!prin,
b = Clexp [i][1];
prin = bnfisprincipal (K, idealinv (K, idealmul (K, K.diff,
idealmul (K, b, complex_conjugate_ideal (K, b))))) [2];
xi0 = nfbasistoalg (K, prin);
units = [1, -1, K.fu [1].pol, -K.fu [1].pol];
for (j = 1, 4,
xi = (units [j] * xi0) % K.pol;
if (is_totally_positive_imaginary (K0, Krel, Phi, xi),
return ([b, xi, bdlog]);
);
);
);
);
};
/* {{{ Shimura group: basic stuff */
shimura_group_identity_element(cm) =
/* trivial */
{
my (K,K0);
K0 = cm[1];
K= cm[2];
return([idealhnf(K,1),Mod(1,K0.pol)]);
}
shimura_group_arbitrary_lift_of_ideal(cm, idl) =
/* cm: as output by init_cmfield
idl: an ideal of K (= cm[2])
returns an arbitrary lift of a as an element of the Shimura group
*/
{
my (K, K0, Krel, idlnorm, index, g, good);
K = cm [2];
K0 = cm[1];
Krel = cm [3];
eps0 = cm [4];
idlnorm=rnfidealnormrel(Krel,rnfidealabstorel(Krel,K.zk*idl));
index=bnfisprincipal(K0,idlnorm);
g=K0.zk*index[2];
good=select(x->is_totally_positive(K0,x*g),[1,-1,eps0,-eps0]);
g*=good[1];
return ([idl,Mod(g,K0.pol)]);
}
shimura_group_is_element(cm, a) =
/* check what its name suggests */
{
my (K0,idl,g,idlnorm);
K0 = cm[1];
K0narrow = bnrinit (K0, [1, [1, 1]], 1);
Krel = cm [3];
idl=a[1];
g=a[2];
idlnorm=rnfidealnormrel(Krel,rnfidealabstorel(Krel,K.zk*idl));
return (is_totally_positive(K0, g) && bnrisprincipal(K0narrow, idealdiv(K0,idlnorm,g), 0) == 0);
}
shimura_group_element_mul(cm, a, b) =
{
my (K,aidl,ag,bidl,bg);
K=cm[2];
aidl=a[1];
ag=a[2];
bidl=b[1];
bg=b[2];
return ([idealmul(K,aidl,bidl),ag*bg]);
}
shimura_group_element_pow(cm, a, n) =
{
my (K0,K,aidl,ag,nx,i,np);
K0=cm[1];
K=cm[2];
if(n==0, return(shimura_group_identity_element(cm)););
if(type(a)=="t_VEC",
if(n*sign(n)>100,
warning("Prefer precomputations using shimura_group_element_powtable for raising to the power ", n);
error("no, really, it's not serious");
);
/* normal code -- we actually compute the power using the
* generator, and perhaps later we might consider reducing.
* Note that doing this is inherently dangerous. We might fall
* short of precision when handling a beast such as ag^n !!
*/
aidl=a[1];
ag=a[2];
ax=[idealpow(K,aidl,n),ag^n];
if(n*sign(n)>100,
warning("resulting log-trace of g^n: ", round(log(abs(nfelttrace(K0, ax[2])))));
);
,
/* Otherwise a is in fact a list of powers. */
np=#a;
/* Elements of this list correspond to reduced versions of:
[1] a^1
[2] a^2
[3] a^4
[i] a^(2^(i-1))
[np] a^(2^(np-1))
*/
nx=n*sign(n);
if (2^np<=nx,error("*** precomputed powers too short"););
ax=shimura_group_identity_element(cm);
for(i=1,np,
if (bittest(nx, i-1),
ax=shimura_group_element_mul(cm, ax, a[i]););
);
if(n<0,ax=[idealpow(K,ax[1],-1),ax[2]^-1];);
ax=shimura_group_element_reduce(cm, ax);
);
return (ax);
}
shimura_group_element_eq(cm, a, b) =
/* cm: as output by init_cmfield
tell whether a and b are representatives of the same class in the
Shimura group.
*/
{
my (K0, K,Krel,aidl,ag,bidl,bg,v,s,t);
K0 = cm[1];
K=cm[2];
Krel = cm [3];
aidl=a[1];
ag=a[2];
bidl=b[1];
bg=b[2];
v=bnfisprincipal(K,idealdiv(K, aidl, bidl),3);
if(v[1],return(0));
s=K.zk*v[2];
t=Mod(s*complex_conjugate_element(s),K.pol);
t=rnfeltdown(Krel,rnfeltabstorel(Krel,t));
if(!is_totally_positive(K0,t), error("argh"););
v=bnfisunit(K0,ag/bg/t);
if(#v!=2,error("argh"););
if(type(v[2])!="t_INTMOD",error("argh"););
/* Hmm. I'm having problems with the above. weird.
if(v[2]!=0,error("argh");); */
return(v[1]%2==0);
}
shimura_group_element_is_identity(cm, a) =
{
return(shimura_group_element_eq(cm, a, shimura_group_identity_element(cm)));
}
quadratic_field_unit_reduce(K0, g, m)=
/* K0: a real quadratic field.
* g: an element of K0 of type t_POLMOD
* m: any integer, indicating a power of the fundamental unit in K0.
*/
{
my (e, ltrace_g, ltrace_gx, x, k);
e = K0.fu[1];
e = sign(nfelttrace(K0,e))*e;
ltrace_g=log(abs(nfelttrace(K0, g)));
/* Use this very rough approximation only when the trace of the
* element to be reduced dominates.
*/
if(ltrace_g > 4*m*log(abs(nfelttrace(K0, e))),
ltrace_gx=ltrace_g;
x=0;
while(abs(ltrace_gx-ltrace_g) < 4,
if(x,x*=2,x=1);
ltrace_gx=log(abs(nfelttrace(K0, g*e^(m*x))));
);
k=round(x*ltrace_g/(ltrace_gx-ltrace_g));
g1=(g/e^(m*k));
, /* else */
g1=g;
);
g1=Mod(g1, K0.pol);
my (vg,ve);
vg=vector(2);
for(i=1,2,
/* example 35,48 reaches here (in case default precision has been
* set to 32) */
while(abs(subst(g1.pol,z,K0.roots[i]))==0,
op = default(realprecision);
np = 2*op;
warning("*** Raising precision: ", op, " -> ", np);
default(realprecision, np);
K0 = bnfinit(K0.pol);
);
vg[i] = log(abs(subst(g1.pol,z,K0.roots[i])));
);
ve=vector(2,i,m*log(abs(subst(e.pol,z,K0.roots[i]))));
/* We want to approximate the value x such that (g-xe,e)=0 -- which
* is of course simply x=(g,e)/(e,e) */
x=round((vg[1]*ve[1]+vg[2]*ve[2])/(ve[1]^2+ve[2]^2));
g2=g1/e^(m*x);
return(g2);
}
shimura_group_element_reduce(cm, a) =
/* cm: as output by init_cmfield
gives an equivalent, but smaller representative of the given class in
the Shimura group
*/
{
my (K, Krel, K0, idl, g, ridl, s, t);
K=cm[2];
Krel = cm [3];
K0=cm[1];
idl=a[1];
g=a[2];
ridl = idealred (K, [idl, 1]);
s = nfbasistoalg (K, ridl [2]);
t = (s * complex_conjugate_element (s)).pol;
t = if (poldegree (t) > 0, polcoeff (t, 2) * z, 0) + polcoeff (t, 0);
g=Mod (g/t, K0.pol);
g=quadratic_field_unit_reduce(K0, g, 2);
return ([ridl [1], g]);
}
shimura_group_factorback(cm, f, e) =
/* f and e are lists of the same length (well, only the length of f
* counts). f is a list of elements of the Shimura group (*), and e
* contains exponents.
*
* (*) As an alternate syntax, and as per the special handling in
* shimura_group_element_pow, f may also be a list of lists of powers of
* elements of the Shimura group. Such lists may be precomputed using
* shimura_group_element_powtable. This allows dramatically faster
* computations for large components in the class group. OTOH, it's
* useful only for repeated calls to the factorback routine with the same
* set of generators.
*/
{
my (n,x,y);
n=#f;
x=shimura_group_identity_element(cm);
for(i=1,n,
y=shimura_group_element_pow(cm,f[i],e[i]);
x=shimura_group_element_mul(cm,x,y);
);
return (x);
}
shimura_group_element_powtable(cm, a, n) =
/* given an element, and maximum power n, precompute reduced
* representatives for a^(2^k) for all k
* such that 2^k <= n ; this used for later computation of powers in a
* reduced fashion.
*/
{
ax=List();
k=1;
listput(ax,a);
while(2*k<=n,
k*=2;
a=shimura_group_element_mul(cm,a,a);
a=shimura_group_element_reduce(cm,a);
listput(ax,a););
return(ax);
}
/* }}} */
shimura_group_structure(cm) = /* {{{ */
/* cm: as output by init_cmfield
* This returns information identifying the structure of the Shimura
* group. Namely, the following properties of the group are returned:
* [1] Smith invariants d1, ..., dk, with 1|dk|d_{k-1}|...|d1
* [2] generators a1,...,ak, with order(ai)=di
* The rest is auxiliary data used for conversions.
* [3] generators b1,...,bn (with n>=k) (see also below).
* [4] conversion matrix P which gives the bnfisprincipal results
* for the b1[1]...bn[1] (in columns). More precisely, arbitrary
* lifts for the class group generators having been chosen, this
* matrix indicates how a given element of the class group may
* be lifted into she Shimura group in a way which is consistent
* with these arbitrary choices. If q is such that P*q = x, then
* an ideal for which bnfprincipal(I)==x may be consistently
* lifted to (b1...bn)*q.
* [5] unimodular conversion matrix U such that
* (a1...ar,1..1)=(b1..bn)*(U^-1)
*/
{
my (K,Krel,K0,eps0,gens,n,a,m,i,j,v,k,krel);
my (kgens,Cgens,Cgensx,id,nn0,mo,ms,mu,mv,md,CKgens,mui);
K0=cm[1];
K=cm[2];
Krel = cm [3];
eps0 = cm [4];
/* First task: compute the kernel of the norm map from the class
* group of K to the narrow class group of K0 */
K0narrow = bnrinit (K0, [1, [1, 1]], 1);
n=#K.clgp.cyc;
nn0=#K0narrow.clgp.cyc;
m=matrix(nn0,n+nn0);
for(i=1,nn0,m[i,n+i]=K0narrow.clgp.cyc[i];);
for(i=1,n,
idl=K.clgp.gen[i];
ridl=rnfidealabstorel(Krel,K.zk*idl);
idlnorm=rnfidealnormrel(Krel,ridl);
v=bnrisprincipal(K0narrow,idlnorm);
for(j=1,nn0,m[j,i]=v[1][j];);
);
/* Take the part of the hnf which creates the zeroes */
v=mathnf(m,1)[2];
/* Notice that the coefficients of combinations appear in columns */
krn=matrix(n,n,i,j,v[i,j]);
kgens=vector(n,j,idealfactorback(K,K.clgp.gen,krn[,j]));
/* Now all ideals in kgens are known to map via the norm to the
* trivial class in the narrow class group.
* It's not exactly sufficient, because we'd like to know the
* relations among these elements in the ambient group Cl(K). Again
* a hnf computation
*/
m=concat(krn,matdiagonal(K.clgp.cyc));
v=mathnf(m,1)[2];
krel=matrix(n,n,i,j,v[i,j]);
/* We need to expand this presentation based on relations by the
* finer information on how this translates into the unit part. */
Cgens=vector(n,i,shimura_group_arbitrary_lift_of_ideal(cm,kgens[i]));
mo=K.cyc;
if(nfeltnorm(K0,eps0)<0,
/* This is rather trivial if the unit part is trivial, of course. */
m=krel;
,
/* otherwise there's more info to put in */
m=matid(n+1);
for(i=1,n,for(j=1,n,m[i,j]=krel[i,j];););
m[n+1,n+1]=2;
Cgens=concat(Cgens,[[idealhnf(K,1),eps0]]);
mo=concat(mo, [2]);
);
/* Now fix each relation ! We need some equivalent of factorback on
* the Shimura group. Cf related function. */
/* First, we need to precompute tables of each of these. */
Cgensx=vector(#Cgens,i,shimura_group_element_powtable(cm, Cgens[i], mo[i]));
id=shimura_group_identity_element(cm);
for(j=1,n,
v=shimura_group_factorback(cm,Cgensx,m[,j]);
/* There isn't so much choice, so either we have 1, or we don't */
if(!shimura_group_element_eq(cm,v,id),
if(nfeltnorm(K0,eps0)<0,error("argh"););
m[n+1,j]=1;
);
);
ms=matsnf(m,1);
mu=ms[1];
mv=ms[2];
md=ms[3];
r=#select(i->md[i,i]!=1,vector(matsize(md)[1],i,i));
smith_invs=vector(r,i,md[i,i]);
/* There's an additional trick. The entries in mu^-1 may exceed the
* orders of the elements in Cgensx. We have C^m=0, thus
* C^(u^-1*u*m)=C^(u^-1)^d=0, but u^-1 isn't necessarily reduced mod
* m. So we should consider whether there's an equivalent set of
* generators C^(u^-1-m*z), with the entries in u^-1+m*z properly
* reduced. So we need to solve m*z=u^-1
*/
/* Note that the temptation to optimize this based on the fact that m
* is almost diagonal is bogus. m comes from a hnf computation, and
* is in fact upper triangular + some coefficients in the last row.
* Example for which m[1,3] != 0: 352 6257
*/
mui=mu^-1;
my (kk);
kk=matsolve(m,mu^-1);
kk=matrix(#mo,#mo,i,j,round(kk[i,j]));
mui=mu^-1-m*kk; /* apply */
CKgens=vector(r,j,shimura_group_factorback(cm,Cgensx,mui[,j]));
/* Do we need this ?
CKgensx=vector(#CKgens,i,shimura_group_element_powtable(cm, CKgens[i], smith_invs[i]));
*/
return([smith_invs,CKgens,Cgensx,krn,mu]);
}
/* }}} */
/* {{{ Shimura group: complex conjugation action */
shimura_group_element_conjugate(cm, a)=
/* return the complex conjugate (abar, x) of (a,x) */
{
K = cm [2];
return([complex_conjugate_ideal(K, a[1]),a[2]]);
}
shimura_group_precompute_conjugation_action(cm)=
/* return a vector indicating the coefficients of the complex conjugation
* action on the smith generators of the CKdata group. Note that since
* complex conjugation commutes with group morphisms here, we know that its
* matrix is diagonal. As in other cases, we're not returning the members
* as INTMOD's, but it's trivially fixed by the caller if needed.
*/
{
my(CKdata,CKgens,smith_invs,n,i,j);
CKdata = cm [16];
CKgens=CKdata[2];
smith_invs=CKdata[1];
n=#smith_invs;
m=matrix(n,n);
for(i=1,n,
my (a,v);
a=shimura_group_element_conjugate(cm,CKgens[i]);
v=shimura_group_element_dlog(cm,a);
for(j=1,n,m[i,j]=v[j];);
);
if(!matisdiagonal(m),error("argh"););
return(vector(n,i,m[i,i]));
}
/* }}} */
shimura_group_element_dlog(cm, a) = /* {{{ */
/* given a cm and the corresponding Shimura group info as computed
* by the above function, returns the decomposition of a w.r.t the
* elementary divisors of the Shimura group structure.
*/
{
my (K, m,Cgensx,CKgens,krn,mu,idl,q,i,n,ax,qx,CKdata);
/* first step: write the ideal part with respect to the basis of the
* norm map to the narrow class group */
CKdata = cm [16];
K = cm [2];
smith_invs=CKdata[1];
CKgens=CKdata[2];
Cgensx=CKdata[3];
krn=CKdata[4];
mu=CKdata[5];
idl=a[1];
q=matsolve(krn,bnfisprincipal(K,idl)[1]);
/* FIXME: Seems that it's mandatory to append the Smith invariants to
* krn ? */
n=#q;
/* Note that #Cgensx may be n+1. But we know no coordinate
* value for it yet */
qx=vector(#Cgensx);
for(i=1,n,qx[i]=q[i];);
/* recreate the ideal, using Cgensx */
ax=shimura_group_factorback(cm,Cgensx,qx);
if(!shimura_group_element_eq(cm,a,ax),
if(#Cgensx!=n+1,error("argh"););
qx[n+1]=1;
);
qs=mu*Col(qx);
r=#smith_invs;
return(vector(r,i,qs[i] % smith_invs[i]));
/* return(vector(r,i,Mod(qs[i],smith_invs[i]))); */
}
/* }}} */
shimura_group_type_norm_subgroup(cm) = /* {{{ */
/* computes the subgroup of the Shimura subgroup which is spanned by the
* image of the class group of the reflex CM field under the reflex type
* norm. Returns:
* [1] the smith invariants of that group,
* [2] corresponding generators (as members of the reflex class group),
* [3] generators of cosets in the ambient group,
* [4] dlog vectors for the subgroup generators as combination of the
* ambient group generators.
* The quotient being a 2-group, all the coset generators have order 2
* (i.e., if x coset generators are returned, there are 2^x cosets in
* total).
*/
{
my (smith_invs, CKgens, T, CKdata, r, Kr, n, m, v, u, tmp, h, pick);
my (ms, mu, mv, md);
my (Hr, Hsmith_invs, mui, Sgens, coset_gens);
CKdata = cm[16];
smith_invs=CKdata[1];
CKgens=CKdata[2];
r = #smith_invs;
Kr = cm [7];
n = #Kr.clgp.cyc;
m0 = matrix(r,n);
T=vector(n,j,type_norm(cm,Kr.clgp.gen[j]));
for(j=1,n, v=shimura_group_element_dlog(cm,T[j]); for(i=1,r,m0[i,j]=v[i];););
m=concat(m0, matdiagonal(smith_invs));
/* These form the subgroup generators. We need to know the relations
* between these generators: hnf. We're only interested in the
* transformation matrix. The hnf itself gives the quotient, but
* we're not too interested at this point (except if we want to
* verify that it's a 2-group). */
my(tmp);
tmp=mathnf(m,1);
h=tmp[1];
u=tmp[2];
/* Still, the info on the coset representatives is there to pick, so
* there's really no reason not to. We have CKgens^m = [primary
* generators of S, and 0], thus CKgens^(m*u) = CKgens^h. The coset
* generators are thus the members of the CKgens array which get
* raised to a non trivial power in order to belong to S.
*/
pick=select(i->h[i,i]==2,vector(r,i,i));
coset_gens=vector(#pick,j,CKgens[j]);
u=matrix(n,n,i,j,u[i,j]);
/* Now this square matrix is the relation matrix among my generators.
* the Smith form eventually gives me the structure. */
ms=matsnf(u,1);
mu=ms[1];
mv=ms[2];
md=ms[3];
Hr=#select(i->md[i,i]!=1,vector(matsize(md)[1],i,i));
Hsmith_invs=vector(Hr,i,md[i,i]);
mui=mu^-1;
/* We give the generators as ideals of the reflex class group. */
Sgens=vector(Hr,j,idealfactorback(Kr,Kr.gen,mui[,j]));
/* We may quite easily provide the matrix which gives the subgroup
* generators as combinations of the ambient group generators (such
* that CKgens^m = Sgens */
my (tmat);
tmat=matrix(r,n,i,j,m[i,j])*mui;
tmat=matrix(r,Hr,i,j,tmat[i,j] % smith_invs[i]);
return ([Hsmith_invs, Sgens, coset_gens, tmat]);
}
/* }}} */
/* {{{ Shimura group: miscellaneous */
shimura_group_allproducts(cm, g, m) =
/* cm: as output by init_cmfield
g: a list of elements of the Shimura group.
m: a list of exponents.
returns the list of all \prod g[j]^x[j], where 0<=x[j] ", np);
default(realprecision, np);
K0 = bnfinit(pol, 1);
);
);
return(K0);
}
abelian_group_name(invs) =
{
my (st);
st="";
for(i=1,#invs,
if(length(st)>0,
st=concat(concat(st," x C"), invs[i]);
,
st=concat("C", invs[i]);
);
);
return(st);
}
/* {{{ CM field basics */
init_cmfield (A, B) =
/* works as init_cmfield_basic, but adds the following to the output vector:
16 CKdata: abstract data for working with the Shimura group
17 Sdata: abstract data for working with the subgroup of Shimura
group which is reached by the reflex typenorm map.
*/
{
my (cm, CKdata, Sdata);
cm = init_cmfield_basic(A, B);
if(cm[11]==0,
print("Galois group = C4");
,
print("Galois group = D4");
);
gettime();
/* create the preliminary list */
CKdata=shimura_group_structure(cm);
print("Shimura group = ", abelian_group_name(CKdata[1]));
cm = concat (cm, [CKdata]);
Sdata=shimura_group_type_norm_subgroup(cm);
print("Type norm subgroup = ", abelian_group_name(Sdata[1]));
cm = concat (cm, [Sdata]);
return (cm);
}
/* }}} */
/* {{{ Enumeration of CM orbits */
triple_new(cm, bp, sh, flag) =
/* cm: 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
* returns a triple [a, xi, flag] = sh*bp
*
* Watch out for a typo in BrGrLa11: The action of sh on bp is by
* multiplying the first components, but _dividing_ the second ones.
*/
{
my (K, Krel, a, xi);
K = cm [2];
Krel = cm [3];
a = idealmul (K, bp [1], sh [1]);
xi = (bp [2] / rnfeltup (Krel, sh [2])).pol; /* .pol to remove the modulus */
return ([a, xi, flag]);
}
reduce_basepoint(cm, Phi, bp)=
{
my (K0, K, Krel, k, unitindex, xi, scale);
K0 = cm [1];
K = cm [2];
Krel = cm [3];
unitindex = cm [5];
xi=bp[2];
k = 2/unitindex;
scale =rnfeltdown(Krel,Mod(xi/Phi,K.pol));
scale = quadratic_field_unit_reduce(K0,scale, k);
xi=(Phi*rnfeltup(Krel,scale)).pol;
return([bp[1], xi]);
}
symplectic_bases(cm) =
/* cm a quartic CM field
If tofile==false, prints on screen the symplectic bases in the format
defined in README.format.
Otherwise, prints them into the file named D_A_B.in, where D is the
discriminant of the real-quadratic subfield. */
{
my (A, B, K0r, Arred, Brred, K0, K, unitindex, Phi,
D, file, cosets, m, noC, bp, trip, basis, den);
K0 = cm [1];
K0r = cm [6];
A = polcoeff (K0.pol, 1);
B = polcoeff (K0.pol, 0);
Arred = cm [14];
Brred = cm [15];
K = cm [2];
unitindex = cm [5];
Phi = cm [9];
noC = cm [10];
D = K0.disc;
if (tofile,
file = concat ([D, "_", A, "_", B, ".in"]);
file = concat (Strexpand("$DATA_PREFIX"), file);
delete_file (file),
/* else */
file = 0;
);
my (CKdata,CKgens,Sdata,coset_gens,Sgens,X,smith_ings,xL,lambda);
CKdata = cm [16];
CKgens = CKdata[2];
Sdata = cm [17];
coset_gens = Sdata[3];
X = shimura_group_allproducts(cm, coset_gens, vector(#coset_gens,i,2));
smith_invs=CKdata[1];
xL = explode_abstractgroup(Sdata[1]);
lambda = shimura_group_precompute_conjugation_action(cm);
bp = find_basepoint (cm, Phi);
bp = reduce_basepoint (cm, Phi, bp);
my (coset_lists,steps,steps_values);
coset_lists=List();
steps=Set();
post_setup = gettime();
/* coset_prep seems to be so ridiculously small that it's not really
* useful to track this measurement...
coset_prep=vector(#X);
*/
for (l = 1, #X,
my (s,h,xh);
/* The triples in this coset are of the form:
* (bp[1],bp[2]) . (s[1],s[2]) * (a[1],a[2])
* where a runs through the subgroup S of the Shimura group.
* Thus the ideal part of a given element is (dropping [1]'s):
* (bp*s*a). The conjugate may thus be written as:
* (bp*s*(bpbar*sbar/bp/s)*abar).
*
* Thus given an abstract representation xa of a, and given a
* (matrix) multiplier lambda which corresponds to the complex
* conjugation ; given again an abstract representation xh for the
* Shimura group element (bpbar*sbar/bp/s,1), we need thus consider
* only one representative of each equivalence class of S modulo
* the relation: xa <--> xh + lambda * xa
*
* Points which are stable under xa <--> xh + lambda * xa
* correspond to real roots of the class polynomial. These receive
* the ``flag'' value of 1.
*
* Other points correspond to conjugate complex roots. Naturally we
* keep only one representative of such pair, and the ``flag''
* value is 2.
*/
/* So we need first to compute xh. Not too difficult. */
s=X[l];
h=idealmul(K, bp[1], s[1]);
h=idealdiv(K, complex_conjugate_ideal(K, h), h);
h=[h, Mod(1, K0.pol)];
xh=shimura_group_element_dlog(cm, h);
/* next, iterate through all abstract representations of the
* subgroup elements. Look for those for which the conjugate is
* also present.
*/
my (S,pick,xa_S,xa,xa2);
S=Set();
pick=List();
for(i=1,#xL,
xa_S = xL[i]~;
/* xL[i] is a vector of exponents w.r.t the generators of the
* subgroup S. Here, we want them as members of CK. The
* conversion matrix is Sdata[4] */
xa = xa_S * mattranspose(Sdata[4]);
if (#xa_S,
xa=vector(#smith_invs,k,Mod(xa[k],smith_invs[k]));
, /* else */
xa=vector(#smith_invs); /* degenerate case */
);
xa2=vector(#xa,i,xh[i] + xa[i]*lambda[i]);
if(xa==xa2, /* a real root */
/* write1("/dev/stderr", "r"); */
listput(pick, [lift(xa_S), 1]);
S = setunion (S, Set ([xa])); /* in fact useless */
, /* else: */ /* complex conjugate pair */
if (!setsearch (S, xa),
/* write1("/dev/stderr", "c"); */ /* not yet met */
S = setunion (S, Set ([xa, xa2]));
listput(pick, [lift(xa_S), 2]);
,
/* write1("/dev/stderr", "-"); */ /* already met */
);
);
);
/* write1("/dev/stderr", "\n"); */
/* Last thing. It's easy to enumerate members of a group when they
* come in increasing lexicographical order. Here, we have gaps.
* Therefore we precompute the steps. */
steps=setunion(steps, Set(vector(#pick-1,i,pick[i+1][1]-pick[i][1])));
steps=setunion(steps, Set([pick[1][1]]));
listput(coset_lists, pick);
/*
coset_prep[l] += gettime();
*/
);
Sgens=vector(#Sdata[2],i,type_norm(cm,Sdata[2][i]));
Sgensx=vector(#Sgens,i,shimura_group_element_powtable(cm, Sgens[i], Sdata[1][i]));
/* print(steps); */
steps_values=vector(#steps,i,shimura_group_factorback(cm, Sgensx, eval(steps[i])));
/* steps_values=vector(#steps,i,shimura_group_factorback(cm, CKgens,
* eval(steps[i]))); */
post_setup += gettime();
printf ("post_setup = %.1f\n", post_setup / 1000.0);
bases=vector(#coset_lists);
/* starting output */
for (l = 1, #coset_lists,
C=coset_lists[l];
b=vector(#Sdata[1]);
a = X[l]; /* starts this coset */
bases[l] = List();
for (k=1, #C,
p = C[k][1] - b;
flag = C[k][2];
if(p != 0,
idx = setsearch(steps, p);
if (!idx, error("argh: ", p, " not found ??"););
a = shimura_group_element_mul(cm, a, steps_values[idx]);
a = shimura_group_element_reduce(cm, a);
b = C[k][1];
);
trip = triple_new (cm, bp, a, flag);
my (af,basis);
af = switch_ideal_basis_for_symplectic_input(K, trip[1]);
basis = symplectic_basis (K, af, trip [2]);
listput(bases[l], [basis, trip[3]]);
);
/*
coset_prep[l] += gettime();
printf ("coset_prep%d = %.1f\n", l, coset_prep[l] / 1000.0);
*/
);
return(bases);
}
/* }}} */