Fix curve vs twist selection in cm.

Fixes #37.

5 files changed, 81 insertions, 66 deletions

diff --git a/src/cm/cm.c b/src/cm/cm.c
index 4690ed8..649df4d 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -202,6 +202,10 @@ static int cm_init(exhaustive_t *setup) {
 			fprintf(err, "Order requested not a number: %s\n", cfg->cm_order);
 			return 1;
 		}
+		if (cmpiu(order, 5) < 0) {
+			pari_fprintf(err, "Order requested too small: %Pi\n", order);
+			return 1;
+		}
 		long ord_log = logint0(order, gen_2, NULL);
 		if (ord_log > cfg->bits) {
 			pari_fprintf(err,
diff --git a/src/cm/cm_any.c b/src/cm/cm_any.c
index 5e4bab2..0e5ed7c 100644
--- a/src/cm/cm_any.c
+++ b/src/cm/cm_any.c
@@ -182,66 +182,24 @@ GEN cm_construct_curve(GEN order, GEN d, GEN p, cm_any_roots_t *roots,
 	for (long i = roots->used; i < roots->total; ++i) {
 		roots->used = i + 1;
 		GEN root = gel(roots->roots, i + 1);
-		if (gequal(root, gen_0)) {
-			debug_log("skipping root = 0");
+		if (gequal(root, gen_0) || equaliu(root, 1728)) {
+			debug_log("skipping root = 0 or 1728");
 			continue;
 		}
 		debug_log("trying root[%i] = %Pi", i + 1, root);
 
-		GEN e = ellinit(ellfromj(mkintmod(root, p)), p, 0);
+		GEN e1 = ellinit(ellfromj(mkintmod(root, p)), p, 0);
 		pari_CATCH(e_TYPE) { continue; }
-		pari_TRY { checkell(e); };
+		pari_TRY { checkell(e1); };
 		pari_ENDCATCH{};
 		debug_log("ellinit done");
 
-		if (ord_prime) {
-			// Easy part, the requested order is prime so
-			// [order]G = 0 iff the curve has exactly order points, for any G on
-			// it. otherwise it is the twist
-			GEN g = genrand(e);
-			if (ell_is_inf(ellmul(e, g, order))) {
-				debug_log("Got curve.");
-				return gerepilecopy(ltop, e);
-			} else {
-				debug_log("Got curve twist.");
-				return gerepilecopy(ltop, ellinit(elltwist(e, NULL), p, 0));
-			}
+		GEN ord = ellff_get_card(e1);
+		if (equalii(ord, order)) {
+			return gerepilecopy(ltop, e1);
 		} else {
-			// Hard part, requested order is composite, so it might share a
-			// factor with the order of the twist, which means [order]G = 0
-			// might be true for a point on the twist as well as a point o the
-			// right curve.
-			//
-			// We calculate what the twist order is, then compute gcd of the
-			// orders which leads to the product of the factors that the orders
-			// do not share. By multiplying a random point by this product on
-			// some curve, we can determine that that curve has that number of
-			// points.
-			GEN twist_order = subii(addis(p, 1), subii(order, addis(p, 1)));
-			GEN twist = ellinit(elltwist(e, NULL), p, 0);
-			GEN gcd = gcdii(order, twist_order);
-			GEN orig_mul = divii(order, gcd);
-			GEN twist_mul = divii(twist_order, gcd);
-			while (true) {
-				GEN orig_point = genrand(e);
-				if (ell_is_inf(ellmul(e, orig_point, orig_mul))) {
-					debug_log("Got curve.");
-					return gerepilecopy(ltop, e);
-				}
-				if (ell_is_inf(ellmul(e, orig_point, twist_mul))) {
-					debug_log("Got curve twist.");
-					return gerepilecopy(ltop, twist);
-				}
-				GEN twist_point = genrand(twist);
-				if (ell_is_inf(ellmul(e, twist_point, twist_mul))) {
-					debug_log("Got curve.");
-					return gerepilecopy(ltop, e);
-				}
-				if (ell_is_inf(ellmul(e, twist_point, orig_mul))) {
-					debug_log("Got curve twist.");
-					return gerepilecopy(ltop, twist);
-				}
-			}
+			GEN e2 = ellinit(elltwist(e1, NULL), p, 0);
+			return gerepilecopy(ltop, e2);
 		}
 	}
 	avma = ltop;
@@ -360,7 +318,8 @@ GENERATOR(cm_gen_curve_unique) {
 	GEN e;
 	if (min_d && min_curve == curve) {
 		debug_log("Reusing min D = %Pi", min_d->d);
-		// We have some discriminant, cannot use the roots because the D with walkdown was too large.
+		// We have some discriminant, cannot use the roots because the D with
+		// walkdown was too large.
 		if (min_d->d && isclone(min_d->d)) {
 			gunclone(min_d->d);
 		}
@@ -428,10 +387,12 @@ GENERATOR(cm_gen_curve_unique) {
 		}
 		long dsize = logint(d, gen_2);
 		if (dsize > 30) {
-            debug_log("Discriminant too large after walkdown (d = %Pi, %li bits).", d, dsize);
-            avma = ltop;
-            return -3;
-        }
+			debug_log(
+			    "Discriminant too large after walkdown (d = %Pi, %li bits).", d,
+			    dsize);
+			avma = ltop;
+			return -3;
+		}
 		if (min_roots) {
 			cm_update_roots(d, min_d->p, min_roots);
 		} else {
diff --git a/src/cm/cm_prime.c b/src/cm/cm_prime.c
index fa4d6f5..bdbf7c1 100644
--- a/src/cm/cm_prime.c
+++ b/src/cm/cm_prime.c
@@ -11,6 +11,11 @@
 #include "obj/subgroup.h"
 #include "util/bits.h"
 
+/*
+ * This file implements Algorithm 2.2 from
+ * "Constructing elliptic curves of prime order"
+ * https://pub.math.leidenuniv.nl/~stevenhagenp/bs.pdf
+ */
 static size_t add_primes(GEN r, GEN order, GEN **primes, size_t nprimes) {
 	debug_log("add_primes r = %Pi, nprimes = %lu", r, nprimes);
 	size_t nalloc = nprimes;
@@ -23,12 +28,16 @@ static size_t add_primes(GEN r, GEN order, GEN **primes, size_t nprimes) {
 
 	GEN rlog = mulii(r, logN);
 	GEN rplog = mulii(addis(r, 1), logN);
+	// debug_log("rlog = %Pi", rlog);
+	// debug_log("rplog = %Pi", rplog);
 
 	forprime_t iter;
 	forprime_init(&iter, rlog, rplog);
 	GEN prime;
 	while ((prime = forprime_next(&iter))) {
+		// debug_log("prime = %Pi", prime);
 		long k = kronecker(order, prime);
+		// debug_log("k = %li", k);
 		if (k == 1) {
 			GEN pstar = prime;
 			GEN ppow = divis(subis(prime, 1), 2);
@@ -91,6 +100,7 @@ static void qdisc_next(cm_prime_qdisc_t *qdisc) {
 				}
 			}
 			bits_free(&ibits);
+			// debug_log("pprod = %Pi", pprod);
 
 			GEN absp = absi(pprod);
 			long m4 = mod4(absp);
diff --git a/test/common.sh b/test/common.sh
index 42978a3..dfbbe94 100644
--- a/test/common.sh
+++ b/test/common.sh
@@ -7,7 +7,6 @@
 ####
 
 ecgen="../ecgen"
-econvert="../econvert"
 ASSERT="lib/assert.sh/assert.sh"
 JSON="lib/JSON.sh/JSON.sh"
 
@@ -41,3 +40,10 @@ canonical_num() {
     num=$(strip_num "$1")
     echo "ibase=16;${num^^}" | bc
 }
+
+get_pari_order() {
+  p=$(canonical_num $(echo $1 | ${JSON} -x field\",\"p | cut -f 2))
+  a=$(canonical_num $(echo $1 | ${JSON} -x \"a | cut -f 2))
+  b=$(canonical_num $(echo $1 | ${JSON} -x \"b | cut -f 2))
+  echo "ellcard(ellinit([${a}, ${b}], ${p}))" | gp -q 2>/dev/null
+}
\ No newline at end of file
diff --git a/test/ecgen.sh b/test/ecgen.sh
index d69f7f7..49a0b1e 100755
--- a/test/ecgen.sh
+++ b/test/ecgen.sh
@@ -101,20 +101,24 @@ function nums() {
 function anomalous() {
     start_test
     assert_raises "${ecgen} --fp --anomalous -r 20"
-    out=$(${ecgen} --fp --anomalous -r 20 2>/dev/null)
-    p=$(echo $out | ${JSON} -x field\",\"p | cut -f 2)
-    order=$(echo $out | ${JSON} -x ^0,\"order\" | cut -f 2)
-    assert "strip_num $p" $(strip_num $order)
+    for i in $(seq 10); do
+      out=$(${ecgen} --fp --anomalous -r 20 2>/dev/null)
+      p=$(echo $out | ${JSON} -x field\",\"p | cut -f 2)
+      order=$(echo $out | ${JSON} -x ^0,\"order\" | cut -f 2)
+      assert "strip_num $p" $(strip_num $order)
+    done
 }
 
 function supersingular() {
     start_test
     assert_raises "${ecgen} --fp --supersingular -r -c 5 20"
-    out=$(${ecgen} --fp --supersingular -r 20 2>/dev/null)
-	  p=$(echo $out | ${JSON} -x field\",\"p | cut -f 2)
-	  order=$(echo $out | ${JSON} -x ^0,\"order\" | cut -f 2)
-	  order_m1=$(echo $(canonical_num $order) - 1 | bc)
-    assert "canonical_num $p" $order_m1
+    for i in $(seq 10); do
+      out=$(${ecgen} --fp --supersingular -r 20 2>/dev/null)
+      p=$(echo $out | ${JSON} -x field\",\"p | cut -f 2)
+      order=$(echo $out | ${JSON} -x ^0,\"order\" | cut -f 2)
+      order_m1=$(echo $(canonical_num $order) - 1 | bc)
+      assert "canonical_num $p" $order_m1
+    done
 
     assert_raises "${ecgen} --fp --supersingular --input=data/prime.in 64"
 
@@ -194,6 +198,35 @@ function cm() {
     assert_raises "${ecgen} --fp --order=0x1000 8" 1
 }
 
+function cm_orders() {
+  start_test
+  for i in $(seq 5 100); do
+    out=$(timeout -k 4 3 ${ecgen} --fp -n $i --points=none 6 2>/dev/null)
+    if [[ -z "$out" || "$out" = "[" ]]; then
+      continue
+    fi
+    order=$(echo $out | ${JSON} -x ^0,\"order\" | cut -f 2)
+    pari_order=$(get_pari_order "$out")
+    assert "canonical_num $order" $pari_order
+  done
+
+  prime_orders=(45678945611413 47889465415131 78246132456157 3879641663983 134537095890397 3790687732807)
+  for ord in "${prime_orders[@]}"; do
+    out=$(${ecgen} --fp -r --order=$ord 64 2>/dev/null)
+    p=$(echo $out | ${JSON} -x field\",\"p | cut -f 2)
+    order=$(echo $out | ${JSON} -x ^0,\"order\" | cut -f 2)
+    assert "canonical_num $order" $ord
+  done
+
+  composite_orders=(106618070007935 32268705670290 78286235471710 93953327960423 17042092126557 43615536370894)
+  for ord in "${composite_orders[@]}"; do
+    out=$(${ecgen} --fp -r --order=$ord 64 2>/dev/null)
+    p=$(echo $out | ${JSON} -x field\",\"p | cut -f 2)
+    order=$(echo $out | ${JSON} -x ^0,\"order\" | cut -f 2)
+    assert "canonical_num $order" $ord
+  done
+}
+
 function secg() {
     function test_order() {
         name="${1}"
@@ -236,5 +269,6 @@ twist
 cli
 hex
 cm
+cm_orders
 secg
 end_suite ecgen