avoid duplicate work in computing multiples of basis points

modcurve/torsion-subscheme.c

@@ -160,17 +160,12 @@ multiples(GEN J, GEN P, long l) {
   return mult;
 }
 
-/* assume P, Q != 0 */
 static int
-is_multiple(GEN J, GEN P, GEN Q, long l) {
+point_index(GEN J, GEN V, GEN P) {
   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;
-    }
+  for (i = 1; i < lg(V); i++) {
+    if (jacobian_equal(J, gel(V, i), P))
+      return i;
   }
   return 0;
 }
@@ -181,8 +176,8 @@ is_multiple(GEN J, GEN P, GEN Q, long l) {
 */
 static GEN
 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;
+  GEN bases, P, Q, R, V1, V2;
+  long i, j, k, n = lg(proj) - 1, found = 0;
   pari_sp av = avma;
 
   bases = cgetg(n + 1, t_VEC);
@@ -199,18 +194,21 @@ find_bases(GEN J, unsigned long l, GEN proj, int tries) {
       if (jacobian_is_zero(J, Q))
 	continue;
       if (j == 0) {
+	V1 = multiples(J, Q, l);
 	R = jacobian_Frob(J, Q, 1);
-	if (is_multiple(J, Q, R, l)) {
-	  gel(bases, i) = mkvec(Q);
+	if ((k = point_index(J, V1, R)) != 0) {
+	  gel(bases, i) = mkvec(V1);
 	  found++;
 	} else {
-	  gel(bases, i) = mkvec2(Q, R);
+	  V2 = multiples(J, R, l);
+	  gel(bases, i) = mkvec2(V1, V2);
 	  found += 2;
 	}
       } else {
-	R = gmael(bases, i, 1);
-	if (!is_multiple(J, R, Q, l)) {
-	  gel(bases, i) = mkvec2(R, Q);
+	V1 = gmael(bases, i, 1);
+	if ((k = point_index(J, V1, Q)) == 0) {
+	  V2 = multiples(J, Q, l);
+	  gel(bases, i) = mkvec2(V1, V2);
 	  found++;
 	}
       }
@@ -282,14 +280,10 @@ eval_function(GEN J, GEN D, GEN multiples_O, long *w) {
 }
 
 static GEN
-values_from_basis(GEN J, GEN P, GEN Q, long l, GEN multiples_O) {
+values_from_basis(GEN J, GEN V1, GEN V2, long l, GEN multiples_O) {
   pari_sp av = avma, av1;
   long i, j, w, W = 0;
-  GEN multiples_P, multiples_Q, D, V;
-
-  err_printf("computing multiples of P and Q\n");
-  multiples_P = multiples(J, P, l);
-  multiples_Q = multiples(J, Q, l);
+  GEN D, V;
 
   err_printf("computing function values:");
   av1 = avma;
@@ -298,7 +292,7 @@ values_from_basis(GEN J, GEN P, GEN Q, long l, GEN multiples_O) {
     for (j = 0; j < l; j++) {
       err_printf(" (%li, %li)", i, j);
       err_flush();
-      D = jacobian_addflip(J, gel(multiples_P, i + 1), gel(multiples_Q, j + 1));
+      D = jacobian_addflip(J, gel(V1, i + 1), gel(V2, j + 1));
       gcoeff(V, i + 1, j + 1) = eval_function(J, D, multiples_O, &w);
       W += w;
       if (gc_needed(av1, 1))
@@ -319,7 +313,8 @@ static GEN
 all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
   pari_sp av = avma;
   GEN bases, basis, basis_dual, values, values_dual;
-  GEN O, multiples_O, P, Q, P_dual, Q_dual, z, Z;
+  GEN O, multiples_O, V1, V2, V1_dual, V2_dual;
+  GEN P, Q, P_dual, Q_dual, z, Z;
 
   err_printf("computing multiples of the distinguished point\n");
   O = modular_curve_distinguished_point(J);
@@ -333,9 +328,11 @@ all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
   }
 
   basis = gel(bases, 1);
-  P = gel(basis, 1);
-  Q = gel(basis, 2);
-  values = values_from_basis(J, P, Q, l, multiples_O);
+  V1 = gel(basis, 1);
+  V2 = gel(basis, 2);
+  P = gel(V1, 2);
+  Q = gel(V2, 2);
+  values = values_from_basis(J, V1, V2, l, multiples_O);
 
   if (lg(proj) == 2) {
     err_printf("computing Weil pairing\n");
@@ -346,9 +343,11 @@ all_function_values(GEN J, unsigned long l, GEN proj, int tries) {
   }
   else {
     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);
+    V1_dual = gel(basis_dual, 1);
+    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);
     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)),