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