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); */