From 962e197004b19e0b8e42dc977dc0ec1a85f407c1 Mon Sep 17 00:00:00 2001
From: J08nY
Date: Sat, 23 Sep 2017 22:04:37 +0200
Subject: Rename P1363 form struct.

---
 src/cm/cm.c    |  2 +-
 src/cm/p1363.c | 36 ++++++++++++++++++------------------
 src/cm/p1363.h | 12 ++++++------
 3 files changed, 25 insertions(+), 25 deletions(-)

diff --git a/src/cm/cm.c b/src/cm/cm.c
index ce171d4..f6a593f 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -12,7 +12,7 @@ int cm_do() {
 	fprintf(err, "This is *NOT IMPLEMENTED* currently.\n");
 
 	GEN D = stoi(71);
-	form_t **forms;
+	p1363_form_t **forms;
 	size_t nforms = p1363_forms(D, &forms);
 	for (size_t i = 0; i < nforms; ++i) {
 		p1363_invariant(D, forms[i]);
diff --git a/src/cm/p1363.c b/src/cm/p1363.c
index 24c65f0..3c1bf8c 100644
--- a/src/cm/p1363.c
+++ b/src/cm/p1363.c
@@ -51,13 +51,13 @@ GEN p1363_group(GEN D) {
 
 long p1363_num(GEN group) { return glength(group); }
 
-size_t p1363_forms(GEN D, form_t ***forms) {
+size_t p1363_forms(GEN D, p1363_form_t ***forms) {
 	GEN group = p1363_group(D);
 	size_t nforms = (size_t)p1363_num(group);
-	*forms = try_calloc(nforms * sizeof(form_t *));
+	*forms = try_calloc(nforms * sizeof(p1363_form_t *));
 
 	for (size_t i = 0; i < nforms; ++i) {
-		(*forms)[i] = try_calloc(sizeof(form_t));
+		(*forms)[i] = try_calloc(sizeof(p1363_form_t));
 		(*forms)[i]->A = gel(gel(group, i + 1), 1);
 		(*forms)[i]->B = gel(gel(group, i + 1), 2);
 		(*forms)[i]->C = gel(gel(group, i + 1), 3);
@@ -66,7 +66,7 @@ size_t p1363_forms(GEN D, form_t ***forms) {
 	return nforms;
 }
 
-void p1363_free(form_t ***forms, size_t nforms) {
+void p1363_free(p1363_form_t ***forms, size_t nforms) {
 	if (*forms) {
 		for (size_t i = 0; i < nforms; ++i) {
 			if ((*forms)[i]) {
@@ -117,7 +117,7 @@ GEN p1363_func_F(GEN z) {
 	return gerepilecopy(ltop, sum);
 }
 
-GEN p1363_func_fzero(GEN D, form_t *form) {
+GEN p1363_func_fzero(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 
 	GEN upper = p1363_func_F(gneg(form->theta));
@@ -129,7 +129,7 @@ GEN p1363_func_fzero(GEN D, form_t *form) {
 	return gerepilecopy(ltop, result);
 }
 
-GEN p1363_func_fone(GEN D, form_t *form) {
+GEN p1363_func_fone(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 
 	GEN upper = p1363_func_F(form->theta);
@@ -141,7 +141,7 @@ GEN p1363_func_fone(GEN D, form_t *form) {
 	return gerepilecopy(ltop, result);
 }
 
-GEN p1363_func_ftwo(GEN D, form_t *form) {
+GEN p1363_func_ftwo(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 
 	GEN upper = p1363_func_F(gpowgs(form->theta, 4));
@@ -154,9 +154,9 @@ GEN p1363_func_ftwo(GEN D, form_t *form) {
 	return gerepilecopy(ltop, result);
 }
 
-void p1363_m8(GEN D, form_t *form) { form->m8 = mod8(D); }
+void p1363_m8(GEN D, p1363_form_t *form) { form->m8 = mod8(D); }
 
-void p1363_I(GEN D, form_t *form) {
+void p1363_I(GEN D, p1363_form_t *form) {
 	switch (form->m8) {
 		case 1:
 		case 2:
@@ -179,7 +179,7 @@ void p1363_I(GEN D, form_t *form) {
 	}
 }
 
-void p1363_J(GEN D, form_t *form) {
+void p1363_J(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	GEN ac = mulii(form->A, form->C);
 
@@ -197,7 +197,7 @@ void p1363_J(GEN D, form_t *form) {
 	avma = ltop;
 }
 
-void p1363_K(GEN D, form_t *form) {
+void p1363_K(GEN D, p1363_form_t *form) {
 	switch (form->m8) {
 		case 1:
 		case 2:
@@ -216,7 +216,7 @@ void p1363_K(GEN D, form_t *form) {
 	}
 }
 
-void p1363_L(GEN D, form_t *form) {
+void p1363_L(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	GEN ac = mulii(form->A, form->C);
 	long a2 = mod2(form->A);
@@ -243,7 +243,7 @@ void p1363_L(GEN D, form_t *form) {
 	form->L = gerepileupto(ltop, form->L);
 }
 
-void p1363_M(GEN D, form_t *form) {
+void p1363_M(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	GEN quot;
 	if (mod2(form->A) == 0) {
@@ -254,7 +254,7 @@ void p1363_M(GEN D, form_t *form) {
 	form->M = gerepileupto(ltop, powii(gen_m1, quot));
 }
 
-void p1363_N(GEN D, form_t *form) {
+void p1363_N(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	long ac2 = mod2(mulii(form->A, form->C));
 
@@ -287,14 +287,14 @@ void p1363_N(GEN D, form_t *form) {
 	form->N = gerepileupto(ltop, form->N);
 }
 
-void p1363_lambda(GEN D, form_t *form) {
+void p1363_lambda(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	GEN pik = mulri(mppi(BIGDEFAULTPREC), stoi(form->K));
 	GEN quot = divrs(pik, 24);
 	form->lambda = gerepileupto(ltop, expIr(quot));
 }
 
-void p1363_theta(GEN D, form_t *form) {
+void p1363_theta(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 
 	GEN upper = gadd(gneg(gsqrt(D, BIGDEFAULTPREC)), gmul(form->B, gen_I()));
@@ -303,7 +303,7 @@ void p1363_theta(GEN D, form_t *form) {
 	form->theta = gerepileupto(ltop, gexp(quot, BIGDEFAULTPREC));
 }
 
-GEN p1363_invariant(GEN D, form_t *form) {
+GEN p1363_invariant(GEN D, p1363_form_t *form) {
 	pari_printf("[A,B,C] = %Pi %Pi %Pi\n", form->A, form->B, form->C);
 	pari_sp ltop = avma;
 
@@ -368,7 +368,7 @@ GEN p1363_invariant(GEN D, form_t *form) {
 	return gerepilecopy(ltop, result);
 }
 
-GEN p1363_poly(GEN D, form_t **forms, size_t nforms) {
+GEN p1363_poly(GEN D, p1363_form_t **forms, size_t nforms) {
 	pari_sp ltop = avma;
 	long v = fetch_var();
 	name_var(v, "t");
diff --git a/src/cm/p1363.h b/src/cm/p1363.h
index d80d70f..ef6f7a2 100644
--- a/src/cm/p1363.h
+++ b/src/cm/p1363.h
@@ -10,7 +10,7 @@
 
 #include <pari/pari.h>
 
-typedef struct form_t {
+typedef struct {
 	GEN A;
 	GEN B;
 	GEN C;
@@ -26,14 +26,14 @@ typedef struct form_t {
 
 	GEN lambda;
 	GEN theta;
-} form_t;
+} p1363_form_t;
 
-size_t p1363_forms(GEN D, form_t ***forms);
+size_t p1363_forms(GEN D, p1363_form_t ***forms);
 
-void p1363_free(form_t ***forms, size_t nforms);
+void p1363_free(p1363_form_t ***forms, size_t nforms);
 
-GEN p1363_invariant(GEN D, form_t *form);
+GEN p1363_invariant(GEN D, p1363_form_t *form);
 
-GEN p1363_poly(GEN D, form_t **forms, size_t nforms);
+GEN p1363_poly(GEN D, p1363_form_t **forms, size_t nforms);
 
 #endif  // ECGEN_CM_P1363_H
-- 
cgit v1.2.3-70-g09d2


From 99e192422d95240adbe806a53caabfa9ddb258c3 Mon Sep 17 00:00:00 2001
From: J08nY
Date: Fri, 29 Sep 2017 16:50:13 +0200
Subject: Added custom CM generation skeleton.

---
 src/cm/cm.c               |   7 ++-
 src/cm/custom.c           | 129 ++++++++++++++++++++++++++++++++++++++++++++++
 src/cm/custom.h           |  24 +++++++++
 test/src/cm/test_custom.c |  42 +++++++++++++++
 4 files changed, 200 insertions(+), 2 deletions(-)
 create mode 100644 src/cm/custom.c
 create mode 100644 src/cm/custom.h
 create mode 100644 test/src/cm/test_custom.c

diff --git a/src/cm/cm.c b/src/cm/cm.c
index f6a593f..04fbd9b 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -4,21 +4,24 @@
  */
 #include "cm.h"
 #include "io/output.h"
-#include "p1363.h"
+#include "custom.h"
 
 int cm_do() {
 	debug_log_start("Starting Complex Multiplication method");
 
 	fprintf(err, "This is *NOT IMPLEMENTED* currently.\n");
 
+	/*
 	GEN D = stoi(71);
 	p1363_form_t **forms;
 	size_t nforms = p1363_forms(D, &forms);
 	for (size_t i = 0; i < nforms; ++i) {
 		p1363_invariant(D, forms[i]);
 	}
-
 	p1363_free(&forms, nforms);
+	*/
+	custom_curve();
+
 	debug_log_start("Finished Complex Multiplication method");
 	return 0;
 }
diff --git a/src/cm/custom.c b/src/cm/custom.c
new file mode 100644
index 0000000..fc403e5
--- /dev/null
+++ b/src/cm/custom.c
@@ -0,0 +1,129 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017 J08nY
+ */
+#include "custom.h"
+#include "io/input.h"
+
+static bool custom_is_disc(GEN test) {
+	bool result = true;
+	long md4;
+	pari_CATCH(e_DOMAIN)
+			{
+				result = false;
+			}
+		pari_TRY{ check_quaddisc_imag(test, &md4, ""); };
+	pari_ENDCATCH
+	return result;
+}
+
+static GEN custom_disc(GEN order, GEN prime) {
+	pari_sp ltop = avma;
+
+	GEN t = subii(order, addii(prime, gen_1));
+	GEN v2D = subii(mulis(prime, 4), sqri(t));
+
+	GEN max_divisor = gen_1;
+	GEN min_D = v2D;
+
+	// Optional part that sometimes makes the discriminant smaller, at the cost
+	// factoring v2D.
+	GEN fctr = factor(v2D);
+	long len = glength(gel(fctr, 1));
+	for (long i = 1; i <= len; ++i) {
+		GEN value = gmael(fctr, 1, i);
+		GEN mul = gmael(fctr, 2, i);
+		GEN pow;
+		if (mod2(mul) == 0) {
+			pow = mul;
+		} else {
+			pow = subis(mul, 1);
+		}
+		if (equalis(pow, 0)) {
+			continue;
+		}
+		GEN divi = powii(value, pow);
+		if (mpcmp(divi, max_divisor) > 0) {
+			GEN candidate_D = divii(v2D, divi);
+			if (signe(candidate_D) == 1) {
+				setsigne(candidate_D, -1);
+			}
+			if (custom_is_disc(candidate_D)) {
+				max_divisor = divi;
+				min_D = candidate_D;
+			}
+		}
+	}
+
+	if (signe(min_D) == 1) {
+		setsigne(min_D, -1);
+	}
+	if (custom_is_disc(min_D)) {
+		return gerepileupto(ltop, min_D);
+	} else {
+		avma = ltop;
+		return NULL;
+	}
+}
+
+static custom_quadr_t custom_prime_input(GEN order) {
+	pari_sp ltop = avma;
+	custom_quadr_t result = {0};
+	GEN p;
+	GEN D;
+	do {
+		p = input_prime("p:", cfg->bits);
+		if (gequalm1(p)) {
+			continue;
+		}
+		D = custom_disc(order, p);
+	} while (D == NULL);
+
+	gerepileall(ltop, 2, &p, &D);
+	result.p = p;
+	result.D = D;
+	return result;
+}
+
+static custom_quadr_t custom_prime_random(GEN order) {
+	pari_sp ltop = avma;
+	custom_quadr_t result = {0};
+	GEN p;
+	GEN t = gen_0;
+	GEN v = gen_0;
+	GEN D;
+	long i = 2;
+	//TODO: this is wrong, order = p + 1 +- t  is not kept.
+	while (gequal0(t)) {
+		D = stoi(-4 * (i / 2) + (i % 2));
+		for (long j = 0; j < 100; ++j) {
+			p = random_prime(cfg->bits);
+			if (cornacchia2(negi(D), p, &t, &v)) {
+				break;
+			}
+		}
+	};
+
+	gerepileall(ltop, 4, &p, &t, &v, &D);
+	result.p = p;
+	result.t = t;
+	result.v = v;
+	result.D = D;
+	return result;
+}
+
+curve_t *custom_curve() {
+	GEN order = strtoi(cfg->cm_order);
+
+	custom_quadr_t quadr;
+	if (cfg->random) {
+		quadr = custom_prime_random(order);
+	} else {
+		quadr = custom_prime_input(order);
+	}
+	pari_printf("p = %Pi, D = %Pi, ", quadr.p, quadr.D);
+	GEN H = polclass(quadr.D, 0, 0);
+	pari_printf("H = %Ps\n", H);
+
+	return NULL;
+}
\ No newline at end of file
diff --git a/src/cm/custom.h b/src/cm/custom.h
new file mode 100644
index 0000000..ed6f3d3
--- /dev/null
+++ b/src/cm/custom.h
@@ -0,0 +1,24 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017 J08nY
+ */
+#ifndef ECGEN_CUSTOM_H
+#define ECGEN_CUSTOM_H
+
+#include "misc/types.h"
+#include "misc/config.h"
+
+typedef struct {
+	GEN p;
+	GEN t;
+	GEN v;
+	GEN D;
+} custom_quadr_t;
+
+/**
+ * @brief
+ * @return
+ */
+curve_t * custom_curve();
+
+#endif //ECGEN_CUSTOM_H
diff --git a/test/src/cm/test_custom.c b/test/src/cm/test_custom.c
new file mode 100644
index 0000000..3309de7
--- /dev/null
+++ b/test/src/cm/test_custom.c
@@ -0,0 +1,42 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017 J08nY
+ */
+
+#include <criterion/criterion.h>
+#include "misc/config.h"
+#include "util/random.h"
+#include "test/default.h"
+#include "test/input.h"
+#include "cm/custom.h"
+
+void custom_setup() {
+	default_setup();
+	input_setup();
+	random_init();
+}
+
+void custom_teardown() {
+	input_teardown();
+}
+
+TestSuite(custom, .init = custom_setup, .fini = custom_teardown);
+
+Test(custom, test_curve_one) {
+	cr_skip("Doesnt work yet.");
+	cfg->bits = 128;
+	cfg->cm_order = "263473633827487324648193013259296339349";
+	cfg->random = false;
+
+	fprintf(write_in, "191345683451069598953886857691544477637\n");
+	custom_curve();
+}
+
+Test(custom, test_curve_other) {
+	cr_skip("Doesnt work yet.");
+	cfg->bits = 128;
+	cfg->cm_order = "263473633827487324648193013259296339349";
+	cfg->random = true;
+
+	custom_curve();
+}
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2


From 0e4cfed57e2e23c7854b7a3f85efbc952deabcda Mon Sep 17 00:00:00 2001
From: J08nY
Date: Thu, 5 Apr 2018 21:31:49 +0200
Subject: Add test for p1363_forms.

---
 src/cm/cm.c              |  3 +--
 src/cm/p1363.c           | 52 ++++++++++++++++++++++++++++++++++--------------
 test/src/cm/test_p1363.c | 33 ++++++++++++++++++++++++++++++
 3 files changed, 71 insertions(+), 17 deletions(-)
 create mode 100644 test/src/cm/test_p1363.c

diff --git a/src/cm/cm.c b/src/cm/cm.c
index 04fbd9b..3d684d9 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -5,13 +5,13 @@
 #include "cm.h"
 #include "io/output.h"
 #include "custom.h"
+#include "p1363.h"
 
 int cm_do() {
 	debug_log_start("Starting Complex Multiplication method");
 
 	fprintf(err, "This is *NOT IMPLEMENTED* currently.\n");
 
-	/*
 	GEN D = stoi(71);
 	p1363_form_t **forms;
 	size_t nforms = p1363_forms(D, &forms);
@@ -19,7 +19,6 @@ int cm_do() {
 		p1363_invariant(D, forms[i]);
 	}
 	p1363_free(&forms, nforms);
-	*/
 	custom_curve();
 
 	debug_log_start("Finished Complex Multiplication method");
diff --git a/src/cm/p1363.c b/src/cm/p1363.c
index 3c1bf8c..119b0cb 100644
--- a/src/cm/p1363.c
+++ b/src/cm/p1363.c
@@ -5,7 +5,7 @@
 #include "p1363.h"
 #include "util/memory.h"
 
-GEN p1363_group(GEN D) {
+static GEN p1363_group(GEN D) {
 	pari_sp ltop = avma;
 	GEN s = mpfloor(sqrtr(rdivis(D, 3, BIGDEFAULTPREC)));
 	long llen = itos(s) * 2;
@@ -49,11 +49,11 @@ GEN p1363_group(GEN D) {
 	return gerepilecopy(ltop, vec_shorten(l, j - 1));
 }
 
-long p1363_num(GEN group) { return glength(group); }
+static size_t p1363_num(GEN group) { return (size_t)glength(group); }
 
 size_t p1363_forms(GEN D, p1363_form_t ***forms) {
 	GEN group = p1363_group(D);
-	size_t nforms = (size_t)p1363_num(group);
+	size_t nforms = p1363_num(group);
 	*forms = try_calloc(nforms * sizeof(p1363_form_t *));
 
 	for (size_t i = 0; i < nforms; ++i) {
@@ -79,7 +79,7 @@ void p1363_free(p1363_form_t ***forms, size_t nforms) {
 	}
 }
 
-GEN p1363_func_F(GEN z) {
+static GEN p1363_func_F(GEN z) {
 	pari_sp ltop = avma;
 
 	if (gequal0(z)) {
@@ -117,7 +117,7 @@ GEN p1363_func_F(GEN z) {
 	return gerepilecopy(ltop, sum);
 }
 
-GEN p1363_func_fzero(GEN D, p1363_form_t *form) {
+static GEN p1363_func_fzero(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 
 	GEN upper = p1363_func_F(gneg(form->theta));
@@ -129,7 +129,7 @@ GEN p1363_func_fzero(GEN D, p1363_form_t *form) {
 	return gerepilecopy(ltop, result);
 }
 
-GEN p1363_func_fone(GEN D, p1363_form_t *form) {
+static GEN p1363_func_fone(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 
 	GEN upper = p1363_func_F(form->theta);
@@ -141,7 +141,7 @@ GEN p1363_func_fone(GEN D, p1363_form_t *form) {
 	return gerepilecopy(ltop, result);
 }
 
-GEN p1363_func_ftwo(GEN D, p1363_form_t *form) {
+static GEN p1363_func_ftwo(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 
 	GEN upper = p1363_func_F(gpowgs(form->theta, 4));
@@ -154,9 +154,9 @@ GEN p1363_func_ftwo(GEN D, p1363_form_t *form) {
 	return gerepilecopy(ltop, result);
 }
 
-void p1363_m8(GEN D, p1363_form_t *form) { form->m8 = mod8(D); }
+static void p1363_m8(GEN D, p1363_form_t *form) { form->m8 = mod8(D); }
 
-void p1363_I(GEN D, p1363_form_t *form) {
+static void p1363_I(GEN D, p1363_form_t *form) {
 	switch (form->m8) {
 		case 1:
 		case 2:
@@ -179,7 +179,7 @@ void p1363_I(GEN D, p1363_form_t *form) {
 	}
 }
 
-void p1363_J(GEN D, p1363_form_t *form) {
+static void p1363_J(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	GEN ac = mulii(form->A, form->C);
 
@@ -197,7 +197,7 @@ void p1363_J(GEN D, p1363_form_t *form) {
 	avma = ltop;
 }
 
-void p1363_K(GEN D, p1363_form_t *form) {
+static void p1363_K(GEN D, p1363_form_t *form) {
 	switch (form->m8) {
 		case 1:
 		case 2:
@@ -243,7 +243,7 @@ void p1363_L(GEN D, p1363_form_t *form) {
 	form->L = gerepileupto(ltop, form->L);
 }
 
-void p1363_M(GEN D, p1363_form_t *form) {
+static void p1363_M(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	GEN quot;
 	if (mod2(form->A) == 0) {
@@ -254,7 +254,7 @@ void p1363_M(GEN D, p1363_form_t *form) {
 	form->M = gerepileupto(ltop, powii(gen_m1, quot));
 }
 
-void p1363_N(GEN D, p1363_form_t *form) {
+static void p1363_N(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	long ac2 = mod2(mulii(form->A, form->C));
 
@@ -287,14 +287,14 @@ void p1363_N(GEN D, p1363_form_t *form) {
 	form->N = gerepileupto(ltop, form->N);
 }
 
-void p1363_lambda(GEN D, p1363_form_t *form) {
+static void p1363_lambda(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	GEN pik = mulri(mppi(BIGDEFAULTPREC), stoi(form->K));
 	GEN quot = divrs(pik, 24);
 	form->lambda = gerepileupto(ltop, expIr(quot));
 }
 
-void p1363_theta(GEN D, p1363_form_t *form) {
+static void p1363_theta(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 
 	GEN upper = gadd(gneg(gsqrt(D, BIGDEFAULTPREC)), gmul(form->B, gen_I()));
@@ -303,6 +303,28 @@ void p1363_theta(GEN D, p1363_form_t *form) {
 	form->theta = gerepileupto(ltop, gexp(quot, BIGDEFAULTPREC));
 }
 
+/**
+ * Bit-precision computation for a Weber class polynomial from:
+ * 		On the Efficient Generation of Elliptic Curves,
+ * 		Elisavet Konstantinou, Yiannis C. Stamatiou, Christos Zaroliagis
+ * @param D
+ * @param forms
+ * @param nforms
+ * @return
+ */
+static long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) {
+	pari_sp ltop = avma;
+	long v0 = 32;
+	GEN mul = divrr(mulrr(mppi(BIGDEFAULTPREC), sqrti(D)), mplog2(BIGDEFAULTPREC));
+	GEN sum = gen_0;
+	for (size_t i = 0; i < nforms; ++i) {
+		sum = gadd(sum, ginv(forms[i]->A));
+	}
+	long result = v0 + itos(mpceil(mulrr(mul, sum)));
+	avma = ltop;
+	return result;
+}
+
 GEN p1363_invariant(GEN D, p1363_form_t *form) {
 	pari_printf("[A,B,C] = %Pi %Pi %Pi\n", form->A, form->B, form->C);
 	pari_sp ltop = avma;
diff --git a/test/src/cm/test_p1363.c b/test/src/cm/test_p1363.c
new file mode 100644
index 0000000..fdc1823
--- /dev/null
+++ b/test/src/cm/test_p1363.c
@@ -0,0 +1,33 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017 J08nY
+ */
+
+#include <criterion/criterion.h>
+#include "test/default.h"
+#include "cm/p1363.h"
+
+
+TestSuite(p1363, .init = default_setup, .fini = default_teardown);
+
+Test(p1363, test_p1363_forms) {
+	GEN D = stoi(71);
+	p1363_form_t **forms;
+	size_t nforms = p1363_forms(D, &forms);
+
+	p1363_form_t expected[7] = {
+			{gen_1,   gen_0,    stoi(71)},
+			{stoi(3), gen_1,    stoi(24)},
+			{stoi(3), gen_m1,   stoi(24)},
+			{stoi(8), gen_1,    stoi(9)},
+			{stoi(8), gen_m1,   stoi(9)},
+			{stoi(5), stoi(2),  stoi(15)},
+			{stoi(5), stoi(-2), stoi(15)}
+	};
+	cr_assert_eq(nforms, 7,);
+	for (size_t i = 0; i < nforms; ++i) {
+		cr_assert(equalii(forms[i]->A, expected[i].A), );
+		cr_assert(equalii(forms[i]->B, expected[i].B), );
+		cr_assert(equalii(forms[i]->C, expected[i].C), );
+	}
+}
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2


From d9e0c33e375905fcec15f4aea23ac61798415edb Mon Sep 17 00:00:00 2001
From: J08nY
Date: Thu, 5 Apr 2018 22:41:05 +0200
Subject: Use computed precision in p1363.

---
 src/cm/cm.c               |  8 ++---
 src/cm/custom.c           | 12 +++----
 src/cm/custom.h           |  6 ++--
 src/cm/p1363.c            | 87 ++++++++++++++++++++++++-----------------------
 src/cm/p1363.h            |  4 ++-
 test/src/cm/test_custom.c |  8 ++---
 test/src/cm/test_p1363.c  | 25 ++++++++------
 7 files changed, 75 insertions(+), 75 deletions(-)

diff --git a/src/cm/cm.c b/src/cm/cm.c
index 3d684d9..48a9c5f 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -3,8 +3,8 @@
  * Copyright (C) 2017-2018 J08nY
  */
 #include "cm.h"
-#include "io/output.h"
 #include "custom.h"
+#include "io/output.h"
 #include "p1363.h"
 
 int cm_do() {
@@ -15,11 +15,9 @@ int cm_do() {
 	GEN D = stoi(71);
 	p1363_form_t **forms;
 	size_t nforms = p1363_forms(D, &forms);
-	for (size_t i = 0; i < nforms; ++i) {
-		p1363_invariant(D, forms[i]);
-	}
+	GEN WD = p1363_poly(D, forms, nforms);
+	pari_printf("%Ps\n", WD);
 	p1363_free(&forms, nforms);
-	custom_curve();
 
 	debug_log_start("Finished Complex Multiplication method");
 	return 0;
diff --git a/src/cm/custom.c b/src/cm/custom.c
index fc403e5..9309ed9 100644
--- a/src/cm/custom.c
+++ b/src/cm/custom.c
@@ -8,13 +8,9 @@
 static bool custom_is_disc(GEN test) {
 	bool result = true;
 	long md4;
-	pari_CATCH(e_DOMAIN)
-			{
-				result = false;
-			}
-		pari_TRY{ check_quaddisc_imag(test, &md4, ""); };
-	pari_ENDCATCH
-	return result;
+	pari_CATCH(e_DOMAIN) { result = false; }
+	pari_TRY { check_quaddisc_imag(test, &md4, ""); };
+	pari_ENDCATCH return result;
 }
 
 static GEN custom_disc(GEN order, GEN prime) {
@@ -93,7 +89,7 @@ static custom_quadr_t custom_prime_random(GEN order) {
 	GEN v = gen_0;
 	GEN D;
 	long i = 2;
-	//TODO: this is wrong, order = p + 1 +- t  is not kept.
+	// TODO: this is wrong, order = p + 1 +- t  is not kept.
 	while (gequal0(t)) {
 		D = stoi(-4 * (i / 2) + (i % 2));
 		for (long j = 0; j < 100; ++j) {
diff --git a/src/cm/custom.h b/src/cm/custom.h
index ed6f3d3..da5007a 100644
--- a/src/cm/custom.h
+++ b/src/cm/custom.h
@@ -5,8 +5,8 @@
 #ifndef ECGEN_CUSTOM_H
 #define ECGEN_CUSTOM_H
 
-#include "misc/types.h"
 #include "misc/config.h"
+#include "misc/types.h"
 
 typedef struct {
 	GEN p;
@@ -19,6 +19,6 @@ typedef struct {
  * @brief
  * @return
  */
-curve_t * custom_curve();
+curve_t* custom_curve();
 
-#endif //ECGEN_CUSTOM_H
+#endif  // ECGEN_CUSTOM_H
diff --git a/src/cm/p1363.c b/src/cm/p1363.c
index 119b0cb..5ead675 100644
--- a/src/cm/p1363.c
+++ b/src/cm/p1363.c
@@ -79,15 +79,15 @@ void p1363_free(p1363_form_t ***forms, size_t nforms) {
 	}
 }
 
-static GEN p1363_func_F(GEN z) {
+static GEN p1363_func_F(GEN z, long precision) {
 	pari_sp ltop = avma;
 
 	if (gequal0(z)) {
 		return gcopy(gen_1);
 	}
-	GEN sum = cgetc(BIGDEFAULTPREC * 5);
-	gel(sum, 1) = real_1(BIGDEFAULTPREC * 5);
-	gel(sum, 2) = real_0(BIGDEFAULTPREC * 5);
+	GEN sum = cgetc(precision * 2);
+	gel(sum, 1) = real_1(precision);
+	gel(sum, 2) = real_0(precision);
 	pari_printf("initial sum = %Ps\n", sum);
 	GEN last;
 
@@ -99,8 +99,7 @@ static GEN p1363_func_F(GEN z) {
 		GEN quota = divis(subii(j3, j), 2);
 		GEN quotb = divis(addii(j3, j), 2);
 
-		GEN term = gadd(gpow(z, quota, BIGDEFAULTPREC),
-		                gpow(z, quotb, BIGDEFAULTPREC));
+		GEN term = gadd(gpow(z, quota, precision), gpow(z, quotb, precision));
 		pari_printf("%Ps %Ps \n", sum, term);
 		if (mod2(j) == 0) {
 			sum = gadd(sum, term);
@@ -117,37 +116,37 @@ static GEN p1363_func_F(GEN z) {
 	return gerepilecopy(ltop, sum);
 }
 
-static GEN p1363_func_fzero(GEN D, p1363_form_t *form) {
+static GEN p1363_func_fzero(GEN D, p1363_form_t *form, long precision) {
 	pari_sp ltop = avma;
 
-	GEN upper = p1363_func_F(gneg(form->theta));
-	GEN lower = p1363_func_F(gsqr(form->theta));
-	GEN front = gpow(form->theta, gdivgs(gen_m1, 24), BIGDEFAULTPREC * 5);
+	GEN upper = p1363_func_F(gneg(form->theta), precision);
+	GEN lower = p1363_func_F(gsqr(form->theta), precision);
+	GEN front = gpow(form->theta, gdivgs(gen_m1, 24), precision);
 
 	GEN result = gmul(gdiv(upper, lower), front);
 
 	return gerepilecopy(ltop, result);
 }
 
-static GEN p1363_func_fone(GEN D, p1363_form_t *form) {
+static GEN p1363_func_fone(GEN D, p1363_form_t *form, long precision) {
 	pari_sp ltop = avma;
 
-	GEN upper = p1363_func_F(form->theta);
-	GEN lower = p1363_func_F(gsqr(form->theta));
-	GEN front = gpow(form->theta, gdivgs(gen_m1, 24), BIGDEFAULTPREC * 5);
+	GEN upper = p1363_func_F(form->theta, precision);
+	GEN lower = p1363_func_F(gsqr(form->theta), precision);
+	GEN front = gpow(form->theta, gdivgs(gen_m1, 24), precision);
 
 	GEN result = gmul(gdiv(upper, lower), front);
 
 	return gerepilecopy(ltop, result);
 }
 
-static GEN p1363_func_ftwo(GEN D, p1363_form_t *form) {
+static GEN p1363_func_ftwo(GEN D, p1363_form_t *form, long precision) {
 	pari_sp ltop = avma;
 
-	GEN upper = p1363_func_F(gpowgs(form->theta, 4));
-	GEN lower = p1363_func_F(gsqr(form->theta));
-	GEN front = gmul(sqrti(gen_2),
-	                 gpow(form->theta, gdivgs(gen_1, 12), BIGDEFAULTPREC * 5));
+	GEN upper = p1363_func_F(gpowgs(form->theta, 4), precision);
+	GEN lower = p1363_func_F(gsqr(form->theta), precision);
+	GEN front =
+	    gmul(sqrti(gen_2), gpow(form->theta, gdivgs(gen_1, 12), precision));
 
 	GEN result = gmul(gdiv(upper, lower), front);
 
@@ -287,20 +286,20 @@ static void p1363_N(GEN D, p1363_form_t *form) {
 	form->N = gerepileupto(ltop, form->N);
 }
 
-static void p1363_lambda(GEN D, p1363_form_t *form) {
+static void p1363_lambda(GEN D, p1363_form_t *form, long precision) {
 	pari_sp ltop = avma;
-	GEN pik = mulri(mppi(BIGDEFAULTPREC), stoi(form->K));
+	GEN pik = mulri(mppi(precision), stoi(form->K));
 	GEN quot = divrs(pik, 24);
 	form->lambda = gerepileupto(ltop, expIr(quot));
 }
 
-static void p1363_theta(GEN D, p1363_form_t *form) {
+static void p1363_theta(GEN D, p1363_form_t *form, long precision) {
 	pari_sp ltop = avma;
 
-	GEN upper = gadd(gneg(gsqrt(D, BIGDEFAULTPREC)), gmul(form->B, gen_I()));
-	GEN quot = gmul(gdiv(upper, form->A), mppi(BIGDEFAULTPREC));
+	GEN upper = gadd(gneg(gsqrt(D, precision)), gmul(form->B, gen_I()));
+	GEN quot = gmul(gdiv(upper, form->A), mppi(precision));
 
-	form->theta = gerepileupto(ltop, gexp(quot, BIGDEFAULTPREC));
+	form->theta = gerepileupto(ltop, gexp(quot, precision));
 }
 
 /**
@@ -310,22 +309,23 @@ static void p1363_theta(GEN D, p1363_form_t *form) {
  * @param D
  * @param forms
  * @param nforms
- * @return
+ * @return The pari precision required for W_D.
  */
-static long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) {
+long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) {
 	pari_sp ltop = avma;
 	long v0 = 32;
-	GEN mul = divrr(mulrr(mppi(BIGDEFAULTPREC), sqrti(D)), mplog2(BIGDEFAULTPREC));
+	GEN mul =
+	    divrr(mulrr(mppi(BIGDEFAULTPREC), sqrti(D)), mplog2(BIGDEFAULTPREC));
 	GEN sum = gen_0;
 	for (size_t i = 0; i < nforms; ++i) {
 		sum = gadd(sum, ginv(forms[i]->A));
 	}
-	long result = v0 + itos(mpceil(mulrr(mul, sum)));
+	long result = v0 + itos(gceil(gmul(mul, sum)));
 	avma = ltop;
-	return result;
+	return nbits2prec(result);
 }
 
-GEN p1363_invariant(GEN D, p1363_form_t *form) {
+GEN p1363_invariant(GEN D, p1363_form_t *form, long precision) {
 	pari_printf("[A,B,C] = %Pi %Pi %Pi\n", form->A, form->B, form->C);
 	pari_sp ltop = avma;
 
@@ -342,9 +342,9 @@ GEN p1363_invariant(GEN D, p1363_form_t *form) {
 	pari_printf("M = %Pi\n", form->M);
 	p1363_N(D, form);
 	pari_printf("N = %Pi\n", form->N);
-	p1363_lambda(D, form);
+	p1363_lambda(D, form, precision);
 	pari_printf("lambda = %Ps\n", form->lambda);
-	p1363_theta(D, form);
+	p1363_theta(D, form, precision);
 	pari_printf("theta = %Ps\n", form->theta);
 
 	GEN G = gcdii(D, stoi(3));
@@ -352,23 +352,23 @@ GEN p1363_invariant(GEN D, p1363_form_t *form) {
 
 	GEN bl = mulii(form->B, form->L);
 
-	GEN lmbl = gpow(form->lambda, negi(bl), BIGDEFAULTPREC);
+	GEN lmbl = gpow(form->lambda, negi(bl), precision);
 	pari_printf("lmbl = %Ps\n", lmbl);
 
 	GEN mi6 = gneg(gdiv(stoi(form->I), stoi(6)));
-	GEN powmi6 = gpow(gen_2, mi6, BIGDEFAULTPREC);
+	GEN powmi6 = gpow(gen_2, mi6, precision);
 	pari_printf("powmi6 = %Ps\n", powmi6);
 
 	GEN f = gen_0;
 	switch (form->J) {
 		case 0:
-			f = p1363_func_fzero(D, form);
+			f = p1363_func_fzero(D, form, precision);
 			break;
 		case 1:
-			f = p1363_func_fone(D, form);
+			f = p1363_func_fone(D, form, precision);
 			break;
 		case 2:
-			f = p1363_func_ftwo(D, form);
+			f = p1363_func_ftwo(D, form, precision);
 			break;
 		default:
 			pari_err_DOMAIN("p1363_invariant", "J", ">", stoi(2),
@@ -377,14 +377,14 @@ GEN p1363_invariant(GEN D, p1363_form_t *form) {
 
 	pari_printf("f = %Ps\n", f);
 
-	GEN fK = gpow(f, stoi(form->K), BIGDEFAULTPREC);
+	GEN fK = gpow(f, stoi(form->K), precision);
 	pari_printf("fK = %Ps\n", fK);
 
 	GEN in = gmul(lmbl, powmi6);
 	pari_printf("in = %Ps\n", in);
 	GEN inner = gmul(gmul(form->N, in), fK);
 	pari_printf("inner = %Ps\n", inner);
-	GEN result = gpow(inner, G, BIGDEFAULTPREC);
+	GEN result = gpow(inner, G, precision);
 	pari_printf("result = %Ps\n", result);
 
 	return gerepilecopy(ltop, result);
@@ -395,14 +395,17 @@ GEN p1363_poly(GEN D, p1363_form_t **forms, size_t nforms) {
 	long v = fetch_var();
 	name_var(v, "t");
 
+	long precision = p1363_bit_precision(D, forms, nforms);
+
 	GEN terms = gtovec0(gen_0, nforms);
 	for (size_t i = 0; i < nforms; ++i) {
-		gel(terms, i + 1) = gsub(pol_x(v), p1363_invariant(D, forms[i]));
+		gel(terms, i + 1) =
+		    gsub(pol_x(v), p1363_invariant(D, forms[i], precision));
 	}
 
 	GEN result = gen_1;
 	for (size_t i = 0; i < nforms; ++i) {
-		gmulz(result, gel(terms, i + 1), result);
+		result = gmul(result, gel(terms, i + 1));
 	}
 	return gerepilecopy(ltop, result);
 }
diff --git a/src/cm/p1363.h b/src/cm/p1363.h
index ef6f7a2..7201d80 100644
--- a/src/cm/p1363.h
+++ b/src/cm/p1363.h
@@ -32,7 +32,9 @@ size_t p1363_forms(GEN D, p1363_form_t ***forms);
 
 void p1363_free(p1363_form_t ***forms, size_t nforms);
 
-GEN p1363_invariant(GEN D, p1363_form_t *form);
+GEN p1363_invariant(GEN D, p1363_form_t *form, long precision);
+
+long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms);
 
 GEN p1363_poly(GEN D, p1363_form_t **forms, size_t nforms);
 
diff --git a/test/src/cm/test_custom.c b/test/src/cm/test_custom.c
index 3309de7..dd331e6 100644
--- a/test/src/cm/test_custom.c
+++ b/test/src/cm/test_custom.c
@@ -4,11 +4,11 @@
  */
 
 #include <criterion/criterion.h>
+#include "cm/custom.h"
 #include "misc/config.h"
-#include "util/random.h"
 #include "test/default.h"
 #include "test/input.h"
-#include "cm/custom.h"
+#include "util/random.h"
 
 void custom_setup() {
 	default_setup();
@@ -16,9 +16,7 @@ void custom_setup() {
 	random_init();
 }
 
-void custom_teardown() {
-	input_teardown();
-}
+void custom_teardown() { input_teardown(); }
 
 TestSuite(custom, .init = custom_setup, .fini = custom_teardown);
 
diff --git a/test/src/cm/test_p1363.c b/test/src/cm/test_p1363.c
index fdc1823..7bdc6be 100644
--- a/test/src/cm/test_p1363.c
+++ b/test/src/cm/test_p1363.c
@@ -4,9 +4,8 @@
  */
 
 #include <criterion/criterion.h>
-#include "test/default.h"
 #include "cm/p1363.h"
-
+#include "test/default.h"
 
 TestSuite(p1363, .init = default_setup, .fini = default_teardown);
 
@@ -16,18 +15,22 @@ Test(p1363, test_p1363_forms) {
 	size_t nforms = p1363_forms(D, &forms);
 
 	p1363_form_t expected[7] = {
-			{gen_1,   gen_0,    stoi(71)},
-			{stoi(3), gen_1,    stoi(24)},
-			{stoi(3), gen_m1,   stoi(24)},
-			{stoi(8), gen_1,    stoi(9)},
-			{stoi(8), gen_m1,   stoi(9)},
-			{stoi(5), stoi(2),  stoi(15)},
-			{stoi(5), stoi(-2), stoi(15)}
-	};
-	cr_assert_eq(nforms, 7,);
+	    {gen_1, gen_0, stoi(71)},     {stoi(3), gen_1, stoi(24)},
+	    {stoi(3), gen_m1, stoi(24)},  {stoi(8), gen_1, stoi(9)},
+	    {stoi(8), gen_m1, stoi(9)},   {stoi(5), stoi(2), stoi(15)},
+	    {stoi(5), stoi(-2), stoi(15)}};
+	cr_assert_eq(nforms, 7, );
 	for (size_t i = 0; i < nforms; ++i) {
 		cr_assert(equalii(forms[i]->A, expected[i].A), );
 		cr_assert(equalii(forms[i]->B, expected[i].B), );
 		cr_assert(equalii(forms[i]->C, expected[i].C), );
 	}
+}
+
+Test(p1363, test_p1363_poly) {
+	GEN D = stoi(71);
+	p1363_form_t **forms;
+	size_t nforms = p1363_forms(D, &forms);
+	GEN WD = p1363_poly(D, forms, nforms);
+	pari_printf("%Ps\n", WD);
 }
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2


From 07430c61d6a2a1895bbcec729ded99cb69ce1866 Mon Sep 17 00:00:00 2001
From: J08nY
Date: Sat, 7 Apr 2018 14:37:39 +0200
Subject: Fix P1363 poly generation.

---
 src/cm/cm.c    |   3 +-
 src/cm/p1363.c | 118 ++++++++++++++++++++++++++++++++++++++-------------------
 2 files changed, 80 insertions(+), 41 deletions(-)

diff --git a/src/cm/cm.c b/src/cm/cm.c
index 48a9c5f..28815d7 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -15,8 +15,7 @@ int cm_do() {
 	GEN D = stoi(71);
 	p1363_form_t **forms;
 	size_t nforms = p1363_forms(D, &forms);
-	GEN WD = p1363_poly(D, forms, nforms);
-	pari_printf("%Ps\n", WD);
+	p1363_poly(D, forms, nforms);
 	p1363_free(&forms, nforms);
 
 	debug_log_start("Finished Complex Multiplication method");
diff --git a/src/cm/p1363.c b/src/cm/p1363.c
index 5ead675..0485849 100644
--- a/src/cm/p1363.c
+++ b/src/cm/p1363.c
@@ -3,6 +3,7 @@
  * Copyright (C) 2017-2018 J08nY
  */
 #include "p1363.h"
+#include "io/output.h"
 #include "util/memory.h"
 
 static GEN p1363_group(GEN D) {
@@ -85,10 +86,10 @@ static GEN p1363_func_F(GEN z, long precision) {
 	if (gequal0(z)) {
 		return gcopy(gen_1);
 	}
-	GEN sum = cgetc(precision * 2);
+	GEN sum = cgetc(precision);
 	gel(sum, 1) = real_1(precision);
 	gel(sum, 2) = real_0(precision);
-	pari_printf("initial sum = %Ps\n", sum);
+	debug_log("initial sum = %Ps\n", sum);
 	GEN last;
 
 	GEN j = gcopy(gen_1);
@@ -100,10 +101,12 @@ static GEN p1363_func_F(GEN z, long precision) {
 		GEN quotb = divis(addii(j3, j), 2);
 
 		GEN term = gadd(gpow(z, quota, precision), gpow(z, quotb, precision));
-		pari_printf("%Ps %Ps \n", sum, term);
+
 		if (mod2(j) == 0) {
+			debug_log("%Ps (add) %Ps \n", sum, term);
 			sum = gadd(sum, term);
 		} else {
+			debug_log("%Ps (sub) %Ps \n", sum, term);
 			sum = gsub(sum, term);
 		}
 
@@ -111,7 +114,7 @@ static GEN p1363_func_F(GEN z, long precision) {
 		if (gc_needed(ltop, 1)) gerepileall(ltop, 3, &j, &sum, &last);
 	} while (!gequal(sum, last));
 
-	pari_printf("end sum = %Ps, last = %Ps\n", sum, last);
+	debug_log("end sum = %Ps, last = %Ps\n", sum, last);
 
 	return gerepilecopy(ltop, sum);
 }
@@ -121,9 +124,17 @@ static GEN p1363_func_fzero(GEN D, p1363_form_t *form, long precision) {
 
 	GEN upper = p1363_func_F(gneg(form->theta), precision);
 	GEN lower = p1363_func_F(gsqr(form->theta), precision);
-	GEN front = gpow(form->theta, gdivgs(gen_m1, 24), precision);
+	GEN front = gpow(form->theta, mkfrac(gen_m1, stoi(24)), precision);
+
+	debug_log("F(-theta) = %Ps\n", upper);
+	debug_log("F(theta^2) = %Ps\n", lower);
+
+	GEN divd = gdiv(upper, lower);
+	debug_log("F(-theta) / F(theta^2) = %Ps\n", divd);
+
+	GEN result = gmul(divd, front);
 
-	GEN result = gmul(gdiv(upper, lower), front);
+	gel(result, 2) = gneg(gel(result, 2));
 
 	return gerepilecopy(ltop, result);
 }
@@ -133,9 +144,18 @@ static GEN p1363_func_fone(GEN D, p1363_form_t *form, long precision) {
 
 	GEN upper = p1363_func_F(form->theta, precision);
 	GEN lower = p1363_func_F(gsqr(form->theta), precision);
-	GEN front = gpow(form->theta, gdivgs(gen_m1, 24), precision);
+	GEN front = gpow(form->theta, mkfrac(gen_m1, stoi(24)), precision);
+
+	debug_log("F(theta) = %Ps\n", upper);
+	debug_log("F(theta^2) = %Ps\n", lower);
+
+	GEN divd = gdiv(upper, lower);
+	debug_log("F(theta) / F(theta^2) = %Ps\n", divd);
+
+	GEN result = gmul(front, divd);
 
-	GEN result = gmul(gdiv(upper, lower), front);
+	// TODO: WHY????
+	gel(result, 2) = gneg(gel(result, 2));
 
 	return gerepilecopy(ltop, result);
 }
@@ -145,10 +165,19 @@ static GEN p1363_func_ftwo(GEN D, p1363_form_t *form, long precision) {
 
 	GEN upper = p1363_func_F(gpowgs(form->theta, 4), precision);
 	GEN lower = p1363_func_F(gsqr(form->theta), precision);
-	GEN front =
-	    gmul(sqrti(gen_2), gpow(form->theta, gdivgs(gen_1, 12), precision));
+	GEN front = gmul(gsqrt(gen_2, precision),
+	                 gpow(form->theta, mkfrac(gen_1, stoi(12)), precision));
 
-	GEN result = gmul(gdiv(upper, lower), front);
+	debug_log("F(theta^4) = %Ps\n", upper);
+	debug_log("F(theta^2) = %Ps\n", lower);
+
+	GEN divd = gdiv(upper, lower);
+	debug_log("F(theta^4) / F(theta^2) = %Ps\n", divd);
+
+	GEN result = gmul(divd, front);
+
+	// TODO: WHY????
+	gel(result, 2) = gneg(gel(result, 2));
 
 	return gerepilecopy(ltop, result);
 }
@@ -299,7 +328,7 @@ static void p1363_theta(GEN D, p1363_form_t *form, long precision) {
 	GEN upper = gadd(gneg(gsqrt(D, precision)), gmul(form->B, gen_I()));
 	GEN quot = gmul(gdiv(upper, form->A), mppi(precision));
 
-	form->theta = gerepileupto(ltop, gexp(quot, precision));
+	form->theta = gerepilecopy(ltop, gexp(quot, precision));
 }
 
 /**
@@ -313,9 +342,9 @@ static void p1363_theta(GEN D, p1363_form_t *form, long precision) {
  */
 long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) {
 	pari_sp ltop = avma;
-	long v0 = 32;
-	GEN mul =
-	    divrr(mulrr(mppi(BIGDEFAULTPREC), sqrti(D)), mplog2(BIGDEFAULTPREC));
+	long v0 = 64;
+	GEN mul = divrr(mulrr(mppi(BIGDEFAULTPREC), gsqrt(D, BIGDEFAULTPREC)),
+	                mplog2(BIGDEFAULTPREC));
 	GEN sum = gen_0;
 	for (size_t i = 0; i < nforms; ++i) {
 		sum = gadd(sum, ginv(forms[i]->A));
@@ -326,38 +355,38 @@ long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) {
 }
 
 GEN p1363_invariant(GEN D, p1363_form_t *form, long precision) {
-	pari_printf("[A,B,C] = %Pi %Pi %Pi\n", form->A, form->B, form->C);
+	debug_log("[A,B,C] = %Pi %Pi %Pi\n", form->A, form->B, form->C);
 	pari_sp ltop = avma;
 
 	p1363_m8(D, form);
 	p1363_I(D, form);
-	printf("I = %li\n", form->I);
+	debug_log("I = %li\n", form->I);
 	p1363_J(D, form);
-	printf("J = %li\n", form->J);
+	debug_log("J = %li\n", form->J);
 	p1363_K(D, form);
-	printf("K = %li\n", form->K);
+	debug_log("K = %li\n", form->K);
 	p1363_L(D, form);
-	pari_printf("L = %Pi\n", form->L);
+	debug_log("L = %Pi\n", form->L);
 	p1363_M(D, form);
-	pari_printf("M = %Pi\n", form->M);
+	debug_log("M = %Pi\n", form->M);
 	p1363_N(D, form);
-	pari_printf("N = %Pi\n", form->N);
+	debug_log("N = %Pi\n", form->N);
 	p1363_lambda(D, form, precision);
-	pari_printf("lambda = %Ps\n", form->lambda);
+	debug_log("lambda = %Ps\n", form->lambda);
 	p1363_theta(D, form, precision);
-	pari_printf("theta = %Ps\n", form->theta);
+	debug_log("theta = %Ps\n", form->theta);
 
 	GEN G = gcdii(D, stoi(3));
-	pari_printf("G = %Pi\n", G);
+	debug_log("G = %Pi\n", G);
 
 	GEN bl = mulii(form->B, form->L);
 
 	GEN lmbl = gpow(form->lambda, negi(bl), precision);
-	pari_printf("lmbl = %Ps\n", lmbl);
+	debug_log("lmbl = %Ps\n", lmbl);
 
-	GEN mi6 = gneg(gdiv(stoi(form->I), stoi(6)));
+	GEN mi6 = gneg(mkfrac(stoi(form->I), stoi(6)));
 	GEN powmi6 = gpow(gen_2, mi6, precision);
-	pari_printf("powmi6 = %Ps\n", powmi6);
+	debug_log("powmi6 = %Ps\n", powmi6);
 
 	GEN f = gen_0;
 	switch (form->J) {
@@ -375,17 +404,17 @@ GEN p1363_invariant(GEN D, p1363_form_t *form, long precision) {
 			                stoi(form->J));
 	}
 
-	pari_printf("f = %Ps\n", f);
+	debug_log("f = %Ps\n", f);
 
 	GEN fK = gpow(f, stoi(form->K), precision);
-	pari_printf("fK = %Ps\n", fK);
+	debug_log("f^K = %Ps\n", fK);
 
 	GEN in = gmul(lmbl, powmi6);
-	pari_printf("in = %Ps\n", in);
-	GEN inner = gmul(gmul(form->N, in), fK);
-	pari_printf("inner = %Ps\n", inner);
-	GEN result = gpow(inner, G, precision);
-	pari_printf("result = %Ps\n", result);
+	debug_log("l^-BL * 2^-I/6 = %Ps\n", in);
+	GEN inner = gmul(form->N, in);
+	debug_log("N * l^-BL * 2^-I/6 = %Ps\n", inner);
+	GEN result = gpow(gmul(inner, fK), G, precision);
+	debug_log("result = %Ps\n", result);
 
 	return gerepilecopy(ltop, result);
 }
@@ -399,13 +428,24 @@ GEN p1363_poly(GEN D, p1363_form_t **forms, size_t nforms) {
 
 	GEN terms = gtovec0(gen_0, nforms);
 	for (size_t i = 0; i < nforms; ++i) {
-		gel(terms, i + 1) =
-		    gsub(pol_x(v), p1363_invariant(D, forms[i], precision));
+		GEN invariant = p1363_invariant(D, forms[i], precision);
+		gel(terms, i + 1) = gsub(pol_x(v), invariant);
 	}
 
-	GEN result = gen_1;
+	GEN approximate = gen_1;
 	for (size_t i = 0; i < nforms; ++i) {
-		result = gmul(result, gel(terms, i + 1));
+		debug_log("%Ps\n", gel(terms, i + 1));
+		approximate = gmul(approximate, gel(terms, i + 1));
+	}
+	GEN vec = gtovec(approximate);
+	long veclen = glength(vec);
+	for (long i = 1; i <= veclen; ++i) {
+		gel(vec, i) = ground(gel(vec, i));
+	}
+
+	GEN result = pol_0(v);
+	for (long i = 1; i <= veclen; ++i) {
+		result = gadd(result, gmul(gel(vec, i), pol_xn(veclen - i, v)));
 	}
 	return gerepilecopy(ltop, result);
 }
-- 
cgit v1.2.3-70-g09d2


From d47eb30cf7928093e5677f2ab53f1ec4f02e6e06 Mon Sep 17 00:00:00 2001
From: J08nY
Date: Sat, 7 Apr 2018 15:11:46 +0200
Subject: Add test for P1363.

---
 test/src/cm/test_p1363.c | 21 +++++++++++++++++++--
 1 file changed, 19 insertions(+), 2 deletions(-)

diff --git a/test/src/cm/test_p1363.c b/test/src/cm/test_p1363.c
index 7bdc6be..6df19db 100644
--- a/test/src/cm/test_p1363.c
+++ b/test/src/cm/test_p1363.c
@@ -6,8 +6,19 @@
 #include <criterion/criterion.h>
 #include "cm/p1363.h"
 #include "test/default.h"
+#include "test/output.h"
 
-TestSuite(p1363, .init = default_setup, .fini = default_teardown);
+void p1363_setup() {
+	default_setup();
+	output_setup();
+}
+
+void p1363_teardown() {
+	default_teardown();
+	output_teardown();
+}
+
+TestSuite(p1363, .init = p1363_setup, .fini = p1363_teardown);
 
 Test(p1363, test_p1363_forms) {
 	GEN D = stoi(71);
@@ -25,6 +36,7 @@ Test(p1363, test_p1363_forms) {
 		cr_assert(equalii(forms[i]->B, expected[i].B), );
 		cr_assert(equalii(forms[i]->C, expected[i].C), );
 	}
+	p1363_free(&forms, nforms);
 }
 
 Test(p1363, test_p1363_poly) {
@@ -32,5 +44,10 @@ Test(p1363, test_p1363_poly) {
 	p1363_form_t **forms;
 	size_t nforms = p1363_forms(D, &forms);
 	GEN WD = p1363_poly(D, forms, nforms);
-	pari_printf("%Ps\n", WD);
+	GEN coeffs = gtovec(WD);
+
+	GEN expected =
+	    mkvecn(8, gen_1, gen_m2, gen_m1, gen_1, gen_1, gen_1, gen_m1, gen_m1);
+	cr_assert(gequal(coeffs, expected), );
+	p1363_free(&forms, nforms);
 }
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2


From 13ad6fdecf1067a72c5dd7ae995d890792fda31d Mon Sep 17 00:00:00 2001
From: J08nY
Date: Sun, 8 Apr 2018 17:20:18 +0200
Subject: Add p1363_polclass.

---
 src/cm/p1363.c           | 19 ++++++++++---------
 src/cm/p1363.h           | 37 +++++++++++++++++++++++++++++++++++++
 test/src/cm/test_p1363.c | 10 ++++++++++
 3 files changed, 57 insertions(+), 9 deletions(-)

diff --git a/src/cm/p1363.c b/src/cm/p1363.c
index 0485849..fad2a05 100644
--- a/src/cm/p1363.c
+++ b/src/cm/p1363.c
@@ -134,6 +134,7 @@ static GEN p1363_func_fzero(GEN D, p1363_form_t *form, long precision) {
 
 	GEN result = gmul(divd, front);
 
+	// TODO: WHY????
 	gel(result, 2) = gneg(gel(result, 2));
 
 	return gerepilecopy(ltop, result);
@@ -331,15 +332,6 @@ static void p1363_theta(GEN D, p1363_form_t *form, long precision) {
 	form->theta = gerepilecopy(ltop, gexp(quot, precision));
 }
 
-/**
- * Bit-precision computation for a Weber class polynomial from:
- * 		On the Efficient Generation of Elliptic Curves,
- * 		Elisavet Konstantinou, Yiannis C. Stamatiou, Christos Zaroliagis
- * @param D
- * @param forms
- * @param nforms
- * @return The pari precision required for W_D.
- */
 long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) {
 	pari_sp ltop = avma;
 	long v0 = 64;
@@ -449,3 +441,12 @@ GEN p1363_poly(GEN D, p1363_form_t **forms, size_t nforms) {
 	}
 	return gerepilecopy(ltop, result);
 }
+
+GEN p1363_polclass(GEN D) {
+	pari_sp ltop = avma;
+	p1363_form_t **forms;
+	size_t nforms = p1363_forms(D, &forms);
+	GEN WD = p1363_poly(D, forms, nforms);
+	p1363_free(&forms, nforms);
+	return gerepileupto(ltop, WD);
+}
\ No newline at end of file
diff --git a/src/cm/p1363.h b/src/cm/p1363.h
index 7201d80..1064eda 100644
--- a/src/cm/p1363.h
+++ b/src/cm/p1363.h
@@ -28,14 +28,51 @@ typedef struct {
 	GEN theta;
 } p1363_form_t;
 
+/**
+ * Compute all the primitive reduced quadratic forms for a given discriminant.
+ * @param D
+ * @param forms
+ * @return
+ */
 size_t p1363_forms(GEN D, p1363_form_t ***forms);
 
+/**
+ * Free the computed quadratic forms.
+ * @param forms
+ * @param nforms
+ */
 void p1363_free(p1363_form_t ***forms, size_t nforms);
 
+/**
+ * Compute the class invariant for discriminant D and form.
+ * @param D
+ * @param form
+ * @param precision
+ * @return
+ */
 GEN p1363_invariant(GEN D, p1363_form_t *form, long precision);
 
+/**
+ * Bit-precision computation for a Weber class polynomial from:
+ * 		On the Efficient Generation of Elliptic Curves,
+ * 		Elisavet Konstantinou, Yiannis C. Stamatiou, Christos Zaroliagis
+ * @param D
+ * @param forms
+ * @param nforms
+ * @return The pari precision required for W_D.
+ */
 long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms);
 
+/**
+ * Compute the reduced Webe class polynomial for discriminant D and quadratic
+ * forms.
+ * @param D
+ * @param forms
+ * @param nforms
+ * @return
+ */
 GEN p1363_poly(GEN D, p1363_form_t **forms, size_t nforms);
 
+GEN p1363_polclass(GEN D);
+
 #endif  // ECGEN_CM_P1363_H
diff --git a/test/src/cm/test_p1363.c b/test/src/cm/test_p1363.c
index 6df19db..6ff3c10 100644
--- a/test/src/cm/test_p1363.c
+++ b/test/src/cm/test_p1363.c
@@ -50,4 +50,14 @@ Test(p1363, test_p1363_poly) {
 	    mkvecn(8, gen_1, gen_m2, gen_m1, gen_1, gen_1, gen_1, gen_m1, gen_m1);
 	cr_assert(gequal(coeffs, expected), );
 	p1363_free(&forms, nforms);
+}
+
+Test(p1363, test_p1363_polclass) {
+	GEN WD = p1363_polclass(stoi(71));
+
+	GEN coeffs = gtovec(WD);
+
+	GEN expected =
+	    mkvecn(8, gen_1, gen_m2, gen_m1, gen_1, gen_1, gen_1, gen_m1, gen_m1);
+	cr_assert(gequal(coeffs, expected), );
 }
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2


From 65ec6b9789294a5bf8319b5eb14d2cb65b1b8fe8 Mon Sep 17 00:00:00 2001
From: J08nY
Date: Mon, 9 Apr 2018 18:33:15 +0200
Subject: Add the custom CM method, finally working.

---
 src/cm/cm.c               |  12 +--
 src/cm/custom.c           | 198 +++++++++++++++++++++++++++++++++-------------
 src/gen/point.c           |   2 +-
 src/gen/point.h           |   9 +++
 src/util/memory.c         |  17 ++++
 test/src/cm/test_custom.c |  32 ++++----
 6 files changed, 194 insertions(+), 76 deletions(-)

diff --git a/src/cm/cm.c b/src/cm/cm.c
index 28815d7..587f497 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -5,18 +5,18 @@
 #include "cm.h"
 #include "custom.h"
 #include "io/output.h"
+#include "obj/curve.h"
 #include "p1363.h"
 
 int cm_do() {
 	debug_log_start("Starting Complex Multiplication method");
 
-	fprintf(err, "This is *NOT IMPLEMENTED* currently.\n");
+	curve_t *curve = custom_curve();
+	output_o_begin();
+	output_o(curve);
+	output_o_end();
 
-	GEN D = stoi(71);
-	p1363_form_t **forms;
-	size_t nforms = p1363_forms(D, &forms);
-	p1363_poly(D, forms, nforms);
-	p1363_free(&forms, nforms);
+	curve_free(&curve);
 
 	debug_log_start("Finished Complex Multiplication method");
 	return 0;
diff --git a/src/cm/custom.c b/src/cm/custom.c
index 9309ed9..c8d6d8b 100644
--- a/src/cm/custom.c
+++ b/src/cm/custom.c
@@ -3,7 +3,13 @@
  * Copyright (C) 2017 J08nY
  */
 #include "custom.h"
+#include <obj/obj.h>
 #include "io/input.h"
+#include "io/output.h"
+#include "obj/curve.h"
+#include "obj/point.h"
+#include "obj/subgroup.h"
+#include "util/bits.h"
 
 static bool custom_is_disc(GEN test) {
 	bool result = true;
@@ -19,37 +25,7 @@ static GEN custom_disc(GEN order, GEN prime) {
 	GEN t = subii(order, addii(prime, gen_1));
 	GEN v2D = subii(mulis(prime, 4), sqri(t));
 
-	GEN max_divisor = gen_1;
-	GEN min_D = v2D;
-
-	// Optional part that sometimes makes the discriminant smaller, at the cost
-	// factoring v2D.
-	GEN fctr = factor(v2D);
-	long len = glength(gel(fctr, 1));
-	for (long i = 1; i <= len; ++i) {
-		GEN value = gmael(fctr, 1, i);
-		GEN mul = gmael(fctr, 2, i);
-		GEN pow;
-		if (mod2(mul) == 0) {
-			pow = mul;
-		} else {
-			pow = subis(mul, 1);
-		}
-		if (equalis(pow, 0)) {
-			continue;
-		}
-		GEN divi = powii(value, pow);
-		if (mpcmp(divi, max_divisor) > 0) {
-			GEN candidate_D = divii(v2D, divi);
-			if (signe(candidate_D) == 1) {
-				setsigne(candidate_D, -1);
-			}
-			if (custom_is_disc(candidate_D)) {
-				max_divisor = divi;
-				min_D = candidate_D;
-			}
-		}
-	}
+	GEN min_D = core(v2D);
 
 	if (signe(min_D) == 1) {
 		setsigne(min_D, -1);
@@ -81,31 +57,107 @@ static custom_quadr_t custom_prime_input(GEN order) {
 	return result;
 }
 
+static size_t custom_add_primes(GEN r, GEN order, GEN **primes,
+                                size_t nprimes) {
+	size_t nalloc = nprimes;
+	if (nprimes == 0) {
+		nalloc = 10;
+		*primes = try_calloc(sizeof(GEN) * nalloc);
+		debug_log("calloc %lu", sizeof(GEN) * nalloc);
+	}
+
+	GEN logN = ground(glog(order, BIGDEFAULTPREC));
+
+	GEN rlog = mulii(r, logN);
+	GEN rplog = mulii(addis(r, 1), logN);
+
+	forprime_t iter;
+	forprime_init(&iter, rlog, rplog);
+	GEN prime;
+	while ((prime = forprime_next(&iter))) {
+		long k = kronecker(order, prime);
+		if (k == 1) {
+			GEN pstar = prime;
+			GEN ppow = divis(subis(prime, 1), 2);
+			if (mod2(ppow) == 1) {
+				pstar = negi(prime);
+			} else {
+				pstar = gcopy(pstar);
+			}
+			(*primes)[nprimes++] = pstar;
+			if (nprimes == nalloc) {
+				nalloc *= 2;
+				*primes = try_realloc(*primes, sizeof(GEN) * nalloc);
+			}
+		}
+	}
+
+	return nprimes;
+}
+
 static custom_quadr_t custom_prime_random(GEN order) {
 	pari_sp ltop = avma;
 	custom_quadr_t result = {0};
-	GEN p;
-	GEN t = gen_0;
-	GEN v = gen_0;
-	GEN D;
-	long i = 2;
-	// TODO: this is wrong, order = p + 1 +- t  is not kept.
-	while (gequal0(t)) {
-		D = stoi(-4 * (i / 2) + (i % 2));
-		for (long j = 0; j < 100; ++j) {
-			p = random_prime(cfg->bits);
-			if (cornacchia2(negi(D), p, &t, &v)) {
-				break;
+
+	GEN r = gen_0;
+	GEN *Sp;
+	size_t nprimes = custom_add_primes(r, order, &Sp, 0);
+
+	GEN logN = ground(glog(order, BIGDEFAULTPREC));
+	GEN rlog2 = sqri(mulii(r, logN));
+
+	GEN i = gen_0;
+
+	while (true) {
+		GEN imax = int2n(nprimes);
+
+		while (cmpii(i, imax) < 0) {
+			// debug_log("i %Pi", i);
+			pari_sp btop = avma;
+			GEN pprod = gen_1;
+			bits_t *ibits = bits_from_i_len(i, nprimes);
+			for (size_t j = 0; j < nprimes; ++j) {
+				if (GET_BIT(ibits->bits, j) == 1) {
+					// debug_log("multiplying %Pi", Sp[j]);
+					pprod = mulii(pprod, Sp[j]);
+				}
+			}
+			bits_free(&ibits);
+			if (cmpii(pprod, rlog2) < 0 && equalii(modis(pprod, 8), stoi(5))) {
+				// debug_log("candidate D = %Pi", pprod);
+				GEN x;
+				GEN y;
+				cornacchia2(negi(pprod), order, &x, &y);
+				GEN pp1 = addii(addis(order, 1), x);
+				GEN pp2 = subii(addis(order, 1), x);
+				if (isprime(pp1)) {
+					result.p = pp1;
+					result.D = pprod;
+					result.t = x;
+					result.v = gen_0;
+					gerepileall(ltop, 4, &result.p, &result.t, &result.v,
+					            &result.D);
+					try_free(Sp);
+					return result;
+				}
+				if (isprime(pp2)) {
+					result.p = pp2;
+					result.D = pprod;
+					result.t = x;
+					result.v = gen_0;
+					gerepileall(ltop, 4, &result.p, &result.t, &result.v,
+					            &result.D);
+					try_free(Sp);
+					return result;
+				}
 			}
+			avma = btop;
+			i = addis(i, 1);
 		}
-	};
 
-	gerepileall(ltop, 4, &p, &t, &v, &D);
-	result.p = p;
-	result.t = t;
-	result.v = v;
-	result.D = D;
-	return result;
+		r = addis(r, 1);
+		nprimes = custom_add_primes(r, order, &Sp, nprimes);
+	}
 }
 
 curve_t *custom_curve() {
@@ -117,9 +169,49 @@ curve_t *custom_curve() {
 	} else {
 		quadr = custom_prime_input(order);
 	}
-	pari_printf("p = %Pi, D = %Pi, ", quadr.p, quadr.D);
+	debug_log("order = %Pi", order);
+	debug_log("p = %Pi, t = %Pi, v = %Pi, D = %Pi, ", quadr.p, quadr.t, quadr.v,
+	          quadr.D);
 	GEN H = polclass(quadr.D, 0, 0);
-	pari_printf("H = %Ps\n", H);
 
-	return NULL;
+	debug_log("H = %Ps", H);
+
+	GEN r = FpX_roots(H, quadr.p);
+	debug_log("roots = %Ps", r);
+
+	GEN a = NULL;
+	GEN e = NULL;
+	GEN g = NULL;
+
+	long rlen = glength(r);
+	for (long i = 1; i <= rlen; ++i) {
+		GEN root = gel(r, i);
+		a = Fp_div(Fp_mul(stoi(27), root, quadr.p),
+		           Fp_mul(stoi(4), Fp_sub(stoi(1728), root, quadr.p), quadr.p),
+		           quadr.p);
+		e = ellinit(mkvec2(a, negi(a)), quadr.p, 0);
+		g = genrand(e);
+		GEN gmul = ellmul(e, g, order);
+		if (gequal(gmul, gtovec(gen_0))) {
+			debug_log("YES %Ps", e);
+			break;
+		}
+	}
+
+	curve_t *result = curve_new();
+	result->field = quadr.p;
+	result->a = a;
+	result->b = negi(a);
+	result->curve = e;
+	result->order = order;
+	result->generators = subgroups_new(1);
+	result->generators[0] = subgroup_new();
+	result->generators[0]->generator = point_new();
+	result->generators[0]->generator->point = g;
+	result->generators[0]->generator->order = order;
+	result->generators[0]->generator->cofactor = stoi(1);
+	result->generators[0]->npoints = 0;
+	result->ngens = 1;
+
+	return result;
 }
\ No newline at end of file
diff --git a/src/gen/point.c b/src/gen/point.c
index 00e7c3c..98422db 100644
--- a/src/gen/point.c
+++ b/src/gen/point.c
@@ -55,7 +55,7 @@ GENERATOR(points_gen_random) {
 	return 1;
 }
 
-static point_t **points_from_orders(GEN curve, point_t *generator, GEN orders) {
+point_t **points_from_orders(GEN curve, point_t *generator, GEN orders) {
 	size_t norders = (size_t)glength(orders);
 	point_t **result = points_new(norders);
 
diff --git a/src/gen/point.h b/src/gen/point.h
index c10b742..43d162d 100644
--- a/src/gen/point.h
+++ b/src/gen/point.h
@@ -28,6 +28,15 @@ GENERATOR(point_gen_random);
  */
 GENERATOR(points_gen_random);
 
+/**
+ *
+ * @param curve
+ * @param generator
+ * @param orders
+ * @return
+ */
+point_t **points_from_orders(GEN curve, point_t *generator, GEN orders);
+
 /**
  * GENERATOR(gen_f)
  * Generates prime order points using trial division.
diff --git a/src/util/memory.c b/src/util/memory.c
index 439385c..d18c01f 100644
--- a/src/util/memory.c
+++ b/src/util/memory.c
@@ -5,6 +5,9 @@
 #include "memory.h"
 #include <pari/pari.h>
 
+#define ECGEN_PARI_MEM
+#ifdef ECGEN_PARI_MEM
+
 static void *(*malloc_func)(size_t) = pari_malloc;
 
 static void *(*calloc_func)(size_t) = pari_calloc;
@@ -13,6 +16,20 @@ static void *(*realloc_func)(void *, size_t) = pari_realloc;
 
 static void (*free_func)(void *) = pari_free;
 
+#else
+
+static void *(*malloc_func)(size_t) = malloc;
+
+void *calloc_wrapper(size_t n) { return calloc(1, n); }
+
+static void *(*calloc_func)(size_t) = calloc_wrapper;
+
+static void *(*realloc_func)(void *, size_t) = realloc;
+
+static void (*free_func)(void *) = free;
+
+#endif
+
 void *alloc(void *(*fun)(size_t), size_t size) {
 	void *result = fun(size);
 	if (!result) {
diff --git a/test/src/cm/test_custom.c b/test/src/cm/test_custom.c
index dd331e6..c94657d 100644
--- a/test/src/cm/test_custom.c
+++ b/test/src/cm/test_custom.c
@@ -5,36 +5,36 @@
 
 #include <criterion/criterion.h>
 #include "cm/custom.h"
-#include "misc/config.h"
+#include "obj/curve.h"
 #include "test/default.h"
 #include "test/input.h"
+#include "test/output.h"
 #include "util/random.h"
 
 void custom_setup() {
 	default_setup();
 	input_setup();
+	output_setup();
 	random_init();
 }
 
-void custom_teardown() { input_teardown(); }
+void custom_teardown() {
+	default_teardown();
+	input_teardown();
+	output_teardown();
+}
 
 TestSuite(custom, .init = custom_setup, .fini = custom_teardown);
 
 Test(custom, test_curve_one) {
-	cr_skip("Doesnt work yet.");
 	cfg->bits = 128;
 	cfg->cm_order = "263473633827487324648193013259296339349";
-	cfg->random = false;
-
-	fprintf(write_in, "191345683451069598953886857691544477637\n");
-	custom_curve();
-}
-
-Test(custom, test_curve_other) {
-	cr_skip("Doesnt work yet.");
-	cfg->bits = 128;
-	cfg->cm_order = "263473633827487324648193013259296339349";
-	cfg->random = true;
-
-	custom_curve();
+	GEN order = strtoi(cfg->cm_order);
+	cfg->random = RANDOM_ALL;
+
+	curve_t *curve = custom_curve();
+	cr_assert_not_null(curve, );
+	cr_assert(equalii(curve->order, order), );
+	cr_assert(equalii(ellcard(curve->curve, NULL), order), );
+	curve_free(&curve);
 }
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2


From 4662f2e2977e400925b3816d32782d88ccb39504 Mon Sep 17 00:00:00 2001
From: J08nY
Date: Mon, 9 Apr 2018 19:22:54 +0200
Subject: Remove unused code in custom CM method.

---
 src/cm/custom.c | 54 +-----------------------------------------------------
 1 file changed, 1 insertion(+), 53 deletions(-)

diff --git a/src/cm/custom.c b/src/cm/custom.c
index c8d6d8b..33ebde3 100644
--- a/src/cm/custom.c
+++ b/src/cm/custom.c
@@ -3,7 +3,6 @@
  * Copyright (C) 2017 J08nY
  */
 #include "custom.h"
-#include <obj/obj.h>
 #include "io/input.h"
 #include "io/output.h"
 #include "obj/curve.h"
@@ -11,52 +10,6 @@
 #include "obj/subgroup.h"
 #include "util/bits.h"
 
-static bool custom_is_disc(GEN test) {
-	bool result = true;
-	long md4;
-	pari_CATCH(e_DOMAIN) { result = false; }
-	pari_TRY { check_quaddisc_imag(test, &md4, ""); };
-	pari_ENDCATCH return result;
-}
-
-static GEN custom_disc(GEN order, GEN prime) {
-	pari_sp ltop = avma;
-
-	GEN t = subii(order, addii(prime, gen_1));
-	GEN v2D = subii(mulis(prime, 4), sqri(t));
-
-	GEN min_D = core(v2D);
-
-	if (signe(min_D) == 1) {
-		setsigne(min_D, -1);
-	}
-	if (custom_is_disc(min_D)) {
-		return gerepileupto(ltop, min_D);
-	} else {
-		avma = ltop;
-		return NULL;
-	}
-}
-
-static custom_quadr_t custom_prime_input(GEN order) {
-	pari_sp ltop = avma;
-	custom_quadr_t result = {0};
-	GEN p;
-	GEN D;
-	do {
-		p = input_prime("p:", cfg->bits);
-		if (gequalm1(p)) {
-			continue;
-		}
-		D = custom_disc(order, p);
-	} while (D == NULL);
-
-	gerepileall(ltop, 2, &p, &D);
-	result.p = p;
-	result.D = D;
-	return result;
-}
-
 static size_t custom_add_primes(GEN r, GEN order, GEN **primes,
                                 size_t nprimes) {
 	size_t nalloc = nprimes;
@@ -163,12 +116,7 @@ static custom_quadr_t custom_prime_random(GEN order) {
 curve_t *custom_curve() {
 	GEN order = strtoi(cfg->cm_order);
 
-	custom_quadr_t quadr;
-	if (cfg->random) {
-		quadr = custom_prime_random(order);
-	} else {
-		quadr = custom_prime_input(order);
-	}
+	custom_quadr_t quadr = custom_prime_random(order);
 	debug_log("order = %Pi", order);
 	debug_log("p = %Pi, t = %Pi, v = %Pi, D = %Pi, ", quadr.p, quadr.t, quadr.v,
 	          quadr.D);
-- 
cgit v1.2.3-70-g09d2


From 376ff06154b6364c5983a3f67244c8f2d822a282 Mon Sep 17 00:00:00 2001
From: J08nY
Date: Mon, 9 Apr 2018 20:48:04 +0200
Subject: Properly handle edge-cases for CM method.

---
 src/cm/cm.c               | 15 ++++++++++-----
 src/cm/custom.c           | 10 ++++++----
 src/io/cli.c              |  4 ++++
 test/ecgen.sh             |  7 +++++++
 test/src/cm/test_custom.c |  9 ++++++++-
 5 files changed, 35 insertions(+), 10 deletions(-)

diff --git a/src/cm/cm.c b/src/cm/cm.c
index 587f497..8fa174d 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -11,13 +11,18 @@
 int cm_do() {
 	debug_log_start("Starting Complex Multiplication method");
 
+	int result = 0;
 	curve_t *curve = custom_curve();
-	output_o_begin();
-	output_o(curve);
-	output_o_end();
+	if (curve) {
+		output_o_begin();
+		output_o(curve);
+		output_o_end();
 
-	curve_free(&curve);
+		curve_free(&curve);
+	} else {
+		result = 1;
+	}
 
 	debug_log_start("Finished Complex Multiplication method");
-	return 0;
+	return result;
 }
diff --git a/src/cm/custom.c b/src/cm/custom.c
index 33ebde3..71d6625 100644
--- a/src/cm/custom.c
+++ b/src/cm/custom.c
@@ -3,7 +3,6 @@
  * Copyright (C) 2017 J08nY
  */
 #include "custom.h"
-#include "io/input.h"
 #include "io/output.h"
 #include "obj/curve.h"
 #include "obj/point.h"
@@ -16,7 +15,6 @@ static size_t custom_add_primes(GEN r, GEN order, GEN **primes,
 	if (nprimes == 0) {
 		nalloc = 10;
 		*primes = try_calloc(sizeof(GEN) * nalloc);
-		debug_log("calloc %lu", sizeof(GEN) * nalloc);
 	}
 
 	GEN logN = ground(glog(order, BIGDEFAULTPREC));
@@ -48,7 +46,7 @@ static size_t custom_add_primes(GEN r, GEN order, GEN **primes,
 	return nprimes;
 }
 
-static custom_quadr_t custom_prime_random(GEN order) {
+static custom_quadr_t custom_quadr(GEN order) {
 	pari_sp ltop = avma;
 	custom_quadr_t result = {0};
 
@@ -115,8 +113,12 @@ static custom_quadr_t custom_prime_random(GEN order) {
 
 curve_t *custom_curve() {
 	GEN order = strtoi(cfg->cm_order);
+	if (!isprime(order)) {
+		fprintf(err, "Currently, order must be prime for CM to work.\n");
+		return NULL;
+	}
 
-	custom_quadr_t quadr = custom_prime_random(order);
+	custom_quadr_t quadr = custom_quadr(order);
 	debug_log("order = %Pi", order);
 	debug_log("p = %Pi, t = %Pi, v = %Pi, D = %Pi, ", quadr.p, quadr.t, quadr.v,
 	          quadr.D);
diff --git a/src/io/cli.c b/src/io/cli.c
index bc5764f..37c1a8e 100644
--- a/src/io/cli.c
+++ b/src/io/cli.c
@@ -149,6 +149,10 @@ static void cli_end(struct argp_state *state) {
 		argp_failure(state, 1, 0,
 		             "Brainpool algorithm only creates prime field curves.");
 	}
+	if (cfg->method == METHOD_CM && cfg->field == FIELD_BINARY) {
+		argp_failure(state, 1, 0,
+		             "Complex multiplication only creates prime field curves.");
+	}
 	// default values
 	if (!cfg->count) {
 		cfg->count = 1;
diff --git a/test/ecgen.sh b/test/ecgen.sh
index 2383e34..46ee998 100755
--- a/test/ecgen.sh
+++ b/test/ecgen.sh
@@ -136,6 +136,12 @@ function hex() {
 	assert_raises "${ecgen} --fp -r --hex-check=\"abc\" 32 | grep \"abc\""
 }
 
+function cm() {
+    start_test
+    assert_raises "${ecgen} --fp --order=2147483723 32" 1
+    assert_raises "${ecgen} --fp --order=2147483783 32"
+}
+
 . ${ASSERT} -v
 start_suite
 runs
@@ -148,4 +154,5 @@ invalid
 twist
 cli
 hex
+cm
 end_suite ecgen
diff --git a/test/src/cm/test_custom.c b/test/src/cm/test_custom.c
index c94657d..df1ada8 100644
--- a/test/src/cm/test_custom.c
+++ b/test/src/cm/test_custom.c
@@ -30,11 +30,18 @@ Test(custom, test_curve_one) {
 	cfg->bits = 128;
 	cfg->cm_order = "263473633827487324648193013259296339349";
 	GEN order = strtoi(cfg->cm_order);
-	cfg->random = RANDOM_ALL;
 
 	curve_t *curve = custom_curve();
 	cr_assert_not_null(curve, );
 	cr_assert(equalii(curve->order, order), );
 	cr_assert(equalii(ellcard(curve->curve, NULL), order), );
 	curve_free(&curve);
+}
+
+Test(custom, test_curve_other) {
+	cfg->bits = 32;
+	cfg->cm_order = "2147483723";
+
+	curve_t *curve = custom_curve();
+	cr_assert_null(curve, );
 }
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2


From ff8cac25c04afb3a1ecbe617e6bb58122f4a824d Mon Sep 17 00:00:00 2001
From: J08nY
Date: Tue, 10 Apr 2018 01:12:34 +0200
Subject: Retry looking for the correct discriminant in CM.

---
 src/cm/custom.c | 228 +++++++++++++++++++++++++++++++++++++++++---------------
 src/cm/custom.h |   9 ++-
 2 files changed, 174 insertions(+), 63 deletions(-)

diff --git a/src/cm/custom.c b/src/cm/custom.c
index 71d6625..0ebd708 100644
--- a/src/cm/custom.c
+++ b/src/cm/custom.c
@@ -11,6 +11,7 @@
 
 static size_t custom_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;
 	if (nprimes == 0) {
 		nalloc = 10;
@@ -46,71 +47,155 @@ static size_t custom_add_primes(GEN r, GEN order, GEN **primes,
 	return nprimes;
 }
 
-static custom_quadr_t custom_quadr(GEN order) {
-	pari_sp ltop = avma;
-	custom_quadr_t result = {0};
+static void custom_quadr_init(custom_quadr_t *quadr, GEN order) {
+	quadr->D = gen_0;
+	quadr->p = gen_0;
+	quadr->t = gen_0;
+
+	quadr->order = order;
+	quadr->r = gen_0;
+	quadr->i = gen_0;
+	quadr->Sp = NULL;
+	quadr->nprimes = 0;
+}
 
-	GEN r = gen_0;
-	GEN *Sp;
-	size_t nprimes = custom_add_primes(r, order, &Sp, 0);
+static void custom_quadr_next(custom_quadr_t *quadr) {
+	// Ok I get some state in r, i, Sp and nprimes.
+	// Here I want to check if I want to generate more primes into Sp
+	// Then continue with i
 
-	GEN logN = ground(glog(order, BIGDEFAULTPREC));
-	GEN rlog2 = sqri(mulii(r, logN));
+	GEN logN = ground(glog(quadr->order, BIGDEFAULTPREC));
+	GEN rlog2 = sqri(mulii(quadr->r, logN));
 
-	GEN i = gen_0;
+	// When Do I want more primes? If i == imax, or nprimes == 0
+	GEN imax = int2n(quadr->nprimes);
+	if (equalii(quadr->i, imax) || quadr->nprimes == 0) {
+		quadr->nprimes = custom_add_primes(quadr->r, quadr->order, &(quadr->Sp),
+		                                   quadr->nprimes);
+	}
 
 	while (true) {
-		GEN imax = int2n(nprimes);
+		imax = int2n(quadr->nprimes);
 
-		while (cmpii(i, imax) < 0) {
-			// debug_log("i %Pi", i);
+		while (cmpii(quadr->i, imax) < 0) {
+			debug_log("i %Pi", quadr->i);
 			pari_sp btop = avma;
 			GEN pprod = gen_1;
-			bits_t *ibits = bits_from_i_len(i, nprimes);
-			for (size_t j = 0; j < nprimes; ++j) {
+			bits_t *ibits = bits_from_i_len(quadr->i, quadr->nprimes);
+			for (size_t j = 0; j < quadr->nprimes; ++j) {
 				if (GET_BIT(ibits->bits, j) == 1) {
-					// debug_log("multiplying %Pi", Sp[j]);
-					pprod = mulii(pprod, Sp[j]);
+					// debug_log("multiplying %Pi", quadr->Sp[j]);
+					pprod = mulii(pprod, quadr->Sp[j]);
 				}
 			}
 			bits_free(&ibits);
 			if (cmpii(pprod, rlog2) < 0 && equalii(modis(pprod, 8), stoi(5))) {
-				// debug_log("candidate D = %Pi", pprod);
+				debug_log("candidate D = %Pi, rlog2 = %Pi", pprod, rlog2);
 				GEN x;
 				GEN y;
-				cornacchia2(negi(pprod), order, &x, &y);
-				GEN pp1 = addii(addis(order, 1), x);
-				GEN pp2 = subii(addis(order, 1), x);
+				if (!cornacchia2(negi(pprod), quadr->order, &x, &y)) {
+					avma = btop;
+					quadr->i = addis(quadr->i, 1);
+					// debug_log("Cornacchia fail");
+					continue;
+				}
+				GEN pp1 = addii(addis(quadr->order, 1), x);
+				GEN pp2 = subii(addis(quadr->order, 1), x);
 				if (isprime(pp1)) {
-					result.p = pp1;
-					result.D = pprod;
-					result.t = x;
-					result.v = gen_0;
-					gerepileall(ltop, 4, &result.p, &result.t, &result.v,
-					            &result.D);
-					try_free(Sp);
-					return result;
+					quadr->p = pp1;
+					quadr->D = pprod;
+					quadr->t = x;
+					quadr->i = addis(quadr->i, 1);
+					debug_log("good D %Pi", pprod);
+					return;
 				}
 				if (isprime(pp2)) {
-					result.p = pp2;
-					result.D = pprod;
-					result.t = x;
-					result.v = gen_0;
-					gerepileall(ltop, 4, &result.p, &result.t, &result.v,
-					            &result.D);
-					try_free(Sp);
-					return result;
+					quadr->p = pp2;
+					quadr->D = pprod;
+					quadr->t = x;
+					quadr->i = addis(quadr->i, 1);
+					debug_log("good D %Pi", pprod);
+					return;
 				}
 			}
 			avma = btop;
-			i = addis(i, 1);
+			quadr->i = addis(quadr->i, 1);
 		}
 
-		r = addis(r, 1);
-		nprimes = custom_add_primes(r, order, &Sp, nprimes);
+		quadr->r = addis(quadr->r, 1);
+		quadr->nprimes = custom_add_primes(quadr->r, quadr->order, &(quadr->Sp),
+		                                   quadr->nprimes);
+		rlog2 = sqri(mulii(quadr->r, logN));
 	}
 }
 
+static void custom_quadr_free(custom_quadr_t *quadr) { try_free(quadr->Sp); }
+
+/*
+static custom_quadr_t custom_quadr(GEN order) {
+    pari_sp ltop = avma;
+    custom_quadr_t result = {0};
+
+    GEN r = gen_0;
+    GEN *Sp;
+    size_t nprimes = custom_add_primes(r, order, &Sp, 0);
+
+    GEN logN = ground(glog(order, BIGDEFAULTPREC));
+    GEN rlog2 = sqri(mulii(r, logN));
+
+    GEN i = gen_0;
+
+    while (true) {
+        GEN imax = int2n(nprimes);
+
+        while (cmpii(i, imax) < 0) {
+            // debug_log("i %Pi", i);
+            pari_sp btop = avma;
+            GEN pprod = gen_1;
+            bits_t *ibits = bits_from_i_len(i, nprimes);
+            for (size_t j = 0; j < nprimes; ++j) {
+                if (GET_BIT(ibits->bits, j) == 1) {
+                    // debug_log("multiplying %Pi", Sp[j]);
+                    pprod = mulii(pprod, Sp[j]);
+                }
+            }
+            bits_free(&ibits);
+            if (cmpii(pprod, rlog2) < 0 && equalii(modis(pprod, 8), stoi(5))) {
+                // debug_log("candidate D = %Pi", pprod);
+                GEN x;
+                GEN y;
+                cornacchia2(negi(pprod), order, &x, &y);
+                GEN pp1 = addii(addis(order, 1), x);
+                GEN pp2 = subii(addis(order, 1), x);
+                if (isprime(pp1)) {
+                    result.p = pp1;
+                    result.D = pprod;
+                    result.t = x;
+                    gerepileall(ltop, 3, &result.p, &result.t,
+                                &result.D);
+                    try_free(Sp);
+                    return result;
+                }
+                if (isprime(pp2)) {
+                    result.p = pp2;
+                    result.D = pprod;
+                    result.t = x;
+                    gerepileall(ltop, 3, &result.p, &result.t,
+                                &result.D);
+                    try_free(Sp);
+                    return result;
+                }
+            }
+            avma = btop;
+            i = addis(i, 1);
+        }
+
+        r = addis(r, 1);
+        nprimes = custom_add_primes(r, order, &Sp, nprimes);
+    }
+}
+*/
+
 curve_t *custom_curve() {
 	GEN order = strtoi(cfg->cm_order);
 	if (!isprime(order)) {
@@ -118,36 +203,55 @@ curve_t *custom_curve() {
 		return NULL;
 	}
 
-	custom_quadr_t quadr = custom_quadr(order);
-	debug_log("order = %Pi", order);
-	debug_log("p = %Pi, t = %Pi, v = %Pi, D = %Pi, ", quadr.p, quadr.t, quadr.v,
-	          quadr.D);
-	GEN H = polclass(quadr.D, 0, 0);
-
-	debug_log("H = %Ps", H);
-
-	GEN r = FpX_roots(H, quadr.p);
-	debug_log("roots = %Ps", r);
-
 	GEN a = NULL;
 	GEN e = NULL;
 	GEN g = NULL;
 
-	long rlen = glength(r);
-	for (long i = 1; i <= rlen; ++i) {
-		GEN root = gel(r, i);
-		a = Fp_div(Fp_mul(stoi(27), root, quadr.p),
-		           Fp_mul(stoi(4), Fp_sub(stoi(1728), root, quadr.p), quadr.p),
-		           quadr.p);
-		e = ellinit(mkvec2(a, negi(a)), quadr.p, 0);
-		g = genrand(e);
-		GEN gmul = ellmul(e, g, order);
-		if (gequal(gmul, gtovec(gen_0))) {
-			debug_log("YES %Ps", e);
-			break;
+	custom_quadr_t quadr;
+	custom_quadr_init(&quadr, order);
+	while (true) {
+		custom_quadr_next(&quadr);
+
+		debug_log("order = %Pi", order);
+		debug_log("p = %Pi, t = %Pi, D = %Pi, ", quadr.p, quadr.t, quadr.D);
+		GEN H = polclass(quadr.D, 0, 0);
+
+		debug_log("H = %Ps", H);
+
+		GEN r = FpX_roots(H, quadr.p);
+		debug_log("roots = %Ps", r);
+		if (gequal(r, gtovec(gen_0))) {
+			continue;
 		}
+
+		bool has_curve = false;
+
+		long rlen = glength(r);
+		for (long i = 1; i <= rlen; ++i) {
+			GEN root = gel(r, i);
+			a = Fp_div(
+			    Fp_mul(stoi(27), root, quadr.p),
+			    Fp_mul(stoi(4), Fp_sub(stoi(1728), root, quadr.p), quadr.p),
+			    quadr.p);
+			e = ellinit(mkvec2(a, negi(a)), quadr.p, 0);
+			pari_CATCH(e_TYPE) { continue; }
+			pari_TRY { checkell(e); };
+			pari_ENDCATCH{};
+
+			g = genrand(e);
+			GEN gmul = ellmul(e, g, order);
+			if (ell_is_inf(gmul)) {
+				debug_log("YES %Ps", e);
+				has_curve = true;
+				break;
+			}
+		}
+
+		if (has_curve) break;
 	}
 
+	custom_quadr_free(&quadr);
+
 	curve_t *result = curve_new();
 	result->field = quadr.p;
 	result->a = a;
diff --git a/src/cm/custom.h b/src/cm/custom.h
index da5007a..221c8be 100644
--- a/src/cm/custom.h
+++ b/src/cm/custom.h
@@ -9,10 +9,17 @@
 #include "misc/types.h"
 
 typedef struct {
+	/* Stuff filled with custom_quadr_next. */
 	GEN p;
 	GEN t;
-	GEN v;
 	GEN D;
+
+	/* Stuff for saving state. */
+	GEN order;
+	GEN r;
+	GEN i;
+	GEN* Sp;
+	size_t nprimes;
 } custom_quadr_t;
 
 /**
-- 
cgit v1.2.3-70-g09d2


From db440c3d3af3e9dff252501c3459cbafa2d2d047 Mon Sep 17 00:00:00 2001
From: J08nY
Date: Tue, 10 Apr 2018 02:14:06 +0200
Subject: Fix custom CM method.

---
 src/cm/custom.c | 87 ++++++++++-----------------------------------------------
 src/cm/custom.h |  4 ++-
 test/ecgen.sh   |  2 +-
 3 files changed, 18 insertions(+), 75 deletions(-)

diff --git a/src/cm/custom.c b/src/cm/custom.c
index 0ebd708..fd58364 100644
--- a/src/cm/custom.c
+++ b/src/cm/custom.c
@@ -36,14 +36,16 @@ static size_t custom_add_primes(GEN r, GEN order, GEN **primes,
 			} else {
 				pstar = gcopy(pstar);
 			}
-			(*primes)[nprimes++] = pstar;
 			if (nprimes == nalloc) {
 				nalloc *= 2;
 				*primes = try_realloc(*primes, sizeof(GEN) * nalloc);
 			}
+			(*primes)[nprimes++] = pstar;
 		}
 	}
 
+	*primes = try_realloc(*primes, sizeof(GEN) * nprimes);
+
 	return nprimes;
 }
 
@@ -65,9 +67,9 @@ static void custom_quadr_next(custom_quadr_t *quadr) {
 	// Then continue with i
 
 	GEN logN = ground(glog(quadr->order, BIGDEFAULTPREC));
-	GEN rlog2 = sqri(mulii(quadr->r, logN));
+	GEN rlog2 = sqri(mulii(addis(quadr->r, 1), logN));
 
-	// When Do I want more primes? If i == imax, or nprimes == 0
+	// When do I want more primes? If i == imax, or nprimes == 0
 	GEN imax = int2n(quadr->nprimes);
 	if (equalii(quadr->i, imax) || quadr->nprimes == 0) {
 		quadr->nprimes = custom_add_primes(quadr->r, quadr->order, &(quadr->Sp),
@@ -78,7 +80,7 @@ static void custom_quadr_next(custom_quadr_t *quadr) {
 		imax = int2n(quadr->nprimes);
 
 		while (cmpii(quadr->i, imax) < 0) {
-			debug_log("i %Pi", quadr->i);
+			// debug_log("i %Pi", quadr->i);
 			pari_sp btop = avma;
 			GEN pprod = gen_1;
 			bits_t *ibits = bits_from_i_len(quadr->i, quadr->nprimes);
@@ -89,11 +91,15 @@ static void custom_quadr_next(custom_quadr_t *quadr) {
 				}
 			}
 			bits_free(&ibits);
-			if (cmpii(pprod, rlog2) < 0 && equalii(modis(pprod, 8), stoi(5))) {
-				debug_log("candidate D = %Pi, rlog2 = %Pi", pprod, rlog2);
+
+			GEN absp = absi(pprod);
+			long m4 = mod4(absp);
+			if (cmpii(absp, rlog2) < 0 && equalii(modis(pprod, 8), stoi(5)) &&
+			    m4 != 1 && m4 != 2) {
+				debug_log("candidate D = %Pi", pprod);
 				GEN x;
 				GEN y;
-				if (!cornacchia2(negi(pprod), quadr->order, &x, &y)) {
+				if (!cornacchia2(absp, quadr->order, &x, &y)) {
 					avma = btop;
 					quadr->i = addis(quadr->i, 1);
 					// debug_log("Cornacchia fail");
@@ -125,77 +131,12 @@ static void custom_quadr_next(custom_quadr_t *quadr) {
 		quadr->r = addis(quadr->r, 1);
 		quadr->nprimes = custom_add_primes(quadr->r, quadr->order, &(quadr->Sp),
 		                                   quadr->nprimes);
-		rlog2 = sqri(mulii(quadr->r, logN));
+		rlog2 = sqri(mulii(addis(quadr->r, 1), logN));
 	}
 }
 
 static void custom_quadr_free(custom_quadr_t *quadr) { try_free(quadr->Sp); }
 
-/*
-static custom_quadr_t custom_quadr(GEN order) {
-    pari_sp ltop = avma;
-    custom_quadr_t result = {0};
-
-    GEN r = gen_0;
-    GEN *Sp;
-    size_t nprimes = custom_add_primes(r, order, &Sp, 0);
-
-    GEN logN = ground(glog(order, BIGDEFAULTPREC));
-    GEN rlog2 = sqri(mulii(r, logN));
-
-    GEN i = gen_0;
-
-    while (true) {
-        GEN imax = int2n(nprimes);
-
-        while (cmpii(i, imax) < 0) {
-            // debug_log("i %Pi", i);
-            pari_sp btop = avma;
-            GEN pprod = gen_1;
-            bits_t *ibits = bits_from_i_len(i, nprimes);
-            for (size_t j = 0; j < nprimes; ++j) {
-                if (GET_BIT(ibits->bits, j) == 1) {
-                    // debug_log("multiplying %Pi", Sp[j]);
-                    pprod = mulii(pprod, Sp[j]);
-                }
-            }
-            bits_free(&ibits);
-            if (cmpii(pprod, rlog2) < 0 && equalii(modis(pprod, 8), stoi(5))) {
-                // debug_log("candidate D = %Pi", pprod);
-                GEN x;
-                GEN y;
-                cornacchia2(negi(pprod), order, &x, &y);
-                GEN pp1 = addii(addis(order, 1), x);
-                GEN pp2 = subii(addis(order, 1), x);
-                if (isprime(pp1)) {
-                    result.p = pp1;
-                    result.D = pprod;
-                    result.t = x;
-                    gerepileall(ltop, 3, &result.p, &result.t,
-                                &result.D);
-                    try_free(Sp);
-                    return result;
-                }
-                if (isprime(pp2)) {
-                    result.p = pp2;
-                    result.D = pprod;
-                    result.t = x;
-                    gerepileall(ltop, 3, &result.p, &result.t,
-                                &result.D);
-                    try_free(Sp);
-                    return result;
-                }
-            }
-            avma = btop;
-            i = addis(i, 1);
-        }
-
-        r = addis(r, 1);
-        nprimes = custom_add_primes(r, order, &Sp, nprimes);
-    }
-}
-*/
-
 curve_t *custom_curve() {
 	GEN order = strtoi(cfg->cm_order);
 	if (!isprime(order)) {
diff --git a/src/cm/custom.h b/src/cm/custom.h
index 221c8be..ddb89fe 100644
--- a/src/cm/custom.h
+++ b/src/cm/custom.h
@@ -23,7 +23,9 @@ typedef struct {
 } custom_quadr_t;
 
 /**
- * @brief
+ * Algorithm mostly from:
+ *     Constructing elliptic curves of prime order
+ *     by Reinier Broker and Peter Stevenhagen
  * @return
  */
 curve_t* custom_curve();
diff --git a/test/ecgen.sh b/test/ecgen.sh
index 46ee998..319128d 100755
--- a/test/ecgen.sh
+++ b/test/ecgen.sh
@@ -139,7 +139,7 @@ function hex() {
 function cm() {
     start_test
     assert_raises "${ecgen} --fp --order=2147483723 32" 1
-    assert_raises "${ecgen} --fp --order=2147483783 32"
+    assert_raises "${ecgen} --fp --order=123456789012345678901234567890123456789012345678901234568197 137"
 }
 
 . ${ASSERT} -v
-- 
cgit v1.2.3-70-g09d2