compute Frobenius matrices together with the bases

modcurve/torsion-subscheme.c

 #include <pari/pari.h>
+#include <pari/paripriv.h>  /* ucoeff */
 
 #include "curve.h"
 #include "divisor.h"
@@ -79,11 +80,11 @@ check_multiplicity(GEN group, long l, long p, GEN f_primary,
 
 static GEN
 find_extension(GEN group, unsigned long l, unsigned long p,
-	       GEN f, GEN *projector, GEN *projector_dual,
+	       GEN f, GEN *f_dual, GEN *projector, GEN *projector_dual,
 	       unsigned long max_degree) {
   pari_sp av = avma;
   GEN charpoly_Frob = modular_curve_characteristic_polynomial_frobenius(group, p);
-  GEN f_dual, f_primary, extension, x;
+  GEN f_primary, extension, x;
   long extension_degree;
 
   x = pol_x(gvar(f));
@@ -97,24 +98,24 @@ find_extension(GEN group, unsigned long l, unsigned long p,
   extension = finite_field_extension(extension_degree, p, NULL);
 
   /* f = x^2 + a*x + b  ==>  g = x^2 + (p*a/b)*x + p^2/b */
-  f_dual = shallowcopy(f);
-  gel(f_dual, 3) = gmulsg(p, gdiv(gel(f, 3), gel(f, 2)));
-  gel(f_dual, 2) = gdiv(sqru(p), gel(f, 2));
+  *f_dual = shallowcopy(f);
+  gel(*f_dual, 3) = gmulsg(p, gdiv(gel(f, 3), gel(f, 2)));
+  gel(*f_dual, 2) = gdiv(sqru(p), gel(f, 2));
 
   factor_coprime(charpoly_Frob, f, &f_primary, projector, l);
   check_multiplicity(group, l, p, f_primary, extension_degree);
 
   err_printf("characteristic polynomial of Frobenius: %Ps\n", f);
-  if(gequal(f_dual, f))
+  if(gequal(*f_dual, f))
     *projector_dual = *projector;
   else {
     err_printf("characteristic polynomial of Frobenius for the"
-	       " dual representation: %Ps\n", f_dual);
-    factor_coprime(charpoly_Frob, f_dual, &f_primary, projector_dual, l);
+	       " dual representation: %Ps\n", *f_dual);
+    factor_coprime(charpoly_Frob, *f_dual, &f_primary, projector_dual, l);
     check_multiplicity(group, l, p, f_primary, extension_degree);
   }
 
-  gerepileall(av, 3, &extension, projector, projector_dual);
+  gerepileall(av, 4, &extension, f_dual, projector, projector_dual);
   return extension;
 }
 
@@ -174,19 +175,21 @@ point_index(GEN J, GEN V, GEN P) {
   Naïve algorithm for finding bases of the subspaces of the
   l-torsion defined by the elements of  proj.
 */
-static GEN
-find_bases(GEN J, unsigned long l, GEN proj, int tries) {
-  GEN bases, P, Q, R, V1, V2;
-  long i, j, k, n = lg(proj) - 1, found = 0;
+static void
+find_bases(GEN J, unsigned long l, GEN f, GEN proj, int tries,
+	   GEN *bases, GEN *matrices) {
+  GEN P, Q, R, V1, V2, M;
+  long i, j, k, n = lg(proj) - 1, found = 0, d;
   pari_sp av = avma;
 
-  bases = cgetg(n + 1, t_VEC);
+  *bases = cgetg(n + 1, t_VEC);
+  *matrices = cgetg(n + 1, t_VEC);
   for (i = 1; i <= n; i++)
-    gel(bases, i) = cgetg(1, t_VEC);
+    gel(*bases, i) = cgetg(1, t_VEC);
   while (tries-- > 0 && found < 2*n) {
     P = random_torsion_point(J, l);
     for(i = 1; i <= n; i++) {
-      j = lg(gel(bases, i)) - 1;
+      j = lg(gel(*bases, i)) - 1;
       if (j == 2)
 	continue;
       err_printf("applying Frobenius polynomial %Ps\n", gel(proj, i));
@@ -197,59 +200,43 @@ find_bases(GEN J, unsigned long l, GEN proj, int tries) {
 	V1 = multiples(J, Q, l);
 	R = jacobian_Frob(J, Q, 1);
 	if ((k = point_index(J, V1, R)) != 0) {
-	  gel(bases, i) = mkvec(V1);
+	  gel(*bases, i) = mkvec(V1);
+	  gel(*matrices, i) = mkmat2(mkvecsmall2(k - 1, 0), zero_Flv(2));
 	  found++;
 	} else {
 	  V2 = multiples(J, R, l);
-	  gel(bases, i) = mkvec2(V1, V2);
+	  gel(*bases, i) = mkvec2(V1, V2);
+	  gel(*matrices, i) = RgM_to_Flm(matcompanion(gel(f, i)), l);
 	  found += 2;
 	}
       } else {
-	V1 = gmael(bases, i, 1);
+	V1 = gmael(*bases, i, 1);
 	if ((k = point_index(J, V1, Q)) == 0) {
 	  V2 = multiples(J, Q, l);
-	  gel(bases, i) = mkvec2(V1, V2);
+	  gel(*bases, i) = mkvec2(V1, V2);
+	  /* with S = V1[1]:  Frob(S) = a*S,  Frob(Q) = b*S + d*Q */
+	  M = gel(*matrices, i);
+	  d = Fl_neg(Fl_add(ucoeff(M, 1, 1),
+			    Rg_to_Fl(gmael(f, i, 3), l), l), l);
+	  ucoeff(M, 2, 2) = d;
+	  R = jacobian_subtract(J, jacobian_Frob(J, Q, 1),
+				jacobian_multiply(J, Q, centerlift(gmodulss(d, l))));
+	  if ((k = point_index(J, V1, R)) == 0)
+	    pari_err(e_MISC, "inconsistent Frobenius action");
+	  ucoeff(M, 1, 2) = k - 1;
 	  found++;
 	}
       }
     }
     if (gc_needed(av, 1))
-      bases = gerepileupto(av, bases);
+      gerepileall(av, 2, bases, matrices);
   }
   if (found < 2*n) {
     if (GIVEUP)
       pari_err(e_MISC, "too many tries to find bases");
-    return NULL;
-  }
-  return gerepilecopy(av, bases);
-}
-
-static GEN
-Frob_matrix(GEN J, unsigned long l, GEN basis, unsigned long m, GEN f) {
-  GEN P = gel(basis, 1), Q = gel(basis, 2), Frob;
-  pari_sp av = avma;
-
-  err_printf("computing Frobenius matrix\n");
-  if(jacobian_equal(J, Q, jacobian_Frob(J, P, m))) {
-    /*
-      If we have taken  Q = Frob(P), then the matrix of  Frob
-      with respect to (P, Q) is the companion matrix of the
-      characteristic polynomial of  Frob, which is  f.
-      We just do a consistency check.
-    */
-    GEN lincomb = jacobian_add(J,
-			       jacobian_multiply(J, P, gneg(lift(polcoeff0(f, 0, -1)))),
-			       jacobian_multiply(J, Q, gneg(lift(polcoeff0(f, 1, -1)))));
-    if(!jacobian_equal(J, jacobian_Frob(J, Q, m), lincomb))
-      pari_err(e_MISC, "inconsistent matrix of Frobenius");
-    Frob = matcompanion(f);
-  } else {
-    /* Q != Frob(P) */
-    Frob = jacobian_l_torsion_Frob_matrix(J, l, basis, m);
+    *bases = *matrices = NULL;
   }
-
-  err_printf("Frobenius matrix = %Ps\n", Frob);
-  return gerepileupto(av, Frob);
+  err_printf("Frobenius matrices: %Ps\n", *matrices);
 }
 
 static GEN
@@ -280,7 +267,8 @@ eval_function(GEN J, GEN D, GEN multiples_O, long *w) {
 }
 
 static GEN
-values_from_basis(GEN J, GEN V1, GEN V2, long l, GEN multiples_O) {
+values_from_basis(GEN J, GEN V1, GEN V2, long l, GEN matrix,
+		  GEN multiples_O) {
   pari_sp av = avma;
   long i, j, w, W = 0;
   GEN D, V;
@@ -309,9 +297,10 @@ values_from_basis(GEN J, GEN V1, GEN V2, long l, GEN multiples_O) {
   the Weil pairing on  V  (resp. between  V  and  V_dual).
 */
 static GEN
-all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
+all_function_values(GEN J, unsigned long l, GEN f, GEN proj, int tries) {
   pari_sp av = avma;
-  GEN bases, basis, basis_dual, values, values_dual;
+  GEN bases, matrices, basis, basis_dual;
+  GEN matrix, matrix_dual, values, values_dual;
   GEN O, multiples_O, V1, V2, V1_dual, V2_dual;
   GEN P, Q, P_dual, Q_dual, z, Z;
 
@@ -320,7 +309,7 @@ all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
   multiples_O = curve_point_multiples(J, O);
 
   err_printf("computing bases for the desired subspaces of the %li-torsion\n", l);
-  bases = find_bases(J, l, proj, tries);
+  find_bases(J, l, f, proj, tries, &bases, &matrices);
   if (bases == NULL) {
     basis = jacobian_l_torsion_basis(J, l, multiples_O);
     /* TODO: project */
@@ -331,7 +320,8 @@ all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
   V2 = gel(basis, 2);
   P = gel(V1, 2);
   Q = gel(V2, 2);
-  values = values_from_basis(J, V1, V2, l, multiples_O);
+  matrix = gel(matrices, 1);
+  values = values_from_basis(J, V1, V2, l, matrix, multiples_O);
 
   if (lg(proj) == 2) {
     err_printf("computing Weil pairing\n");
@@ -346,7 +336,9 @@ all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
     V2_dual = gel(basis_dual, 2);
     P_dual = gel(V1_dual, 2);
     Q_dual = gel(V2_dual, 2);
-    values_dual = values_from_basis(J, V1_dual, V2_dual, l, multiples_O);
+    matrix_dual = gel(matrices, 2);
+    values_dual = values_from_basis(J, V1_dual, V2_dual, l, matrix_dual,
+				    multiples_O);
     err_printf("computing Weil pairings\n");
     Z = mkmat2(mkcol2(jacobian_weil_pairing(J, P, P_dual, l),
 		      jacobian_weil_pairing(J, Q, P_dual, l)),
@@ -367,13 +359,13 @@ all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
 GEN
 torsion_subscheme(GEN Gamma, unsigned long l, unsigned long p,
 		  GEN f, unsigned long max_degree, unsigned long tries) {
-  GEN projector, projector_dual, proj;
+  GEN f_dual, projector, projector_dual, proj;
   GEN extension, J, J_k, order_J, values;
   long v;
   pari_sp av = avma;
 
   /* Find the field of definition of the representation space.  */
-  extension = find_extension(Gamma, l, p, f, &projector,
+  extension = find_extension(Gamma, l, p, f, &f_dual, &projector,
 			     &projector_dual, max_degree);
 
   /*
@@ -391,10 +383,13 @@ torsion_subscheme(GEN Gamma, unsigned long l, unsigned long p,
   pari_printf("Jacobian has order %Pi = %li^%li * %Pi\n",
 	      order_J, l, v, gdiv(order_J, powuu(l, v)));
 
-  if (RgX_equal(projector, projector_dual))
+  if (RgX_equal(f, f_dual)) {
+    f = mkvec(f);
     proj = mkvec(projector);
-  else
+  } else {
+    f = mkvec2(f, f_dual);
     proj = mkvec2(projector, projector_dual);
-  values = all_function_values(J_k, l, proj, tries);
+  }
+  values = all_function_values(J_k, l, f, proj, tries);
   return gerepileupto(av, values);
 }