/* Computes Weil pairing, Tate pairing using Miller's algorithm * Ben Lynn * * For speed, point_random assumes the curve is y^2 = x^3 + 1 and p = 2 mod 3 *//*Copyright (C) 2001 Benjamin Lynn (blynn@cs.stanford.edu)See LICENSE for license*/#include <stdlib.h>#include "curve.h"#include "benchmark.h"#include "mm.h"#include "crypto.h" //for random functions#include <assert.h>enum {    //constants for sliding window algorithms    windowsize = 5,    windowsizepower = 15,	    //this is 2^(windowsize-1) - 1};void point_init(point_ptr P)//allocates memory for a point{    fp2_init(P->x);    fp2_init(P->y);    P->infinity = 0;    mm_tally("point", 1, "init");}void point_clear(point_ptr P)//deallocates memory for a point{    fp2_clear(P->x);    fp2_clear(P->y);    mm_tally("point", -1, "clear");}size_t point_out_str(FILE *stream, int base, point_ptr P){    FILE *fp;    size_t s, status;    if (!stream) fp = stdout;    else fp = stream;    status = fp2_out_str(fp, base, P->x);    if (status < 0) return status;    s = status;    status = fprintf(fp, " ");    if (status < 0) return status;    s += status;    status = fp2_out_str(fp, base, P->y);    if (status < 0) return status;    s += status;    return s;}void point_set(point_ptr src, point_ptr dst)//src = dst{    fp2_set(src->x, dst->x);    fp2_set(src->y, dst->y);    src->infinity = dst->infinity;}void point_set_O(point_ptr P)//P = O{    fp2_set_0(P->x);    fp2_set_1(P->y);    P->infinity = 1;}int point_equal(point_t P, point_t Q)//P == Q?{    if (P->infinity) return Q->infinity;    return fp2_equal(P->x, Q->x) && fp2_equal(P->y, Q->y);}void curve_init(curve_t curve, mpz_t prime, mpz_t qprime)//initializes system parameters//not thread-safe{    int i;    int m;    int c0 = 0, c1;    int count = 0;    int j;    mpz_init(curve->p);    mpz_init(curve->q);    mpz_init(curve->p1onq);    mpz_init(curve->cbrtpwr);    mpz_init(curve->tatepwr);    mpz_set(curve->p, prime);    mpz_set(curve->q, qprime);    m = mpz_sizeinbase(curve->q, 2);    //(uses NAF algorithm)    for (i=0; i<=m; i++) {	c1 = (mpz_tstbit(curve->q, i) + mpz_tstbit(curve->q, i+1) + c0) >> 1;	j = mpz_tstbit(curve->q, i) + c0 - 2 * c1;	if (j != 0) {	    if (count >= 3) {		curve->solinasa = 0;		curve->solinasb = 0;		break;	    }	    if (i == 0) {		curve->solinasa = j;	    } else if (count == 1) {		curve->solinasb = i * j;	    } else {		curve->solinasa *= i;	    }	    count++;	}	c0 = c1;    }    if (count == 2) {	curve->solinasa *= curve->solinasb;	curve->solinasb = 0;    }    //printf("Solinas: a, b : %d %d \n", curve->solinasa, curve->solinasb);    //(p + 1) / q    mpz_add_ui(curve->p1onq, curve->p, 1);    mpz_div(curve->p1onq, curve->p1onq, curve->q);    //(2*p - 1)/3;    mpz_mul_ui(curve->cbrtpwr, curve->p, 2);    mpz_sub_ui(curve->cbrtpwr, curve->cbrtpwr, 1);    mpz_div_ui(curve->cbrtpwr, curve->cbrtpwr, 3);    //(p^2-1)/q    // = (p-1)p1onq    mpz_sub_ui(curve->tatepwr, curve->p, 1);    mpz_mul(curve->tatepwr, curve->tatepwr, curve->p1onq);    curve->pre_x = (mpz_t *) malloc(sizeof(mpz_t) * (m + 1));    curve->pre_y = (mpz_t *) malloc(sizeof(mpz_t) * (m + 1));    for (i=0; i<=m; i++) {	mpz_init(curve->pre_x[i]);	mpz_init(curve->pre_y[i]);    }}void curve_clear(curve_t curve){    int i;    int m = mpz_sizeinbase(curve->q, 2);    mpz_clear(curve->p);    mpz_clear(curve->q);    mpz_clear(curve->p1onq);    mpz_clear(curve->cbrtpwr);    mpz_clear(curve->tatepwr);    for (i=0; i<=m; i++) {	mpz_clear(curve->pre_x[i]);	mpz_clear(curve->pre_y[i]);    }    free(curve->pre_x);    free(curve->pre_y);}void miller_cache_init(miller_cache_t mc, curve_t curve){    int i;    int m = mpz_sizeinbase(curve->q, 2);    mc->numa = (mpz_t *) malloc(sizeof(mpz_t) * m);    mc->numc = (mpz_t *) malloc(sizeof(mpz_t) * m);    mc->denomc = (mpz_t *) malloc(sizeof(mpz_t) * m);    mpz_init(mc->denoms1);    mpz_init(mc->denomsb);    mpz_init(mc->numl1a);    mpz_init(mc->numl1c);    mpz_init(mc->denoml1c);    mpz_init(mc->numl2c);    for (i=0; i<m; i++) {	mpz_init(mc->numa[i]);	mpz_init(mc->numc[i]);	mpz_init(mc->denomc[i]);    }    mc->count = m;}void miller_cache_clear(miller_cache_t mc){    int i;    int m = mc->count;    mpz_clear(mc->denoms1);    mpz_clear(mc->denomsb);    mpz_clear(mc->numl1a);    mpz_clear(mc->numl1c);    mpz_clear(mc->denoml1c);    mpz_clear(mc->numl2c);    for (i=0; i<m; i++) {	mpz_clear(mc->numa[i]);	mpz_clear(mc->numc[i]);	mpz_clear(mc->denomc[i]);    }    free(mc->numa);    free(mc->numc);    free(mc->denomc);}void x_from_y(mpz_t x, mpz_t y, curve_t curve){    //x = cuberoot(y^2 - 1)    mpz_mul(x, y, y);    mpz_sub_ui(x, x, 1);    mpz_mod(x, x, curve->p);    mpz_powm(x, x, curve->cbrtpwr, curve->p);}void fp2_random(fp2_t r, mpz_t p)//r = random element of F_p^2{    mympz_randomm(r->a, p);    mympz_randomm(r->b, p);}void point_random(point_ptr P, curve_t curve)//P = random point on E/F_p{    //this only works for p = 2 mod 3    //and y^2 = x^3 + 1    fp2_t x, y;    fp2_init(x);    fp2_init(y);    mpz_set_ui(x->b, 0);    mpz_set_ui(y->b, 0);    mympz_randomm(y->a, curve->p);    x_from_y(x->a, y->a, curve);    fp2_set(P->x, x);    fp2_set(P->y, y);    fp2_clear(x);    fp2_clear(y);}void general_point_random(point_ptr P, curve_t curve)//P = random point on E/F_p^2{    fp2_t zeta;    point_t P2;    point_init(P2);    point_random(P, curve);    point_random(P2, curve);    fp2_init(zeta);    fp2_set_cbrt_unity(zeta, curve->p);        fp2_mul(P2->x, P2->x, zeta, curve->p);    point_add(P, P, P2, curve);    point_clear(P2);    fp2_clear(zeta);}void point_add(point_ptr R, point_ptr P, point_ptr Q, curve_t curve)//R = P + Q{    mpz_ptr p = curve->p;    fp2_t lambda, temp, temp2;    if (P->infinity) {	point_set(R, Q);	return;    }    if (Q->infinity) {	point_set(R, P);	return;    }    R->infinity = 0;    fp2_init(lambda);    fp2_init(temp);    fp2_init(temp2);    if (fp2_equal(P->x, Q->x)) { // Px == Py	fp2_neg(temp, Q->y, p);	if (fp2_equal(P->y, temp)) { // Py == -Qy	    point_set_O(R);	} else { //Py == Qy	    //line: Y - (lambda X + mu)	    //lambda = (x * (x + x + x + *twicea2) + *a4) / (y + y);	    //we assume twicea2 = 0, a4 = 0	    fp2_add(lambda, P->x, P->x, p);	    fp2_add(lambda, lambda, P->x, p);	    fp2_mul(lambda, lambda, P->x, p);	    fp2_add(temp, P->y, P->y, p);	    fp2_div(lambda, lambda, temp, p);	    //Rx = lambda^2 - 2Px	    fp2_set(temp, P->x); //in case &P = &R	    fp2_sqr(R->x, lambda, p);	    fp2_add(temp2, temp, temp, p);	    fp2_sub(R->x, R->x, temp2, p);	    //Ry = (Px - Rx) * lambda - Py	    fp2_sub(temp, temp, R->x, p);	    fp2_mul(temp, temp, lambda, p);	    fp2_sub(R->y, temp, P->y, p);	}    } else {	//line: Y - (lambda X + mu)	//lambda = (Qy - Py) / (Qx - Px);	fp2_sub(lambda, Q->y, P->y, p);	fp2_sub(temp, Q->x, P->x, p);	fp2_div(lambda, lambda, temp, p);	//Rx = lambda^2 - Px - Qx	fp2_set(temp, P->x); //in case &P = &R	fp2_sqr(temp2, lambda, p);	fp2_sub(temp2, temp2, temp, p);	fp2_sub(R->x, temp2, Q->x, p);	//Ry = (Px - Rx) * lambda - Py	fp2_sub(temp, temp, R->x, p);	fp2_mul(temp, temp, lambda, p);	fp2_sub(R->y, temp, P->y, p);    }    fp2_clear(lambda);    fp2_clear(temp);    fp2_clear(temp2);}static void proj_double(mpz_t x, mpz_t y, mpz_t z, mpz_t p)//(x, y, z) *= 2//see Blake, Seroussi & Smart, Fig IV.2//assumes (x, y, z) is not O, or a point of order 2 (i.e. y != 0)//we have a = 0 in our curve{    mpz_t t1, t2, t3, t4, t5;    mpz_init(t1);    mpz_init(t2);    mpz_init(t3);    mpz_init(t4);    mpz_init(t5);    //t1 = 3x^2    mpz_mul(t1, x, x);    mpz_add(t2, t1, t1);    mpz_add(t1, t1, t2);    mpz_mod(t1, t1, p);    //z' = 2yz    mpz_mul(z, z, y);    mpz_add(z, z, z);    mpz_mod(z, z, p);    //t2 = 4xy^2, t5 holds y^2    mpz_mul(t5, y, y);    mpz_mod(t5, t5, p);    mpz_mul(t2, t5, x);    mpz_mul_2exp(t2, t2, 2);    mpz_mod(t2, t2, p);    //x' = t1^2 - 2t2    mpz_mul(t3, t1, t1);    //mpz_mod(t3, t3, p);    mpz_add(t4, t2, t2);    mpz_sub(x, t3, t4);    mpz_mod(x, x, p);    //t3 = 8y^2 (recall t5 holds y^2)    mpz_mul(t3, t5, t5);    //mpz_mod(t3, t3, p);    mpz_mul_2exp(t3, t3, 3);    mpz_mod(t3, t3, p);    //y' = t1(t2 - x) - t3    mpz_sub(t4, t2, x);    mpz_mul(t4, t4, t1);    mpz_sub(y, t4, t3);    mpz_mod(y, y, p);    mpz_clear(t1);    mpz_clear(t2);    mpz_clear(t3);    mpz_clear(t4);    mpz_clear(t5);}static void proj_mix_in(mpz_t x, mpz_t y, mpz_t z, mpz_t a, mpz_t b, mpz_t p)//(x, y, z) += (a, b, 1)//assumes neither is O, and they are distinct points//for now also assume their sum is not O//see Blake, Seroussi & Smart, Fig IV.1{    //we take z_2 = 1, so t1 = x, t4 = y    mpz_t t2, t3, t5, t6, t7, t8;    mpz_init(t2);    mpz_init(t3);    mpz_init(t5);    mpz_init(t6);    mpz_init(t7);    mpz_init(t8);    //lambda_2 = x_2 * z_1^2    //t8 holds z^2 until t5 has been computed    mpz_mul(t8, z, z);    mpz_mod(t8, t8, p);    mpz_mul(t2, t8, a);    mpz_mod(t2, t2, p);    //lambda_3 = lambda_1 - lambda_2    mpz_sub(t3, x, t2);    //if (!mpz_size(t3)) {	//answer is O    //}    //lambda_5 = y_2 * z_1^3    mpz_mul(t5, t8, z);    mpz_mod(t5, t5, p);    mpz_mul(t5, t5, b);    mpz_mod(t5, t5, p);    //lambda_6 = lambda_4 - lambda_5    mpz_sub(t6, y, t5);    //lambda_7 = lambda_1 + lambda_2    mpz_add(t7, x, t2);    //lambda_8 = lambda_4 + lambda_5    mpz_add(t8, y, t5);    //z_3 = z_1 z_2 lambda_3    mpz_mul(z, z, t3);    mpz_mod(z, z, p);    //x_3 = lambda_6^2 - lambda_7 lambda_3^2    //t2, t5 no longer needed    //t2 holds t3^2    mpz_mul(t5, t6, t6);    mpz_mul(t2, t3, t3);    mpz_mod(t2, t2, p);    mpz_mul(x, t2, t7);    mpz_sub(x, t5, x);    mpz_mod(x, x, p);    //lambda_9 = lambda_7 lambda_3^2 - 2 x_3    //t5 doubles as t9    //t7 no longer needed after first line    mpz_mul(t5, t7, t2);    mpz_add(t7, x, x);    mpz_sub(t5, t5, t7);    mpz_mod(t5, t5, p);    //y_3 = (lambda_9 lambda_6 - lambda_8 lambda_3^3)/2    //t8 no longer needed after second line    mpz_mul(t7, t5, t6);    mpz_mul(t8, t8, t2);    mpz_mod(t8, t8, p);    mpz_mul(t8, t8, t3);    mpz_sub(y, t7, t8);    mpz_mod(y, y, p);    if (mpz_odd_p(y)) {	mpz_add(y, y, p);    }    mpz_fdiv_q_2exp(y, y, 1);    //is divexact better here?    mpz_clear(t2);    mpz_clear(t3);    mpz_clear(t5);    mpz_clear(t6);    mpz_clear(t7);    mpz_clear(t8);}static void tate_power(fp2_t res, curve_t curve){    fp2_t t0;    fp2_init(t0);    fp2_pow(t0, res, curve->p1onq, curve->p);    mpz_set(res->a, t0->a);    mpz_sub(res->b, curve->p, t0->b);    fp2_div(res, res, t0, curve->p);    fp2_clear(t0);}static void pts_get_vertical(fp2_ptr v,	point_ptr A, point_ptr P, mpz_ptr z, mpz_t p){    mpz_t z2;    fp2_t temp;    mpz_init(z2);    fp2_init(temp);    assert(!A->infinity);    assert(!P->infinity); // (could handle with a = b = 0; c = 1;)    //a = 1; b = 0; c = -P.x;    zp_mul(z2, z, z, p);    fp2_mul_mpz(temp, A->x, z2, p);    zp_sub(temp->a, temp->a, P->x->a, p);    fp2_mul(v, v, temp, p);    mpz_clear(z2);    fp2_clear(temp);}