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]);