diff options
author | wukong <rebirthmonkey@gmail.com> | 2015-11-23 17:48:48 +0100 |
---|---|---|
committer | wukong <rebirthmonkey@gmail.com> | 2015-11-23 17:48:48 +0100 |
commit | fca74d4bc3569506a6659880a89aa009dc11f552 (patch) | |
tree | 4cefd06af989608ea8ebd3bc6306889e2a1ad175 /moon-abe/pbc-0.5.14/ecc/mnt.c | |
parent | 840ac3ebca7af381132bf7e93c1e4c0430d6b16a (diff) |
moon-abe cleanup
Change-Id: Ie1259856db03f0b9e80de3e967ec6bd1f03191b3
Diffstat (limited to 'moon-abe/pbc-0.5.14/ecc/mnt.c')
-rw-r--r-- | moon-abe/pbc-0.5.14/ecc/mnt.c | 496 |
1 files changed, 0 insertions, 496 deletions
diff --git a/moon-abe/pbc-0.5.14/ecc/mnt.c b/moon-abe/pbc-0.5.14/ecc/mnt.c deleted file mode 100644 index 230442fc..00000000 --- a/moon-abe/pbc-0.5.14/ecc/mnt.c +++ /dev/null @@ -1,496 +0,0 @@ -// Routines for finding: -// * MNT curves with embedding degree 6 -// * Freeman curves (which have embedding degree 10) - -#include <stdio.h> -#include <stdlib.h> -#include <stdint.h> // for intptr_t -#include <gmp.h> -#include "pbc_mnt.h" -#include "pbc_memory.h" -#include "pbc_utils.h" -#include "misc/darray.h" - -struct pell_solution_s { - int count; - mpz_t minx; //minimal solution of x^2 - Dy^2 = 1 - mpz_t miny; - mpz_t *x; - mpz_t *y; -}; -typedef struct pell_solution_s pell_solution_t[1]; -typedef struct pell_solution_s *pell_solution_ptr; - -static void freempz(void *data) { - mpz_clear(data); - pbc_free(data); -} - -// Solves x^2 - Dy^2 = N where D not a square. -// For square D, we have (x+Dy)(x-Dy) = N so we look at the factors of N. -static void general_pell(pell_solution_t ps, mpz_t D, int N) { - // TODO: Use brute force for small D. - int i, sgnN = N > 0 ? 1 : -1; - intptr_t f, n; - - // Find square factors of N. - darray_t listf; - darray_init(listf); - - f = 1; - for (;;) { - n = f * f; - if (n > abs(N)) break; - if (!(abs(N) % n)) { - darray_append(listf, int_to_voidp(f)); - } - f++; - } - - //a0, twice_a0 don't change once initialized - //a1 is a_i every iteration - //P0, P1 become P_{i-1}, P_i every iteration - //similarly for Q0, Q1 - mpz_t a0, twice_a0, a1; - mpz_t P0, P1; - mpz_t Q0, Q1; - //variables to compute the convergents - mpz_t p0, p1, pnext; - mpz_t q0, q1, qnext; - - int d; - - darray_t listp, listq; - mpz_ptr zptr; - - mpz_init(a0); - mpz_init(twice_a0); - mpz_init(a1); - mpz_init(P0); mpz_init(P1); - mpz_init(Q0); mpz_init(Q1); - mpz_init(p0); mpz_init(p1); mpz_init(pnext); - mpz_init(q0); mpz_init(q1); mpz_init(qnext); - - darray_init(listp); - darray_init(listq); - - mpz_sqrt(a0, D); - mpz_set_ui(P0, 0); - mpz_set_ui(Q0, 1); - - mpz_set(P1, a0); - mpz_mul(Q1, a0, a0); - mpz_sub(Q1, D, Q1); - mpz_add(a1, a0, P1); - mpz_tdiv_q(a1, a1, Q1); - - mpz_add(twice_a0, a0, a0); - - mpz_set(p0, a0); - mpz_set_ui(q0, 1); - mpz_mul(p1, a0, a1); - mpz_add_ui(p1, p1, 1); - mpz_set(q1, a1); - - d = -1; - for(;;) { - if (d == sgnN) { - for (i=0; i<listf->count; i++) { - f = (intptr_t) listf->item[i]; - if (!mpz_cmp_ui(Q1, abs(N) / (f * f))) { -//element_printf("found %Zd, %Zd, %d\n", p0, q0, f); - zptr = (mpz_ptr) pbc_malloc(sizeof(mpz_t)); - mpz_init(zptr); - mpz_set(zptr, p0); - mpz_mul_ui(zptr, p0, f); - darray_append(listp, zptr); - zptr = (mpz_ptr) pbc_malloc(sizeof(mpz_t)); - mpz_init(zptr); - mpz_set(zptr, q0); - mpz_mul_ui(zptr, q0, f); - darray_append(listq, zptr); - } - } - } - - if (!mpz_cmp(twice_a0, a1) && d == 1) break; - //compute more of the continued fraction expansion - mpz_set(P0, P1); - mpz_mul(P1, a1, Q1); - mpz_sub(P1, P1, P0); - mpz_set(Q0, Q1); - mpz_mul(Q1, P1, P1); - mpz_sub(Q1, D, Q1); - mpz_divexact(Q1, Q1, Q0); - mpz_add(a1, a0, P1); - mpz_tdiv_q(a1, a1, Q1); - - //compute next convergent - mpz_mul(pnext, a1, p1); - mpz_add(pnext, pnext, p0); - mpz_set(p0, p1); - mpz_set(p1, pnext); - - mpz_mul(qnext, a1, q1); - mpz_add(qnext, qnext, q0); - mpz_set(q0, q1); - mpz_set(q1, qnext); - d = -d; - } - darray_clear(listf); - - mpz_init(ps->minx); - mpz_init(ps->miny); - mpz_set(ps->minx, p0); - mpz_set(ps->miny, q0); - n = listp->count; - ps->count = n; - if (n) { - ps->x = (mpz_t *) pbc_malloc(sizeof(mpz_t) * n); - ps->y = (mpz_t *) pbc_malloc(sizeof(mpz_t) * n); - for (i = 0; i < n; i++) { - mpz_init(ps->x[i]); - mpz_init(ps->y[i]); - mpz_set(ps->x[i], (mpz_ptr) listp->item[i]); - mpz_set(ps->y[i], (mpz_ptr) listq->item[i]); - } - } - - mpz_clear(a0); - mpz_clear(twice_a0); - mpz_clear(a1); - mpz_clear(P0); mpz_clear(P1); - mpz_clear(Q0); mpz_clear(Q1); - mpz_clear(p0); mpz_clear(p1); mpz_clear(pnext); - mpz_clear(q0); mpz_clear(q1); mpz_clear(qnext); - - darray_forall(listp, freempz); - darray_forall(listq, freempz); - darray_clear(listp); - darray_clear(listq); -} - -static void pell_solution_clear(pell_solution_t ps) { - int i, n = ps->count; - - if (n) { - for (i=0; i<n; i++) { - mpz_clear(ps->x[i]); - mpz_clear(ps->y[i]); - } - pbc_free(ps->x); - pbc_free(ps->y); - } - mpz_clear(ps->minx); - mpz_clear(ps->miny); -} - -void pbc_cm_init(pbc_cm_t cm) { - mpz_init(cm->q); - mpz_init(cm->r); - mpz_init(cm->h); - mpz_init(cm->n); -} - -void pbc_cm_clear(pbc_cm_t cm) { - mpz_clear(cm->q); - mpz_clear(cm->r); - mpz_clear(cm->h); - mpz_clear(cm->n); -} - -static int mnt_step2(int (*callback)(pbc_cm_t, void *), void *data, - unsigned int D, mpz_t U) { - int d; - mpz_t n, l, q; - mpz_t p; - mpz_t r, cofac; - - mpz_init(l); - mpz_mod_ui(l, U, 6); - if (!mpz_cmp_ui(l, 1)) { - mpz_sub_ui(l, U, 1); - d = 1; - } else if (!mpz_cmp_ui(l, 5)) { - mpz_add_ui(l, U, 1); - d = -1; - } else { - mpz_clear(l); - return 0; - } - - mpz_divexact_ui(l, l, 3); - mpz_init(q); - - mpz_mul(q, l, l); - mpz_add_ui(q, q, 1); - if (!mpz_probab_prime_p(q, 10)) { - mpz_clear(q); - mpz_clear(l); - return 0; - } - - mpz_init(n); - if (d < 0) { - mpz_sub(n, q, l); - } else { - mpz_add(n, q, l); - } - - mpz_init(p); - mpz_init(r); - mpz_init(cofac); - { - mpz_set_ui(cofac, 1); - mpz_set(r, n); - mpz_set_ui(p, 2); - if (!mpz_probab_prime_p(r, 10)) for(;;) { - if (mpz_divisible_p(r, p)) do { - mpz_mul(cofac, cofac, p); - mpz_divexact(r, r, p); - } while (mpz_divisible_p(r, p)); - if (mpz_probab_prime_p(r, 10)) break; - //TODO: use a table of primes instead? - mpz_nextprime(p, p); - if (mpz_sizeinbase(p, 2) > 16) { - //printf("has 16+ bit factor\n"); - mpz_clear(r); - mpz_clear(p); - mpz_clear(cofac); - mpz_clear(q); - mpz_clear(l); - mpz_clear(n); - return 0; - } - } - } - - pbc_cm_t cm; - pbc_cm_init(cm); - cm->k = 6; - cm->D = D; - mpz_set(cm->q, q); - mpz_set(cm->r, r); - mpz_set(cm->h, cofac); - mpz_set(cm->n, n); - int res = callback(cm, data); - pbc_cm_clear(cm); - - mpz_clear(cofac); - mpz_clear(r); - mpz_clear(p); - mpz_clear(q); - mpz_clear(l); - mpz_clear(n); - return res; -} - -int pbc_cm_search_d(int (*callback)(pbc_cm_t, void *), void *data, - unsigned int D, unsigned int bitlimit) { - mpz_t D3; - mpz_t t0, t1, t2; - - mpz_init(D3); - mpz_set_ui(D3, D * 3); - - if (mpz_perfect_square_p(D3)) { - // The only squares that differ by 8 are 1 and 9, - // which we get if U=V=1, D=3, but then l is not an integer. - mpz_clear(D3); - return 0; - } - - mpz_init(t0); - mpz_init(t1); - mpz_init(t2); - - pell_solution_t ps; - general_pell(ps, D3, -8); - - int i, n; - int res = 0; - n = ps->count; - if (n) for (;;) { - for (i=0; i<n; i++) { - //element_printf("%Zd, %Zd\n", ps->x[i], ps->y[i]); - res = mnt_step2(callback, data, D, ps->x[i]); - if (res) goto toobig; - //compute next solution as follows - //if p, q is current solution - //compute new solution p', q' via - //(p + q sqrt{3D})(t + u sqrt{3D}) = p' + q' sqrt(3D) - //where t, u is min. solution to Pell equation - mpz_mul(t0, ps->minx, ps->x[i]); - mpz_mul(t1, ps->miny, ps->y[i]); - mpz_mul(t1, t1, D3); - mpz_add(t0, t0, t1); - if (2 * mpz_sizeinbase(t0, 2) > bitlimit + 10) goto toobig; - mpz_mul(t2, ps->minx, ps->y[i]); - mpz_mul(t1, ps->miny, ps->x[i]); - mpz_add(t2, t2, t1); - mpz_set(ps->x[i], t0); - mpz_set(ps->y[i], t2); - } - } -toobig: - - pell_solution_clear(ps); - mpz_clear(t0); - mpz_clear(t1); - mpz_clear(t2); - mpz_clear(D3); - return res; -} - -static int freeman_step2(int (*callback)(pbc_cm_t, void *), void *data, - unsigned int D, mpz_t U) { - mpz_t n, x, q; - mpz_t p; - mpz_t r, cofac; - pbc_cm_t cm; - - mpz_init(x); - mpz_mod_ui(x, U, 15); - if (!mpz_cmp_ui(x, 5)) { - mpz_sub_ui(x, U, 5); - } else if (!mpz_cmp_ui(x, 10)) { - mpz_add_ui(x, U, 5); - } else { - pbc_die("should never reach here"); - mpz_clear(x); - return 0; - } - - mpz_divexact_ui(x, x, 15); - mpz_init(q); - mpz_init(r); - - //q = 25x^4 + 25x^3 + 25x^2 + 10x + 3 - mpz_mul(r, x, x); - mpz_add(q, x, x); - mpz_mul_ui(r, r, 5); - mpz_add(q, q, r); - mpz_mul(r, r, x); - mpz_add(q, q, r); - mpz_mul(r, r, x); - mpz_add(q, q, r); - mpz_mul_ui(q, q, 5); - mpz_add_ui(q, q, 3); - - if (!mpz_probab_prime_p(q, 10)) { - mpz_clear(q); - mpz_clear(r); - mpz_clear(x); - return 0; - } - - //t = 10x^2 + 5x + 3 - //n = q - t + 1 - mpz_init(n); - - mpz_mul_ui(n, x, 5); - mpz_mul(r, n, x); - mpz_add(r, r, r); - mpz_add(n, n, r); - mpz_sub(n, q, n); - mpz_sub_ui(n, n, 2); - - mpz_init(p); - mpz_init(cofac); - { - mpz_set_ui(cofac, 1); - mpz_set(r, n); - mpz_set_ui(p, 2); - if (!mpz_probab_prime_p(r, 10)) for(;;) { - if (mpz_divisible_p(r, p)) do { - mpz_mul(cofac, cofac, p); - mpz_divexact(r, r, p); - } while (mpz_divisible_p(r, p)); - if (mpz_probab_prime_p(r, 10)) break; - //TODO: use a table of primes instead? - mpz_nextprime(p, p); - if (mpz_sizeinbase(p, 2) > 16) { - //printf("has 16+ bit factor\n"); - mpz_clear(r); - mpz_clear(p); - mpz_clear(cofac); - mpz_clear(q); - mpz_clear(x); - mpz_clear(n); - return 0; - } - } - } - - pbc_cm_init(cm); - cm->k = 10; - cm->D = D; - mpz_set(cm->q, q); - mpz_set(cm->r, r); - mpz_set(cm->h, cofac); - mpz_set(cm->n, n); - int res = callback(cm, data); - pbc_cm_clear(cm); - - mpz_clear(cofac); - mpz_clear(r); - mpz_clear(p); - mpz_clear(q); - mpz_clear(x); - mpz_clear(n); - return res; -} - -int pbc_cm_search_g(int (*callback)(pbc_cm_t, void *), void *data, - unsigned int D, unsigned int bitlimit) { - int res = 0; - mpz_t D15; - mpz_t t0, t1, t2; - - mpz_init(D15); - mpz_set_ui(D15, D); - mpz_mul_ui(D15, D15, 15); - if (mpz_perfect_square_p(D15)) { - mpz_clear(D15); - return 0; - } - - mpz_init(t0); - mpz_init(t1); - mpz_init(t2); - - pell_solution_t ps; - general_pell(ps, D15, -20); - - int i, n; - n = ps->count; - if (n) for (;;) { - for (i=0; i<n; i++) { - res = freeman_step2(callback, data, D, ps->x[i]); - if (res) goto toobig; - // Compute next solution as follows: - // If p, q is current solution - // then compute new solution p', q' via - // (p + q sqrt{15D})(t + u sqrt{15D}) = p' + q' sqrt(15D) - // where t, u is min. solution to Pell equation - mpz_mul(t0, ps->minx, ps->x[i]); - mpz_mul(t1, ps->miny, ps->y[i]); - mpz_mul(t1, t1, D15); - mpz_add(t0, t0, t1); - if (2 * mpz_sizeinbase(t0, 2) > bitlimit + 10) goto toobig; - mpz_mul(t2, ps->minx, ps->y[i]); - mpz_mul(t1, ps->miny, ps->x[i]); - mpz_add(t2, t2, t1); - mpz_set(ps->x[i], t0); - mpz_set(ps->y[i], t2); - } - } -toobig: - - pell_solution_clear(ps); - mpz_clear(t0); - mpz_clear(t1); - mpz_clear(t2); - mpz_clear(D15); - return res; -} |