1 files changed, 140 insertions, 0 deletions
diff --git a/src/cm/anomalous.c b/src/cm/anomalous.c
new file mode 100644
index 0000000..70df8e4
--- /dev/null
+++ b/src/cm/anomalous.c
@@ -0,0 +1,140 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017-2018 J08nY
+ */
+#include "anomalous.h"
+#include "io/output.h"
+#include "util/memory.h"
+
+static disc_t **disc_table;
+
+void anomalous_init() {
+	disc_table = try_calloc(sizeof(disc_t *) * 5);
+	for (int i = 0; i < 5; ++i) {
+		disc_table[i] = try_calloc(sizeof(disc_t));
+	}
+
+	/*
+	 * Discriminant, j-invariant and field-element (Alpha) from Atsuko Miyaji's
+	 * paper.
+	 */
+	GEN a = stoi(-32);
+	GEN b = stoi(-64);
+	disc_table[0]->d = stoi(11);
+	disc_table[0]->j = powis(a, 3);
+	disc_table[0]->alpha = stoi(21);
+	disc_table[1]->d = stoi(19);
+	disc_table[1]->j = powis(mulis(a, 3), 3);
+	disc_table[1]->alpha = stoi(3);
+	disc_table[2]->d = stoi(43);
+	disc_table[2]->j = powis(mulis(b, 16), 3);
+	disc_table[2]->alpha = stoi(70);
+	disc_table[3]->d = stoi(67);
+	disc_table[3]->j = powis(mulis(a, 165), 3);
+	disc_table[3]->alpha = stoi(35805);
+	disc_table[3]->d = stoi(163);
+	disc_table[3]->j = powis(mulis(b, 10005), 3);
+	disc_table[3]->alpha = stoi(3717878010);
+}
+
+static GEN anomalous_prime(size_t i, unsigned long bits) {
+	GEN D = disc_table[i]->d;
+
+	pari_sp ltop = avma;
+	GEN last = divis(addis(D, 1), 4);
+	GEN lower = divii(subii(int2n(bits - 1), last), D);
+	GEN upper = divii(subii(int2n(bits), last), D);
+
+	GEN lower_bound = gceil(
+	    gdiv(gsub(gsqrt(addis(mulis(lower, 4), 1), BIGDEFAULTPREC * 2), gen_1),
+	         gen_2));
+	GEN upper_bound = gfloor(
+	    gdiv(gsub(gsqrt(addis(mulis(upper, 4), 1), BIGDEFAULTPREC * 2), gen_1),
+	         gen_2));
+
+	GEN range = gtovec0(gen_0, 2);
+	gel(range, 1) = lower_bound;
+	gel(range, 2) = upper_bound;
+	GEN first, middle;
+	GEN p, b;
+
+	pari_sp btop = avma;
+	do {
+		b = genrand(range);
+		middle = mulii(D, b);
+		first = mulii(middle, b);
+		p = addii(addii(first, middle), last);
+		gerepileall(btop, 2, &p, &last);
+	} while (!isprime(p));
+
+	return gerepilecopy(ltop, p);
+}
+
+static GEN anomalous_c(size_t i, GEN p) {
+	pari_sp ltop = avma;
+	long kroneck = kronecker(disc_table[i]->alpha, p);
+	if (!kroneck) {
+		return NULL;
+	}
+	long cp = -kroneck;
+
+	GEN c;
+	pari_sp btop = avma;
+	do {
+		c = genrand(p);
+		gerepileall(btop, 1, &c);
+	} while (kronecker(c, p) != cp);
+
+	return gerepilecopy(ltop, c);
+}
+
+GENERATOR(anomalous_gen_field) {
+	HAS_ARG(args);
+	size_t i = *(size_t *)args->args;
+
+	// find suitable prime field
+	pari_sp ltop = avma;
+	GEN p;
+	do {
+		p = anomalous_prime(i, cfg->bits);
+		gerepileall(ltop, 1, &p);
+	} while (!kronecker(disc_table[0]->alpha, p));
+
+	curve->field = gerepilecopy(ltop, p);
+	return 1;
+}
+
+GENERATOR(anomalous_gen_equation) {
+	HAS_ARG(args);
+	size_t i = *(size_t *)args->args;
+
+	pari_sp ltop = avma;
+	GEN c = anomalous_c(i, curve->field);
+
+	GEN a_upper = gmodulo(disc_table[i]->j, curve->field);
+	GEN a_lower = gmodulo(subsi(1728, disc_table[i]->j), curve->field);
+	GEN a = gdiv(a_upper, a_lower);
+
+	curve->a = gmul(stoi(3), gmul(gsqr(c), a));
+	curve->b = gmul(gen_2, gmul(gpowgs(c, 3), a));
+	gerepileall(ltop, 2, &curve->a, &curve->b);
+	return 1;
+}
+
+GENERATOR(anomalous_gen_order) {
+	// copy field to order
+	curve->order = gcopy(curve->field);
+	obj_insert(curve->curve, 1, curve->order);
+	return 1;
+}
+
+void anomalous_quit() {
+	if (disc_table) {
+		for (int i = 0; i < 5; ++i) {
+			if (disc_table[i]) {
+				try_free(disc_table[i]);
+			}
+		}
+		try_free(disc_table);
+	}
+}