[Ecm-commits] r2389 - trunk

cvs commits ecm-commits at lists.gforge.inria.fr
Mon Feb 11 18:49:33 CET 2013


Author: morain
Date: 2013-02-11 18:49:33 +0100 (Mon, 11 Feb 2013)
New Revision: 2389

Log:
* fixed compiler warnings
* improved management of squareroots in CMECM. Not finished.


Modified:
   trunk/addlaws.c
   trunk/cmecm.c
   trunk/manyecm.c
Modified: trunk/addlaws.c
===================================================================
--- trunk/addlaws.c	2013-02-11 13:31:29 UTC (rev 2388)
+++ trunk/addlaws.c	2013-02-11 17:49:33 UTC (rev 2389)
@@ -903,9 +903,12 @@
 {
     /* first use automata */
     size_t le = mpz_sizeinbase(e, 2), iT = 0;
+#if DEBUG_ADD_LAWS >= 2
     long tp = cputime();
+    int i;
+#endif
     char *T = (char *)malloc((2*le) * sizeof(char)); /* humf */
-    int i, iS;
+    int iS;
 
     MO_automaton(T, &iT, e, le);
 #if DEBUG_ADD_LAWS >= 2
@@ -919,9 +922,9 @@
     else
 	printf("# good check in MO\n");
     printf("# le = %ld, iT = %ld, time = %ldms\n",le,iT,elltime(tp,cputime()));
+    tp = cputime();
 #endif
     /* compact T to fill in S */
-    tp = cputime();
     iS = Split(S, Slen, T, iT, w);
 #if DEBUG_ADD_LAWS >= 2
     printf("# time = %ldms\n", elltime(tp, cputime()));
@@ -1675,7 +1678,10 @@
 ec_point_mul_add_sub (ec_point_t Q, mpz_t e, ec_point_t P,
 		      ec_curve_t E, mpmod_t n)
 {
-    int negated = 0, status = 1, iS = 0, w, Slen, j;
+    int negated = 0, status = 1, iS = 0, w, Slen;
+#if DEBUG_ADD_LAWS >= 2
+    int j;
+#endif
     short *S;
     
     if(ec_point_is_zero(P, E, n)){

Modified: trunk/cmecm.c
===================================================================
--- trunk/cmecm.c	2013-02-11 13:31:29 UTC (rev 2388)
+++ trunk/cmecm.c	2013-02-11 17:49:33 UTC (rev 2389)
@@ -116,6 +116,17 @@
 	}
     }
     else if(disc == -4){
+	mpz_init(tmp);
+	if(sqroots != NULL){
+	    /* sqroots[0] = sqrt(-1) */
+	    mpz_mul(tmp, sqroots[0], sqroots[0]);
+	    mpz_add_si(tmp, tmp, 1);
+	    mpz_mod(tmp, tmp, n->orig_modulus);
+	    if(mpz_sgn(tmp) != 0){
+		gmp_printf("# zeta4^2+1=%Zd\n", tmp);
+		exit(-1);
+	    }
+	}
 #if 0
 	/* Y^2 = X^3 + 9 * X has rank 1 and a 4-torsion point */
 	/* a generator is (4 : 10 : 1) */
@@ -136,14 +147,6 @@
         mpz_set_si(tP[0]->z, 1);
 	if(sqroots != NULL){
 	    /* sqroots[0] = sqrt(-1) */
-	    mpz_init(tmp);
-	    mpz_mul(tmp, sqroots[0], sqroots[0]);
-	    mpz_add_si(tmp, tmp, 1);
-	    mpz_mod(tmp, tmp, n->orig_modulus);
-	    if(mpz_sgn(tmp) != 0){
-		gmp_printf("# zeta4^2+1=%Zd\n", tmp);
-		exit(-1);
-	    }
 	    for(i = 1; i < imax; i++){
 		mpz_mul(tmp, tE[i-1]->A, sqroots[0]);
 		mpz_mod(tE[i]->A, tmp, n->orig_modulus);
@@ -152,19 +155,23 @@
 		ec_force_point(tE[i], tP[i], tmp, &x0, n->orig_modulus);
 	    }
 	    *nE = 4;
-	    mpz_clear(tmp);
 	    for(i = 0; i < imax; i++)
 		mpz_init_set(tE[i]->sq[0], sqroots[0]);
 	}
-#else /* Montgomery form curves, we do not need zeta4 */
-	/* Y^2 = X^3+k^2*X -> (1/k)*y^2 = x^3+x, x=X/k, y=Y/k */
-	/* k=3, gen=(4, 10) => 1/3*y^2 = x^3 + x, gen = (4/3, 10/3) */
-	/* {k, num, den} on E.W; divide by k to obtain E.M */
+#else 
+	/* Montgomery form curves, we do not need zeta4
+	   Y^2 = X^3+k^2*X -> (1/k)*y^2 = x^3+x, x=X/k, y=Y/k
+	   k=3, gen=(4, 10) => 1/3*y^2 = x^3 + x, gen = (4/3, 10/3)
+	   {k, num, den} on E.W; divide by k to obtain E.M.
+	   Remark this imposes 4 | #E, so not all classes are possible...!
+	*/
 	long data4[][3] = {{3, 4, 1}, {7, 16, 9}, {10, 5, 1}, {11, 4900, 9},
+#if 0
 			   {14, 2, 1}, {15, 36, 1}, {17, 144, 1},
 			   {19, 722500, 77841}, {23, 7056, 121},
 			   {26, 13, 9}, {31, 1787598400, 32798529},
 			   {35, 225, 4}, {39, 2025, 4},
+#endif
 			   {0, 0, 0}};
 	for(i = 0; data4[i][0] != 0; i++){
 	    ec_curve_init(tE[i], ECM_EC_TYPE_MONTGOMERY, n);
@@ -185,6 +192,7 @@
 	printf("# using %d curves in Montgomery form for disc=-4\n", i);
 	*nE = i;
 #endif
+	mpz_clear(tmp);
     }
     else if(disc == -7){
 	/* E = y^2 = x^3 - 2222640*x - 1568294784
@@ -405,7 +413,7 @@
 
 /* Testing CM stage2 in a very naive way, for the time being. */
 int stage2_CM(mpz_t f, ec_curve_t E, ec_point_t P, mpmod_t modulus, 
-	      unsigned long dF, mpz_t B2, mpz_t B2min)
+	      unsigned long dF, mpz_t B2, ATTRIBUTE_UNUSED mpz_t B2min)
 {
     int ret = ECM_NO_FACTOR_FOUND;
     ec_point_t Q;
@@ -763,7 +771,8 @@
 
 ecm_roots_state_t *
 ecm_rootsG_init_CM (mpz_t f, curve *X, root_params_t *root_params, 
-		    unsigned long dF, unsigned long blocks, mpmod_t modulus)
+		    ATTRIBUTE_UNUSED unsigned long dF, 
+		    ATTRIBUTE_UNUSED unsigned long blocks, mpmod_t modulus)
 {
     ecm_roots_state_t *state;
     progression_params_t *params; /* for less typing */

Modified: trunk/manyecm.c
===================================================================
--- trunk/manyecm.c	2013-02-11 13:31:29 UTC (rev 2388)
+++ trunk/manyecm.c	2013-02-11 17:49:33 UTC (rev 2389)
@@ -251,7 +251,6 @@
 		if(ret == ECM_PRIME_FAC_PRIME_COFAC 
 		   || ret == ECM_PRIME_FAC_COMP_COFAC){
 		    /* output Magma lines to check #E's mod f */
-		    /*dump_curves(tE, tP, nE, f);*/
 		    printf("infos:=[\"E%d\"];\n", i);
 		    dump_curves(tE+i, tP+i, 1, f);
 		}
@@ -268,6 +267,10 @@
 	}
 
     }
+#if DEBUG_MANY_EC >= 2
+    printf("# let's debug all curves\n");
+    dump_curves(tE, tP, nE, N);
+#endif
     mpz_clear (C);
     mpcandi_t_free(&candi);
     return ret;
@@ -690,13 +693,19 @@
     return ret;
 }
 
-/* Assume b^n = 1 mod N. */
+/* Assume b^n = 1 mod N.
+   status =  0 if the squareroot could not be computed,
+            -1 if the check is bad,
+	     1 otherwise.
+ */
 int
-odd_square_root_mod_N(mpz_t f, mpz_t *sqroots, int b, int n, int q, mpz_t N)
+odd_square_root_mod_N(mpz_t f, int *status, mpz_t *sqroots,
+		      int b, int n, int q, mpz_t N)
 {
     mpz_t zeta, tmp, tmp2;
     int np = n, e = 0, *tab, x, ret = ECM_NO_FACTOR_FOUND;
 
+    *status = 1;
     while(np % q == 0){
 	e++;
 	np /= q;
@@ -706,6 +715,7 @@
     mpz_powm_ui(zeta, zeta, np, N);
     if(mpz_cmp_ui(zeta, 1) == 0){
 	printf("# missed: zeta == 1\n");
+	*status = 0;
     }
     else{
 	/* look for k s.t. zeta^{q^k} = 1 */
@@ -754,7 +764,7 @@
 	else{
 	    gmp_printf("Bad check: %Zd\n", tmp2);
 	    gmp_printf("N:=%Zd;\n", N);
-	    ret = ECM_ERROR;
+	    *status = -1;
 	}
 	mpz_clear(tmp);
 	mpz_clear(tmp2);
@@ -766,10 +776,13 @@
 }
 
 /* b^(2*k) = -1 => (b^k)^2 = -1 mod N */
-static
-int psb_minus_even(mpz_t sqroots[], int b, int k, mpz_t N)
+static void
+psb_minus_even(int *tsq, mpz_t sqroots[], int b, int k, mpz_t N)
 {
+    int isq = 0;
+
     printf("# got sqrt(-1)\n");
+    tsq[isq++] = -1;
     mpz_init_set_si(sqroots[0], b);
     mpz_powm_ui(sqroots[0], sqroots[0], k, N);
     if(k % 2 == 0){
@@ -783,46 +796,43 @@
 	mpz_mul(sqroots[1], sqroots[1], sqroots[0]);
 	mpz_mod(sqroots[1], sqroots[1], N);
 	mpz_sub_si(sqroots[0], sqroots[0], 1);
+	tsq[isq++] = 2;
     }
-    return -4;
+    tsq[isq] = 0;
 }
 
 /* b^(2*k+1) = -1 => (b^(k+1))^2 = -b mod N */
-static
-int psb_minus_odd(mpz_t sqroots[], int b, int k, mpz_t N)
+static void
+psb_minus_odd(int *tsq, mpz_t sqroots[], int b, int k, mpz_t N)
 {
-    int disc1 = 0;
-
-    if(b == 2)
-	disc1 = -8;
-    else if(b == 3 || b == 7 || b == 11)
-	disc1 = -b;
-    if(disc1 != 0){
-	printf("# got sqrt(-%d)\n", b);
-	mpz_init_set_si(sqroots[0], b);
-	mpz_powm_ui(sqroots[0], sqroots[0], k+1, N);
-    }
-    return disc1;
+    printf("# got sqrt(-%d)\n", b);
+    tsq[0] = -b;
+    tsq[1] = 0;
+    mpz_init_set_si(sqroots[0], b);
+    mpz_powm_ui(sqroots[0], sqroots[0], k+1, N);
 }
 
 /* b^(2*k+1) = 1 mod N => (b^(k+1))^2 = b mod N */
-static int
-psb_plus_odd(mpz_t sqroots[], int b, int k, mpz_t N)
+static void
+psb_plus_odd(int *tsq, mpz_t sqroots[], int b, int k, mpz_t N)
 {
     printf("# got sqrt(%d)\n", b);
     mpz_init_set_si(sqroots[0], b);
     mpz_powm_ui(sqroots[0], sqroots[0], k+1, N);
-    return b;
+    tsq[0] = b;
+    tsq[1] = 0;
 }
 
 /* OUTPUT: ECM_NO_FACTOR_FOUND or ECM_FACTOR_FOUND_STEP1 in very rare cases! */
 static int
-prepare_squareroots_from_powers(int *disc1, mpz_t f, mpz_t sqroots[], int b,
+prepare_squareroots_from_powers(mpz_t f, int *tsq, 
+				mpz_t sqroots[], int b,
 				int n, int sgn, mpz_t N)
 {
-    int k, discr = 0, ret = ECM_NO_FACTOR_FOUND;
+    int k, ret = ECM_NO_FACTOR_FOUND;
     mpz_t tmp, tmp2;
 
+    tsq[0] = 0;
     if(sgn == -1){
 	/* b^n = 1 mod N */
 	if(n % 2 == 0){
@@ -843,7 +853,7 @@
 	    gmp_printf("# %d^%d = 1 mod %Zd;\n", b, k, N);
 	    if(k % 2 == 1)
 		/* b^(2*r+1) = 1 mod N => (b^(r+1))^2 = b mod N */
-		discr = psb_plus_odd(sqroots, b, n>>1, N);
+		psb_plus_odd(tsq, sqroots, b, n>>1, N);
 	    else{
 		/* b^(2*r) = 1 */
 		mpz_add_si(tmp2, tmp2, 1);
@@ -852,10 +862,10 @@
 		    printf("# %d^%d = -1 mod N;\n", b, k>>1);
 		    if(k % 4 == 0)
 			/* b^(2*s) = -1 */
-			discr = psb_minus_even(sqroots, b, k>>2, N);
+			psb_minus_even(tsq, sqroots, b, k>>2, N);
 		    else
 			/* b^(2*s+1) = -1 */
-			discr = psb_minus_odd(sqroots, b, k>>2, N);
+			psb_minus_odd(tsq, sqroots, b, k>>2, N);
 		}
 		else{
 		    /* we have a factor, since tmp2^2 = 1, tmp2 != -1, +1 */
@@ -870,42 +880,51 @@
 	}
 	else
 	    /* b^(2*k+1) = 1 mod N => (b^(k+1))^2 = b mod N */
-	    discr = psb_plus_odd(sqroots, b, n>>1, N);
+	    psb_plus_odd(tsq, sqroots, b, n>>1, N);
     }
     else{
 	/* b^n = -1 mod N */
 	if(n % 2 == 0)
 	    /* b^(2*k) = -1 mod N => (b^k)^2 = -1 mod N */
-	    discr = psb_minus_even(sqroots, b, n>>1, N);
+	    psb_minus_even(tsq, sqroots, b, n>>1, N);
 	else
 	    /* b^(2*k+1) = -1 mod N => (b^(k+1))^2 = -b mod N */
-	    discr = psb_minus_odd(sqroots, b, n>>1, N);
+	    psb_minus_odd(tsq, sqroots, b, n>>1, N);
     }
-    *disc1 = discr;
     return ret;
 }
 
+/* N is a cofactor of b^n+sgn. */
 static int
-prepare_squareroots(int *disc1, mpz_t f, mpz_t sqroots[], int b,
-		    int n, int sgn, mpz_t N)
+prepare_squareroots(mpz_t f, int *tsq, mpz_t sqroots[], 
+		    int b, int n, int sgn, mpz_t N)
 {
-    int ret = prepare_squareroots_from_powers(disc1, f, sqroots, b, n, sgn, N);
-    int tabq[] = {3, 5, 7, 11, 13, 19, 0}, q, iq, qs;
+    int ret, tabq[] = {3, 5, 7, 11, 13, 19, 0}, q, iq, qs, isq, nn, status;
 
+    tsq[0] = 0;
+    ret = prepare_squareroots_from_powers(f,tsq,sqroots,b,n,sgn,N);
     if(ret != ECM_NO_FACTOR_FOUND)
 	return ret;
-#if 0
-    if(*disc1 != 0)
-	/* for the time being, stop if some discriminant was found */
-	return ret;
-#endif
-    /* let's find some squareroots using small powers of n */
+    printf("# I already found squareroots for:");
+    for(isq = 0; tsq[isq] != 0; isq++)
+	printf(" %d", tsq[isq]);
+    printf("\n");
+    /* let's find some squareroots using small odd prime divisors of n */
     for(iq = 0; tabq[iq] != 0; iq++){
 	q = tabq[iq];
 	qs = (q % 4 == 1 ? q : -q);
-	if(n % q == 0)
+	if(n % q == 0){
 	    printf("# I can find sqrt(%d)\n", qs);
+	    /* make sure that b^nn = 1 */
+	    nn = (sgn == -1 ? n : 2*n);
+	    ret = odd_square_root_mod_N(f, &status, sqroots+isq, b, nn, q, N);
+	    if(ret != ECM_NO_FACTOR_FOUND)
+		break;
+	    if(status == 1)
+		tsq[isq++] = qs;
+	}
     }
+    tsq[isq] = 0;
     return ret;
 }
 
@@ -920,10 +939,10 @@
 		      double B1, mpz_t B2, 
 		      ecm_params params, char *savefilename)
 {
-    int sgn = 1, disc1 = 0, q, nn, disc, i, j, k;
+    int sgn = 1, disc1 = 0, i;
     int ret = ECM_NO_FACTOR_FOUND;
+#if 0
     int tabd[][4] = {{-15, -3, 5, 0}, 
-#if 0
 		     {-20, -4, 5, 0}, {-24, 8, -3, 0},
 		     {-35, 5, -7, 0}, {-40, -8, 5, 0}, {-51, -3, 17, 0},
 		     {-52, -4, 13, 0}, {-88, 8, -11, 0},
@@ -940,34 +959,45 @@
 		  1995, 3003, 3315,
 		     /* h = g = 16 */
 		  5460
-#endif
 		     {0, 0, 0}
     };
+#endif
     mpz_t sqroots[10];
+    int tsq[10], isq;
 
     if(n < 0){
 	sgn = -1;
 	n = -n;
     }
-    /* try discriminants of class number 1, hence rational CM curves */
-    ret = prepare_squareroots(&disc1, tf[0], sqroots, b, n, sgn, N);
+    tsq[0] = 0;
+    ret = prepare_squareroots(tf[0], tsq, sqroots, b, n, sgn, N);
     if(ret != ECM_NO_FACTOR_FOUND)
 	return ret;
-    if(disc1 == discref){
-	gmp_printf("# Let us use disc=%d with B1=%1.0f B2=%Zd\n",disc1,B1,B2);
-	*tried = 1;
-	ret = process_many_curves_loop(tf, nf, N, B1, B2, params,
-				       NULL, NULL, 0, 0, 1,
-				       disc1, sqroots,
-				       savefilename);
-	if(ret != ECM_NO_FACTOR_FOUND){
-	    if(disc1 != 0)
-		mpz_clear(sqroots[0]); /* really? */
-	    return ret;
+    /* try discriminants of class number 1, hence rational CM curves */
+    for(isq = 0; tsq[isq] != 0; isq++){
+	if(tsq[isq] == -1)
+	    disc1 = -4;
+	else if(tsq[isq] == -2)
+	    disc1 = -8;
+	else
+	    disc1 = tsq[isq];
+	if(disc1 == discref){
+	    gmp_printf("# Let us use disc=%d with B1=%1.0f B2=%Zd\n",
+		       disc1, B1, B2);
+	    *tried = 1;
+	    ret = process_many_curves_loop(tf, nf, N, B1, B2, params,
+					   NULL, NULL, 0, 0, 1,
+					   disc1, sqroots+isq,
+					   savefilename);
+	    if(ret != ECM_NO_FACTOR_FOUND)
+		break;
 	}
     }
-    if(disc1 != 0)
-	mpz_clear(sqroots[0]); /* really? */
+    if(ret != ECM_NO_FACTOR_FOUND){
+	for(i = 0; tsq[i] != 0; i++)
+	    mpz_clear(sqroots[i]);
+	return ret;
+    }
 #if 0 /* not ready yet */
 	nn = 2 * n; /* used below */
     /* each odd prime factor q of n can lead to sqrt(q*) */




More information about the Ecm-commits mailing list