aboutsummaryrefslogtreecommitdiff
path: root/src/cm/custom.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/cm/custom.c')
-rw-r--r--src/cm/custom.c198
1 files changed, 145 insertions, 53 deletions
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
generated by cgit v1.2.3-70-g09d2 (git 2.51.2) at 2025-11-08 17:39:28 +0000