Add anomalous curve generation

5 files changed, 102 insertions, 27 deletions

diff --git a/src/exhaustive/anomalous.c b/src/exhaustive/anomalous.c
index fea27ed..89cd55b 100644
--- a/src/exhaustive/anomalous.c
+++ b/src/exhaustive/anomalous.c
@@ -6,7 +6,8 @@
  * @file anomalous.c
  */
 #include "anomalous.h"
-#include "math/random.h"
+#include <io/config.h>
+#include <math/types.h>
 
 static disc_t **disc_table;
 
@@ -30,56 +31,110 @@ void anomalous_init() {
 	 */
 	GEN a = stoi(-32);
 	GEN b = stoi(-64);
-	disc_table[0]->D = stoi(11);
+	disc_table[0]->d = stoi(11);
 	disc_table[0]->j = powis(a, 3);
-	disc_table[0]->a = stoi(21);
-	disc_table[1]->D = stoi(19);
+	disc_table[0]->alpha = stoi(21);
+	disc_table[1]->d = stoi(19);
 	disc_table[1]->j = powis(mulis(a, 3), 3);
-	disc_table[1]->a = stoi(3);
-	disc_table[2]->D = stoi(43);
+	disc_table[1]->alpha = stoi(3);
+	disc_table[2]->d = stoi(43);
 	disc_table[2]->j = powis(mulis(b, 16), 3);
-	disc_table[2]->a = stoi(70);
-	disc_table[3]->D = stoi(67);
+	disc_table[2]->alpha = stoi(70);
+	disc_table[3]->d = stoi(67);
 	disc_table[3]->j = powis(mulis(a, 165), 3);
-	disc_table[3]->a = stoi(35805);
-	disc_table[3]->D = stoi(163);
+	disc_table[3]->alpha = stoi(35805);
+	disc_table[3]->d = stoi(163);
 	disc_table[3]->j = powis(mulis(b, 10005), 3);
-	disc_table[3]->a = stoi(3717878010);
+	disc_table[3]->alpha = stoi(3717878010);
 }
 
-static GEN anomalous_prime(GEN D, const config_t *cfg) {
+static GEN anomalous_prime(size_t i, unsigned long bits) {
+	GEN D = disc_table[i]->d;
+
 	pari_sp ltop = avma;
 	GEN last = divis(addis(D, 1), 4);
-	GEN middle;
-	GEN first;
-	GEN p;
-	GEN b;
+	GEN lower = divii(subii(int2n(bits - 1), last), D);
+	GEN upper = divii(subii(int2n(bits), last), D);
+
+	GEN lower_bound = gceil(gdiv(
+	    gsub(gsqrt(addis(mulis(lower, 4), 1), BIGDEFAULTPREC), gen_1), gen_2));
+	GEN upper_bound = gfloor(gdiv(
+	    gsub(gsqrt(addis(mulis(upper, 4), 1), BIGDEFAULTPREC), gen_1), gen_2));
 
+	GEN range = gtovec0(gen_0, 2);
+	gel(range, 1) = lower_bound;
+	gel(range, 2) = upper_bound;
+	GEN first, middle;
+	GEN p, b;
+
+	pari_sp btop = avma;
 	do {
-		b = random_int(cfg->bits);
+		b = genrand(range);
 		middle = mulii(D, b);
 		first = mulii(middle, b);
 		p = addii(addii(first, middle), last);
-		if (gc_needed(ltop, 1)) {
-			gerepileall(ltop, 2, &p, &last);
-		}
+		gerepileall(btop, 2, &p, &last);
 	} while (!isprime(p));
 
-	// TODO: this does not give me a bits sized prime..
-	// but quite a larger one
 	return gerepilecopy(ltop, p);
 }
 
+static GEN anomalous_c(size_t i, GEN p) {
+	pari_sp ltop = avma;
+	long kroneck = kronecker(disc_table[i]->alpha, p);
+	if (!kroneck) {
+		return NULL;
+	}
+	long cp = -kroneck;
+
+	GEN c;
+	pari_sp btop = avma;
+	do {
+		c = genrand(p);
+		gerepileall(btop, 1, &c);
+	} while (kronecker(c, p) != cp);
+
+	return gerepilecopy(ltop, c);
+}
+
 int anomalous_field(curve_t *curve, const config_t *cfg, arg_t *args) {
+	if (!args) {
+		fprintf(stderr, "No args to an arged function. anomalous_field\n");
+		return INT_MIN;
+	}
+	size_t i = *(size_t *)args->args;
+
 	// find suitable prime field
-	curve->field = anomalous_prime(disc_table[0]->D, cfg);
+	pari_sp ltop = avma;
+	GEN p;
+	do {
+		p = anomalous_prime(i, cfg->bits);
+		gerepileall(ltop, 1, &p);
+	} while (!kronecker(disc_table[0]->alpha, p));
+
+	curve->field = gerepilecopy(ltop, p);
 	return 1;
 }
 
 int anomalous_equation(curve_t *curve, const config_t *cfg, arg_t *args) {
+	if (!args) {
+		fprintf(stderr, "No args to an arged function. anomalous_equation\n");
+		return INT_MIN;
+	}
+	size_t i = *(size_t *)args->args;
+
 	// 1 on success, -2 on failure
 	pari_sp ltop = avma;
-	return INT_MIN;  // TODO this
+	GEN c = anomalous_c(i, curve->field);
+
+	GEN a_upper = gmodulo(disc_table[i]->j, curve->field);
+	GEN a_lower = gmodulo(subsi(1728, disc_table[i]->j), curve->field);
+	GEN a = gdiv(a_upper, a_lower);
+
+	curve->a = gmul(stoi(3), gmul(gsqr(c), a));
+	curve->b = gmul(gen_2, gmul(gpowgs(c, 3), a));
+	gerepileall(ltop, 2, &curve->a, &curve->b);
+	return 1;
 }
 
 int anomalous_order(curve_t *curve, const config_t *cfg, arg_t *args) {
diff --git a/src/exhaustive/anomalous.h b/src/exhaustive/anomalous.h
index eff4686..ac86018 100644
--- a/src/exhaustive/anomalous.h
+++ b/src/exhaustive/anomalous.h
@@ -14,9 +14,9 @@
 #include "math/types.h"
 
 typedef struct disc_t {
-	GEN D;
+	GEN d;
 	GEN j;
-	GEN a;
+	GEN alpha;
 } disc_t;
 
 int anomalous_field(curve_t *curve, const config_t *cfg, arg_t *args);
diff --git a/src/exhaustive/exhaustive.c b/src/exhaustive/exhaustive.c
index 5e2fcc8..8fab389 100644
--- a/src/exhaustive/exhaustive.c
+++ b/src/exhaustive/exhaustive.c
@@ -3,6 +3,8 @@
  * Copyright (C) 2017 J08nY
  */
 #include "exhaustive.h"
+#include <io/config.h>
+#include <math/types.h>
 #include "anomalous.h"
 #include "io/output.h"
 #include "math/arg.h"
@@ -93,6 +95,19 @@ static void exhaustive_ainit(arg_t **argss, const config_t *cfg) {
 	for (size_t i = 0; i < OFFSET_END; ++i) {
 		argss[i] = NULL;
 	}
+	if (cfg->anomalous) {
+		arg_t *field_arg = arg_new();
+		arg_t *eq_arg = arg_new();
+		size_t *i = pari_malloc(sizeof(size_t));
+		*i = 3;
+		field_arg->args = i;
+		field_arg->nargs = 1;
+		eq_arg->args = i;
+		eq_arg->nargs = 1;
+		eq_arg->mallocd = i;
+		argss[OFFSET_FIELD] = field_arg;
+		argss[OFFSET_B] = eq_arg;
+	}
 	if (cfg->points.type == POINTS_RANDOM) {
 		arg_t *points_arg = arg_new();
 		points_arg->args = &cfg->points.amount;
@@ -159,7 +174,7 @@ int exhaustive_gen_retry(curve_t *curve, const config_t *cfg,
 			avma = tops[new_state - start_offset];
 		}
 
-		if (diff == 0) {
+		if (diff <= 0) {
 			int tried = ++tries[state - start_offset];
 			if (retry && tried >= retry) {
 				fprintf(stderr, "Reached retry limit: %i, state = %i\n", retry,
diff --git a/src/math/arg.c b/src/math/arg.c
index e313adb..d2eba9b 100644
--- a/src/math/arg.c
+++ b/src/math/arg.c
@@ -16,6 +16,10 @@ arg_t *arg_new(void) {
 
 void arg_free(arg_t **arg) {
 	if (*arg) {
+		if ((*arg)->mallocd) {
+			pari_free((*arg)->mallocd);
+			(*arg)->mallocd = NULL;
+		}
 		pari_free(*arg);
 		*arg = NULL;
 	}
diff --git a/src/math/types.h b/src/math/types.h
index 6793453..8a7d004 100644
--- a/src/math/types.h
+++ b/src/math/types.h
@@ -67,6 +67,7 @@ typedef enum {
 typedef struct {
 	const void *args;
 	size_t nargs;
+	void *mallocd;
 } arg_t;
 
 /**