a_param.c

这是一个C的源代码

#include <assert.h>
#include <stdio.h>
#include <stdlib.h> //for rand, pbc_malloc, pbc_free
#include <string.h> //for strcmp
#include <gmp.h>
#include "pbc_symtab.h"
#include "pbc_fops.h"
#include "pbc_field.h"
#include "pbc_fp.h"
#include "pbc_fieldquadratic.h"
#include "pbc_pairing.h"
#include "pbc_a_param.h"
#include "pbc_a1_param.h"
#include "pbc_curve.h"
#include "pbc_param.h"
#include "pbc_random.h"
#include "pbc_tracker.h"
#include "pbc_memory.h"
#include "pbc_utils.h"

struct a_pairing_data_s {
    field_t Fq, Fq2, Eq;
    int exp2, exp1;
    int sign1;
};
typedef struct a_pairing_data_s a_pairing_data_t[1];
typedef struct a_pairing_data_s *a_pairing_data_ptr;

void a_param_init(a_param_t sp)
{
    mpz_init(sp->r);
    mpz_init(sp->q);
    mpz_init(sp->h);
}

void a_param_clear(a_param_t sp)
{
    mpz_clear(sp->r);
    mpz_clear(sp->q);
    mpz_clear(sp->h);
}

void a_param_gen(a_param_t sp, int rbits, int qbits)
{
    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);
}

void a_param_out_str(FILE *stream, a_param_ptr p)
{
    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);
}

void a_param_inp_generic (a_param_ptr p, fetch_ops_t fops, void *ctx)
{
    assert (fops);
    assert (ctx);
    symtab_t tab;

    symtab_init(tab);
    param_read_generic (tab, fops, ctx);

    lookup_mpz(p->q, tab, "q");
    lookup_mpz(p->r, tab, "r");
    lookup_mpz(p->h, tab, "h");
    p->exp2 = lookup_int(tab, "exp2");
    p->exp1 = lookup_int(tab, "exp1");
    p->sign1 = lookup_int(tab, "sign1");
    p->sign0 = lookup_int(tab, "sign0");

    param_clear_tab(tab);
    symtab_clear(tab);
}

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;

    void do_tangent(void) {
	compute_abc_tangent(a, b, c, Vx, Vy, e0);
	pp_coeff_set(coeff[i], a, b, c);
    }

    void do_line(void) {
	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);
}

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);
}

static void lucas_odd(element_ptr out, element_ptr in, element_ptr temp, mpz_t cofactor)
//assumes cofactor is odd
//overwrites in and temp, out must not be in
//luckily this touchy routine is only used internally
{
    element_ptr in0 = fi_re(in);
    element_ptr in1 = fi_im(in);
    element_ptr v0 = fi_re(out);
    element_ptr v1 = fi_im(out);
    element_ptr t0 = fi_re(temp);
    element_ptr t1 = fi_im(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 = fi_im(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(fi_im(out), a, Qx);
    element_sub(fi_re(out), c, fi_im(out));
    element_mul(fi_im(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);
}

static void a_pairing_ellnet(element_ptr out, element_ptr in1, element_ptr in2,
	pairing_t pairing)
//in1, in2 are from E(F_q), out from F_q^2
//uses elliptic nets (see Stange)
{
    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);