aboutsummaryrefslogtreecommitdiffstats
path: root/moon-abe/pbc-0.5.14/ecc/a_param.c
diff options
context:
space:
mode:
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.c2315
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);
+}