curve.c

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

static void pts_get_tangent(fp2_ptr v,
	point_ptr A, point_ptr P, mpz_t z, mpz_t p)
{
    //note: z is in F_p so we can multiply by arbitrary
    //powers of z because it gets killed by Tate exp.
    //to derive what I have: x --> x/z^2, y/z^3 and then multiply
    //by appropriate power of z
    mpz_t a, b, c;
    mpz_t temp, temp2;
    fp2_t f1, f2;

    fp2_init(f1);
    fp2_init(f2);
    //it should be
    //a = -slope_tangent(P.x, P.y);
    //b = 1;
    //c = -(P.y + a * P.x);
    //but we multiply by 2*P.y to avoid division

    assert(!P->infinity); //(could handle with a = b = 0; c = 1;)
    assert(!fp2_is_0(P->y)); // need to be wary of points of order 2
    //could handle with pts_get_vertical(v, A, P, z);

    mpz_init(a);
    mpz_init(b);
    mpz_init(c);
    mpz_init(temp);
    mpz_init(temp2);

    //a = -3x^2 * z^2 (assume a2 = a4 = 0)
#if 1
    zp_mul(temp2, P->x->a, P->x->a, p);
    zp_add(temp, temp2, temp2, p);
    zp_add(temp, temp, temp2, p);
    zp_neg(a, temp, p);
#else
    mpz_mul(temp2, P->x->a, P->x->a);
    mpz_add(temp, temp2, temp2);
    mpz_add(temp, temp, temp2);
    mpz_neg(a, temp);
    mpz_mod(a, a, p);
#endif
    zp_mul(temp2, z, z, p);
    zp_mul(a, a, temp2, p);

    zp_mul(temp2, temp2, z, p);

    //b = (y + y)z^3;
    zp_add(b, P->y->a, P->y->a, p);
    zp_mul(b, b, temp2, p);

    //c = - (2y) * y - (-3x^2) * x;
    zp_mul(c, temp, P->x->a, p);
    zp_mul(temp2, P->y->a, P->y->a, p);
    zp_add(temp2, temp2, temp2, p);
    zp_sub(c, c, temp2, p);

    assert(!A->infinity);

    fp2_mul_mpz(f1, A->x, a, p);
    fp2_mul_mpz(f2, A->y, b, p);
    fp2_add(f1, f1, f2, p);
    //zp_mul(b, A->y->a, b, p);
    //zp_add(f1->a, f1->a, b, p);
    zp_add(f1->a, f1->a, c, p);
    fp2_mul(v, v, f1, p);

    fp2_clear(f1);
    fp2_clear(f2);

    mpz_clear(a);
    mpz_clear(b);
    mpz_clear(c);
    mpz_clear(temp);
    mpz_clear(temp2);
    return;
}

static void tate_get_vertical(fp2_ptr v, point_ptr A, point_ptr P, mpz_t p)
{
    fp2_t temp;
    fp2_init(temp);

    if (A->infinity) {
	fprintf(stderr, "can't evaluate at O!\n");
	//fp2_set_0(v);
	return;
    }
    if (P->infinity) {
	//a = b = 0; c = 1;
	return;
    }
    //a = 1; b = 0; c = -P.x;
    fp2_sub(temp, A->x, P->x, p);
    fp2_mul(v, v, temp, p);

    fp2_clear(temp);
}

static void tate_get_tangent(fp2_ptr v, point_ptr A, point_ptr P, mpz_t p)
{
    fp2_t a, b, c;
    fp2_t temp, temp2;
    //it should be
    //a = -slope_tangent(P.x, P.y);
    //b = 1;
    //c = -(P.y + a * P.x);
    //but we multiply by 2*P.y to avoid division

    //a = -Px * (Px + Px + Px + *twicea2) - *a4;
    //assume a2 = a4 = 0

    if (P->infinity) {
	//a = b = 0; c = 1;
	return;
    }

    fp2_init(temp);

    //need to be wary of points of order 2
    if (fp2_is_0(P->y)) {
	//a = 1; b = 0; c = -P.x;
	fp2_sub(temp, A->x, P->x, p);
	fp2_mul(v, v, temp, p);

	fp2_clear(temp);
	return;
    }

    fp2_init(a);
    fp2_init(b);
    fp2_init(c);
    fp2_init(temp2);

    fp2_add(temp, P->x, P->x, p);
    fp2_add(temp, temp, P->x, p);
    fp2_mul(temp, temp, P->x, p);
    fp2_neg(a, temp, p);

    //b = Py + Py;
    fp2_add(b, P->y, P->y, p);

    //c = - b * Py - a * Px;
    fp2_mul(temp, b, P->y, p);
    fp2_mul(c, a, P->x, p);
    fp2_add(c, c, temp, p);
    fp2_neg(c, c, p);

    if (A->infinity) {
	fprintf(stderr, "can't evaluate at O!\n");
	//fp2_mul(v, v, b, p);
    } else {
	fp2_mul(temp, a, A->x, p);
	fp2_mul(temp2, b, A->y, p);
	fp2_add(temp, temp, temp2, p);
	fp2_add(temp, temp, c, p);
	fp2_mul(v, v, temp, p);
    }

    fp2_clear(a);
    fp2_clear(b);
    fp2_clear(c);
    fp2_clear(temp);
    fp2_clear(temp2);
    return;
}

static void tate_get_line(fp2_ptr v, point_ptr A, point_ptr P, point_ptr Q, mpz_t p)
{
    fp2_t a, b, c;
    fp2_t temp, temp2;

    //cases involving O
    if (P->infinity) {
	tate_get_vertical(v, A, Q, p);
	return;
    }
    if (Q->infinity) {
	tate_get_vertical(v, A, P, p);
	return;
    }

    fp2_init(temp);

    //check if we need a tangent or vertical
    if (fp2_equal(P->x, Q->x)) {
	fp2_neg(temp, Q->y, p);
	if (fp2_equal(P->y, temp)) {
	    //a = 1; b = 0; c = -P.x;
	    fp2_sub(temp, A->x, P->x, p);
	    fp2_mul(v, v, temp, p);

	    fp2_clear(temp);
	    return;
	}
	//othewise P = Q
	fp2_clear(temp);
	tate_get_tangent(v, A, P, p);
	return;

    }
    //normal case (P != Q)
    fp2_init(a);
    fp2_init(b);
    fp2_init(c);
    fp2_init(temp2);

    //it should be
    //a = -(Q.y - P.y) / (Q.x - P.x);
    //b = 1;
    //c = -(P.y + a * P.x);
    //but we'll multiply by Q.x - P.x to avoid division
    fp2_sub(b, Q->x, P->x, p);
    fp2_sub(a, P->y, Q->y, p);
    fp2_mul(temp, b, P->y, p);
    fp2_mul(c, a, P->x, p);
    fp2_add(c, c, temp, p);
    fp2_neg(c, c, p);

    if (A->infinity) {
	fprintf(stderr, "can't evaluate at O!\n");
    } else {
	fp2_mul(temp, a, A->x, p);
	fp2_mul(temp2, b, A->y, p);
	fp2_add(temp, temp, temp2, p);
	fp2_add(temp, temp, c, p);
	fp2_mul(v, v, temp, p);
    }

    fp2_clear(a);
    fp2_clear(b);
    fp2_clear(c);
    fp2_clear(temp);
    fp2_clear(temp2);
}

static void get_vertical(fp2_ptr v, fp2_ptr vdenom,
	point_ptr A, point_ptr B, point_ptr P, mpz_t p)
{
    fp2_t temp;
    fp2_init(temp);

    if (A->infinity || B->infinity) {
	fprintf(stderr, "can't evaluate at O!\n");
	//fp2_set_0(v);
	//fp2_set_0(vdenom);
	return;
    }
    if (P->infinity) {
	//a = b = 0; c = 1;
	return;
    }
    //a = 1; b = 0; c = -P.x;
    fp2_sub(temp, A->x, P->x, p);
    fp2_mul(v, v, temp, p);

    fp2_sub(temp, B->x, P->x, p);
    fp2_mul(vdenom, vdenom, temp, p);

    fp2_clear(temp);
}

static void get_tangent(fp2_ptr v, fp2_ptr vdenom,
	point_ptr A, point_ptr B, point_ptr P, mpz_t p)
{
    fp2_t a, b, c;
    fp2_t temp, temp2;
    //it should be
    //a = -slope_tangent(P.x, P.y);
    //b = 1;
    //c = -(P.y + a * P.x);
    //but we multiply by 2*P.y to avoid division

    //a = -Px * (Px + Px + Px + *twicea2) - *a4;
    //assume a2 = a4 = 0

    if (P->infinity) {
	//a = b = 0; c = 1;
	return;
    }

    fp2_init(temp);

    //need to be wary of points of order 2
    if (fp2_is_0(P->y)) {
	//a = 1; b = 0; c = -P.x;
	fp2_sub(temp, A->x, P->x, p);
	fp2_mul(v, v, temp, p);

	fp2_sub(temp, B->x, P->x, p);
	fp2_mul(vdenom, vdenom, temp, p);

	fp2_clear(temp);
	return;
    }

    fp2_init(a);
    fp2_init(b);
    fp2_init(c);
    fp2_init(temp2);

    fp2_add(temp, P->x, P->x, p);
    fp2_add(temp, temp, P->x, p);
    fp2_mul(temp, temp, P->x, p);
    fp2_neg(a, temp, p);

    //b = Py + Py;
    fp2_add(b, P->y, P->y, p);

    //c = - b * Py - a * Px;
    fp2_mul(temp, b, P->y, p);
    fp2_mul(c, a, P->x, p);
    fp2_add(c, c, temp, p);
    fp2_neg(c, c, p);

    if (A->infinity) {
	fprintf(stderr, "can't evaluate at O!\n");
    } else {
	fp2_mul(temp, a, A->x, p);
	fp2_mul(temp2, b, A->y, p);
	fp2_add(temp, temp, temp2, p);
	fp2_add(temp, temp, c, p);
	fp2_mul(v, v, temp, p);
    }
    if (B->infinity) {
	fprintf(stderr, "can't evaluate at O!\n");
    } else {
	fp2_mul(temp, a, B->x, p);
	fp2_mul(temp2, b, B->y, p);
	fp2_add(temp, temp, temp2, p);
	fp2_add(temp, temp, c, p);
	fp2_mul(vdenom, vdenom, temp, p);
    }

    fp2_clear(a);
    fp2_clear(b);
    fp2_clear(c);
    fp2_clear(temp);
    fp2_clear(temp2);
    return;
}

static void get_line(fp2_ptr v, fp2_ptr vdenom,
	point_ptr A, point_ptr B, point_ptr P, point_ptr Q, mpz_t p)
{
    fp2_t a, b, c;
    fp2_t temp, temp2;

    //cases involving O
    if (P->infinity) {
	get_vertical(v, vdenom, A, B, Q, p);
	return;
    }
    if (Q->infinity) {
	get_vertical(v, vdenom, A, B, P, p);
	return;
    }

    fp2_init(temp);

    //check if we need a tangent or vertical
    if (fp2_equal(P->x, Q->x)) {
	fp2_neg(temp, Q->y, p);
	if (fp2_equal(P->y, temp)) {
	    //a = 1; b = 0; c = -P.x;
	    fp2_sub(temp, A->x, P->x, p);
	    fp2_mul(v, v, temp, p);

	    fp2_sub(temp, B->x, P->x, p);
	    fp2_mul(vdenom, vdenom, temp, p);

	    fp2_clear(temp);
	    return;
	}
	//othewise P = Q
	fp2_clear(temp);
	get_tangent(v, vdenom, A, B, P, p);
	return;

    }
    //normal case (P != Q)
    fp2_init(a);
    fp2_init(b);
    fp2_init(c);
    fp2_init(temp2);

    //it should be
    //a = -(Q.y - P.y) / (Q.x - P.x);
    //b = 1;
    //c = -(P.y + a * P.x);
    //but we'll multiply by Q.x - P.x to avoid division
    fp2_sub(b, Q->x, P->x, p);
    fp2_sub(a, P->y, Q->y, p);
    fp2_mul(temp, b, P->y, p);
    fp2_mul(c, a, P->x, p);
    fp2_add(c, c, temp, p);
    fp2_neg(c, c, p);

    if (A->infinity) {
	fprintf(stderr, "can't evaluate at O!\n");
	//fp2_mul(v, v, b, p);
    } else {
	fp2_mul(temp, a, A->x, p);
	fp2_mul(temp2, b, A->y, p);
	fp2_add(temp, temp, temp2, p);
	fp2_add(temp, temp, c, p);
	fp2_mul(v, v, temp, p);
    }
    if (B->infinity) {
	fprintf(stderr, "can't evaluate at O!\n");
	//fp2_mul(vdenom, vdenom, b, p);
    } else {
	fp2_mul(temp, a, B->x, p);
	fp2_mul(temp2, b, B->y, p);
	fp2_add(temp, temp, temp2, p);
	fp2_add(temp, temp, c, p);
	fp2_mul(vdenom, vdenom, temp, p);
    }

    fp2_clear(a);
    fp2_clear(b);
    fp2_clear(c);
    fp2_clear(temp);
    fp2_clear(temp2);
}

int simple_miller(fp2_ptr res, point_ptr P, point_ptr Phat,
	point_ptr Qhat, point_ptr R1, point_ptr R2, curve_t curve)
{
    int m;
    fp2_t f1;
    fp2_t v, vdenom;
    point_t Z;
    mpz_ptr p = curve->p;

    point_init(Z);
    point_set_O(Z);

    fp2_init(f1);
    fp2_init(v);
    fp2_init(vdenom);

    fp2_set_1(v);
    fp2_set_1(vdenom);

    get_vertical(v, vdenom, Qhat, R2, Phat, p);
    get_line(v, vdenom, R2, Qhat, P, R1, p);

    fp2_div(f1, v, vdenom, p);

    fp2_set_1(v);
    fp2_set_1(vdenom);

    m = mpz_sizeinbase(curve->q, 2) - 1;
    for(;;) {
	if (mpz_tstbit(curve->q, m)) {
	    fp2_mul(v, v, f1, p);
	    //g
	    get_line(v, vdenom, Qhat, R2, Z, P, p);
	    //h
	    point_add(Z, Z, P, curve);
	    get_vertical(v, vdenom, R2, Qhat, Z, p);
	}
	if (m == 0) break;
	m--;
	fp2_sqr(v, v, p);
	fp2_sqr(vdenom, vdenom, p);
	//g
	get_line(v, vdenom, Qhat, R2, Z, Z, p);
	//h
	point_add(Z, Z, Z, curve);
	get_vertical(v, vdenom, R2, Qhat, Z, p);
    }
    fp2_div(res, v, vdenom, p);

    fp2_clear(f1);
    fp2_clear(v);
    fp2_clear(vdenom);

    point_clear(Z);

    return 1;
}

int general_miller(fp2_ptr res, point_ptr P, point_ptr Phat,
	point_ptr Qhat, point_ptr R1, point_ptr R2, curve_t curve)
{
    //TODO: use NAF
    //use sliding-window method
    fp2_t g[2*windowsizepower+2];
    point_t bP[2*windowsizepower+2];

    fp2_t vdenom, v;
    point_t Z;
    int i, j, k, l, m;
    mpz_ptr p = curve->p;

    fp2_init(vdenom);
    fp2_init(v);

    fp2_set_1(v);
    fp2_set_1(vdenom);

    //work out f_1
    get_vertical(v, vdenom, Qhat, R2, Phat, p);
    get_line(v, vdenom, R2, Qhat, P, R1, p);
    if (fp2_is_0(v) || fp2_is_0(vdenom)) goto bad_R;
    fp2_init(g[1]);
    fp2_div(g[1], v, vdenom, p);
    point_init(bP[1]);
    point_set(bP[1], P);
    fp2_set_1(v);
    fp2_set_1(vdenom);

    //work out f_2
    point_init(bP[2]);
    point_add(bP[2], P, P, curve);
    //get_line(v, vdenom, Qhat, R2, P, P);
    get_tangent(v, vdenom, Qhat, R2, P, p);
    get_vertical(v, vdenom, R2, Qhat, bP[2], p);
    if (fp2_is_0(v) || fp2_is_0(vdenom)) goto bad_R;
    fp2_init(g[2]);
    fp2_sqr(g[2], g[1], p);
    fp2_mul(g[2], g[2], v, p);
    fp2_div(g[2], g[2], vdenom, p);
    fp2_set_1(v);
    fp2_set_1(vdenom);

    //work out rest of g[], bP[]
    for (i=1; i<=windowsizepower; i++) {
	j = 2 * i + 1;
	k = j - 2;
	point_init(bP[j]);
	point_add(bP[j], bP[k], bP[2], curve);
	get_line(v, vdenom, Qhat, R2, bP[k], bP[2], p);
	get_vertical(v, vdenom, R2, Qhat, bP[j], p);
	if (fp2_is_0(v) || fp2_is_0(vdenom)) goto bad_R;
	fp2_init(g[j]);