find a basis for the dual subspace at the same time

modcurve/torsion-subscheme.c

@@ -160,64 +160,68 @@ multiples(GEN J, GEN P, long l) {
   return mult;
 }
 
-static GEN
-random_representation_space_point(GEN J, unsigned long l, GEN projector, int tries) {
-  GEN P;
-  pari_sp av = avma;
-  while (tries-- > 0) {
-    avma = av;
-    P = random_torsion_point(J, l);
-    err_printf("applying Frobenius polynomial %Ps\n", projector);
-    P = jacobian_Frob_polynomial(J, P, centerlift(projector), 1);
-    if(!jacobian_is_zero(J, P))
-      return gerepileupto(av, P);
+/* assume P, Q != 0 */
+static int
+is_multiple(GEN J, GEN P, GEN Q, long l) {
+  long i;
+  GEN iP;
+  for (i = 1; i < l; i++) {
+    iP = (i == 1) ? P : jacobian_add(J, iP, P);
+    if (jacobian_equal(J, Q, iP)) {
+      err_printf("Q = %li*P\n", i - 1);
+      return 1;
+    }
   }
-  if (GIVEUP)
-    pari_err(e_MISC, "too many tries to find a random point");
-  return NULL;
+  return 0;
 }
 
-/* Naïve algorithm for finding a basis of the  l-torsion.  */
+/*
+  Naïve algorithm for finding bases of the subspaces of the
+  l-torsion defined by the elements of  proj.
+*/
 static GEN
-find_basis(GEN J, unsigned long l, GEN projector, int tries) {
-  GEN P, Q, multiples_P;
-  long i;
+find_bases(GEN J, unsigned long l, GEN proj, int tries) {
+  GEN bases, P, Q, R;
+  long i, j, n = lg(proj) - 1, found = 0;
   pari_sp av = avma;
-  int tries0 = tries;
-
-  /* Find a point  P  in the representation space.  */
-  P = random_representation_space_point(J, l, projector, tries0);
-  if (P == NULL)
-    return NULL;
 
-  /* Create the vector [0, P, 2*P, ..., (l - 1)*P].  */
-  multiples_P = cgetg(l + 1, t_VEC);
-  gel(multiples_P, 1) = jacobian_zero(J);
-  for(i = 1; i <= l - 1; i++)
-    gel(multiples_P, i + 1) = jacobian_add(J, gel(multiples_P, i), P);
-
-  err_printf("looking for a linearly independent point\n");
-  Q = jacobian_Frob(J, P, 1);
-
-  while(tries-- > 0) {
-    int Q_multiple_of_P = 0;
-    for(i = 1; i <= l; i++) {
-      if(jacobian_equal(J, Q, gel(multiples_P, i))) {
-	err_printf("Q = %li*P\n", i - 1);
-	Q_multiple_of_P = 1;
-	break;
+  bases = cgetg(n + 1, t_VEC);
+  for (i = 1; i <= n; i++)
+    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;
+      if (j == 2)
+	continue;
+      err_printf("applying Frobenius polynomial %Ps\n", gel(proj, i));
+      Q = jacobian_Frob_polynomial(J, P, centerlift(gel(proj, i)), 1);
+      if (jacobian_is_zero(J, Q))
+	continue;
+      found++;
+      if (j == 0) {
+	R = jacobian_Frob(J, Q, 1);
+	if (is_multiple(J, Q, R, l))
+	  gel(bases, i) = mkvec(Q);
+	else {
+	  gel(bases, i) = mkvec2(Q, R);
+	  found++;
+	}
+      } else {
+	R = gmael(bases, i, 1);
+	if (!is_multiple(J, R, Q, l))
+	  gel(bases, i) = mkvec2(R, Q);
       }
     }
-    if(!Q_multiple_of_P)
-      return gerepilecopy(av, mkvec2(P, Q));;
-    err_printf("generating new torsion point\n");
-    Q = random_representation_space_point(J, l, projector, tries0);
-    if (Q == NULL)
-      break;
+    if (gc_needed(av, 1))
+      bases = gerepileupto(av, bases);
+  }
+  if (found < 2*n) {
+    if (GIVEUP)
+      pari_err(e_MISC, "too many tries to find bases");
+    return NULL;
   }
-  if (GIVEUP)
-    pari_err(e_MISC, "too many tries to find a linearly independent point");
-  return NULL;
+  return gerepilecopy(av, bases);
 }
 
 static GEN
@@ -304,14 +308,15 @@ values_from_basis(GEN J, GEN P, GEN Q, long l, GEN multiples_O) {
 }
 
 /*
-  Return a pair [V, W], where  V  is the matrix of all
-  function values and  W  is the total weight.
+  Return either [V, W, Z] or [V, W, V_dual, W_dual, Z], where
+  V  is the matrix of all function values,  W  is the total
+  weight, similarly for  V_dual  and  W_dual, and  Z  encodes
+  the Weil pairing on  V  (resp. between  V  and  V_dual).
 */
 static GEN
-all_function_values(GEN J, unsigned long l, GEN projector,
-		    GEN projector_dual, int tries) {
+all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
   pari_sp av = avma;
-  GEN basis, basis_dual, values, values_dual;
+  GEN bases, basis, basis_dual, values, values_dual;
   GEN O, multiples_O, P, Q, P_dual, Q_dual, z, Z;
 
   err_printf("computing multiples of the distinguished point\n");
@@ -319,19 +324,18 @@ all_function_values(GEN J, unsigned long l, GEN projector,
   multiples_O = curve_point_multiples(J, O);
 
   err_printf("computing bases for the desired subspaces of the %li-torsion\n", l);
-  basis = find_basis(J, l, projector, tries);
-  if (!basis) {
+  bases = find_bases(J, l, proj, tries);
+  if (bases == NULL) {
     basis = jacobian_l_torsion_basis(J, l, multiples_O);
     /* TODO: project */
   }
-  /* TODO: find the dual subspace at the same time */
 
+  basis = gel(bases, 1);
   P = gel(basis, 1);
   Q = gel(basis, 2);
   values = values_from_basis(J, P, Q, l, multiples_O);
 
-  if (gequal(projector, projector_dual)) {
-    basis_dual = basis;
+  if (lg(proj) == 2) {
     err_printf("computing Weil pairing\n");
     z = jacobian_weil_pairing(J, P, Q, l);
     if(gequal1(z) || !gequal1(gpowgs(z, l)))
@@ -339,7 +343,7 @@ all_function_values(GEN J, unsigned long l, GEN projector,
     return gerepilecopy(av, shallowconcat(values, mkvec(z)));
   }
   else {
-    basis_dual = find_basis(J, l, projector_dual, tries);
+    basis_dual = gel(bases, 2);
     P_dual = gel(basis_dual, 1);
     Q_dual = gel(basis_dual, 2);
     values_dual = values_from_basis(J, P_dual, Q_dual, l, multiples_O);
@@ -363,7 +367,7 @@ all_function_values(GEN J, unsigned long l, GEN projector,
 GEN
 torsion_subscheme(GEN Gamma, unsigned long l, unsigned long p,
 		  GEN f, unsigned long max_degree, unsigned long tries) {
-  GEN projector, projector_dual;
+  GEN projector, projector_dual, proj;
   GEN extension, J, J_k, order_J, values;
   long v;
   pari_sp av = avma;
@@ -387,6 +391,10 @@ 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)));
 
-  values = all_function_values(J_k, l, projector, projector_dual, tries);
+  if (RgX_equal(projector, projector_dual))
+    proj = mkvec(projector);
+  else
+    proj = mkvec2(projector, projector_dual);
+  values = all_function_values(J_k, l, proj, tries);
   return gerepileupto(av, values);
 }