g_param.c

这是一个C的源代码

	    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^2 - q + 1) / r is odd
    //thus v1 = V_k, v0 = V_{k-1}
    //     U = (P v1 - 2 v0) / (P^2 - 4)

    element_double(v0, v0);
    element_mul(in0, t1, v1);
    element_sub(in0, in0, v0);

    element_square(t1, t1);
    element_sub(t1, t1, t0);
    element_sub(t1, t1, t0);

    element_halve(v0, v1);
    element_div(v1, in0, t1);
    element_mul(v1, v1, in1);
    element_clear(temp);
}

static void tatepower10(element_ptr out, element_ptr in, pairing_t pairing)
{
    mnt_pairing_data_ptr p = pairing->data;
    element_t e0, e1, e2, e3;
    element_init(e0, p->Fqk);
    element_init(e1, p->Fqd);
    element_init(e2, p->Fqd);
    element_init(e3, p->Fqk);
    element_ptr e0re = fi_re(e0);
    element_ptr e0im = fi_im(e0);
    element_ptr e0re0 = ((element_t *) e0re->data)[0];
    element_ptr e0im0 = ((element_t *) e0im->data)[0];
    element_t *inre = fi_re(in)->data;
    element_t *inim = fi_im(in)->data;
    //see thesis
    void qpower(int sign) {
	polymod_const_mul(e2, inre[1], p->xpowq);
	element_set(e0re, e2);
	polymod_const_mul(e2, inre[2], p->xpowq2);
	element_add(e0re, e0re, e2);
	polymod_const_mul(e2, inre[3], p->xpowq3);
	element_add(e0re, e0re, e2);
	polymod_const_mul(e2, inre[4], p->xpowq4);
	element_add(e0re, e0re, e2);
	element_add(e0re0, e0re0, inre[0]);

	if (sign > 0) {
	    polymod_const_mul(e2, inim[1], p->xpowq);
	    element_set(e0im, e2);
	    polymod_const_mul(e2, inim[2], p->xpowq2);
	    element_add(e0im, e0im, e2);
	    polymod_const_mul(e2, inim[3], p->xpowq3);
	    element_add(e0im, e0im, e2);
	    polymod_const_mul(e2, inim[4], p->xpowq4);
	    element_add(e0im, e0im, e2);
	    element_add(e0im0, e0im0, inim[0]);
	} else {
	    polymod_const_mul(e2, inim[1], p->xpowq);
	    element_neg(e0im, e2);
	    polymod_const_mul(e2, inim[2], p->xpowq2);
	    element_sub(e0im, e0im, e2);
	    polymod_const_mul(e2, inim[3], p->xpowq3);
	    element_sub(e0im, e0im, e2);
	    polymod_const_mul(e2, inim[4], p->xpowq4);
	    element_sub(e0im, e0im, e2);
	    element_sub(e0im0, e0im0, inim[0]);
	}
    }
    qpower(1);
    element_set(e3, e0);
    element_set(e0re, fi_re(in));
    element_neg(e0im, fi_im(in));
    element_mul(e3, e3, e0);
    qpower(-1);
    element_mul(e0, e0, in);
    element_invert(e0, e0);
    element_mul(in, e3, e0);

    element_set(e0, in);
    lucas_even(out, e0, pairing->phikonr);

    element_clear(e0);
    element_clear(e1);
    element_clear(e2);
    element_clear(e3);
}

static void (*cc_miller_no_denom_fn)(element_t res, mpz_t q, element_t P,
	element_ptr Qx, element_ptr Qy);

static void cc_pairing(element_ptr out, element_ptr in1, element_ptr in2,
	pairing_t pairing)
{
    element_ptr Qbase = in2;
    element_t Qx, Qy;
    mnt_pairing_data_ptr p = pairing->data;

    element_init(Qx, p->Fqd);
    element_init(Qy, p->Fqd);
    //map from twist: (x, y) --> (v^-1 x, v^-(3/2) y)
    //where v is the quadratic nonresidue used to construct the twist
    element_mul(Qx, curve_x_coord(Qbase), p->nqrinv);
    //v^-3/2 = v^-2 * v^1/2
    element_mul(Qy, curve_y_coord(Qbase), p->nqrinv2);
    cc_miller_no_denom_fn(out, pairing->r, in1, Qx, Qy);
    tatepower10(out, out, pairing);
    element_clear(Qx);
    element_clear(Qy);
}

static int cc_is_almost_coddh(element_ptr a, element_ptr b,
	element_ptr c, element_ptr d,
	pairing_t pairing)
{
    int res = 0;
    element_t t0, t1, t2;
    element_t cx, cy;
    element_t dx, dy;
    mnt_pairing_data_ptr p = pairing->data;

    element_init(cx, p->Fqd);
    element_init(cy, p->Fqd);
    element_init(dx, p->Fqd);
    element_init(dy, p->Fqd);

    element_init(t0, pairing->GT);
    element_init(t1, pairing->GT);
    element_init(t2, pairing->GT);
    //map from twist: (x, y) --> (v^-1 x, v^-(3/2) y)
    //where v is the quadratic nonresidue used to construct the twist
    element_mul(cx, curve_x_coord(c), p->nqrinv);
    element_mul(dx, curve_x_coord(d), p->nqrinv);
    //v^-3/2 = v^-2 * v^1/2
    element_mul(cy, curve_y_coord(c), p->nqrinv2);
    element_mul(dy, curve_y_coord(d), p->nqrinv2);

    cc_miller_no_denom_fn(t0, pairing->r, a->data, dx, dy);
    cc_miller_no_denom_fn(t1, pairing->r, b->data, cx, cy);
    tatepower10(t0, t0, pairing);
    tatepower10(t1, t1, pairing);
    element_mul(t2, t0, t1);
    if (element_is1(t2)) {
	//g, g^x, h, h^-x case
	res = 1;
    } else {
	element_invert(t1, t1);
	element_mul(t2, t0, t1);
	if (element_is1(t2)) {
	    //g, g^x, h, h^x case
	    res = 1;
	}
    }
    element_clear(cx);
    element_clear(cy);
    element_clear(dx);
    element_clear(dy);
    element_clear(t0);
    element_clear(t1);
    element_clear(t2);
    return res;
}

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 g_pairing_pp_init(pairing_pp_t p, element_ptr in1, pairing_t pairing)
{
    element_ptr P = in1;
    const element_ptr Px = curve_x_coord(P);
    const element_ptr Py = curve_y_coord(P);
    element_t Z;
    int m;
    mnt_pairing_data_ptr info = pairing->data;
    element_t t0;
    element_t a, b, c;
    field_ptr Fq = info->Fq;
    pp_coeff_t *coeff;
    mpz_ptr q = pairing->r;
    pp_coeff_ptr pp;
    const element_ptr cca = curve_a_coeff(P);
    element_ptr Zx;
    element_ptr Zy;

    void store_abc(void)
    {
	element_init(pp->a, Fq);
	element_init(pp->b, Fq);
	element_init(pp->c, Fq);
	element_set(pp->a, a);
	element_set(pp->b, b);
	element_set(pp->c, c);
	pp++;
    }

    void do_tangent(void)
    {
	//a = -slope_tangent(Z.x, Z.y);
	//b = 1;
	//c = -(Z.y + a * Z.x);
	//but we multiply by 2*Z.y to avoid division

	//a = -Zx * (3 Zx + twicea_2) - a_4;
	//Common curves: a2 = 0 (and cc->a is a_4), so
	//a = -(3 Zx^2 + cc->a)
	//b = 2 * Zy
	//c = -(2 Zy^2 + a Zx);

	element_square(a, Zx);
	element_double(t0, a);
	element_add(a, a, t0);
	element_add(a, a, cca);
	element_neg(a, a);

	element_add(b, Zy, Zy);

	element_mul(t0, b, Zy);
	element_mul(c, a, Zx);
	element_add(c, c, t0);
	element_neg(c, c);

	store_abc();
    }

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

	element_sub(b, Px, Zx);
	element_sub(a, Zy, Py);
	element_mul(t0, b, Zy);
	element_mul(c, a, Zx);
	element_add(c, c, t0);
	element_neg(c, c);

	store_abc();
    }

    element_init(Z, P->field);
    element_set(Z, P);
    Zx = curve_x_coord(Z);
    Zy = curve_y_coord(Z);

    element_init(t0, Fq);
    element_init(a, Fq);
    element_init(b, Fq);
    element_init(c, Fq);

    m = mpz_sizeinbase(q, 2) - 2;
    p->data = pbc_malloc(sizeof(pp_coeff_t) * 2 * m);
    coeff = (pp_coeff_t *) p->data;
    pp = coeff[0];

    for(;;) {
	do_tangent();

	if (!m) break;

	element_double(Z, Z);
	if (mpz_tstbit(q, m)) {
	    do_line();
	    element_add(Z, Z, P);
	}
	m--;
    }

    element_clear(t0);
    element_clear(a);
    element_clear(b);
    element_clear(c);
    element_clear(Z);
}

static void g_pairing_pp_clear(pairing_pp_t p)
{
    //TODO: better to store a sentinel value in p->data?
    mpz_ptr q = p->pairing->r;
    int m = mpz_sizeinbase(q, 2) + mpz_popcount(q) - 3;
    int i;
    pp_coeff_t *coeff = (pp_coeff_t *) p->data;
    pp_coeff_ptr pp;
    for (i=0; i<m; i++) {
	pp = coeff[i];
	element_clear(pp->a);
	element_clear(pp->b);
	element_clear(pp->c);
    }
    pbc_free(p->data);
}

static void g_pairing_pp_apply(element_ptr out, element_ptr in2, pairing_pp_t p)
{
    mpz_ptr q = p->pairing->r;
    mnt_pairing_data_ptr info = p->pairing->data;
    int m = mpz_sizeinbase(q, 2) - 2;
    pp_coeff_t *coeff = (pp_coeff_t *) p->data;
    pp_coeff_ptr pp = coeff[0];
    element_ptr Qbase = in2;
    element_t e0;
    element_t Qx, Qy;
    element_t v;
    element_init_GT(e0, p->pairing);
    element_init_GT(v, p->pairing);
    element_init(Qx, info->Fqd);
    element_init(Qy, info->Fqd);

    //map from twist: (x, y) --> (v^-1 x, v^-(3/2) y)
    //where v is the quadratic nonresidue used to construct the twist
    element_mul(Qx, curve_x_coord(Qbase), info->nqrinv);
    //v^-3/2 = v^-2 * v^1/2
    element_mul(Qy, curve_y_coord(Qbase), info->nqrinv2);

    element_set1(out);
    for(;;) {
	d_miller_evalfn(e0, pp->a, pp->b, pp->c, Qx, Qy);
	element_mul(out, out, e0);
	pp++;

	if (!m) break;

	if (mpz_tstbit(q, m)) {
	    d_miller_evalfn(e0, pp->a, pp->b, pp->c, Qx, Qy);
	    element_mul(out, out, e0);
	    pp++;
	}
	m--;
	element_square(out, out);
    }
    tatepower10(out, out, p->pairing);

    element_clear(e0);
    element_clear(Qx);
    element_clear(Qy);
    element_clear(v);
}

static void g_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)
{
    mnt_pairing_data_ptr p = pairing->data;

    const element_ptr a = curve_a_coeff(in1);
    const element_ptr b = curve_b_coeff(in1);

    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, out);
    element_init_same_as(B, out);

    // c1 = 2y
    // cm3 = -2y
    element_double(c1, y);
    element_neg(cm3, c1);

    //use c0, cm1, cm2, C, c4 as temp variables for now
    //compute c3, c2
    element_square(cm2, x);
    element_square(C, cm2);
    element_mul(cm1, b, x);
    element_double(cm1, cm1);
    element_square(c4, a);

    element_mul(c2, cm1, cm2);
    element_double(c2, c2);
    element_mul(c0, a, C);
    element_add(c2, c2, c0);
    element_mul(c0, c4, cm2);
    element_sub(c2, c2, c0);
    element_double(c0, c2);
    element_double(c0, c0);
    element_add(c2, c2, c0);

    element_mul(c0, cm1, a);
    element_square(c3, b);
    element_double(c3, c3);
    element_double(c3, c3);
    element_add(c0, c0, c3);
    element_double(c0, c0);
    element_mul(c3, a, c4);
    element_add(c0, c0, c3);
    element_sub(c2, c2, c0);
    element_mul(c0, cm2, C);
    element_add(c3, c0, c2);
    element_mul(c3, c3, c1);
    element_double(c3, c3);

    element_mul(c0, a, cm2);
    element_add(c0, c0, cm1);
    element_double(c0, c0);
    element_add(c0, c0, C);
    element_double(c2, c0);
    element_add(c0, c0, c2);
    element_sub(c2, c0, c4);

    // c0 = 1
    // cm2 = -1
    element_set1(c0);
    element_neg(cm2, c0);

    // 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

    element_mul(fi_re(d0), x2, p->nqrinv);
    element_neg(A, d0);
    element_add(polymod_coeff(fi_re(A), 0), polymod_coeff(fi_re(A), 0), x);

    element_double(C, x);
    element_add(polymod_coeff(fi_re(d0), 0), polymod_coeff(fi_re(d0), 0), C);

    element_square(dm1, A);
    element_mul(dm1, d0, dm1);

    element_mul(fi_im(d1), y2, p->nqrinv2);
    element_set(polymod_coeff(fi_re(d1), 0), y);

    element_square(d1, d1);
    element_sub(d1, dm1, d1);
    element_invert(B, d1);

    element_invert(A, A);

    element_mul(fi_im(d1), y2, p->nqrinv2);
    element_set0(fi_re(d1));
    element_neg(polymod_coeff(fi_re(d1), 0), y);
    element_mul(d1, d1, A);
    element_square(d1, d1);
    element_sub(d1, d0, d1);