9 files changed, 477 insertions, 93 deletions

diff --git a/src/cm/cm.c b/src/cm/cm.c
index ce171d4..8fa174d 100644
--- a/src/cm/cm.c
+++ b/src/cm/cm.c
@@ -3,22 +3,26 @@
  * Copyright (C) 2017-2018 J08nY
  */
 #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");
+	int result = 0;
+	curve_t *curve = custom_curve();
+	if (curve) {
+		output_o_begin();
+		output_o(curve);
+		output_o_end();
 
-	GEN D = stoi(71);
-	form_t **forms;
-	size_t nforms = p1363_forms(D, &forms);
-	for (size_t i = 0; i < nforms; ++i) {
-		p1363_invariant(D, forms[i]);
+		curve_free(&curve);
+	} else {
+		result = 1;
 	}
 
-	p1363_free(&forms, nforms);
 	debug_log_start("Finished Complex Multiplication method");
-	return 0;
+	return result;
 }
diff --git a/src/cm/custom.c b/src/cm/custom.c
new file mode 100644
index 0000000..fd58364
--- /dev/null
+++ b/src/cm/custom.c
@@ -0,0 +1,212 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017 J08nY
+ */
+#include "custom.h"
+#include "io/output.h"
+#include "obj/curve.h"
+#include "obj/point.h"
+#include "obj/subgroup.h"
+#include "util/bits.h"
+
+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;
+		*primes = try_calloc(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);
+			}
+			if (nprimes == nalloc) {
+				nalloc *= 2;
+				*primes = try_realloc(*primes, sizeof(GEN) * nalloc);
+			}
+			(*primes)[nprimes++] = pstar;
+		}
+	}
+
+	*primes = try_realloc(*primes, sizeof(GEN) * nprimes);
+
+	return nprimes;
+}
+
+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;
+}
+
+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(quadr->order, BIGDEFAULTPREC));
+	GEN rlog2 = sqri(mulii(addis(quadr->r, 1), logN));
+
+	// 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) {
+		imax = int2n(quadr->nprimes);
+
+		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(quadr->i, quadr->nprimes);
+			for (size_t j = 0; j < quadr->nprimes; ++j) {
+				if (GET_BIT(ibits->bits, j) == 1) {
+					// debug_log("multiplying %Pi", quadr->Sp[j]);
+					pprod = mulii(pprod, quadr->Sp[j]);
+				}
+			}
+			bits_free(&ibits);
+
+			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(absp, 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)) {
+					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)) {
+					quadr->p = pp2;
+					quadr->D = pprod;
+					quadr->t = x;
+					quadr->i = addis(quadr->i, 1);
+					debug_log("good D %Pi", pprod);
+					return;
+				}
+			}
+			avma = btop;
+			quadr->i = addis(quadr->i, 1);
+		}
+
+		quadr->r = addis(quadr->r, 1);
+		quadr->nprimes = custom_add_primes(quadr->r, quadr->order, &(quadr->Sp),
+		                                   quadr->nprimes);
+		rlog2 = sqri(mulii(addis(quadr->r, 1), logN));
+	}
+}
+
+static void custom_quadr_free(custom_quadr_t *quadr) { try_free(quadr->Sp); }
+
+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;
+	}
+
+	GEN a = NULL;
+	GEN e = NULL;
+	GEN g = NULL;
+
+	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;
+	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/cm/custom.h b/src/cm/custom.h
new file mode 100644
index 0000000..ddb89fe
--- /dev/null
+++ b/src/cm/custom.h
@@ -0,0 +1,33 @@
+/*
+ * ecgen, tool for generating Elliptic curve domain parameters
+ * Copyright (C) 2017 J08nY
+ */
+#ifndef ECGEN_CUSTOM_H
+#define ECGEN_CUSTOM_H
+
+#include "misc/config.h"
+#include "misc/types.h"
+
+typedef struct {
+	/* Stuff filled with custom_quadr_next. */
+	GEN p;
+	GEN t;
+	GEN D;
+
+	/* Stuff for saving state. */
+	GEN order;
+	GEN r;
+	GEN i;
+	GEN* Sp;
+	size_t nprimes;
+} custom_quadr_t;
+
+/**
+ * Algorithm mostly from:
+ *     Constructing elliptic curves of prime order
+ *     by Reinier Broker and Peter Stevenhagen
+ * @return
+ */
+curve_t* custom_curve();
+
+#endif  // ECGEN_CUSTOM_H
diff --git a/src/cm/p1363.c b/src/cm/p1363.c
index 24c65f0..fad2a05 100644
--- a/src/cm/p1363.c
+++ b/src/cm/p1363.c
@@ -3,9 +3,10 @@
  * Copyright (C) 2017-2018 J08nY
  */
 #include "p1363.h"
+#include "io/output.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,15 +50,15 @@ 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, 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 *));
+	size_t nforms = p1363_num(group);
+	*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 +67,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]) {
@@ -79,16 +80,16 @@ void p1363_free(form_t ***forms, size_t nforms) {
 	}
 }
 
-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);
-	pari_printf("initial sum = %Ps\n", sum);
+	GEN sum = cgetc(precision);
+	gel(sum, 1) = real_1(precision);
+	gel(sum, 2) = real_0(precision);
+	debug_log("initial sum = %Ps\n", sum);
 	GEN last;
 
 	GEN j = gcopy(gen_1);
@@ -99,12 +100,13 @@ 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));
-		pari_printf("%Ps %Ps \n", sum, term);
+		GEN term = gadd(gpow(z, quota, precision), gpow(z, quotb, precision));
+
 		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);
 		}
 
@@ -112,51 +114,78 @@ GEN p1363_func_F(GEN z) {
 		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);
 }
 
-GEN p1363_func_fzero(GEN D, 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, 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(gdiv(upper, lower), front);
+	GEN result = gmul(divd, front);
+
+	// TODO: WHY????
+	gel(result, 2) = gneg(gel(result, 2));
 
 	return gerepilecopy(ltop, result);
 }
 
-GEN p1363_func_fone(GEN D, 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, 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(gdiv(upper, lower), front);
+	GEN result = gmul(front, divd);
+
+	// TODO: WHY????
+	gel(result, 2) = gneg(gel(result, 2));
 
 	return gerepilecopy(ltop, result);
 }
 
-GEN p1363_func_ftwo(GEN D, 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(gsqrt(gen_2, precision),
+	                 gpow(form->theta, mkfrac(gen_1, stoi(12)), precision));
+
+	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(gdiv(upper, lower), front);
+	GEN result = gmul(divd, front);
+
+	// TODO: WHY????
+	gel(result, 2) = gneg(gel(result, 2));
 
 	return gerepilecopy(ltop, result);
 }
 
-void p1363_m8(GEN D, 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, form_t *form) {
+static void p1363_I(GEN D, p1363_form_t *form) {
 	switch (form->m8) {
 		case 1:
 		case 2:
@@ -179,7 +208,7 @@ void p1363_I(GEN D, form_t *form) {
 	}
 }
 
-void p1363_J(GEN D, 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 +226,7 @@ void p1363_J(GEN D, form_t *form) {
 	avma = ltop;
 }
 
-void p1363_K(GEN D, form_t *form) {
+static void p1363_K(GEN D, p1363_form_t *form) {
 	switch (form->m8) {
 		case 1:
 		case 2:
@@ -216,7 +245,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 +272,7 @@ void p1363_L(GEN D, form_t *form) {
 	form->L = gerepileupto(ltop, form->L);
 }
 
-void p1363_M(GEN D, 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 +283,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) {
+static void p1363_N(GEN D, p1363_form_t *form) {
 	pari_sp ltop = avma;
 	long ac2 = mod2(mulii(form->A, form->C));
 
@@ -287,100 +316,137 @@ void p1363_N(GEN D, form_t *form) {
 	form->N = gerepileupto(ltop, form->N);
 }
 
-void p1363_lambda(GEN D, 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));
 }
 
-void p1363_theta(GEN D, 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 = gerepilecopy(ltop, gexp(quot, precision));
+}
 
-	form->theta = gerepileupto(ltop, gexp(quot, BIGDEFAULTPREC));
+long p1363_bit_precision(GEN D, p1363_form_t **forms, size_t nforms) {
+	pari_sp ltop = avma;
+	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));
+	}
+	long result = v0 + itos(gceil(gmul(mul, sum)));
+	avma = ltop;
+	return nbits2prec(result);
 }
 
-GEN p1363_invariant(GEN D, form_t *form) {
-	pari_printf("[A,B,C] = %Pi %Pi %Pi\n", form->A, form->B, form->C);
+GEN p1363_invariant(GEN D, p1363_form_t *form, long precision) {
+	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);
-	p1363_lambda(D, form);
-	pari_printf("lambda = %Ps\n", form->lambda);
-	p1363_theta(D, form);
-	pari_printf("theta = %Ps\n", form->theta);
+	debug_log("N = %Pi\n", form->N);
+	p1363_lambda(D, form, precision);
+	debug_log("lambda = %Ps\n", form->lambda);
+	p1363_theta(D, form, precision);
+	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), BIGDEFAULTPREC);
-	pari_printf("lmbl = %Ps\n", lmbl);
+	GEN lmbl = gpow(form->lambda, negi(bl), precision);
+	debug_log("lmbl = %Ps\n", lmbl);
 
-	GEN mi6 = gneg(gdiv(stoi(form->I), stoi(6)));
-	GEN powmi6 = gpow(gen_2, mi6, BIGDEFAULTPREC);
-	pari_printf("powmi6 = %Ps\n", powmi6);
+	GEN mi6 = gneg(mkfrac(stoi(form->I), stoi(6)));
+	GEN powmi6 = gpow(gen_2, mi6, precision);
+	debug_log("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),
 			                stoi(form->J));
 	}
 
-	pari_printf("f = %Ps\n", f);
+	debug_log("f = %Ps\n", f);
 
-	GEN fK = gpow(f, stoi(form->K), BIGDEFAULTPREC);
-	pari_printf("fK = %Ps\n", fK);
+	GEN fK = gpow(f, stoi(form->K), precision);
+	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, BIGDEFAULTPREC);
-	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);
 }
 
-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");
 
+	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]));
+		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) {
-		gmulz(result, gel(terms, i + 1), result);
+		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);
 }
+
+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 d80d70f..1064eda 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,53 @@ typedef struct form_t {
 
 	GEN lambda;
 	GEN theta;
-} form_t;
+} p1363_form_t;
 
-size_t p1363_forms(GEN D, form_t ***forms);
+/**
+ * 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);
 
-void p1363_free(form_t ***forms, size_t nforms);
+/**
+ * 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);
 
-GEN p1363_invariant(GEN D, form_t *form);
+/**
+ * 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_poly(GEN D, form_t **forms, size_t nforms);
+GEN p1363_polclass(GEN D);
 
 #endif  // ECGEN_CM_P1363_H
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
@@ -29,6 +29,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/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/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) {