d_param.c

这是一个C的源代码

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>
#include "pbc_fops.h"
#include "pbc_symtab.h"
#include "pbc_darray.h"
#include "pbc_field.h"
#include "pbc_poly.h"
#include "pbc_fp.h"
#include "pbc_fieldquadratic.h"
#include "pbc_hilbert.h"
#include "pbc_mnt.h"
#include "pbc_curve.h"
#include "pbc_pairing.h"
#include "pbc_d_param.h"
#include "pbc_param.h"
#include "pbc_tracker.h"
#include "pbc_memory.h"
#include "pbc_utils.h"

struct mnt_pairing_data_s {
    field_t Fq, Fqx, Fqd, Fqk;
    field_t Eq, Etwist;
    element_t nqrinv, nqrinv2;
    mpz_t tateexp;
    int k;
    element_t xpowq, xpowq2;
};
typedef struct mnt_pairing_data_s mnt_pairing_data_t[1];
typedef struct mnt_pairing_data_s *mnt_pairing_data_ptr;

void d_param_init(d_param_ptr param)
{
    mpz_init(param->q);
    mpz_init(param->n);
    mpz_init(param->h);
    mpz_init(param->r);
    mpz_init(param->a);
    mpz_init(param->b);
    mpz_init(param->nk);
    mpz_init(param->hk);
    param->k = 0;
    param->coeff = NULL;
    mpz_init(param->nqr);
}

void d_param_clear(d_param_ptr param)
{
    int d = param->k / 2;
    int i;
    mpz_clear(param->q);
    mpz_clear(param->n);
    mpz_clear(param->h);
    mpz_clear(param->r);
    mpz_clear(param->a);
    mpz_clear(param->b);
    mpz_clear(param->nk);
    mpz_clear(param->hk);
    mpz_clear(param->nqr);
    for (i=0; i<d; i++) {
	mpz_clear(param->coeff[i]);
    }
    pbc_free(param->coeff);
}

void d_param_out_str(FILE *stream, d_param_ptr p)
{
    int d = p->k / 2;
    int i;
    char s[80];
    param_out_type(stream, "d");
    param_out_mpz(stream, "q", p->q);
    param_out_mpz(stream, "n", p->n);
    param_out_mpz(stream, "h", p->h);
    param_out_mpz(stream, "r", p->r);
    param_out_mpz(stream, "a", p->a);
    param_out_mpz(stream, "b", p->b);
    param_out_int(stream, "k", p->k);
    param_out_mpz(stream, "nk", p->nk);
    param_out_mpz(stream, "hk", p->hk);
    for (i=0; i<d; i++) {
	sprintf(s, "coeff%d", i);
	param_out_mpz(stream, s, p->coeff[i]);
    }
    param_out_mpz(stream, "nqr", p->nqr);
}

void d_param_inp_generic (d_param_ptr p, fetch_ops_t fops, void *ctx)
{
    assert (fops);
    assert (ctx);
    symtab_t tab;
    char s[80];
    int i, d;

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

    lookup_mpz(p->q, tab, "q");
    lookup_mpz(p->n, tab, "n");
    lookup_mpz(p->h, tab, "h");
    lookup_mpz(p->r, tab, "r");
    lookup_mpz(p->a, tab, "a");
    lookup_mpz(p->b, tab, "b");
    p->k = lookup_int(tab, "k");
    lookup_mpz(p->nk, tab, "nk");
    lookup_mpz(p->hk, tab, "hk");
    lookup_mpz(p->nqr, tab, "nqr");

    d = p->k / 2;
    p->coeff = pbc_realloc(p->coeff, sizeof(mpz_t) * d);
    for (i=0; i<d; i++) {
	sprintf(s, "coeff%d", i);
	mpz_init(p->coeff[i]);
	lookup_mpz(p->coeff[i], tab, s);
    }

    param_clear_tab(tab);
    symtab_clear(tab);
}

static inline void d_miller_evalfn(element_t e0,
	element_t a, element_t b, element_t c,
	element_t Qx, element_t Qy)
{
    //a, b, c are in Fq
    //point Q is (Qx, Qy * sqrt(nqr)) where nqr is used to construct
    //the quadratic field extension Fqk of Fqd
    element_ptr re_out = fi_re(e0);
    element_ptr im_out = fi_im(e0);

    int i;
    int d = polymod_field_degree(re_out->field);
    for (i=0; i<d; i++) {
	element_mul(polymod_coeff(re_out, i), polymod_coeff(Qx, i), a);
	element_mul(polymod_coeff(im_out, i), polymod_coeff(Qy, i), b);
    }
    element_add(polymod_coeff(re_out, 0), polymod_coeff(re_out, 0), c);
}

static void cc_miller_no_denom_proj(element_t res, mpz_t q, element_t P,
	element_ptr Qx, element_ptr Qy)
{
    int m;
    element_t v;
    element_t Z;
    element_t a, b, c;
    element_t t0, t1;
    element_ptr t2 = a, t3 = b, t4 = c;
    element_t e0;
    element_t z, z2;
    element_ptr Zx, Zy;
    const element_ptr curve_a = curve_a_coeff(P);
    const element_ptr Px = curve_x_coord(P);
    const element_ptr Py = curve_y_coord(P);

    void proj_double(void)
    {
	//t0 = 3x^2 + (curve_a) z^4
	element_square(t0, Zx);
	//element_mul_si(t0, t0, 3);
	element_double(t1, t0);
	element_add(t0, t0, t1);
	element_square(t1, z2);
	element_mul(t1, t1, curve_a);
	element_add(t0, t0, t1);

	//z_out = 2 y z
	element_mul(z, Zy, z);
	//element_mul_si(z, z, 2);
	element_double(z, z);
	element_square(z2, z);

	//t1 = 4 x y^2
	element_square(t2, Zy);
	element_mul(t1, Zx, t2);
	//element_mul_si(t1, t1, 4);
	element_double(t1, t1);
	element_double(t1, t1);

	//x_out = t0^2 - 2 t1
	//element_mul_si(t3, t1, 2);
	element_double(t3, t1);
	element_square(Zx, t0);
	element_sub(Zx, Zx, t3);

	//t2 = 8y^4
	element_square(t2, t2);
	//element_mul_si(t2, t2, 8);
	element_double(t2, t2);
	element_double(t2, t2);
	element_double(t2, t2);

	//y_out = t0(t1 - x_out) - t2
	element_sub(t1, t1, Zx);
	element_mul(t0, t0, t1);
	element_sub(Zy, t0, t2);
    }

    void proj_mixin(void)
    {
	//t2 = Px z^2
	element_mul(t2, z2, Px);

	//t3 = Zx - t2
	element_sub(t3, Zx, t2);

	//t0 = Py z^3
	element_mul(t0, z2, Py);
	element_mul(t0, t0, z);

	//t1 = Zy - t0
	element_sub(t1, Zy, t0);

	//e7 = Zx + t2, use t2 to double for e7
	element_add(t2, Zx, t2);

	//e8 = Zy + t0, use t0 to double for e8
	element_add(t0, Zy, t0);

	//z = z t3
	element_mul(z, z, t3);
	element_square(z2, z);

	//Zx = t1^2 - e7 t3^2
	//t3 now holds t3^3,
	//t4 holds e7 t3^2
	element_square(t4, t3);
	element_mul(t3, t4, t3);
	element_square(Zx, t1);
	element_mul(t4, t2, t4);
	element_sub(Zx, Zx, t4);

	//t4 = e7 t3^2 - 2 Zx
	element_sub(t4, t4, Zx);
	element_sub(t4, t4, Zx);

	//Zy = (t4 t1 - e8 t3^3)/2
	element_mul(t4, t4, t1);
	element_mul(t0, t0, t3);
	element_sub(t4, t4, t0);
	element_halve(Zy, t4);
    }

    void do_tangent(void)
    {
	//a = -(3x^2 + cca z^4)
	//b = 2 y z^3
	//c = -(2 y^2 + x a)
	//a = z^2 a
	element_square(a, z2);
	element_mul(a, a, curve_a);
	element_square(b, Zx);
	//element_mul_si(b, b, 3);
	element_double(t0, b);
	element_add(b, b, t0);
	element_add(a, a, b);
	element_neg(a, a);

	element_mul(b, z, z2);
	element_mul(b, b, Zy);
	element_mul_si(b, b, 2);

	element_mul(c, Zx, a);
	element_mul(a, a, z2);
	element_square(t0, Zy);
	element_mul_si(t0, t0, 2);
	element_add(c, c, t0);
	element_neg(c, c);

	d_miller_evalfn(e0, a, b, c, Qx, Qy);
	element_mul(v, v, e0);
    }

    void do_line(void)
    {
	//a = -(Py z^3 - Zy)
	//b = Px z^3 - Zx z
	//c = Zx z Py - Zy Px;

	element_mul(t0, Zx, z);
	element_mul(t1, z2, z);

	element_mul(a, Py, t1);
	element_sub(a, Zy, a);

	element_mul(b, Px, t1);
	element_sub(b, b, t0);

	element_mul(t0, t0, Py);
	element_mul(c, Zy, Px);
	element_sub(c, t0, c);

	d_miller_evalfn(e0, a, b, c, Qx, Qy);
	element_mul(v, v, e0);
    }

    element_init(a, Px->field);
    element_init(b, a->field);
    element_init(c, a->field);
    element_init(t0, a->field);
    element_init(t1, a->field);
    element_init(e0, res->field);
    element_init(z, a->field);
    element_init(z2, a->field);
    element_set1(z);
    element_set1(z2);

    element_init(v, res->field);
    element_init(Z, P->field);

    element_set(Z, P);
    Zx = curve_x_coord(Z);
    Zy = curve_x_coord(Z);

    element_set1(v);
    m = mpz_sizeinbase(q, 2) - 2;

    for(;;) {
	do_tangent();
	if (!m) break;
	proj_double();
	if (mpz_tstbit(q, m)) {
	    do_line();
	    proj_mixin();
	}
	m--;
	element_square(v, v);
    }

    element_set(res, v);

    element_clear(v);
    element_clear(Z);
    element_clear(a);
    element_clear(b);
    element_clear(c);
    element_clear(t0);
    element_clear(t1);
    element_clear(e0);
    element_clear(z);
    element_clear(z2);
}

static void cc_miller_no_denom_affine(element_t res, mpz_t q, element_t P,
	element_ptr Qx, element_ptr Qy)
{
    int m;
    element_t v;
    element_t Z;
    element_t a, b, c;
    element_t t0;
    element_t e0;
    const element_ptr cca = curve_a_coeff(P);
    const element_ptr Px = curve_x_coord(P);
    const element_ptr Py = curve_y_coord(P);
    element_ptr Zx, Zy;

    /* TODO: when exactly is this not needed?
    void do_vertical(void)
    {
	mapbase(e0, Z->x);
	element_sub(e0, Qx, e0);
	element_mul(v, v, e0);
    }
    */

    void do_tangent(void)
    {
	//a = -(3 Zx^2 + cc->a)
	//b = 2 * Zy
	//c = -(2 Zy^2 + a Zx);

	element_square(a, Zx);
	element_mul_si(a, a, 3);
	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);

	d_miller_evalfn(e0, a, b, c, Qx, Qy);
	element_mul(v, v, e0);
    }

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

	d_miller_evalfn(e0, a, b, c, Qx, Qy);
	element_mul(v, v, e0);
    }

    element_init(a, Px->field);
    element_init(b, a->field);
    element_init(c, a->field);
    element_init(t0, a->field);
    element_init(e0, res->field);

    element_init(v, res->field);
    element_init(Z, P->field);

    element_set(Z, P);
    Zx = curve_x_coord(Z);
    Zy = curve_y_coord(Z);

    element_set1(v);
    m = mpz_sizeinbase(q, 2) - 2;

    for(;;) {
	do_tangent();

	if (!m) break;

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

    element_set(res, v);

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

static void lucas_even(element_ptr out, element_ptr in, mpz_t cofactor)
//assumes cofactor is even
//mangles in
//in cannot be out
{
    element_t temp;
    element_init_same_as(temp, out);
    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^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 cc_tatepower(element_ptr out, element_ptr in, pairing_t pairing)
{
    mnt_pairing_data_ptr p = pairing->data;
    if (p->k == 6) {
	element_t e0, e2, e3;
	element_init(e0, p->Fqk);
	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);
	    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);
		element_add(e0im0, e0im0, inim[0]);
	    } else {