curve.c

斯坦福大学密码学家Boneh的基于身份的公钥密码系统

	fp2_mul(g[j], g[k], g[2], p);
	fp2_mul(g[j], g[j], v, p);
	fp2_div(g[j], g[j], vdenom, p);

	fp2_set_1(v);
	fp2_set_1(vdenom);
    }

    point_init(Z);
    point_set_O(Z);

    m = mpz_sizeinbase(curve->q, 2) - 1;
    while(m>=0) {
	if (!mpz_tstbit(curve->q, m)) {
	    fp2_sqr(v, v, p);
	    fp2_sqr(vdenom, vdenom, p);
	    get_tangent(v, vdenom, Qhat, R2, Z, p);
	    point_add(Z, Z, Z, curve);
	    get_vertical(v, vdenom, R2, Qhat, Z, p);
	    m--;
	} else {
	    l = m - windowsize + 1;
	    if (l < 0) l = 0;
	    l = mpz_scan1(curve->q, l);
	    j = 1;
	    fp2_sqr(v, v, p);
	    fp2_sqr(vdenom, vdenom, p);
	    get_tangent(v, vdenom, Qhat, R2, Z, p);
	    point_add(Z, Z, Z, curve);
	    get_vertical(v, vdenom, R2, Qhat, Z, p);
	    for (k = m - 1; k>=l; k--) {
		j = j << 1;
		if (mpz_tstbit(curve->q, k)) j++;
		fp2_sqr(v, v, p);
		fp2_sqr(vdenom, vdenom, p);
		get_tangent(v, vdenom, Qhat, R2, Z, p);
		point_add(Z, Z, Z, curve);
		get_vertical(v, vdenom, R2, Qhat, Z, p);
	    }
	    fp2_mul(v, v, g[j], p);
	    get_line(v, vdenom, Qhat, R2, Z, bP[j], p);
	    point_add(Z, Z, bP[j], curve);
	    get_vertical(v, vdenom, R2, Qhat, Z, p);
	    m = l-1;
	}
    }

    if (fp2_is_0(v) || fp2_is_0(vdenom)) goto bad_R;
    else fp2_div(res, v, vdenom, p);

    // check with simple algorithm
    if (0) {
	fp2_t check;

	fp2_init(check);

	simple_miller(check, P, Phat, Qhat, R1, R2, curve);

	if (!fp2_equal(check, res)) {
	    printf("res check = ");
	    fp2_out_str(NULL, 0, res);
	    printf(" ");
	    fp2_out_str(NULL, 0, check);
	    printf("\n");
	    fprintf(stderr, "BUG: miller() is incorrectly implemented\n");
	    //fp2_set(res, check);
	}

	fp2_clear(check);
    }

    fp2_clear(v);
    fp2_clear(vdenom);
    point_clear(Z);

    fp2_clear(g[2]);
    point_clear(bP[2]);
    for (i=0; i<=windowsizepower; i++) {
	j = 2*i + 1;
	point_clear(bP[j]);
	fp2_clear(g[j]);
    }

    return 1;

bad_R:
    printf("bad R\n");

    fp2_clear(v);
    fp2_clear(vdenom);

    /* TODO: memory leak if jumped here after g[] inited
    point_clear(Z);
    fp2_clear(g[2]);
    point_clear(bP[2]);
    for (i=0; i<=windowsizepower; i++) {
	j = 2*i + 1;
	point_clear(bP[j]);
	fp2_clear(g[j]);
    }
    */

    return 0;
}

static void pts_preprocess_vertical(miller_cache_t mc,
	int i, point_t P, mpz_t z, mpz_t p)
{
    mpz_t t0;

    mpz_init(t0);

    mpz_mul(t0, z, z);
    mpz_invert(t0, t0, p);

    mpz_mul(t0, t0, P->x->a);
    mpz_neg(t0, t0);
    mpz_mod(mc->denomc[i], t0, p);

    mpz_clear(t0);
}

void pts_preprocess_tangent(miller_cache_t mc,
	int i, point_t P, mpz_t z, mpz_t p)
{
    mpz_t t0;
    mpz_t *numa, *numc;

    mpz_init(t0);
    numa = mc->numa;
    numc = mc->numc;

    assert(mpz_cmp_ui(P->y->a, 0)); //assume P not of order 2
    //could handle with:
    //{
    //	mpz_set_ui(numa[i], 1);
    //	mpz_set(numc[i], P->x->a);
    //	return;
    //}

    //a = -3x^2/2yz
    //c = -(y + a * xz)/z^3;

    mpz_mul_2exp(t0, P->y->a, 1);
    mpz_mul(t0, t0, z);
    mpz_invert(t0, t0, p);
    mpz_mul(t0, t0, P->x->a);
    mpz_mul(t0, t0, P->x->a);
    mpz_mul_si(t0, t0, -3);
    mpz_mod(numa[i], t0, p);

    mpz_mul(t0, z, z);
    mpz_mul(t0, t0, z);
    mpz_invert(t0, t0, p);
    mpz_mul(numc[i], numa[i], P->x->a);
    mpz_mul(numc[i], numc[i], z);
    mpz_add(numc[i], numc[i], P->y->a);
    mpz_neg(numc[i], numc[i]);
    mpz_mul(t0, numc[i], t0);
    mpz_mod(numc[i], t0, p);

    mpz_clear(t0);
}

void pts_preprocess_line(mpz_t a, mpz_t c, point_t P, point_t Q, mpz_t p)
{
    mpz_t t0;

    mpz_init(t0);

    //assume P.x != Q.x

    //a = -(Q.y - P.y) / (Q.x - P.x);
    //c = -(P.y + a * P.x);

    mpz_sub(c, Q->x->a, P->x->a);
    mpz_invert(c, c, p);
    mpz_sub(a, P->y->a, Q->y->a);
    mpz_mul(a, a, c);
    mpz_mod(a, a, p);

    mpz_mul(c, a, P->x->a);
    mpz_add(c, c, P->y->a);
    mpz_neg(c, c);
    mpz_mod(c, c, p);

    mpz_clear(t0);
}

void tate_preprocess(miller_cache_t mc, point_ptr P, curve_t curve)
//for primes of the form 2^a +- 2^b +- 1
//uses proj. coords, assumes P is a point over F_p
//and that order of group = Solinas prime
{
    //specialized for Solinas primes
    int a, b;
    point_t bP;

    point_t Z;
    int i;
    mpz_t z, temp;
    mpz_ptr p = curve->p;

    mpz_init(z);
    mpz_init(temp);
    mpz_set_ui(z, 1);

    a = abs(curve->solinasa);
    b = abs(curve->solinasb);

    point_init(Z);
    point_init(bP);

    point_set(Z, P);
    i = 0;
    if (b != 0) {
	//work out f_2^b
	for(;i<b; i++) {
	    //g
	    //tate_get_tangent(v, Qhat, Z, p);
	    //pts_get_tangent(v, Qhat, Z, z);
	    pts_preprocess_tangent(mc, i, Z, z, p);
	    //h
	    //point_add(Z, Z, Z, curve);
	    proj_double(Z->x->a, Z->y->a, z, p);
	    //tate_get_vertical(vdenom, Qhat, Z, p);
	    //pts_get_vertical(vdenom, Qhat, Z, z, p);
	    pts_preprocess_vertical(mc, i, Z, z, p);
	}

	mpz_invert(z, z, p);
	zp_mul(temp, z, z, p);
	fp2_mul_mpz(Z->x, Z->x, temp, p);
	zp_mul(temp, temp, z, p);
	fp2_mul_mpz(Z->y, Z->y, temp, p);
	mpz_set_ui(z, 1);

	point_set(bP, Z);

	if (curve->solinasb < 0) {
	    //tate_get_vertical(fbdenom, Qhat, Z);
	    mpz_neg(mc->denomsb, Z->x->a);
	    fp2_neg(bP->y, bP->y, p);
	}
    }

    //work out f_2^a
    for(; i<a; i++) {
	//g
	//pts_get_tangent(v, Qhat, Z, z);
	pts_preprocess_tangent(mc, i, Z, z, p);
	//h
	proj_double(Z->x->a, Z->y->a, z, p);
	//pts_get_vertical(vdenom, Qhat, Z, z);
	pts_preprocess_vertical(mc, i, Z, z, p);
    }

    mpz_invert(z, z, p);
    zp_mul(temp, z, z, p);
    fp2_mul_mpz(Z->x, Z->x, temp, p);
    zp_mul(temp, temp, z, p);
    fp2_mul_mpz(Z->y, Z->y, temp, p);
    mpz_set_ui(z, 1);

    //work out f_(2^a +- 2^b +- 1)
    if (b != 0) {
	//g
	//tate_get_line(v, Qhat, Z, bP);
	pts_preprocess_line(mc->numl1a, mc->numl1c, Z, bP, p);
	//h
	point_add(Z, Z, bP, curve);
	//tate_get_vertical(vdenom, Qhat, Z);
	mpz_sub(mc->denoml1c, p, Z->x->a);
    }
    //the sign of solinasa records whether it's +1 or -1
    if (curve->solinasa < 0) {
	//tate_get_vertical(vdenom, Qhat, Z);
	mpz_sub(mc->denoms1, p, Z->x->a);
    }

    //g
    //tate_get_line(v, Qhat, Z, cP);
    //now Z = -cP so g is vertical
    mpz_sub(mc->numl2c, p, Z->x->a);
    //h
    //point_add(Z, Z, cP, curve);
    //now Z = O, so h = 1
    //tate_get_vertical(vdenom, Qhat, Z);

    point_clear(Z);
    point_clear(bP);

    mpz_clear(z);
    mpz_clear(temp);
}

void miller_postprocess(fp2_ptr res, miller_cache_t mc,
	point_ptr Q, curve_t curve)
//for primes of the form 2^a +- 2^b +- 1
//uses proj. coords, assumes P is a point over F_p
//and that order of group = Solinas prime
{
    //specialized for Solinas primes
    int a, b;
    fp2_t fb, fbdenom;
    fp2_t t0;
    fp2_t vdenom, v;
    int i;
    mpz_t *numa, *numc;
    mpz_t *denomc;
    mpz_ptr p = curve->p;

    numa = mc->numa;
    numc = mc->numc;
    denomc = mc->denomc;

    a = abs(curve->solinasa);
    b = abs(curve->solinasb);

    fp2_init(t0);
    fp2_init(vdenom);
    fp2_init(v);
    fp2_init(fb);
    fp2_init(fbdenom);

    fp2_set_1(v);
    fp2_set_1(vdenom);

    //f_1 = 1

    i = 0;
    if (b != 0) {
	//work out f_2^b
	for(;i<b; i++) {
	    fp2_sqr(v, v, p);
	    fp2_sqr(vdenom, vdenom, p);
	    //g
	    //tate_get_tangent(v, Qhat, Z, p);
	    //pts_get_tangent(v, Qhat, Z, z, p);
	    {
		fp2_mul_mpz(t0, Q->x, numa[i], p);
		fp2_add(t0, t0, Q->y, p);
		mpz_add(t0->a, t0->a, numc[i]);
		fp2_mul(v, v, t0, p);
	    }
	    //h
	    //tate_get_vertical(vdenom, Qhat, Z, p);
	    //pts_get_vertical(vdenom, Qhat, Z, z, p);
	    {
		fp2_set(t0, Q->x);
		mpz_add(t0->a, t0->a, denomc[i]);
		fp2_mul(vdenom, vdenom, t0, p);
	    }
	}

	if (curve->solinasb < 0) {
	    fp2_set(fbdenom, v);
	    fp2_set(fb, vdenom);
	    //tate_get_vertical(fbdenom, Qhat, Z, p);
	    {
		fp2_set(t0, Q->x);
		mpz_add(t0->a, t0->a, mc->denomsb);
		fp2_mul(fbdenom, fbdenom, t0, p);
	    }
	} else {
	    fp2_set(fb, v);
	    fp2_set(fbdenom, vdenom);
	}
    }

    //work out f_2^a
    for(; i<a; i++) {
	fp2_sqr(v, v, p);
	fp2_sqr(vdenom, vdenom, p);
	//g
	//pts_get_tangent(v, Qhat, Z, z, p);
	{
	    fp2_mul_mpz(t0, Q->x, numa[i], p);
	    fp2_add(t0, t0, Q->y, p);
	    mpz_add(t0->a, t0->a, numc[i]);
	    fp2_mul(v, v, t0, p);
	}
	//h
	//pts_get_vertical(vdenom, Qhat, Z, z, p);
	{
	    fp2_set(t0, Q->x);
	    mpz_add(t0->a, t0->a, denomc[i]);
	    fp2_mul(vdenom, vdenom, t0, p);
	}
    }

    //work out f_(2^a +- 2^b +- 1)
    if (b != 0) {
	fp2_mul(v, v, fb, p);
	fp2_mul(vdenom, vdenom, fbdenom, p);
	//g
	//tate_get_line(v, Qhat, Z, bP, p);
	{
	    fp2_mul_mpz(t0, Q->x, mc->numl1a, p);
	    fp2_add(t0, t0, Q->y, p);
	    mpz_add(t0->a, t0->a, mc->numl1c);
	    fp2_mul(v, v, t0, p);
	}
	//h
	//tate_get_vertical(vdenom, Qhat, Z, p);
	{
	    fp2_set(t0, Q->x);
	    mpz_add(t0->a, t0->a, mc->denoml1c);
	    fp2_mul(vdenom, vdenom, t0, p);
	}
    }

    if (curve->solinasa < 0) {
    //the sign of solinasa records whether it's +1 or -1
	//tate_get_vertical(vdenom, Qhat, P);
	{
	    fp2_set(t0, Q->x);
	    mpz_add(t0->a, t0->a, mc->denoms1);
	    fp2_mul(vdenom, vdenom, t0, p);
	}
    }

    //g
    //tate_get_line(v, Qhat, Z, cP);
    {
	fp2_set(t0, Q->x);
	mpz_add(t0->a, t0->a, mc->numl2c);
	fp2_mul(v, v, t0, p);
    }
    //h
    //tate_get_vertical(vdenom, Qhat, Z);
    fp2_div(res, v, vdenom, p);

    fp2_clear(v);
    fp2_clear(vdenom);
    fp2_clear(fb);
    fp2_clear(fbdenom);
    fp2_clear(t0);
}

void tate_postprocess(fp2_ptr res, miller_cache_t mc, point_ptr Q, curve_t curve)
{
    bm_put(bm_get_time(), "miller0");
    miller_postprocess(res, mc, Q, curve);
    bm_put(bm_get_time(), "miller1");
    tate_power(res, curve);
}

void tate_solinas_miller(fp2_ptr res, point_ptr P, point_ptr Qhat, curve_t curve)
//for primes of the form 2^a +- 2^b +- 1
//uses proj. coords, assumes P is a point over F_p
//and that order of group = Solinas prime
{
    //specialized for Solinas primes
    int a, b;
    fp2_t fb, fbdenom;
    point_t bP;

    fp2_t vdenom, v;
    point_t Z;
    int i;
    mpz_t z, temp;
    mpz_ptr p = curve->p;

    mpz_init(z);
    mpz_init(temp);
    mpz_set_ui(z, 1);

    a = abs(curve->solinasa);
    b = abs(curve->solinasb);

    point_init(Z);
    point_init(bP);
    fp2_init(vdenom);
    fp2_init(v);
    fp2_init(fb);
    fp2_init(fbdenom);

    fp2_set_1(v);
    fp2_set_1(vdenom);

    //f_1 = 1

    point_set(Z, P);
    i = 0;
    if (b != 0) {
	//work out f_2^b
	for(;i<b; i++) {
	    fp2_sqr(v, v, p);
	    fp2_sqr(vdenom, vdenom, p);
	    //g
	    //tate_get_tangent(v, Qhat, Z);
	    pts_get_tangent(v, Qhat, Z, z, p);
	    //h
	    //point_add(Z, Z, Z, curve);
	    proj_double(Z->x->a, Z->y->a, z, p);
	    //tate_get_vertical(vdenom, Qhat, Z);
	    pts_get_vertical(vdenom, Qhat, Z, z, p);
	}

	mpz_invert(z, z, p);
	zp_mul(temp, z, z, p);
	fp2_mul_mpz(Z->x, Z->x, temp, p);
	zp_mul(temp, temp, z, p);
	fp2_mul_mpz(Z->y, Z->y, temp, p);
	mpz_set_ui(z, 1);

	point_set(bP, Z);

	if (curve->solinasb < 0) {
	    fp2_set(fbdenom, v);
	    fp2_set(fb, vdenom);
	    tate_get_vertical(fbdenom, Qhat, Z, p);
	    fp2_neg(bP->y, bP->y, p);
	} else {
	    fp2_set(fb, v);
	    fp2_set(fbdenom, vdenom);
	}
    }

    //work out f_2^a
    for(; i<a; i++) {
	fp2_sqr(v, v, p);
	fp2_sqr(vdenom, vdenom, p);
	//g
	pts_get_tangent(v, Qhat, Z, z, p);
	//h
	proj_double(Z->x->a, Z->y->a, z, p);
	pts_get_vertical(vdenom, Qhat, Z, z, p);
    }

    mpz_invert(z, z, p);
    zp_mul(temp, z, z, p);
    fp2_mul_mpz(Z->x, Z->x, temp, p);
    zp_mul(temp, temp, z, p);
    fp2_mul_mpz(Z->y, Z->y, temp, p);
    mpz_set_ui(z, 1);

    //work out f_(2^a +- 2^b +- 1)