diff options
Diffstat (limited to 'moon-abe/pbc-0.5.14/ecc/a_param.c')
| -rw-r--r-- | moon-abe/pbc-0.5.14/ecc/a_param.c | 2315 |
1 files changed, 2315 insertions, 0 deletions
diff --git a/moon-abe/pbc-0.5.14/ecc/a_param.c b/moon-abe/pbc-0.5.14/ecc/a_param.c new file mode 100644 index 00000000..6cf8dd96 --- /dev/null +++ b/moon-abe/pbc-0.5.14/ecc/a_param.c @@ -0,0 +1,2315 @@ +#include <stdarg.h> +#include <stdio.h> +#include <stdint.h> // for intptr_t +#include <stdlib.h> //for rand, pbc_malloc, pbc_free +#include <string.h> //for strcmp +#include <gmp.h> +#include "pbc_utils.h" +#include "pbc_field.h" +#include "pbc_fp.h" +#include "pbc_fieldquadratic.h" +#include "pbc_param.h" +#include "pbc_pairing.h" +#include "pbc_curve.h" +#include "pbc_random.h" +#include "pbc_memory.h" +#include "ecc/param.h" +#include "pbc_a_param.h" +#include "pbc_a1_param.h" + +typedef struct { + int exp2; + int exp1; + int sign1; + int sign0; + mpz_t r; // r = 2^exp2 + sign1 * 2^exp1 + sign0 * 1 + mpz_t q; // we work in E(F_q) (and E(F_q^2)) + mpz_t h; // r * h = q + 1 +} *a_param_ptr; + +typedef struct { + field_t Fq, Fq2, Eq; + int exp2, exp1; + int sign1; +} *a_pairing_data_ptr; + +static void a_out_str(FILE *stream, void *data) { + a_param_ptr p = data; + param_out_type(stream, "a"); + param_out_mpz(stream, "q", p->q); + param_out_mpz(stream, "h", p->h); + param_out_mpz(stream, "r", p->r); + param_out_int(stream, "exp2", p->exp2); + param_out_int(stream, "exp1", p->exp1); + param_out_int(stream, "sign1", p->sign1); + param_out_int(stream, "sign0", p->sign0); +} + +static void a_clear(void *data) { + a_param_ptr sp = data; + mpz_clear(sp->r); + mpz_clear(sp->q); + mpz_clear(sp->h); + pbc_free(data); +} + +static void phi_identity(element_ptr out, element_ptr in, pairing_ptr pairing) { + UNUSED_VAR(pairing); + element_set(out, in); +} + +static void compute_abc_tangent(element_ptr a, element_ptr b, element_ptr c, + element_ptr Vx, element_ptr Vy, element_ptr e0) { + //a = -slope_tangent(V.x, V.y); + //b = 1; + //c = -(V.y + aV.x); + //but we multiply by -2*V.y to avoid division so: + //a = -(3 Vx^2 + cc->a) + //b = 2 * Vy + //c = -(2 Vy^2 + a Vx); + element_square(a, Vx); + //element_mul_si(a, a, 3); + element_add(e0, a, a); + element_add(a, e0, a); + element_set1(b); + element_add(a, a, b); + element_neg(a, a); + + element_double(b, Vy); + + element_mul(e0, b, Vy); + element_mul(c, a, Vx); + element_add(c, c, e0); + element_neg(c, c); +} + +static void compute_abc_tangent_proj(element_ptr a, element_ptr b, element_ptr c, + element_ptr Vx, element_ptr Vy, + element_ptr z, element_ptr z2, element_ptr e0) { + //a = -(3x^2 + cca z^4) + //for this case cca = 1 + //b = 2 y z^3 + //c = -(2 y^2 + x a) + //a = z^2 a + element_square(a, z2); + element_square(b, Vx); + ////element_mul_si(b, b, 3); + element_double(e0, b); + element_add(b, e0, b); + element_add(a, a, b); + element_neg(a, a); + + ////element_mul_si(e0, Vy, 2); + element_double(e0, Vy); + element_mul(b, e0, z2); + element_mul(b, b, z); + + element_mul(c, Vx, a); + element_mul(a, a, z2); + element_mul(e0, e0, Vy); + element_add(c, c, e0); + element_neg(c, c); +} + +static void compute_abc_line(element_ptr a, element_ptr b, element_ptr c, + element_ptr Vx, element_ptr Vy, + element_ptr V1x, element_ptr V1y, + element_ptr e0) { + //a = -(B.y - A.y) / (B.x - A.x); + //b = 1; + //c = -(A.y + a * A.x); + //but we'll multiply by B.x - A.x to avoid division, so + //a = -(By - Ay) + //b = Bx - Ax + //c = -(Ay b + a Ax); + element_sub(a, Vy, V1y); + element_sub(b, V1x, Vx); + element_mul(c, Vx, V1y); + element_mul(e0, Vy, V1x); + element_sub(c, c, e0); +} + +struct pp_coeff_s { + element_t a; + element_t b; + element_t c; +}; +typedef struct pp_coeff_s pp_coeff_t[1]; +typedef struct pp_coeff_s *pp_coeff_ptr; + +static void pp_coeff_set(pp_coeff_ptr p, element_t a, element_t b, element_t c) { + element_init(p->a, a->field); + element_init(p->b, b->field); + element_init(p->c, c->field); + element_set(p->a, a); + element_set(p->b, b); + element_set(p->c, c); +} + +static void a_pairing_pp_init(pairing_pp_t p, element_ptr in1, pairing_t pairing) { + int i, n; + a_pairing_data_ptr ainfo = pairing->data; + p->data = pbc_malloc(sizeof(pp_coeff_t) * (ainfo->exp2 + 1)); + pp_coeff_t *coeff = (pp_coeff_t *) p->data; + element_t V, V1; + element_t a, b, c; + element_t e0; + element_ptr Vx, Vy; + element_ptr V1x, V1y; + + #define do_tangent() \ + compute_abc_tangent(a, b, c, Vx, Vy, e0); \ + pp_coeff_set(coeff[i], a, b, c); + + #define do_line() \ + compute_abc_line(a, b, c, Vx, Vy, V1x, V1y, e0); \ + pp_coeff_set(coeff[i], a, b, c); + + element_init(V, ainfo->Eq); + element_init(V1, ainfo->Eq); + element_set(V, in1); + Vx = curve_x_coord(V); + Vy = curve_y_coord(V); + V1x = curve_x_coord(V1); + V1y = curve_y_coord(V1); + element_init(e0, ainfo->Fq); + element_init(a, ainfo->Fq); + element_init(b, ainfo->Fq); + element_init(c, ainfo->Fq); + + n = ainfo->exp1; + for (i=0; i<n; i++) { + do_tangent(); + element_double(V, V); + } + + if (ainfo->sign1 < 0) { + element_neg(V1, V); + } else { + element_set(V1, V); + } + n = ainfo->exp2; + for (; i<n; i++) { + do_tangent(); + element_double(V, V); + } + + do_line(); + + element_clear(e0); + element_clear(a); + element_clear(b); + element_clear(c); + element_clear(V); + element_clear(V1); + #undef do_tangent + #undef do_line +} + +static void a_pairing_pp_clear(pairing_pp_t p) { + a_pairing_data_ptr ainfo = p->pairing->data; + pp_coeff_t *coeff = (pp_coeff_t *) p->data; + int i, n = ainfo->exp2 + 1; + for (i=0; i<n; i++) { + pp_coeff_ptr pp = coeff[i]; + element_clear(pp->a); + element_clear(pp->b); + element_clear(pp->c); + } + pbc_free(p->data); +} + +// Requires cofactor to be odd. +// Overwrites in and temp, out != in. +// Luckily this touchy routine is only used internally. +// TODO: rewrite to allow (out == in)? would simplify a_finalpow() +static void lucas_odd(element_ptr out, element_ptr in, element_ptr temp, mpz_t cofactor) { + element_ptr in0 = element_x(in); + element_ptr in1 = element_y(in); + element_ptr v0 = element_x(out); + element_ptr v1 = element_y(out); + element_ptr t0 = element_x(temp); + element_ptr t1 = element_y(temp); + int j; + + element_set_si(t0, 2); + element_double(t1, in0); + + element_set(v0, t0); + element_set(v1, t1); + + j = mpz_sizeinbase(cofactor, 2) - 1; + for (;;) { + if (!j) { + element_mul(v1, v0, v1); + element_sub(v1, v1, t1); + element_square(v0, v0); + element_sub(v0, v0, t0); + break; + } + if (mpz_tstbit(cofactor, j)) { + element_mul(v0, v0, v1); + element_sub(v0, v0, t1); + element_square(v1, v1); + element_sub(v1, v1, t0); + } else { + element_mul(v1, v0, v1); + element_sub(v1, v1, t1); + element_square(v0, v0); + element_sub(v0, v0, t0); + } + j--; + } + + //assume cofactor = (q + 1) / r is even + //(r should be odd and q + 1 is always even) + //thus v0 = V_k, v1 = V_{k+1} + //and V_{k-1} = P v0 - v1 + + //so U_k = (P V_k - 2 V_{k-1}) / (P^2 - 4) + // = (2 v1 - P v0) / (P^2 - 4) + + element_mul(in0, v0, t1); + element_double(v1, v1); + element_sub(v1, v1, in0); + + element_square(t1, t1); + element_sub(t1, t1, t0); + element_sub(t1, t1, t0); + element_div(v1, v1, t1); + + element_halve(v0, v0); + element_mul(v1, v1, in1); +} + +static inline void a_tateexp(element_ptr out, element_ptr in, element_ptr temp, mpz_t cofactor) { + element_ptr in1 = element_y(in); + //simpler but slower: + //element_pow_mpz(out, f, tateexp); + + //1. Exponentiate by q-1 + //which is equivalent to the following + + element_invert(temp, in); + element_neg(in1, in1); + element_mul(in, in, temp); + + //2. Exponentiate by (q+1)/r + + //Instead of: + // element_pow_mpz(out, in, cofactor); + //we use Lucas sequences (see "Compressed Pairings", Scott and Barreto) + lucas_odd(out, in, temp, cofactor); +} + +//computes a Qx + b Qy + c for type A pairing +static inline void a_miller_evalfn(element_ptr out, + element_ptr a, element_ptr b, element_ptr c, + element_ptr Qx, element_ptr Qy) { + //we'll map Q via (x,y) --> (-x, iy) + //hence Re(a Qx + b Qy + c) = -a Q'x + c and + //Im(a Qx + b Qy + c) = b Q'y + element_mul(element_y(out), a, Qx); + element_sub(element_x(out), c, element_y(out)); + element_mul(element_y(out), b, Qy); +} + +static void a_pairing_pp_apply(element_ptr out, element_ptr in2, pairing_pp_t p) { + //TODO: use proj coords here too to shave off a little time + element_ptr Qx = curve_x_coord(in2); + element_ptr Qy = curve_y_coord(in2); + element_t f, f0; + int i, n; + a_pairing_data_ptr ainfo = p->pairing->data; + pp_coeff_t *coeff = p->data; + element_init(f, ainfo->Fq2); + element_init(f0, ainfo->Fq2); + + element_set1(f); + n = ainfo->exp1; + for (i=0; i<n; i++) { + pp_coeff_ptr pp = coeff[i]; + element_square(f, f); + a_miller_evalfn(f0, pp->a, pp->b, pp->c, Qx, Qy); + element_mul(f, f, f0); + } + if (ainfo->sign1 < 0) { + element_invert(out, f); + } else { + element_set(out, f); + } + n = ainfo->exp2; + for (; i<n; i++) { + element_square(f, f); + pp_coeff_ptr pp = coeff[i]; + a_miller_evalfn(f0, pp->a, pp->b, pp->c, Qx, Qy); + element_mul(f, f, f0); + } + + element_mul(f, f, out); + { + pp_coeff_ptr pp = coeff[i]; + a_miller_evalfn(f0, pp->a, pp->b, pp->c, Qx, Qy); + element_mul(f, f, f0); + } + + a_tateexp(out, f, f0, p->pairing->phikonr); + + element_clear(f); + element_clear(f0); +} + +// in1, in2 are from E(F_q), out from F_q^2. +// Pairing via elliptic nets (see Stange). +static void a_pairing_ellnet(element_ptr out, element_ptr in1, element_ptr in2, + pairing_t pairing) { + element_ptr x = curve_x_coord(in1); + element_ptr y = curve_y_coord(in1); + + element_ptr x2 = curve_x_coord(in2); + element_ptr y2 = curve_y_coord(in2); + + //we map (x2,y2) to (-x2, i y2) before pairing + //notation: cmi means c_{k-i}, ci means c_{k+i} + element_t cm3, cm2, cm1, c0, c1, c2, c3, c4; + element_t dm1, d0, d1; + element_t A, B, C; + + element_init_same_as(cm3, x); + element_init_same_as(cm2, x); + element_init_same_as(cm1, x); + element_init_same_as(c0, x); + element_init_same_as(c1, x); + element_init_same_as(c2, x); + element_init_same_as(c3, x); + element_init_same_as(c4, x); + element_init_same_as(C, x); + + element_init_same_as(dm1, out); + element_init_same_as(d0, out); + element_init_same_as(d1, out); + element_init_same_as(A, x); + element_init_same_as(B, out); + + // c1 = 2y + // c0 = 1 + // cm2 = -1 + // cm3 = -2y + element_double(c1, y); + element_set1(c0); + element_neg(cm3, c1); + element_neg(cm2, c0); + + // a = 1, b = 0 for Y^2 = X^3 + X + //hence c3 = c_{k+3} = c_4 = 4y(x^6 + 5(x^4 - x^2) - 1) + //use cm1, C, c2 as temp variables for now + element_square(cm1, x); + element_square(C, cm1); + element_sub(c2, C, cm1); + element_double(c3, c2); + element_double(c3, c3); + element_add(c3, c3, c2); + element_mul(c2, C, cm1); + element_add(c3, c3, c2); + element_add(c3, c3, cm2); + element_mul(c3, c3, c1); + element_double(c3, c3); + + // c2 = c_3 = 3x^4 + 6x^2 - 1 + element_double(cm1, cm1); + element_add(cm1, cm1, C); + element_double(C, cm1); + element_add(C, C, cm1); + element_add(c2, C, cm2); + + // c4 = c_5 = c_2^3 c_4 - c_3^3 = c1^3 c3 - c2^3 + element_square(C, c1); + element_mul(c4, C, c1); + element_mul(c4, c4, c3); + element_square(C, c2); + element_mul(C, C, c2); + element_sub(c4, c4, C); + + //compute A, B, d1 (which is d_2 since k = 1) + //(recall phi takes x2 to -x2, y2 to i y2) + element_add(A, x, x2); + element_double(C, x); + element_sub(C, C, x2); + element_square(cm1, A); + element_mul(cm1, C, cm1); + element_set(element_x(d1), y); + element_set(element_y(d1), y2); + element_square(d1, d1); + element_sub(element_x(d1), element_x(d1), cm1); + element_neg(B, d1); + element_invert(B, B); + element_invert(A, A); + element_mul(element_x(d1), y, A); + element_neg(element_x(d1), element_x(d1)); + element_mul(element_y(d1), y2, A); + element_square(d1, d1); + element_sub(element_x(d1), C, element_x(d1)); + element_neg(element_y(d1), element_y(d1)); + + // cm1 = 0 + // C = (2y)^-1 + element_set0(cm1); + element_invert(C, c1); + + element_set1(dm1); + element_set1(d0); + + element_t sm2, sm1; + element_t s0, s1, s2, s3; + element_t tm2, tm1; + element_t t0, t1, t2, t3; + element_t e0, e1; + element_t u, v; + + element_init_same_as(sm2, x); + element_init_same_as(sm1, x); + element_init_same_as(s0, x); + element_init_same_as(s1, x); + element_init_same_as(s2, x); + element_init_same_as(s3, x); + + element_init_same_as(tm2, x); + element_init_same_as(tm1, x); + element_init_same_as(t0, x); + element_init_same_as(t1, x); + element_init_same_as(t2, x); + element_init_same_as(t3, x); + + element_init_same_as(e0, x); + element_init_same_as(e1, x); + + element_init_same_as(u, d0); + element_init_same_as(v, d0); + + int m = mpz_sizeinbase(pairing->r, 2) - 2; + for (;;) { + element_square(sm2, cm2); + element_square(sm1, cm1); + element_square(s0, c0); + element_square(s1, c1); + element_square(s2, c2); + element_square(s3, c3); + + element_mul(tm2, cm3, cm1); + element_mul(tm1, cm2, c0); + element_mul(t0, cm1, c1); + element_mul(t1, c0, c2); + element_mul(t2, c1, c3); + element_mul(t3, c2, c4); + + element_square(u, d0); + element_mul(v, dm1, d1); + + if (mpz_tstbit(pairing->r, m)) { + //double-and-add + element_mul(e0, t0, sm2); + element_mul(e1, tm2, s0); + element_sub(cm3, e0, e1); + element_mul(cm3, cm3, C); + + element_mul(e0, t0, sm1); + element_mul(e1, tm1, s0); + element_sub(cm2, e0, e1); + + element_mul(e0, t1, sm1); + element_mul(e1, tm1, s1); + element_sub(cm1, e0, e1); + element_mul(cm1, cm1, C); + + element_mul(e0, t1, s0); + element_mul(e1, t0, s1); + element_sub(c0, e0, e1); + + element_mul(e0, t2, s0); + element_mul(e1, t0, s2); + element_sub(c1, e0, e1); + element_mul(c1, c1, C); + + element_mul(e0, t2, s1); + element_mul(e1, t1, s2); + element_sub(c2, e0, e1); + + element_mul(e0, t3, s1); + element_mul(e1, t1, s3); + element_sub(c3, e0, e1); + element_mul(c3, c3, C); + + element_mul(e0, t3, s2); + element_mul(e1, t2, s3); + element_sub(c4, e0, e1); + + element_mul(element_x(out), element_x(u), t0); + element_mul(element_y(out), element_y(u), t0); + element_mul(element_x(dm1), element_x(v), s0); + element_mul(element_y(dm1), element_y(v), s0); + element_sub(dm1, dm1, out); + + element_mul(element_x(out), element_x(u), t1); + element_mul(element_y(out), element_y(u), t1); + element_mul(element_x(d0), element_x(v), s1); + element_mul(element_y(d0), element_y(v), s1); + element_sub(d0, d0, out); + element_mul(element_x(d0), element_x(d0), A); + element_mul(element_y(d0), element_y(d0), A); + + element_mul(element_x(out), element_x(u), t2); + element_mul(element_y(out), element_y(u), t2); + element_mul(element_x(d1), element_x(v), s2); + element_mul(element_y(d1), element_y(v), s2); + element_sub(d1, d1, out); + element_mul(d1, d1, B); + } else { + //double + element_mul(e0, tm1, sm2); + element_mul(e1, tm2, sm1); + element_sub(cm3, e0, e1); + + element_mul(e0, t0, sm2); + element_mul(e1, tm2, s0); + element_sub(cm2, e0, e1); + element_mul(cm2, cm2, C); + + element_mul(e0, t0, sm1); + element_mul(e1, tm1, s0); + element_sub(cm1, e0, e1); + + element_mul(e0, t1, sm1); + element_mul(e1, tm1, s1); + element_sub(c0, e0, e1); + element_mul(c0, c0, C); + + element_mul(e0, t1, s0); + element_mul(e1, t0, s1); + element_sub(c1, e0, e1); + + element_mul(e0, t2, s0); + element_mul(e1, t0, s2); + element_sub(c2, e0, e1); + element_mul(c2, c2, C); + + element_mul(e0, t2, s1); + element_mul(e1, t1, s2); + element_sub(c3, e0, e1); + + element_mul(e0, t3, s1); + element_mul(e1, t1, s3); + element_sub(c4, e0, e1); + element_mul(c4, c4, C); + + element_mul(element_x(out), element_x(u), tm1); + element_mul(element_y(out), element_y(u), tm1); + element_mul(element_x(dm1), element_x(v), sm1); + element_mul(element_y(dm1), element_y(v), sm1); + element_sub(dm1, dm1, out); + + element_mul(element_x(out), element_x(u), t0); + element_mul(element_y(out), element_y(u), t0); + element_mul(element_x(d0), element_x(v), s0); + element_mul(element_y(d0), element_y(v), s0); + element_sub(d0, d0, out); + + element_mul(element_x(out), element_x(u), t1); + element_mul(element_y(out), element_y(u), t1); + element_mul(element_x(d1), element_x(v), s1); + element_mul(element_y(d1), element_y(v), s1); + element_sub(d1, d1, out); + element_mul(element_x(d1), element_x(d1), A); + element_mul(element_y(d1), element_y(d1), A); + } + if (!m) break; + m--; + } + // since c_k lies base field + // it gets killed by the final powering + //element_invert(c1, c1); + //element_mul(element_x(d1), element_x(d1), c1); + //element_mul(element_y(d1), element_y(d1), c1); + + a_tateexp(out, d1, d0, pairing->phikonr); + + element_clear(dm1); + element_clear(d0); + element_clear(d1); + + element_clear(cm3); + element_clear(cm2); + element_clear(cm1); + element_clear(c0); + element_clear(c1); + element_clear(c2); + element_clear(c3); + element_clear(c4); + + element_clear(sm2); + element_clear(sm1); + element_clear(s0); + element_clear(s1); + element_clear(s2); + element_clear(s3); + + element_clear(tm2); + element_clear(tm1); + element_clear(t0); + element_clear(t1); + element_clear(t2); + element_clear(t3); + + element_clear(e0); + element_clear(e1); + element_clear(A); + element_clear(B); + element_clear(C); + element_clear(u); + element_clear(v); +} + +struct ellnet_pp_st_s { + element_t sm1, s0, s1, s2; + element_t tm1, t0, t1, t2; +}; +typedef struct ellnet_pp_st_s ellnet_pp_st_t[1]; +typedef struct ellnet_pp_st_s *ellnet_pp_st_ptr; + +struct ellnet_pp_s { + element_t x; + element_t y; + ellnet_pp_st_t *seq; +}; +typedef struct ellnet_pp_s ellnet_pp_t[1]; +typedef struct ellnet_pp_s *ellnet_pp_ptr; + +static void a_pairing_ellnet_pp_init(pairing_pp_t p, element_ptr in1, pairing_t pairing) { + element_ptr x = curve_x_coord(in1); + element_ptr y = curve_y_coord(in1); + int i, rbits = mpz_sizeinbase(pairing->r, 2); + ellnet_pp_ptr pp = p->data = pbc_malloc(sizeof(ellnet_pp_t)); + pp->seq = pbc_malloc(sizeof(ellnet_pp_st_t) * rbits); + element_init_same_as(pp->x, x); + element_init_same_as(pp->y, y); + element_set(pp->x, x); + element_set(pp->y, y); + for (i=0; i<rbits; i++) { + ellnet_pp_st_ptr seq = pp->seq[i]; + element_init_same_as(seq->sm1, x); + element_init_same_as(seq->s0, x); + element_init_same_as(seq->s1, x); + element_init_same_as(seq->s2, x); + element_init_same_as(seq->tm1, x); + element_init_same_as(seq->t0, x); + element_init_same_as(seq->t1, x); + element_init_same_as(seq->t2, x); + } + + //we map (x2,y2) to (-x2, i y2) before pairing + //notation: cmi means c_{k-i}, ci means c_{k+i} + element_t cm3, cm2, cm1, c0, c1, c2, c3, c4; + element_t C; + + element_init_same_as(cm3, x); + element_init_same_as(cm2, x); + element_init_same_as(cm1, x); + element_init_same_as(c0, x); + element_init_same_as(c1, x); + element_init_same_as(c2, x); + element_init_same_as(c3, x); + element_init_same_as(c4, x); + element_init_same_as(C, x); + + // c1 = 2y + // c0 = 1 + // cm2 = -1 + // cm3 = -2y + element_double(c1, y); + element_set1(c0); + element_neg(cm3, c1); + element_neg(cm2, c0); + + // a = 1, b = 0 for Y^2 = X^3 + X + //hence c3 = c_{k+3} = c_4 = 4y(x^6 + 5(x^4 - x^2) - 1) + //use cm1, C, c2 as temp variables for now + element_square(cm1, x); + element_square(C, cm1); + element_sub(c2, C, cm1); + element_double(c3, c2); + element_double(c3, c3); + element_add(c3, c3, c2); + element_mul(c2, C, cm1); + element_add(c3, c3, c2); + element_add(c3, c3, cm2); + element_mul(c3, c3, c1); + element_double(c3, c3); + + // c2 = c_3 = 3x^4 + 6x^2 - 1 + element_double(cm1, cm1); + element_add(cm1, cm1, C); + element_double(C, cm1); + element_add(C, C, cm1); + element_add(c2, C, cm2); + + // c4 = c_5 = c_2^3 c_4 - c_3^3 = c1^3 c3 - c2^3 + element_square(C, c1); + element_mul(c4, C, c1); + element_mul(c4, c4, c3); + element_square(C, c2); + element_mul(C, C, c2); + element_sub(c4, c4, C); + + // cm1 = 0 + // C = (2y)^-1 + element_set0(cm1); + element_invert(C, c1); + + int k = 0; + element_t sm2, s3; + element_t tm2, t3; + element_ptr sm1, s0, s1, s2; + element_ptr tm1, t0, t1, t2; + element_t e0, e1; + + element_init_same_as(sm2, x); + element_init_same_as(s3, x); + + element_init_same_as(tm2, x); + element_init_same_as(t3, x); + + element_init_same_as(e0, x); + element_init_same_as(e1, x); + + int m = rbits - 2; + for (;;) { + ellnet_pp_st_ptr seq = pp->seq[k]; + sm1 = seq->sm1; + s0 = seq->s0; + s1 = seq->s1; + s2 = seq->s2; + tm1 = seq->tm1; + t0 = seq->t0; + t1 = seq->t1; + t2 = seq->t2; + + element_square(sm2, cm2); + element_square(sm1, cm1); + element_square(s0, c0); + element_square(s1, c1); + element_square(s2, c2); + element_square(s3, c3); + + element_mul(tm2, cm3, cm1); + element_mul(tm1, cm2, c0); + element_mul(t0, cm1, c1); + element_mul(t1, c0, c2); + element_mul(t2, c1, c3); + element_mul(t3, c2, c4); + + if (!m) break; + k++; + + if (mpz_tstbit(pairing->r, m)) { + //double-and-add + element_mul(e0, t0, sm2); + element_mul(e1, tm2, s0); + element_sub(cm3, e0, e1); + element_mul(cm3, cm3, C); + + element_mul(e0, t0, sm1); + element_mul(e1, tm1, s0); + element_sub(cm2, e0, e1); + + element_mul(e0, t1, sm1); + element_mul(e1, tm1, s1); + element_sub(cm1, e0, e1); + element_mul(cm1, cm1, C); + + element_mul(e0, t1, s0); + element_mul(e1, t0, s1); + element_sub(c0, e0, e1); + + element_mul(e0, t2, s0); + element_mul(e1, t0, s2); + element_sub(c1, e0, e1); + element_mul(c1, c1, C); + + element_mul(e0, t2, s1); + element_mul(e1, t1, s2); + element_sub(c2, e0, e1); + + element_mul(e0, t3, s1); + element_mul(e1, t1, s3); + element_sub(c3, e0, e1); + element_mul(c3, c3, C); + + element_mul(e0, t3, s2); + element_mul(e1, t2, s3); + element_sub(c4, e0, e1); + + } else { + //double + element_mul(e0, tm1, sm2); + element_mul(e1, tm2, sm1); + element_sub(cm3, e0, e1); + + element_mul(e0, t0, sm2); + element_mul(e1, tm2, s0); + element_sub(cm2, e0, e1); + element_mul(cm2, cm2, C); + + element_mul(e0, t0, sm1); + element_mul(e1, tm1, s0); + element_sub(cm1, e0, e1); + + element_mul(e0, t1, sm1); + element_mul(e1, tm1, s1); + element_sub(c0, e0, e1); + element_mul(c0, c0, C); + + element_mul(e0, t1, s0); + element_mul(e1, t0, s1); + element_sub(c1, e0, e1); + + element_mul(e0, t2, s0); + element_mul(e1, t0, s2); + element_sub(c2, e0, e1); + element_mul(c2, c2, C); + + element_mul(e0, t2, s1); + element_mul(e1, t1, s2); + element_sub(c3, e0, e1); + + element_mul(e0, t3, s1); + element_mul(e1, t1, s3); + element_sub(c4, e0, e1); + element_mul(c4, c4, C); + } + m--; + } + + element_clear(cm3); + element_clear(cm2); + element_clear(cm1); + element_clear(c0); + element_clear(c1); + element_clear(c2); + element_clear(c3); + element_clear(c4); + + element_clear(sm2); + element_clear(s3); + + element_clear(tm2); + element_clear(t3); + + element_clear(e0); + element_clear(e1); + element_clear(C); +} + +static void a_pairing_ellnet_pp_clear(pairing_pp_t p) { + ellnet_pp_ptr pp = p->data; + int i, rbits = mpz_sizeinbase(p->pairing->r, 2); + for (i=0; i<rbits; i++) { + ellnet_pp_st_ptr seq = pp->seq[i]; + element_clear(seq->sm1); + element_clear(seq->s0); + element_clear(seq->s1); + element_clear(seq->s2); + element_clear(seq->tm1); + element_clear(seq->t0); + element_clear(seq->t1); + element_clear(seq->t2); + } + element_clear(pp->x); + element_clear(pp->y); + pbc_free(pp->seq); + pbc_free(p->data); +} + +static void a_pairing_ellnet_pp_apply(element_ptr out, element_ptr in2, pairing_pp_t p) { + element_ptr x2 = curve_x_coord(in2); + element_ptr y2 = curve_y_coord(in2); + ellnet_pp_ptr pp = p->data; + int rbits = mpz_sizeinbase(p->pairing->r, 2); + int k = 0; + int m = rbits - 2; + element_t A, B; + element_t e0, e1; + element_t dm1, d0, d1; + element_t u, v; + + element_init_same_as(A, x2); + element_init_same_as(B, out); + element_init_same_as(e0, x2); + element_init_same_as(e1, x2); + element_init_same_as(dm1, out); + element_init_same_as(d0, out); + element_init_same_as(d1, out); + element_init_same_as(u, out); + element_init_same_as(v, out); + + element_add(A, pp->x, x2); + element_double(e0, pp->x); + element_sub(e0, e0, x2); + element_square(e1, A); + element_mul(e1, e0, e1); + element_set(element_x(d1), pp->y); + element_set(element_y(d1), y2); + element_square(d1, d1); + element_sub(element_x(d1), element_x(d1), e1); + element_neg(B, d1); + element_invert(B, B); + element_invert(A, A); + element_mul(element_x(d1), pp->y, A); + element_neg(element_x(d1), element_x(d1)); + element_mul(element_y(d1), y2, A); + element_square(d1, d1); + element_sub(element_x(d1), e0, element_x(d1)); + element_neg(element_y(d1), element_y(d1)); + + element_set1(dm1); + element_set1(d0); + for (;;) { + element_ptr sm1, s0, s1, s2; + element_ptr tm1, t0, t1, t2; + ellnet_pp_st_ptr seq = pp->seq[k]; + sm1 = seq->sm1; + s0 = seq->s0; + s1 = seq->s1; + s2 = seq->s2; + tm1 = seq->tm1; + t0 = seq->t0; + t1 = seq->t1; + t2 = seq->t2; + k++; + + element_square(u, d0); + element_mul(v, dm1, d1); + + if (mpz_tstbit(p->pairing->r, m)) { + //double-and-add + element_mul(element_x(out), element_x(u), t0); + element_mul(element_y(out), element_y(u), t0); + element_mul(element_x(dm1), element_x(v), s0); + element_mul(element_y(dm1), element_y(v), s0); + element_sub(dm1, dm1, out); + + element_mul(element_x(out), element_x(u), t1); + element_mul(element_y(out), element_y(u), t1); + element_mul(element_x(d0), element_x(v), s1); + element_mul(element_y(d0), element_y(v), s1); + element_sub(d0, d0, out); + element_mul(element_x(d0), element_x(d0), A); + element_mul(element_y(d0), element_y(d0), A); + + element_mul(element_x(out), element_x(u), t2); + element_mul(element_y(out), element_y(u), t2); + element_mul(element_x(d1), element_x(v), s2); + element_mul(element_y(d1), element_y(v), s2); + element_sub(d1, d1, out); + element_mul(d1, d1, B); + } else { + //double + element_mul(element_x(out), element_x(u), tm1); + element_mul(element_y(out), element_y(u), tm1); + element_mul(element_x(dm1), element_x(v), sm1); + element_mul(element_y(dm1), element_y(v), sm1); + element_sub(dm1, dm1, out); + + element_mul(element_x(out), element_x(u), t0); + element_mul(element_y(out), element_y(u), t0); + element_mul(element_x(d0), element_x(v), s0); + element_mul(element_y(d0), element_y(v), s0); + element_sub(d0, d0, out); + + element_mul(element_x(out), element_x(u), t1); + element_mul(element_y(out), element_y(u), t1); + element_mul(element_x(d1), element_x(v), s1); + element_mul(element_y(d1), element_y(v), s1); + element_sub(d1, d1, out); + element_mul(element_x(d1), element_x(d1), A); + element_mul(element_y(d1), element_y(d1), A); + } + if (!m) break; + m--; + } + a_tateexp(out, d1, d0, p->pairing->phikonr); + + element_clear(A); + element_clear(B); + element_clear(e0); + element_clear(e1); + element_clear(dm1); + element_clear(d0); + element_clear(d1); + element_clear(u); + element_clear(v); +} + +//in1, in2 are from E(F_q), out from F_q^2 +static void a_pairing_proj(element_ptr out, element_ptr in1, element_ptr in2, + pairing_t pairing) { + a_pairing_data_ptr p = pairing->data; + element_t V, V1; + element_t z, z2; + element_t f, f0, f1; + element_t a, b, c; + element_t e0; + const element_ptr e1 = a, e2 = b, e3 = c; + int i, n; + element_ptr Vx, Vy; + element_ptr V1x, V1y; + element_ptr Qx = curve_x_coord(in2); + element_ptr Qy = curve_y_coord(in2); + + //could save a couple of inversions by avoiding + //this function and rewriting do_line() to handle projective coords + //convert V from weighted projective (Jacobian) to affine + //i.e. (X, Y, Z) --> (X/Z^2, Y/Z^3) + //also sets z to 1 + #define point_to_affine() \ + element_invert(z, z); \ + element_square(e0, z); \ + element_mul(Vx, Vx, e0); \ + element_mul(e0, e0, z); \ + element_mul(Vy, Vy, e0); \ + element_set1(z); \ + element_set1(z2); + + #define proj_double() { \ + /* e0 = 3x^2 + (cc->a) z^4 */ \ + /* for this case a = 1 */ \ + element_square(e0, Vx); \ + /*element_mul_si(e0, e0, 3);*/ \ + element_double(e1, e0); \ + element_add(e0, e1, e0); \ + element_square(e1, z2); \ + element_add(e0, e0, e1); \ + \ + /* z_out = 2 y z */ \ + element_mul(z, Vy, z); \ + /*element_mul_si(z, z, 2);*/ \ + element_double(z, z); \ + element_square(z2, z); \ + \ + /* e1 = 4 x y^2 */ \ + element_square(e2, Vy); \ + element_mul(e1, Vx, e2); \ + /*element_mul_si(e1, e1, 4);*/ \ + element_double(e1, e1); \ + element_double(e1, e1); \ + \ + /* x_out = e0^2 - 2 e1 */ \ + element_double(e3, e1); \ + element_square(Vx, e0); \ + element_sub(Vx, Vx, e3); \ + \ + /* e2 = 8y^4 */ \ + element_square(e2, e2); \ + /*element_mul_si(e2, e2, 8);*/ \ + element_double(e2, e2); \ + element_double(e2, e2); \ + element_double(e2, e2); \ + \ + /*y_out = e0(e1 - x_out) - e2*/\ + element_sub(e1, e1, Vx); \ + element_mul(e0, e0, e1); \ + element_sub(Vy, e0, e2); \ + } + + #define do_tangent() \ + compute_abc_tangent_proj(a, b, c, Vx, Vy, z, z2, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); + + #define do_line() \ + compute_abc_line(a, b, c, Vx, Vy, V1x, V1y, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); + + element_init(V, p->Eq); + element_init(V1, p->Eq); + element_set(V, in1); + + Vx = curve_x_coord(V); + Vy = curve_y_coord(V); + V1x = curve_x_coord(V1); + V1y = curve_y_coord(V1); + + element_init(f, p->Fq2); + element_init(f0, p->Fq2); + element_init(f1, p->Fq2); + element_set1(f); + element_init(a, p->Fq); + element_init(b, p->Fq); + element_init(c, p->Fq); + element_init(e0, p->Fq); + element_init(z, p->Fq); + element_init(z2, p->Fq); + element_set1(z); + element_set1(z2); + n = p->exp1; + for (i=0; i<n; i++) { + //f = f^2 g_V,V(Q) + //where g_V,V = tangent at V + element_square(f, f); + do_tangent(); + proj_double(); + } + point_to_affine(); + if (p->sign1 < 0) { + element_neg(V1, V); + element_invert(f1, f); + } else { + element_set(V1, V); + element_set(f1, f); + } + n = p->exp2; + for (; i<n; i++) { + element_square(f, f); + do_tangent(); + proj_double(); + } + + element_mul(f, f, f1); + point_to_affine(); + do_line(); + + a_tateexp(out, f, f0, pairing->phikonr); + + element_clear(f); + element_clear(f0); + element_clear(f1); + element_clear(z); + element_clear(z2); + element_clear(V); + element_clear(V1); + element_clear(a); + element_clear(b); + element_clear(c); + element_clear(e0); + #undef point_to_affine + #undef proj_double + #undef do_tangent + #undef do_line +} + +//in1, in2 are from E(F_q), out from F_q^2 +static void a_pairing_affine(element_ptr out, element_ptr in1, element_ptr in2, + pairing_t pairing) { + a_pairing_data_ptr p = pairing->data; + element_t V, V1; + element_t f, f0, f1; + element_t a, b, c; + element_t e0; + int i, n; + element_ptr Qx = curve_x_coord(in2); + element_ptr Qy = curve_y_coord(in2); + element_ptr Vx, Vy; + element_ptr V1x, V1y; + + #define do_tangent() \ + compute_abc_tangent(a, b, c, Vx, Vy, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); + + #define do_line() \ + compute_abc_line(a, b, c, Vx, Vy, V1x, V1y, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); + + element_init(V, p->Eq); + element_init(V1, p->Eq); + Vx = curve_x_coord(V); + Vy = curve_y_coord(V); + + V1x = curve_x_coord(V1); + V1y = curve_y_coord(V1); + + element_set(V, in1); + element_init(f, p->Fq2); + element_init(f0, p->Fq2); + element_init(f1, p->Fq2); + element_set1(f); + element_init(a, p->Fq); + element_init(b, p->Fq); + element_init(c, p->Fq); + element_init(e0, p->Fq); + n = p->exp1; + for (i=0; i<n; i++) { + //f = f^2 g_V,V(Q) + //where g_V,V = tangent at V + element_square(f, f); + do_tangent(); + element_double(V, V); + } + if (p->sign1 < 0) { + element_neg(V1, V); + element_invert(f1, f); + } else { + element_set(V1, V); + element_set(f1, f); + } + n = p->exp2; + for (; i<n; i++) { + element_square(f, f); + do_tangent(); + element_double(V, V); + } + + element_mul(f, f, f1); + do_line(); + + a_tateexp(out, f, f0, pairing->phikonr); + + element_clear(f); + element_clear(f0); + element_clear(f1); + element_clear(V); + element_clear(V1); + element_clear(a); + element_clear(b); + element_clear(c); + element_clear(e0); + #undef do_tangent + #undef do_line +} + +// On Computing Products of Pairing +//in1, in2 are from E(F_q), out from F_q^2 +void a_pairings_affine(element_ptr out, element_t in1[], element_t in2[], + int n_prod, pairing_t pairing) { + a_pairing_data_ptr p = pairing->data; + element_t* V = pbc_malloc(sizeof(element_t)*n_prod); + element_t* V1 = pbc_malloc(sizeof(element_t)*n_prod); + element_t f, f0, f1; + element_t a, b, c; + element_t e0; + int i, j, n; + element_ptr Qx, Qy; + element_ptr Vx, Vy; + element_ptr V1x, V1y; + + #define do_tangents() \ + for(j=0; j<n_prod; j++){ \ + Vx = curve_x_coord(V[j]); \ + Vy = curve_y_coord(V[j]); \ + Qx = curve_x_coord(in2[j]); \ + Qy = curve_y_coord(in2[j]); \ + \ + compute_abc_tangent(a, b, c, Vx, Vy, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); \ + } + + #define do_lines() \ + for(j=0;j<n_prod;j++){ \ + Vx = curve_x_coord(V[j]); \ + Vy = curve_y_coord(V[j]); \ + V1x = curve_x_coord(V1[j]); \ + V1y = curve_y_coord(V1[j]); \ + Qx = curve_x_coord(in2[j]); \ + Qy = curve_y_coord(in2[j]); \ + \ + compute_abc_line(a, b, c, Vx, Vy, V1x, V1y, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); \ + } + + for(i=0; i<n_prod; i++){ + element_init(V[i],p->Eq); + element_init(V1[i],p->Eq); + element_set(V[i],in1[i]); + } + + + element_init(f, p->Fq2); + element_init(f0, p->Fq2); + element_init(f1, p->Fq2); + element_set1(f); + element_init(a, p->Fq); + element_init(b, p->Fq); + element_init(c, p->Fq); + element_init(e0, p->Fq); + n = p->exp1; + for (i=0; i<n; i++) { + //f = f^2 g_V,V(Q) + //where g_V,V = tangent at V + element_square(f, f); + do_tangents(); + element_multi_double(V, V, n_prod); //V_i = V_i + V_i for all i at one time. + } + if (p->sign1 < 0) { + for(j=0; j<n_prod; j++){ + element_neg(V1[j], V[j]); + } + element_invert(f1, f); + } else { + for(j=0; j<n_prod; j++){ + element_set(V1[j], V[j]); + } + element_set(f1, f); + } + n = p->exp2; + for (; i<n; i++) { + element_square(f, f); + do_tangents(); + element_multi_double(V, V, n_prod); + } + + element_mul(f, f, f1); + do_lines(); + + a_tateexp(out, f, f0, pairing->phikonr); + + element_clear(f); + element_clear(f0); + element_clear(f1); + for(j=0;j<n_prod;j++){ + element_clear(V[j]); + element_clear(V1[j]); + } + pbc_free(V); + pbc_free(V1); + element_clear(a); + element_clear(b); + element_clear(c); + element_clear(e0); + #undef do_tangents + #undef do_lines +} + +static void a_pairing_clear(pairing_t pairing) { + field_clear(pairing->GT); + + a_pairing_data_ptr p = pairing->data; + field_clear(p->Eq); + field_clear(p->Fq); + field_clear(p->Fq2); + pbc_free(p); + + mpz_clear(pairing->r); + mpz_clear(pairing->phikonr); + field_clear(pairing->Zr); +} + +static void a_pairing_option_set(pairing_t pairing, char *key, char *value) { + if (!strcmp(key, "method")) { + if (!strcmp(value, "miller")) { + pairing->map = a_pairing_proj; + pairing->pp_init = a_pairing_pp_init; + pairing->pp_clear = a_pairing_pp_clear; + pairing->pp_apply = a_pairing_pp_apply; + } else if (!strcmp(value, "miller-affine")) { + pairing->map = a_pairing_affine; + pairing->pp_init = a_pairing_pp_init; + pairing->pp_clear = a_pairing_pp_clear; + pairing->pp_apply = a_pairing_pp_apply; + } else if (!strcmp(value, "shipsey-stange")) { + pairing->map = a_pairing_ellnet; + pairing->pp_init = a_pairing_ellnet_pp_init; + pairing->pp_clear = a_pairing_ellnet_pp_clear; + pairing->pp_apply = a_pairing_ellnet_pp_apply; + } + } +} + +static void a_finalpow(element_t e) { + pairing_ptr pairing = e->field->pairing; + element_t t0, t1; + element_init_same_as(t0, e->data); + element_init_same_as(t1, e->data); + a_tateexp(t0, e->data, t1, pairing->phikonr); + element_set(e->data, t0); + element_clear(t0); + element_clear(t1); +} + +static void a_init_pairing(pairing_ptr pairing, void *data) { + a_param_ptr param = data; + element_t a, b; + a_pairing_data_ptr p; + + p = pairing->data = pbc_malloc(sizeof(*p)); + p->exp2 = param->exp2; + p->exp1 = param->exp1; + p->sign1 = param->sign1; + mpz_init(pairing->r); + mpz_set(pairing->r, param->r); + field_init_fp(pairing->Zr, pairing->r); + pairing->map = a_pairing_proj; + pairing->prod_pairings = a_pairings_affine; + + field_init_fp(p->Fq, param->q); + element_init(a, p->Fq); + element_init(b, p->Fq); + element_set1(a); + element_set0(b); + field_init_curve_ab(p->Eq, a, b, pairing->r, param->h); + element_clear(a); + element_clear(b); + + field_init_fi(p->Fq2, p->Fq); + + //k=2, hence phi_k(q) = q + 1, phikonr = (q+1)/r + mpz_init(pairing->phikonr); + mpz_set(pairing->phikonr, param->h); + + pairing->G1 = p->Eq; + pairing->G2 = pairing->G1; + pairing->phi = phi_identity; + pairing_GT_init(pairing, p->Fq2); + pairing->finalpow = a_finalpow; + + pairing->clear_func = a_pairing_clear; + pairing->option_set = a_pairing_option_set; + pairing->pp_init = a_pairing_pp_init; + pairing->pp_clear = a_pairing_pp_clear; + pairing->pp_apply = a_pairing_pp_apply; +} + +static void a_param_init(pbc_param_ptr par) { + static pbc_param_interface_t interface = {{ + a_clear, + a_init_pairing, + a_out_str, + }}; + par->api = interface; + a_param_ptr p = par->data = pbc_malloc(sizeof(*p)); + mpz_init(p->r); + mpz_init(p->q); + mpz_init(p->h); +} + +// Public interface for type A pairings: + +int pbc_param_init_a(pbc_param_ptr par, struct symtab_s *tab) { + a_param_init(par); + a_param_ptr p = par->data; + + int err = 0; + err += lookup_mpz(p->q, tab, "q"); + err += lookup_mpz(p->r, tab, "r"); + err += lookup_mpz(p->h, tab, "h"); + err += lookup_int(&p->exp2, tab, "exp2"); + err += lookup_int(&p->exp1, tab, "exp1"); + err += lookup_int(&p->sign1, tab, "sign1"); + err += lookup_int(&p->sign0, tab, "sign0"); + return err; +} + +void pbc_param_init_a_gen(pbc_param_ptr par, int rbits, int qbits) { + a_param_init(par); + a_param_ptr sp = par->data; + int found = 0; + + mpz_ptr q = sp->q; + mpz_ptr r = sp->r; + mpz_ptr h = sp->h; + + do { + int i; + mpz_set_ui(r, 0); + + if (rand() % 2) { + sp->exp2 = rbits - 1; + sp->sign1 = 1; + } else { + sp->exp2 = rbits; + sp->sign1 = -1; + } + mpz_setbit(r, sp->exp2); + + //use q as a temp variable + mpz_set_ui(q, 0); + sp->exp1 = (rand() % (sp->exp2 - 1)) + 1; + mpz_setbit(q, sp->exp1); + if (sp->sign1 > 0) { + mpz_add(r, r, q); + } else { + mpz_sub(r, r, q); + } + + if (rand() % 2) { + sp->sign0 = 1; + mpz_add_ui(r, r, 1); + } else { + sp->sign0 = -1; + mpz_sub_ui(r, r, 1); + } + if (!mpz_probab_prime_p(r, 10)) continue; + for (i=0; i<10; i++) { + int bit; + //use q as a temp variable + mpz_set_ui(q, 0); + bit = qbits - rbits - 4 + 1; + if (bit < 3) bit = 3; + mpz_setbit(q, bit); + pbc_mpz_random(h, q); + mpz_mul_ui(h, h, 12); + //finally q takes the value it should + mpz_mul(q, h, r); + mpz_sub_ui(q, q, 1); + if (mpz_probab_prime_p(q, 10)) { + found = 1; + break; + } + } + } while (!found); +} + +// Type A1 pairings: + +struct a1_param_s { + mpz_t p; + mpz_t n; + int l; +}; +typedef struct a1_param_s a1_param_t[1]; +typedef struct a1_param_s *a1_param_ptr; + +struct a1_pairing_data_s { + field_t Fp, Fp2, Ep; +}; +typedef struct a1_pairing_data_s a1_pairing_data_t[1]; +typedef struct a1_pairing_data_s *a1_pairing_data_ptr; + +static void a1_clear(void *data) { + a1_param_ptr param = data; + mpz_clear(param->p); + mpz_clear(param->n); + pbc_free(data); +} + +static void a1_out_str(FILE *stream, void *data) { + a1_param_ptr p = data; + param_out_type(stream, "a1"); + param_out_mpz(stream, "p", p->p); + param_out_mpz(stream, "n", p->n); + param_out_int(stream, "l", p->l); +} + +struct pp2_coeff_s { + element_t cx2; + element_t cy2; + element_t cxy; + element_t cx; + element_t cy; + element_t c; +}; +typedef struct pp2_coeff_s pp2_coeff_t[1]; +typedef struct pp2_coeff_s *pp2_coeff_ptr; + +static void pp2_coeff_set(pp2_coeff_ptr p, + element_t cx2, element_t cy2, element_t cxy, + element_t cx, element_t cy, element_t c) { + element_init(p->cx2, cx2->field); + element_init(p->cy2, cy2->field); + element_init(p->cxy, cxy->field); + element_init(p->cx, cx->field); + element_init(p->cy, cy->field); + element_init(p->c, c->field); + element_set(p->cx2, cx2); + element_set(p->cy2, cy2); + element_set(p->cxy, cxy); + element_set(p->cx, cx); + element_set(p->cy, cy); + element_set(p->c, c); +} + +static void a1_pairing_pp_clear(pairing_pp_t p) { + void **pp = p->data; + while (*pp) { + pbc_free(*pp); + pp++; + } + pbc_free(p->data); +} + +static void a1_pairing_pp_init(pairing_pp_t p, element_ptr in1, pairing_t pairing) { + int m; + element_ptr Px = curve_x_coord(in1); + element_ptr Py = curve_y_coord(in1); + a1_pairing_data_ptr a1info = pairing->data; + p->data = pbc_malloc(sizeof(void *) * mpz_sizeinbase(pairing->r, 2)); + void **pp = p->data; + element_t V; + element_t a, b, c; + element_t a2, b2, c2; + element_t e0, e1, e2; + element_ptr Vx, Vy; + + #define do_tangent() compute_abc_tangent(a, b, c, Vx, Vy, e0); + + #define do_line() compute_abc_line(a2, b2, c2, Vx, Vy, Px, Py, e0); + + element_init(V, a1info->Ep); + element_set(V, in1); + Vx = curve_x_coord(V); + Vy = curve_y_coord(V); + + element_init(a, a1info->Fp); + element_init(b, a1info->Fp); + element_init(c, a1info->Fp); + element_init(e0, a1info->Fp); + element_init(e1, a1info->Fp); + element_init(e2, a1info->Fp); + element_init(a2, a1info->Fp); + element_init(b2, a1info->Fp); + element_init(c2, a1info->Fp); + + m = mpz_sizeinbase(pairing->r, 2) - 2; + + for(;;) { + do_tangent(); + if (!m) break; + element_double(V, V); + + if (mpz_tstbit(pairing->r, m)) { + do_line(); + element_add(V, V, in1); + //preprocess two at once + //e0 = coeff of x + element_mul(e0, a, c2); + element_mul(e1, a2, c); + element_add(e0, e0, e1); + + //e1 = coeff of y + element_mul(e1, b2, c); + element_mul(e2, b, c2); + element_add(e1, e1, e2); + + //c = constant term + element_mul(c, c, c2); + + //c2 = coeff of xy + element_mul(c2, a, b2); + element_mul(e2, a2, b); + element_add(c2, c2, e2); + + //a = coeff of x^2 + element_mul(a, a, a2); + + //b = coeff of y^2 + element_mul(b, b, b2); + + *pp = pbc_malloc(sizeof(pp2_coeff_t)); + pp2_coeff_set(*pp, a, b, c2, e0, e1, c); + } else { + *pp = pbc_malloc(sizeof(pp_coeff_t)); + pp_coeff_set(*pp, a, b, c); + } + pp++; + m--; + } + *pp = pbc_malloc(sizeof(pp_coeff_t)); + pp_coeff_set(*pp, a, b, c); + pp++; + *pp = NULL; + + element_clear(a2); + element_clear(b2); + element_clear(c2); + element_clear(e2); + element_clear(e1); + element_clear(e0); + element_clear(a); + element_clear(b); + element_clear(c); + element_clear(V); + #undef do_tangent + #undef do_line +} + +static void a1_pairing_pp_apply(element_ptr out, element_ptr in2, pairing_pp_t p) { + void **pp = p->data; + a1_pairing_data_ptr a1info = p->pairing->data; + element_t f, f0; + element_t e0, e1; + int m; + element_ptr Qx = curve_x_coord(in2); + element_ptr Qy = curve_y_coord(in2); + element_t Qx2, Qy2, Qxy; + + #define do_tangent() \ + pp_coeff_ptr ppp = *pp; \ + a_miller_evalfn(f0, ppp->a, ppp->b, ppp->c, Qx, Qy); + + #define do_line() { \ + pp2_coeff_ptr ppp = *pp; \ + /*we'll map Q via (x,y) --> (-x, iy) */ \ + /*hence Qx^2 = x^2, Qy^2 = -y^2, Qx Qy = -ixy */\ + /*where x = Q'x, y = Q'y */ \ + \ + /* Re = cx2 x^2 - cy2 y^2 - cx x + c */ \ + /* Im = -cxy xy + cy y */ \ + element_mul(e0, ppp->cx2, Qx2); \ + element_mul(e1, ppp->cy2, Qy2); \ + element_sub(e0, e0, e1); \ + element_mul(e1, ppp->cx, Qx); \ + element_sub(e0, e0, e1); \ + element_add(element_x(f0), e0, ppp->c); \ + \ + element_mul(e0, ppp->cy, Qy); \ + element_mul(e1, ppp->cxy, Qxy); \ + element_sub(element_y(f0), e0, e1); \ + } + + element_init(f, out->field); + element_init(f0, out->field); + + element_set1(f); + + element_init(e0, a1info->Fp); + element_init(e1, a1info->Fp); + element_init(Qx2, a1info->Fp); + element_init(Qy2, a1info->Fp); + element_init(Qxy, a1info->Fp); + + element_square(Qx2, Qx); + element_square(Qy2, Qy); + element_mul(Qxy, Qx, Qy); + + m = mpz_sizeinbase(p->pairing->r, 2) - 2; + + while (m > 0) { + if (mpz_tstbit(p->pairing->r, m)) { + do_line(); + } else { + do_tangent(); + } + element_mul(f, f, f0); + pp++; + m--; + element_square(f, f); + } + do_tangent(); + element_mul(f, f, f0); + + //Tate exponentiation + //simpler but slower: + //element_pow_mpz(out, f, p->tateexp); + //use this trick instead: + element_invert(f0, f); + element_neg(element_y(f), element_y(f)); + element_mul(f, f, f0); + element_pow_mpz(out, f, p->pairing->phikonr); + + /* We could use this instead but p->h is small so this does not help much + a_tateexp(out, f, f0, p->h); + */ + + element_clear(Qx2); + element_clear(Qy2); + element_clear(Qxy); + element_clear(f); + element_clear(f0); + element_clear(e1); + element_clear(e0); + #undef do_tangent + #undef do_line +} + +// e0 is a temp var. +// Mixed coordinates. +static void compute_abc_line_proj(element_ptr a, element_ptr b, element_ptr c, + element_ptr Vx, element_ptr Vy, element_ptr z, element_ptr z2, + element_ptr V1x, element_ptr V1y, element_ptr e0) { + //temporally used to store Z1^3 + element_mul(c,z,z2); + //a = Y1-Y2*Z1^3 + element_mul(e0,V1y,c); + element_sub(a,Vy,e0); + //b = -(X1*Z1-X2*Z1^3) + element_mul(b,c,V1x); + element_mul(e0,Vx,z); + element_sub(b,b,e0); + //c = -(Y2*b+X2*a) + element_mul(c,b,V1y); + element_mul(e0,a,V1x); + element_add(c,c,e0); + element_neg(c,c); +} + +// in1, in2 are from E(F_q), out from F_q^2 +static void a1_pairing_proj(element_ptr out, element_ptr in1, element_ptr in2, + pairing_t pairing) { + a1_pairing_data_ptr p = pairing->data; + element_t V; + element_t z, z2; + element_t f, f0; + element_t a, b, c; + element_t e0; + const element_ptr e1 = a, e2 = b, e3 = c; // used in point_to_affine() etc. + int m; + element_ptr Px = curve_x_coord(in1); + element_ptr Py = curve_y_coord(in1); + element_ptr Qx = curve_x_coord(in2); + element_ptr Qy = curve_y_coord(in2); + element_ptr Vx; + element_ptr Vy; + + #define point_to_affine() \ + element_invert(z, z); \ + element_square(e0, z); \ + element_mul(Vx, Vx, e0); \ + element_mul(e0, e0, z); \ + element_mul(Vy, Vy, e0); \ + element_set1(z); \ + element_set1(z2); + + //TODO: do I need to check if V=-in1? + //Where V=(Vx,Vy,z) and in1=(Px,Py,1), a mixed coordinates. + #define proj_add() { \ + /* H=X2*Z1^2-X1 */ \ + element_mul(e0,Px,z2); \ + element_sub(e0,e0,Vx); \ + /* H^2 */ \ + element_square(e1,e0); \ + /* r=Y2*Z1^3-Y1 */ \ + element_mul(e2,z,z2); \ + element_mul(e2,e2,Py); \ + element_sub(e2,e2,Vy); \ + \ + /* X3=r^2-H^3-2X1*H^2 */ \ + element_set(z2,Vx); /* use z2 to store X1 and update Vx=X3 */ \ + element_square(Vx,e2); \ + element_mul(e3,e0,e1); /* e3=H^3 */ \ + element_sub(Vx,Vx,e3); \ + element_double(e3,z2); \ + element_mul(e3,e3,e1); /* 2X1*H^2 */ \ + element_sub(Vx,Vx,e3); \ + /* Y3=r(X1*H^2-X3)-Y1*H^3 */ \ + element_mul(e3,z2,e1); \ + element_sub(e3,e3,Vx); \ + element_mul(e3,e3,e2); \ + element_mul(e2,e0,e1); /* e2 no longer used. */ \ + element_mul(e2,e2,Vy); \ + element_sub(Vy,e3,e2); \ + /* Z3=Z1*H */ \ + element_mul(z,z,e0); \ + element_square(z2,z); \ + } + + #define proj_double() { \ + /* e0 = 3x^2 + (cc->a) z^4 */ \ + /* for this case a = 1 */ \ + element_square(e0, Vx); \ + /* element_mul_si(e0, e0, 3); */ \ + element_double(e1, e0); \ + element_add(e0, e1, e0); \ + element_square(e1, z2); \ + element_add(e0, e0, e1); \ + \ + /* z_out = 2 y z */ \ + element_mul(z, Vy, z); \ + /* element_mul_si(z, z, 2); */ \ + element_double(z, z); \ + element_square(z2, z); \ + \ + /* e1 = 4 x y^2 */ \ + element_square(e2, Vy); \ + element_mul(e1, Vx, e2); \ + /* element_mul_si(e1, e1, 4); */ \ + element_double(e1, e1); \ + element_double(e1, e1); \ + \ + /* x_out = e0^2 - 2 e1 */ \ + element_double(e3, e1); \ + element_square(Vx, e0); \ + element_sub(Vx, Vx, e3); \ + \ + /* e2 = 8y^4 */ \ + element_square(e2, e2); \ + /* element_mul_si(e2, e2, 8); */ \ + element_double(e2, e2); \ + element_double(e2, e2); \ + element_double(e2, e2); \ + \ + /* y_out = e0(e1 - x_out) - e2 */ \ + element_sub(e1, e1, Vx); \ + element_mul(e0, e0, e1); \ + element_sub(Vy, e0, e2); \ + } + + #define do_tangent() { \ + compute_abc_tangent_proj(a, b, c, Vx, Vy, z, z2, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); \ + } + + #define do_line() { \ + compute_abc_line_proj(a, b, c, Vx, Vy, z, z2, Px, Py, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); \ + } + + element_init(V, p->Ep); + element_set(V, in1); + Vx = curve_x_coord(V); + Vy = curve_y_coord(V); + + element_init(f, p->Fp2); + element_init(f0, p->Fp2); + element_set1(f); + element_init(a, p->Fp); + element_init(b, p->Fp); + element_init(c, p->Fp); + element_init(e0, p->Fp); + element_init(z, p->Fp); + element_init(z2, p->Fp); + element_set1(z); + element_set1(z2); + + m = mpz_sizeinbase(pairing->r, 2) - 2; + //TODO: sliding NAF + for(;;) { + do_tangent(); + if (!m) break; + + proj_double(); //V=2V + if (mpz_tstbit(pairing->r, m)) { + // point_to_affine(); + do_line(); + proj_add(); //V=V+in1 + } + + m--; + element_square(f, f); + } + + // Tate exponentiation. + // Simpler but slower: + // element_pow_mpz(out, f, p->tateexp); + // Use this trick instead: + element_invert(f0, f); + element_neg(element_y(f), element_y(f)); + element_mul(f, f, f0); + element_pow_mpz(out, f, pairing->phikonr); + + /* We could use this instead but p->h is small so this does not help much + a_tateexp(out, f, f0, p->h); + */ + + element_clear(f); + element_clear(f0); + element_clear(z); + element_clear(z2); + element_clear(V); + element_clear(a); + element_clear(b); + element_clear(c); + element_clear(e0); + #undef point_to_affine + #undef proj_add + #undef proj_double + #undef do_tangent + #undef do_line +} + +//in1, in2 are from E(F_q), out from F_q^2 +static void a1_pairing(element_ptr out, element_ptr in1, element_ptr in2, + pairing_t pairing) { + a1_pairing_data_ptr p = pairing->data; + element_t V; + element_t f, f0; + element_t a, b, c; + element_t e0; + int m; + element_ptr Px = curve_x_coord(in1); + element_ptr Py = curve_y_coord(in1); + element_ptr Qx = curve_x_coord(in2); + element_ptr Qy = curve_y_coord(in2); + element_ptr Vx; + element_ptr Vy; + + #define do_tangent() { \ + compute_abc_tangent(a, b, c, Vx, Vy, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); \ + } + + #define do_line() { \ + compute_abc_line(a, b, c, Vx, Vy, Px, Py, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); \ + } + + element_init(V, p->Ep); + element_set(V, in1); + Vx = curve_x_coord(V); + Vy = curve_y_coord(V); + + element_init(f, p->Fp2); + element_init(f0, p->Fp2); + element_set1(f); + element_init(a, p->Fp); + element_init(b, p->Fp); + element_init(c, p->Fp); + element_init(e0, p->Fp); + + m = mpz_sizeinbase(pairing->r, 2) - 2; + + //TODO: sliding NAF + for(;;) { + do_tangent(); + if (!m) break; + + element_double(V, V); + if (mpz_tstbit(pairing->r, m)) { + do_line(); + element_add(V, V, in1); + } + + m--; + element_square(f, f); + } + + // Tate exponentiation. + // Simpler but slower: + // element_pow_mpz(out, f, p->tateexp); + // Use this trick instead: + element_invert(f0, f); + element_neg(element_y(f), element_y(f)); + element_mul(f, f, f0); + element_pow_mpz(out, f, pairing->phikonr); + + /* We could use this instead but p->h is small so this does not help much + a_tateexp(out, f, f0, p->h); + */ + + element_clear(f); + element_clear(f0); + element_clear(V); + element_clear(a); + element_clear(b); + element_clear(c); + element_clear(e0); + #undef do_tangent + #undef do_line +} + +//in1, in2 are from E(F_q), out from F_q^2 +void a1_pairings_affine(element_ptr out, element_t in1[], element_t in2[], + int n_prod, pairing_t pairing) { + a1_pairing_data_ptr p = pairing->data; + element_t* V = pbc_malloc(sizeof(element_t)*n_prod); + element_t f, f0; + element_t a, b, c; + element_t e0; + int m, i; + element_ptr Px, Py; + element_ptr Qx, Qy; + element_ptr Vx, Vy; + + #define do_tangents() { \ + for(i=0; i<n_prod; i++){ \ + Vx = curve_x_coord(V[i]); \ + Vy = curve_y_coord(V[i]); \ + Qx = curve_x_coord(in2[i]); \ + Qy = curve_y_coord(in2[i]); \ + compute_abc_tangent(a, b, c, Vx, Vy, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); \ + } \ + } + + #define do_lines() { \ + for(i=0; i<n_prod; i++){ \ + Vx = curve_x_coord(V[i]); \ + Vy = curve_y_coord(V[i]); \ + Px = curve_x_coord(in1[i]); \ + Py = curve_y_coord(in1[i]); \ + Qx = curve_x_coord(in2[i]); \ + Qy = curve_y_coord(in2[i]); \ + compute_abc_line(a, b, c, Vx, Vy, Px, Py, e0); \ + a_miller_evalfn(f0, a, b, c, Qx, Qy); \ + element_mul(f, f, f0); \ + } \ + } + + for(i=0; i<n_prod; i++){ + element_init(V[i], p->Ep); + element_set(V[i], in1[i]); + } + element_init(f, p->Fp2); + element_init(f0, p->Fp2); + element_set1(f); + element_init(a, p->Fp); + element_init(b, p->Fp); + element_init(c, p->Fp); + element_init(e0, p->Fp); + + m = mpz_sizeinbase(pairing->r, 2) - 2; + + //TODO: sliding NAF + for(;;) { + do_tangents(); + if (!m) break; + element_multi_double(V, V, n_prod); + if (mpz_tstbit(pairing->r, m)) { + do_lines(); + element_multi_add(V, V, in1, n_prod); + } + + m--; + element_square(f, f); + } + + // Tate exponentiation. + // Simpler but slower: + // element_pow_mpz(out, f, p->tateexp); + // Use this trick instead: + element_invert(f0, f); + element_neg(element_y(f), element_y(f)); + element_mul(f, f, f0); + element_pow_mpz(out, f, pairing->phikonr); + + /* We could use this instead but p->h is small so this does not help much + a_tateexp(out, f, f0, p->h); + */ + + element_clear(f); + element_clear(f0); + for(i=0; i<n_prod; i++){ + element_clear(V[i]); + } + pbc_free(V); + element_clear(a); + element_clear(b); + element_clear(c); + element_clear(e0); + #undef do_tangents + #undef do_lines +} + +static void a1_pairing_clear(pairing_t pairing) { + field_clear(pairing->GT); + + a1_pairing_data_ptr p = pairing->data; + field_clear(p->Ep); + field_clear(p->Fp2); + field_clear(p->Fp); + pbc_free(p); + + mpz_clear(pairing->phikonr); + mpz_clear(pairing->r); + field_clear(pairing->Zr); +} + +static void a1_pairing_option_set(pairing_t pairing, char *key, char *value) { + if (!strcmp(key, "method")) { + if (!strcmp(value, "miller")) { + pairing->map = a1_pairing_proj; + pairing->pp_init = a1_pairing_pp_init; + pairing->pp_clear = a1_pairing_pp_clear; + pairing->pp_apply = a1_pairing_pp_apply; + } else if (!strcmp(value, "miller-affine")){ + pairing->map = a1_pairing; + pairing->pp_init = a1_pairing_pp_init; + pairing->pp_clear = a1_pairing_pp_clear; + pairing->pp_apply = a1_pairing_pp_apply; + } else if (!strcmp(value, "shipsey-stange")) { + pairing->map = a_pairing_ellnet; + pairing->pp_init = a_pairing_ellnet_pp_init; + pairing->pp_clear = a_pairing_ellnet_pp_clear; + pairing->pp_apply = a_pairing_ellnet_pp_apply; + } + } +} + +static void a1_init_pairing(pairing_t pairing, void *data) { + a1_param_ptr param = data; + element_t a, b; + mpz_init(pairing->r); + mpz_set(pairing->r, param->n); + field_init_fp(pairing->Zr, pairing->r); + + a1_pairing_data_ptr p; + + p = pairing->data = pbc_malloc(sizeof(a1_pairing_data_t)); + + //k=2, hence phi_k(q) = q + 1, phikonr = (q+1)/r + mpz_init(pairing->phikonr); + mpz_set_ui(pairing->phikonr, param->l); + + field_init_fp(p->Fp, param->p); + element_init(a, p->Fp); + element_init(b, p->Fp); + element_set1(a); + element_set0(b); + field_init_curve_ab(p->Ep, a, b, pairing->r, pairing->phikonr); + + // Turns out to be faster. + field_curve_use_random_solvefory(p->Ep); + + element_clear(a); + element_clear(b); + field_init_fi(p->Fp2, p->Fp); + + pairing->finalpow = a_finalpow; + pairing->G1 = pbc_malloc(sizeof(field_t)); + pairing->G2 = pairing->G1 = p->Ep; + pairing_GT_init(pairing, p->Fp2); + + pairing->map = a1_pairing_proj; //default uses projective coordinates. + pairing->phi = phi_identity; + pairing->prod_pairings = a1_pairings_affine; + + pairing->clear_func = a1_pairing_clear; + + pairing->pp_init = a1_pairing_pp_init; + pairing->pp_clear = a1_pairing_pp_clear; + pairing->pp_apply = a1_pairing_pp_apply; + pairing->option_set = a1_pairing_option_set; +} + +static void a1_init(pbc_param_t p) { + static pbc_param_interface_t interface = {{ + a1_clear, + a1_init_pairing, + a1_out_str, + }}; + p->api = interface; + a1_param_ptr param = p->data = pbc_malloc(sizeof(*param)); + mpz_init(param->p); + mpz_init(param->n); +} + +// Public interface: + +int pbc_param_init_a1(pbc_param_ptr par, struct symtab_s *tab) { + a1_init(par); + a1_param_ptr p = par->data; + + int err = 0; + err += lookup_mpz(p->p, tab, "p"); + err += lookup_mpz(p->n, tab, "n"); + err += lookup_int(&p->l, tab, "l"); + return err; +} + +void pbc_param_init_a1_gen(pbc_param_ptr par, mpz_t order) { + a1_init(par); + a1_param_ptr param = par->data; + // If order is even, ideally check all even l, not just multiples of 4 + // but I don't see a good reason for having an even order. + unsigned int l = 4; + mpz_t n; + mpz_ptr p = param->p; + mpz_init(n); + mpz_mul_ui(n, order, 4); + mpz_sub_ui(p, n, 1); + for (;;) { + if (mpz_probab_prime_p(p, 20)) { + break; + } + mpz_add(p, p, n); + l += 4; + } + param->l = l; + mpz_set(param->n, order); + mpz_clear(n); +} |