enum {    See LICENSE for license    qbits = 10};static mpz_t p, q, p1onq;static mpz_t nq;static mpz_t moq;static mpz_t curve_a, curve_b;static gmp_randstate_t rand_state;static void evaluate_curve(mpz_t z, mpz_t x, mpz_t y){    mpz_t t0, t1;    mpz_init(t0); mpz_init(t1);    mpz_mul(t1, y, y);    mpz_mul(t0, x, x);    mpz_add(t0, t0, curve_a);    mpz_mul(t0, t0, x);    mpz_add(t0, t0, curve_b);    mpz_sub(t1, t1, t0);    mpz_mod(z, t1, p);    mpz_clear(t0); mpz_clear(t1);}static void point_add(mpz_t x3, mpz_t y3,	mpz_t x1, mpz_t y1,	mpz_t x2, mpz_t y2){    mpz_t t0, t1;    mpz_t l;    mpz_init(t0); mpz_init(t1);    mpz_init(l);    mpz_sub(t0, x2, x1);    if (!mpz_invert(t0, t0, p)) {	printf("divide by zero!\n");	exit(1);    }    mpz_sub(l, y2, y1);    mpz_mul(l, l, t0);    mpz_mod(l, l, p);    mpz_mul(t0, l, l);    mpz_sub(t0, t0, x1);    mpz_sub(t0, t0, x2);    mpz_mod(t0, t0, p);    mpz_sub(t1, x1, t0);    mpz_mul(t1, t1, l);    mpz_sub(t1, t1, y1);    mpz_set(x3, t0);    mpz_mod(y3, t1, p);    mpz_clear(l);    mpz_clear(t0); mpz_clear(t1);}static void point_double(mpz_t x3, mpz_t y3, mpz_t x1, mpz_t y1){    mpz_t t0, t1;    mpz_t l;    mpz_init(t0); mpz_init(t1);    mpz_init(l);    mpz_mul(l, x1, x1);    mpz_mul_ui(l, l, 3);    mpz_add(l, l, curve_a);    mpz_add(t0, y1, y1);    if (!mpz_invert(t0, t0, p)) {	printf("/0!\n");	exit(1);    }    mpz_mul(l, l, t0);    mpz_mod(l, l, p);    mpz_mul(t0, l, l);    mpz_sub(t0, t0, x1);    mpz_sub(t0, t0, x1);    mpz_mod(t0, t0, p);    mpz_sub(t1, x1, t0);    mpz_mul(t1, t1, l);    mpz_sub(t1, t1, y1);    mpz_set(x3, t0);    mpz_mod(y3, t1, p);    mpz_clear(l);    mpz_clear(t0); mpz_clear(t1);}static void point_mul(mpz_t x3, mpz_t y3, mpz_t a, mpz_t x1, mpz_t y1){    int j;    mpz_t t0, t1;    mpz_init(t0); mpz_init(t1);    j = mpz_sizeinbase(a, 2) - 1;    mpz_set(t0, x1);    mpz_set(t1, y1);    j--;    while(j >= 0) {	point_double(t0, t1, t0, t1);	if (mpz_tstbit(a, j)) {	    point_add(t0, t1, t0, t1, x1, y1);	}	j--;    }    mpz_set(x3, t0);    mpz_set(y3, t1);    mpz_clear(t0); mpz_clear(t1);}static void tonelli(mpz_t x, mpz_t a){    int i, m;    int r, e;    mpz_t qq;    mpz_t y;    mpz_t b;    mpz_t t;    mpz_init(qq);    mpz_init(y);    mpz_init(b);    mpz_init(t);    mpz_sub_ui(qq, p, 1);    e = mpz_scan1(qq, 0);    mpz_tdiv_q_2exp (qq, qq, e);    mpz_powm(y, nq, qq, p);    r = e;    //(q-1)/2    mpz_tdiv_q_2exp(qq, qq, 1);    mpz_powm(x, a, qq, p);    mpz_mul(b, x, x);    mpz_mul(b, b, a);    mpz_mod(b, b, p);    mpz_mul(x, x, a);    mpz_mod(x, x, p);    while (mpz_cmp_ui(b, 1)) {	m = 1;	mpz_set(t, b);	for (;;) {	    mpz_mul(t, t, t);	    mpz_mod(t, t, p);	    if (!mpz_cmp_ui(t, 1)) break;	    m++;	}	mpz_set(t, y);	for (i=0; i<r-m-1; i++) {	    mpz_mul(t, t, t);	    mpz_mod(t, t, p);	}	mpz_mul(y, t, t);	mpz_mod(y, y, p);	r = m;	mpz_mul(x, x, t);	mpz_mod(x, x, p);	mpz_mul(b, b, y);	mpz_mod(b, b, p);    }    mpz_clear(qq);    mpz_clear(y);    mpz_clear(b);    mpz_clear(t);}static void point_random(mpz_t x, mpz_t y){    mpz_t t0;    mpz_init(t0);    for (;;) {	mpz_urandomm(x, rand_state, p);	mpz_mul(t0, x, x);	mpz_add(t0, t0, curve_a);	mpz_mul(t0, t0, x);	mpz_add(t0, t0, curve_b);	if (mpz_legendre(t0, p) == 1) break;    }    tonelli(y, t0);    mpz_clear(t0);}void make_order_q(mpz_t x, mpz_t y){    point_mul(x, y, moq, x, y);}static void find_curve(){    mpz_t x, y;    mpz_t x2, y2;    mpz_t j, c;    mpz_t n;    mpz_init(j); mpz_init(c);    mpz_init(x); mpz_init(y);    mpz_init(x2); mpz_init(y2);    mpz_init(n);    mpz_set_ui(n, 1);    mpz_mul_2exp(n, n, qbits);    do {	mpz_urandomm(q, rand_state, n);	mpz_nextprime(q, q);	//printf("q: "); mpz_out_str(stdout, 10, q); printf("\n");	mpz_mul(p, q, q);	mpz_mul_ui(p, p, 28);	mpz_add_ui(p, p, 1);	//printf("p: "); mpz_out_str(stdout, 10, p); printf("\n");    } while (!mpz_probab_prime_p(p, 10));    printf("p: "); mpz_out_str(stdout, 10, p); printf("\n");    printf("q: "); mpz_out_str(stdout, 10, q); printf("\n");    mpz_sub_ui(p1onq, p, 1);    mpz_divexact(p1onq, p1onq, q);    mpz_set_ui(moq, 28);    mpz_init(nq);    for(;;) {	mpz_urandomm(nq, rand_state, p);	if (mpz_legendre(nq, p) == -1) break;    }    //printf("nq = "); mpz_out_str(stdout, 10, nq); printf("\n");    mpz_sub_ui(j, p, 3375);    mpz_sub_ui(c, j, 1728);    mpz_invert(c, c, p);    mpz_mul(c, c, j);    mpz_mod(c, c, p);    mpz_mul_si(curve_a, c, -3);    mpz_mod(curve_a, curve_a, p);    mpz_mul_ui(curve_b, c, 2);    mpz_mod(curve_b, curve_b, p);    //mpz_set_ui(p, 59);    //mpz_set_ui(curve_a, 0);    //mpz_set_ui(curve_b, 1);    printf("Y^2 = X^3 + ");    mpz_out_str(stdout, 10, curve_a);    printf("X + ");    mpz_out_str(stdout, 10, curve_b);    printf("\n");    for(;;) {	mpz_t z;	mpz_init(z);	point_random(x2, y2);	make_order_q(x2, y2);	mpz_sub_ui(z, q, 1);	point_mul(x, y, z, x2, y2);//printf("(x, y) = ("); mpz_out_str(stdout, 10, x); printf(", ");//mpz_out_str(stdout, 10, y); printf(")\n");//printf("(x, y) = ("); mpz_out_str(stdout, 10, x2); printf(", ");//mpz_out_str(stdout, 10, y2); printf(")\n");	if (mpz_cmp(x, x2)) {	    //should be order p + 3 instead, unless my conjecture's wrong	    mpz_add_ui(z, p, 2);	    point_mul(x, y, z, x2, y2);	    if (mpz_cmp(x, x2)) {		printf("COUNTEREXAMPLE!\n");		exit(0);	    }	    printf("wrong curve. trying quadratic twist...\n");	    mpz_mul(curve_a, curve_a, nq);	    mpz_mul(curve_a, curve_a, nq);	    mpz_mul(curve_b, curve_b, nq);	    mpz_mul(curve_b, curve_b, nq);	    mpz_mul(curve_b, curve_b, nq);	} else break; 	mpz_clear(z);    }    printf("E[q] lies in F_p\n");    mpz_clear(j); mpz_clear(c);    mpz_clear(x); mpz_clear(y);    mpz_clear(x2); mpz_clear(y2);    mpz_clear(n);}static void get_line(mpz_t v, mpz_t vd, mpz_t Qx, mpz_t Qy,	mpz_t Ax, mpz_t Ay, mpz_t Bx, mpz_t By){    mpz_t a, b, c;    mpz_t t0, t1;    mpz_init(a); mpz_init(b); mpz_init(c);    mpz_init(t0); mpz_init(t1);    //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    mpz_sub(b, Bx, Ax);    mpz_mod(b, b, p);    mpz_sub(a, Ay, By);    mpz_mod(a, a, p);    mpz_mul(t0, b, Ay);    mpz_mul(c, a, Ax);    mpz_add(c, c, t0);    mpz_neg(c, c);    mpz_mod(c, c, p);    mpz_mul(t0, a, Qx);    mpz_mul(t1, b, Qy);    mpz_add(t0, t0, t1);    mpz_add(t0, t0, c);    mpz_mod(t0, t0, p);    mpz_mul(v, v, t0);    mpz_mod(v, v, p);    mpz_mul(vd, vd, b);    mpz_mod(vd, vd, p);    mpz_clear(a); mpz_clear(b); mpz_clear(c);    mpz_clear(t0); mpz_clear(t1);}static void get_tangent(mpz_t v, mpz_t vd, mpz_t Qx, mpz_t Qy,	mpz_t Ax, mpz_t Ay){    mpz_t a, b, c;    mpz_t t0, t1;    mpz_init(a); mpz_init(b); mpz_init(c);    mpz_init(t0); mpz_init(t1);    //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 = 0 (a4 = curve_a)    mpz_add(a, Ax, Ax);    mpz_add(a, a, Ax);    mpz_mul(a, Ax, Ax);    mpz_neg(a, a);    mpz_sub(a, a, curve_a);    mpz_mod(a, a, p);    mpz_add(b, Ay, Ay);    mpz_mod(b, b, p);    mpz_mul(t0, b, Ay);    mpz_mul(c, a, Ax);    mpz_add(c, c, t0);    mpz_neg(c, c);    mpz_mod(c, c, p);    mpz_mul(t0, a, Qx);    mpz_mul(t1, b, Qy);    mpz_add(t0, t0, t1);    mpz_add(t0, t0, c);    mpz_mul(v, v, t0);    mpz_mod(v, v, p);    mpz_mul(vd, vd, b);    mpz_mod(vd, vd, p);    mpz_clear(a); mpz_clear(b); mpz_clear(c);    mpz_clear(t0); mpz_clear(t1);}static void get_vertical(mpz_t v, mpz_t vd, mpz_t Qx, mpz_t Ax){    mpz_t t0;    mpz_init(t0);    mpz_sub(t0, Qx, Ax);    mpz_mul(v, v, t0);    mpz_mod(v, v, p);    mpz_clear(t0);}void miller(mpz_ptr res, mpz_t Px, mpz_t Py, mpz_t Qx, mpz_t Qy)//compute res = f_P(Q)//assumes q is odd{    //f_1 = 1    int m;    mpz_t v, vd;    mpz_t Zx, Zy;    mpz_init(v); mpz_init(vd);    mpz_init(Zx); mpz_init(Zy);    mpz_set_ui(v, 1);    mpz_set_ui(vd, 1);    m = mpz_sizeinbase(q, 2) - 1;    mpz_set(Zx, Px);    mpz_set(Zy, Py);    m--;    get_tangent(v, vd, Qx, Qy, Zx, Zy);    point_double(Zx, Zy, Zx, Zy);    get_vertical(vd, v, Qx, Zx);    for(;;) {	//printf("v: "); mpz_out_str(stdout, 10, v); printf("\n");	//printf("vd: "); mpz_out_str(stdout, 10, vd); printf("\n");	if (mpz_tstbit(q, m)) {	    if (!m) {		get_vertical(v, vd, Qx, Zx);		break;	    }	    get_line(v, vd, Qx, Qy, Zx, Zy, Px, Py);	    point_add(Zx, Zy, Zx, Zy, Px, Py);	    get_vertical(vd, v, Qx, Zx);	}	mpz_mul(v, v, v);	mpz_mod(v, v, p);	mpz_mul(vd, vd, vd);	mpz_mod(vd, vd, p);	get_tangent(v, vd, Qx, Qy, Zx, Zy);	point_double(Zx, Zy, Zx, Zy);	get_vertical(vd, v, Qx, Zx);	m--;    }    //printf("v: "); mpz_out_str(stdout, 10, v); printf("\n");    //printf("vd: "); mpz_out_str(stdout, 10, vd); printf("\n");    mpz_invert(vd, vd, p);    mpz_mul(res, v, vd);    mpz_mod(res, res, p);    mpz_clear(v); mpz_clear(vd);    mpz_clear(Zx); mpz_clear(Zy);}void pairing(mpz_t r, mpz_t Px, mpz_t Py, mpz_t Qx, mpz_t Qy){    miller(r, Px, Py, Qx, Qy);    mpz_powm(r, r, p1onq, p);}void test_pairing(){    mpz_t Px, Py, Qx, Qy, Rx, Ry;    mpz_t r0, r1, r2;    mpz_init(Px); mpz_init(Py);    mpz_init(Qx); mpz_init(Qy);    mpz_init(Rx); mpz_init(Ry);    mpz_init(r0); mpz_init(r1); mpz_init(r2);    point_random(Px, Py);    make_order_q(Px, Py);printf("P = ("); mpz_out_str(stdout, 10, Px); printf(", ");mpz_out_str(stdout, 10, Py); printf(")\n");    point_random(Qx, Qy);    make_order_q(Qx, Qy);printf("Q = ("); mpz_out_str(stdout, 10, Qx); printf(", ");mpz_out_str(stdout, 10, Qy); printf(")\n");    point_random(Rx, Ry);    make_order_q(Rx, Ry);printf("R = ("); mpz_out_str(stdout, 10, Rx); printf(", ");mpz_out_str(stdout, 10, Ry); printf(")\n");    pairing(r0, Px, Py, Rx, Ry);printf("e(P, R) = "); mpz_out_str(stdout, 10, r0); printf("\n");    pairing(r1, Qx, Qy, Rx, Ry);printf("e(Q, R) = "); mpz_out_str(stdout, 10, r1); printf("\n");    point_add(Px, Py, Px, Py, Qx, Qy);    pairing(r2, Px, Py, Rx, Ry);printf("e(P + Q, R) = "); mpz_out_str(stdout, 10, r2); printf("\n");    mpz_mul(r0, r0, r1);    mpz_mod(r0, r0, p);printf("e(P, R)e(Q, R) = "); mpz_out_str(stdout, 10, r0); printf("\n");    mpz_clear(Px); mpz_clear(Py);    mpz_clear(Qx); mpz_clear(Qy);    mpz_clear(Rx); mpz_clear(Ry);    mpz_clear(r0); mpz_clear(r1); mpz_clear(r2);}int main(){    int i;    gmp_randinit_default(rand_state);    mpz_init(p); mpz_init(q);    mpz_init(p1onq);    mpz_init(moq);    mpz_init(curve_a); mpz_init(curve_b);    for (i=0; i<100; i++)    {	find_curve();	test_pairing();    }    mpz_clear(p); mpz_clear(q);    mpz_clear(p1onq);    mpz_clear(curve_a); mpz_clear(curve_b);    mpz_clear(nq);    mpz_clear(moq);    return(0);}