curve.c

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

    if (b != 0) {
	fp2_mul(v, v, fb, p);
	fp2_mul(vdenom, vdenom, fbdenom, p);
	//g
	tate_get_line(v, Qhat, Z, bP, p);
	//h
	point_add(Z, Z, bP, curve);
	tate_get_vertical(vdenom, Qhat, Z, p);
    }

    if (curve->solinasa < 0) {
    //the sign of solinasa records whether it's +1 or -1
	tate_get_vertical(vdenom, Qhat, P, p);
    }

    //g
    tate_get_vertical(v, Qhat, Z, p);
    //h
    fp2_div(res, v, vdenom, p);

    point_clear(Z);
    point_clear(bP);
    fp2_clear(v);
    fp2_clear(vdenom);
    fp2_clear(fb);
    fp2_clear(fbdenom);
    mpz_clear(z);
    mpz_clear(temp);
}

int point_valid_p(point_t P, curve_t curve)
{
    int result = 1;
    mpz_ptr p = curve->p;
    fp2_t error;
    fp2_t temp;

    if (P->infinity) return result;

    fp2_init(error);
    fp2_init(temp);
    fp2_mul(error, P->x, P->x, p);
    fp2_mul(error, error, P->x, p);
    fp2_mul(temp, P->y, P->y, p);
    fp2_sub(error, temp, error, p);
    fp2_sub_ui(error, error, 1, p);

    if (mpz_cmp_ui(error->a, 0) || mpz_cmp_ui(error->b, 0)) {
	//printf("error = ");
	//fp2_out_str(NULL, 0, error);
	//printf("\n");
	result = 0;
    }

    fp2_clear(error);
    fp2_clear(temp);
    return result;
}

int point_special_p(point_t P, curve_t curve)
//check that P is in E/F_p and finite
{
    int result = 1;
    mpz_ptr p = curve->p;
    mpz_t error;
    mpz_t temp;

    if (P->infinity) return 0;

    if (mpz_cmp_ui(P->x->b, 0)) return 0;
    if (mpz_cmp_ui(P->y->b, 0)) return 0;

    mpz_init(error);
    mpz_init(temp);
    mpz_mul(error, P->x->a, P->x->a);
    mpz_mod(error, error, p);
    mpz_mul(error, error, P->x->a);
    mpz_mod(error, error, p);

    mpz_mul(temp, P->y->a, P->y->a);
    mpz_sub(error, temp, error);
    mpz_sub_ui(error, error, 1);
    mpz_mod(error, error, p);

    if (mpz_cmp_ui(error, 0)) {
	//printf("error = ");
	//fp2_out_str(NULL, 0, error);
	//printf("\n");
	result = 0;
    }

    mpz_clear(error);
    mpz_clear(temp);
    return result;
}

void tate_pairing(fp2_ptr res, point_ptr P, point_ptr Q, curve_t curve)
// res = e(P,Q) where e is the Tate pairing
// assume P in E/F_p
{
    bm_put(bm_get_time(), "miller0");
    tate_solinas_miller(res, P, Q, curve);
    bm_put(bm_get_time(), "miller1");

    tate_power(res, curve);
}

#if 0
void weil_pairing(fp2_ptr res, point_ptr P, point_ptr Q)
// res = e(P,Q) where e is the Weil pairing
{
    point_t Phat, Qhat;
    fp2_t num, denom;

    point_init(Phat);
    point_init(Qhat);
    fp2_init(num);
    fp2_init(denom);

    if (!miller_randomized_flag) {
	new_random_Rs();
    }

    for(;;) {
	point_add(Phat, P, curve->R1, curve);
	point_add(Qhat, Q, curve->R2, curve);
	if (!solinas_miller(num, P, Phat, Qhat, curve->R1, curve->R2)) goto millerfail;
	if (!solinas_miller(denom, Q, Qhat, Phat, curve->R2, curve->R1)) goto millerfail;
	break;
millerfail:
	new_random_Rs();
    }

    point_clear(Phat);
    point_clear(Qhat);
    fp2_div(res, num, denom, p);
    fp2_clear(num);
    fp2_clear(denom);
}
#endif

void zzp_point_double(mpz_t x, mpz_t y, mpz_t Px, mpz_t Py, mpz_t p)
//(x,y,1) = 2(Px, Py,1)
//assume points are on E/F_p
//and we don't expect to get O, either as input or output
{
    //line: Y - (lambda X + mu)
    //we assume a2 = a4 = 0
    mpz_t lambda, temp;
    mpz_init(lambda);
    mpz_init(temp);

    mpz_add(lambda, Px, Px);
    mpz_add(lambda, lambda, Px);
    mpz_mul(lambda, lambda, Px);
    mpz_mod(lambda, lambda, p);
    mpz_add(temp, Py, Py);
    mpz_invert(temp, temp, p);
    mpz_mul(lambda, lambda, temp);
    mpz_mod(lambda, lambda, p);

    mpz_set(temp, Px);
    mpz_mul(x, lambda, lambda);
    mpz_sub(x, x, temp);
    mpz_sub(x, x, temp);
    mpz_mod(x, x, p);

    mpz_sub(y, temp, x);
    mpz_mul(y, y, lambda);
    mpz_sub(y, y, Py);
    mpz_mod(y, y, p);

    mpz_clear(temp);
    mpz_clear(lambda);
}

void zzp_point_add(mpz_t x, mpz_t y, mpz_t Px, mpz_t Py,
	mpz_t Qx, mpz_t Qy, mpz_t p)
//(x, y) = (Px, Py) + (Qx, Qy)
//assume points are on E/F_p, assume output != O
{
    mpz_t lambda, temp;
    mpz_init(lambda);
    mpz_init(temp);

    //line: Y - (lambda X + mu)
    mpz_sub(lambda, Qy, Py);
    mpz_sub(temp, Qx, Px);
    mpz_invert(temp, temp, p);
    mpz_mul(lambda, lambda, temp);
    mpz_mod(lambda, lambda, p);

    mpz_set(temp, Px);
    mpz_mul(x, lambda, lambda);
    mpz_sub(x, x, temp);
    mpz_sub(x, x, Qx);
    mpz_mod(x, x, p);

    mpz_sub(y, temp, x);
    mpz_mul(y, y, lambda);
    mpz_sub(y, y, Py);
    mpz_mod(y, y, p);

    mpz_clear(temp);
    mpz_clear(lambda);
}

void point_mul_preprocess(point_ptr P, curve_t curve)
//get ready for multiplications on P
//set up signed sliding-windowS
{
    int i;
    int m;

    assert(point_special_p(P, curve));

    m = mpz_sizeinbase(curve->q, 2);

    //compute x[1], y[1]
    mpz_set(curve->pre_x[0], P->x->a);
    mpz_set(curve->pre_y[0], P->y->a);

    for (i=1; i<=m; i++) {
	zzp_point_double(curve->pre_x[i], curve->pre_y[i], curve->pre_x[i-1], curve->pre_y[i-1], curve->p);
    }
}

void point_mul_postprocess(point_ptr R, mpz_t n, curve_t curve)
//R = nP (P has been preprocessed)
//P must lie on E/F_p, n must be positive
//assumes nP != O (?)
//uses signed sliding-window method
{
    //compute NAF form (see Blake, Seroussi & Smart IV.2.4)
    int j;
    int m = mpz_sizeinbase(n, 2);
    int *s = (int *) malloc(sizeof(int) * (m+1));
    int c0 = 0, c1;

    mpz_t Rx, Ry, Rz;
    mpz_t x0, y0;
    mpz_ptr p = curve->p;

    assert(mpz_cmp_ui(n, 0) > 0);
    assert(mpz_cmp(n, curve->q) < 0);

    for (j=0; j<=m; j++) {
	c1 = (mpz_tstbit(n, j) + mpz_tstbit(n, j+1) + c0) >> 1;
	s[j] = mpz_tstbit(n, j) + c0 - 2 * c1;
	c0 = c1;
    }

    mpz_init(Rx); mpz_init(Ry); mpz_init(Rz);
    mpz_init(x0); mpz_init(y0);

    //first bit might not be used after NAF conversion
    if (!s[m]) m--;

    mpz_set_ui(Rz, 1);
    //now work in projective coordinates

    mpz_set(Rx, curve->pre_x[m]);
    mpz_set(Ry, curve->pre_y[m]);

    m--;
    while(m>=0) {
	if (s[m]) {
	    if (s[m] < 0) {
		mpz_sub(x0, p, curve->pre_y[m]);
		proj_mix_in(Rx, Ry, Rz, curve->pre_x[m], x0, p);
	    } else {
		proj_mix_in(Rx, Ry, Rz, curve->pre_x[m], curve->pre_y[m], p);
	    }
	}
	m--;
    }

    //convert back to affine
    mpz_invert(Rz, Rz, p);

    mpz_mul(x0, Rz, Rz);
    mpz_mod(x0, x0, p);
    mpz_mul(Rx, Rx, x0);
    mpz_mod(Rx, Rx, p);
    mpz_mul(y0, x0, Rz);
    mpz_mod(y0, y0, p);
    mpz_mul(Ry, Ry, y0);
    mpz_mod(Ry, Ry, p);

    mpz_set_ui(R->x->b, 0);
    mpz_set_ui(R->y->b, 0);
    mpz_set(R->x->a, Rx);
    mpz_set(R->y->a, Ry);
    R->infinity = 0;

    mpz_clear(Rx); mpz_clear(Ry); mpz_clear(Rz);
    mpz_clear(x0); mpz_clear(y0);
    free(s);
}

void point_mul(point_ptr R, mpz_t n, point_ptr P, curve_t curve)
//R = nP
//P must lie on E/F_p, 0 < n
//uses signed sliding-window method
//TODO: handle cases when you get O's during the computation
{
    //compute NAF form (see Blake, Seroussi & Smart IV.2.4)
    int i, j, k, l;
    int m = mpz_sizeinbase(n, 2);
    int *s = (int *) malloc(sizeof(int) * (m+1));
    int c0 = 0, c1;

    mpz_t x[2*windowsizepower+2];
    mpz_t y[2*windowsizepower+2];
    mpz_t Rx, Ry, Rz;
    mpz_ptr p = curve->p;

    assert(mpz_cmp_ui(n, 0) > 0);
    assert(point_special_p(P, curve));

    for (j=0; j<=m; j++) {
	c1 = (mpz_tstbit(n, j) + mpz_tstbit(n, j+1) + c0) >> 1;
	s[j] = mpz_tstbit(n, j) + c0 - 2 * c1;
	c0 = c1;
    }

    mpz_init(Rx); mpz_init(Ry); mpz_init(Rz);

    if (0) {
	if (mpz_size(P->x->b) || mpz_size(P->y->b) || n <= 0) {
	    fprintf(stderr, "BUG: point_mul can't handle points not on E/F_p\n");
	    //exit(1);
	}
    }

    mpz_init(y[0]); //we need a temp. variable later

    //compute x[1], y[1]
    mpz_init(x[1]);
    mpz_init(y[1]);
    mpz_set(x[1], P->x->a);
    mpz_set(y[1], P->y->a);

    //compute x[2], y[2]
    mpz_init(x[2]);
    mpz_init(y[2]);
    zzp_point_double(x[2], y[2], x[1], y[1], p);

    for (i=1; i<=windowsizepower; i++) {
	j = 2 * i + 1;
	k = j - 2;
	mpz_init(x[j]);
	mpz_init(y[j]);
	zzp_point_add(x[j], y[j], x[k], y[k], x[2], y[2], p);
    }

    //first bit might not be used after NAF conversion
    if (!s[m]) m--;

    mpz_set_ui(Rz, 1);
    //now work in projective coordinates

    l = m - windowsize + 1;
    if (l < 0) l = 0;
    for (; l<m; l++) {
	if (s[l]) break;
    }
    j = s[l];
    i = 1;
    for (k=l+1; k<=m; k++) {
	i = i << 1;
	j += s[k] * i;
    }
    if (j < 0) {
	j = -j;
	mpz_set(Rx, x[j]);
	//assumes y[j] != 0
	mpz_sub(Ry, p, y[j]);
    } else {
	mpz_set(Rx, x[j]);
	mpz_set(Ry, y[j]);
    }

    m = l-1;

    while(m>=0) {

	if (!s[m]) {
	    proj_double(Rx, Ry, Rz, p);
	    m--;
	} else {
	    l = m - windowsize + 1;
	    if (l < 0) l = 0;
	    for (; l<m; l++) {
		if (s[l]) break;
	    }
	    j = s[l];
	    proj_double(Rx, Ry, Rz, p);
	    i = 1;
	    for (k=l+1; k<=m; k++) {
		i = i << 1;
		j += s[k] * i;
		proj_double(Rx, Ry, Rz, p);
	    }
	    if (j < 0) {
		j = -j;
		mpz_sub(y[0], p, y[j]);
		proj_mix_in(Rx, Ry, Rz, x[j], y[0], p);
	    } else {
		proj_mix_in(Rx, Ry, Rz, x[j], y[j], p);
	    }
	    m = l-1;
	}
    }

    //convert back to affine
    mpz_invert(Rz, Rz, p);

    //x[1], y[1] no longer needed
    mpz_mul(x[1], Rz, Rz);
    mpz_mod(x[1], x[1], p);
    mpz_mul(Rx, Rx, x[1]);
    mpz_mod(Rx, Rx, p);
    mpz_mul(y[1], x[1], Rz);
    mpz_mod(y[1], y[1], p);
    mpz_mul(Ry, Ry, y[1]);
    mpz_mod(Ry, Ry, p);

    mpz_set_ui(R->x->b, 0);
    mpz_set_ui(R->y->b, 0);
    mpz_set(R->x->a, Rx);
    mpz_set(R->y->a, Ry);
    R->infinity = 0;

    mpz_clear(Rx); mpz_clear(Ry); mpz_clear(Rz);

    mpz_clear(x[2]); mpz_clear(y[2]);
    for (i=0; i<=windowsizepower; i++) {
	j = 2 * i + 1;
	mpz_clear(x[j]); mpz_clear(y[j]);
    }
    mpz_clear(y[0]);
    free(s);
}

void general_point_mul(point_t Q, mpz_t a, point_t P, curve_t curve)
//Q = aP
//can handle P on E/F_p^2, any integer a
{
    point_t R;
    int i, j, k, l, m;
    mpz_t n;

    point_t g[2*windowsizepower+2];

    mpz_init(n);

    point_init(R);
    point_set_O(R);

    mpz_mod(n, a, curve->q);

    if (P->infinity) {
	point_set_O(Q);
	return;
    }
    if (!mpz_cmp_ui(a, 0)) {
	point_set_O(Q);
	return;
    }

    assert(point_valid_p(P, curve));

    //use sliding-window method
    point_init(g[1]);
    point_set(g[1], P);
    point_init(g[2]);
    point_add(g[2], P, P, curve);

    for (i=1; i<=windowsizepower; i++) {
	j = 2 * i + 1;
	point_init(g[j]);
	point_add(g[j], g[j-2], g[2], curve);
    }

    m = mpz_sizeinbase(n, 2) - 1;

    while(m>=0) {
	if (!mpz_tstbit(n, m)) {
	    point_add(R, R, R, curve);
	    m--;
	} else {
	    l = m - windowsize + 1;
	    if (l < 0) l = 0;
	    l = mpz_scan1(n, l);
	    j = 1;
	    point_add(R, R, R, curve);
	    for (k = m - 1; k>=l; k--) {
		j = j << 1;
		if (mpz_tstbit(n, k)) j++;
		point_add(R, R, R, curve);
	    }
	    point_add(R, R, g[j], curve);
	    m = l-1;
	}
    }
    /* check (repeated squaring) */
    if (0) {
	point_t Rcheck;
	point_init(Rcheck);
	point_set_O(Rcheck);
	m = mpz_sizeinbase(n, 2) - 1;
	for(;;) {
	    if (mpz_tstbit(n, m)) {
		point_add(Rcheck, Rcheck, P, curve);
	    }
	    if (m == 0) break;
	    point_add(Rcheck, Rcheck, Rcheck, curve);
	    m--;
	}
	if (!point_equal(R, Rcheck)) {
	    fprintf(stderr, "point_mul(): BUG!\n");
	    //exit(1);
	}
    }

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

    mpz_clear(n);

    point_set(Q, R);
    point_clear(R);
}