3 files changed, 158 insertions, 2 deletions

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