fp2_mul(g[j], g[k], g[2], p);	fp2_mul(g[j], g[j], v, p);	fp2_div(g[j], g[j], vdenom, p);	fp2_set_1(v);	fp2_set_1(vdenom);    }    point_init(Z);    point_set_O(Z);    m = mpz_sizeinbase(curve->q, 2) - 1;    while(m>=0) {	if (!mpz_tstbit(curve->q, m)) {	    fp2_sqr(v, v, p);	    fp2_sqr(vdenom, vdenom, p);	    get_tangent(v, vdenom, Qhat, R2, Z, p);	    point_add(Z, Z, Z, curve);	    get_vertical(v, vdenom, R2, Qhat, Z, p);	    m--;	} else {	    l = m - windowsize + 1;	    if (l < 0) l = 0;	    l = mpz_scan1(curve->q, l);	    j = 1;	    fp2_sqr(v, v, p);	    fp2_sqr(vdenom, vdenom, p);	    get_tangent(v, vdenom, Qhat, R2, Z, p);	    point_add(Z, Z, Z, curve);	    get_vertical(v, vdenom, R2, Qhat, Z, p);	    for (k = m - 1; k>=l; k--) {		j = j << 1;		if (mpz_tstbit(curve->q, k)) j++;		fp2_sqr(v, v, p);		fp2_sqr(vdenom, vdenom, p);		get_tangent(v, vdenom, Qhat, R2, Z, p);		point_add(Z, Z, Z, curve);		get_vertical(v, vdenom, R2, Qhat, Z, p);	    }	    fp2_mul(v, v, g[j], p);	    get_line(v, vdenom, Qhat, R2, Z, bP[j], p);	    point_add(Z, Z, bP[j], curve);	    get_vertical(v, vdenom, R2, Qhat, Z, p);	    m = l-1;	}    }    if (fp2_is_0(v) || fp2_is_0(vdenom)) goto bad_R;    else fp2_div(res, v, vdenom, p);    // check with simple algorithm    if (0) {	fp2_t check;	fp2_init(check);	simple_miller(check, P, Phat, Qhat, R1, R2, curve);	if (!fp2_equal(check, res)) {	    printf("res check = ");	    fp2_out_str(NULL, 0, res);	    printf(" ");	    fp2_out_str(NULL, 0, check);	    printf("\n");	    fprintf(stderr, "BUG: miller() is incorrectly implemented\n");	    //fp2_set(res, check);	}	fp2_clear(check);    }    fp2_clear(v);    fp2_clear(vdenom);    point_clear(Z);    fp2_clear(g[2]);    point_clear(bP[2]);    for (i=0; i<=windowsizepower; i++) {	j = 2*i + 1;	point_clear(bP[j]);	fp2_clear(g[j]);    }    return 1;bad_R:    printf("bad R\n");    fp2_clear(v);    fp2_clear(vdenom);    /* TODO: memory leak if jumped here after g[] inited    point_clear(Z);    fp2_clear(g[2]);    point_clear(bP[2]);    for (i=0; i<=windowsizepower; i++) {	j = 2*i + 1;	point_clear(bP[j]);	fp2_clear(g[j]);    }    */    return 0;}static void pts_preprocess_vertical(miller_cache_t mc,	int i, point_t P, mpz_t z, mpz_t p){    mpz_t t0;    mpz_init(t0);    mpz_mul(t0, z, z);    mpz_invert(t0, t0, p);    mpz_mul(t0, t0, P->x->a);    mpz_neg(t0, t0);    mpz_mod(mc->denomc[i], t0, p);    mpz_clear(t0);}void pts_preprocess_tangent(miller_cache_t mc,	int i, point_t P, mpz_t z, mpz_t p){    mpz_t t0;    mpz_t *numa, *numc;    mpz_init(t0);    numa = mc->numa;    numc = mc->numc;    assert(mpz_cmp_ui(P->y->a, 0)); //assume P not of order 2    //could handle with:    //{    //	mpz_set_ui(numa[i], 1);    //	mpz_set(numc[i], P->x->a);    //	return;    //}    //a = -3x^2/2yz    //c = -(y + a * xz)/z^3;    mpz_mul_2exp(t0, P->y->a, 1);    mpz_mul(t0, t0, z);    mpz_invert(t0, t0, p);    mpz_mul(t0, t0, P->x->a);    mpz_mul(t0, t0, P->x->a);    mpz_mul_si(t0, t0, -3);    mpz_mod(numa[i], t0, p);    mpz_mul(t0, z, z);    mpz_mul(t0, t0, z);    mpz_invert(t0, t0, p);    mpz_mul(numc[i], numa[i], P->x->a);    mpz_mul(numc[i], numc[i], z);    mpz_add(numc[i], numc[i], P->y->a);    mpz_neg(numc[i], numc[i]);    mpz_mul(t0, numc[i], t0);    mpz_mod(numc[i], t0, p);    mpz_clear(t0);}void pts_preprocess_line(mpz_t a, mpz_t c, point_t P, point_t Q, mpz_t p){    mpz_t t0;    mpz_init(t0);    //assume P.x != Q.x    //a = -(Q.y - P.y) / (Q.x - P.x);    //c = -(P.y + a * P.x);    mpz_sub(c, Q->x->a, P->x->a);    mpz_invert(c, c, p);    mpz_sub(a, P->y->a, Q->y->a);    mpz_mul(a, a, c);    mpz_mod(a, a, p);    mpz_mul(c, a, P->x->a);    mpz_add(c, c, P->y->a);    mpz_neg(c, c);    mpz_mod(c, c, p);    mpz_clear(t0);}void tate_preprocess(miller_cache_t mc, point_ptr P, curve_t curve)//for primes of the form 2^a +- 2^b +- 1//uses proj. coords, assumes P is a point over F_p//and that order of group = Solinas prime{    //specialized for Solinas primes    int a, b;    point_t bP;    point_t Z;    int i;    mpz_t z, temp;    mpz_ptr p = curve->p;    mpz_init(z);    mpz_init(temp);    mpz_set_ui(z, 1);    a = abs(curve->solinasa);    b = abs(curve->solinasb);    point_init(Z);    point_init(bP);    point_set(Z, P);    i = 0;    if (b != 0) {	//work out f_2^b	for(;i<b; i++) {	    //g	    //tate_get_tangent(v, Qhat, Z, p);	    //pts_get_tangent(v, Qhat, Z, z);	    pts_preprocess_tangent(mc, i, Z, z, p);	    //h	    //point_add(Z, Z, Z, curve);	    proj_double(Z->x->a, Z->y->a, z, p);	    //tate_get_vertical(vdenom, Qhat, Z, p);	    //pts_get_vertical(vdenom, Qhat, Z, z, p);	    pts_preprocess_vertical(mc, i, Z, z, p);	}	mpz_invert(z, z, p);	zp_mul(temp, z, z, p);	fp2_mul_mpz(Z->x, Z->x, temp, p);	zp_mul(temp, temp, z, p);	fp2_mul_mpz(Z->y, Z->y, temp, p);	mpz_set_ui(z, 1);	point_set(bP, Z);	if (curve->solinasb < 0) {	    //tate_get_vertical(fbdenom, Qhat, Z);	    mpz_neg(mc->denomsb, Z->x->a);	    fp2_neg(bP->y, bP->y, p);	}    }    //work out f_2^a    for(; i<a; i++) {	//g	//pts_get_tangent(v, Qhat, Z, z);	pts_preprocess_tangent(mc, i, Z, z, p);	//h	proj_double(Z->x->a, Z->y->a, z, p);	//pts_get_vertical(vdenom, Qhat, Z, z);	pts_preprocess_vertical(mc, i, Z, z, p);    }    mpz_invert(z, z, p);    zp_mul(temp, z, z, p);    fp2_mul_mpz(Z->x, Z->x, temp, p);    zp_mul(temp, temp, z, p);    fp2_mul_mpz(Z->y, Z->y, temp, p);    mpz_set_ui(z, 1);    //work out f_(2^a +- 2^b +- 1)    if (b != 0) {	//g	//tate_get_line(v, Qhat, Z, bP);	pts_preprocess_line(mc->numl1a, mc->numl1c, Z, bP, p);	//h	point_add(Z, Z, bP, curve);	//tate_get_vertical(vdenom, Qhat, Z);	mpz_sub(mc->denoml1c, p, Z->x->a);    }    //the sign of solinasa records whether it's +1 or -1    if (curve->solinasa < 0) {	//tate_get_vertical(vdenom, Qhat, Z);	mpz_sub(mc->denoms1, p, Z->x->a);    }    //g    //tate_get_line(v, Qhat, Z, cP);    //now Z = -cP so g is vertical    mpz_sub(mc->numl2c, p, Z->x->a);    //h    //point_add(Z, Z, cP, curve);    //now Z = O, so h = 1    //tate_get_vertical(vdenom, Qhat, Z);    point_clear(Z);    point_clear(bP);    mpz_clear(z);    mpz_clear(temp);}void miller_postprocess(fp2_ptr res, miller_cache_t mc,	point_ptr Q, curve_t curve)//for primes of the form 2^a +- 2^b +- 1//uses proj. coords, assumes P is a point over F_p//and that order of group = Solinas prime{    //specialized for Solinas primes    int a, b;    fp2_t fb, fbdenom;    fp2_t t0;    fp2_t vdenom, v;    int i;    mpz_t *numa, *numc;    mpz_t *denomc;    mpz_ptr p = curve->p;    numa = mc->numa;    numc = mc->numc;    denomc = mc->denomc;    a = abs(curve->solinasa);    b = abs(curve->solinasb);    fp2_init(t0);    fp2_init(vdenom);    fp2_init(v);    fp2_init(fb);    fp2_init(fbdenom);    fp2_set_1(v);    fp2_set_1(vdenom);    //f_1 = 1    i = 0;    if (b != 0) {	//work out f_2^b	for(;i<b; i++) {	    fp2_sqr(v, v, p);	    fp2_sqr(vdenom, vdenom, p);	    //g	    //tate_get_tangent(v, Qhat, Z, p);	    //pts_get_tangent(v, Qhat, Z, z, p);	    {		fp2_mul_mpz(t0, Q->x, numa[i], p);		fp2_add(t0, t0, Q->y, p);		mpz_add(t0->a, t0->a, numc[i]);		fp2_mul(v, v, t0, p);	    }	    //h	    //tate_get_vertical(vdenom, Qhat, Z, p);	    //pts_get_vertical(vdenom, Qhat, Z, z, p);	    {		fp2_set(t0, Q->x);		mpz_add(t0->a, t0->a, denomc[i]);		fp2_mul(vdenom, vdenom, t0, p);	    }	}	if (curve->solinasb < 0) {	    fp2_set(fbdenom, v);	    fp2_set(fb, vdenom);	    //tate_get_vertical(fbdenom, Qhat, Z, p);	    {		fp2_set(t0, Q->x);		mpz_add(t0->a, t0->a, mc->denomsb);		fp2_mul(fbdenom, fbdenom, t0, p);	    }	} else {	    fp2_set(fb, v);	    fp2_set(fbdenom, vdenom);	}    }    //work out f_2^a    for(; i<a; i++) {	fp2_sqr(v, v, p);	fp2_sqr(vdenom, vdenom, p);	//g	//pts_get_tangent(v, Qhat, Z, z, p);	{	    fp2_mul_mpz(t0, Q->x, numa[i], p);	    fp2_add(t0, t0, Q->y, p);	    mpz_add(t0->a, t0->a, numc[i]);	    fp2_mul(v, v, t0, p);	}	//h	//pts_get_vertical(vdenom, Qhat, Z, z, p);	{	    fp2_set(t0, Q->x);	    mpz_add(t0->a, t0->a, denomc[i]);	    fp2_mul(vdenom, vdenom, t0, p);	}    }    //work out f_(2^a +- 2^b +- 1)    if (b != 0) {	fp2_mul(v, v, fb, p);	fp2_mul(vdenom, vdenom, fbdenom, p);	//g	//tate_get_line(v, Qhat, Z, bP, p);	{	    fp2_mul_mpz(t0, Q->x, mc->numl1a, p);	    fp2_add(t0, t0, Q->y, p);	    mpz_add(t0->a, t0->a, mc->numl1c);	    fp2_mul(v, v, t0, p);	}	//h	//tate_get_vertical(vdenom, Qhat, Z, p);	{	    fp2_set(t0, Q->x);	    mpz_add(t0->a, t0->a, mc->denoml1c);	    fp2_mul(vdenom, vdenom, t0, p);	}    }    if (curve->solinasa < 0) {    //the sign of solinasa records whether it's +1 or -1	//tate_get_vertical(vdenom, Qhat, P);	{	    fp2_set(t0, Q->x);	    mpz_add(t0->a, t0->a, mc->denoms1);	    fp2_mul(vdenom, vdenom, t0, p);	}    }    //g    //tate_get_line(v, Qhat, Z, cP);    {	fp2_set(t0, Q->x);	mpz_add(t0->a, t0->a, mc->numl2c);	fp2_mul(v, v, t0, p);    }    //h    //tate_get_vertical(vdenom, Qhat, Z);    fp2_div(res, v, vdenom, p);    fp2_clear(v);    fp2_clear(vdenom);    fp2_clear(fb);    fp2_clear(fbdenom);    fp2_clear(t0);}void tate_postprocess(fp2_ptr res, miller_cache_t mc, point_ptr Q, curve_t curve){    bm_put(bm_get_time(), "miller0");    miller_postprocess(res, mc, Q, curve);    bm_put(bm_get_time(), "miller1");    tate_power(res, curve);}void tate_solinas_miller(fp2_ptr res, point_ptr P, point_ptr Qhat, curve_t curve)//for primes of the form 2^a +- 2^b +- 1//uses proj. coords, assumes P is a point over F_p//and that order of group = Solinas prime{    //specialized for Solinas primes    int a, b;    fp2_t fb, fbdenom;    point_t bP;    fp2_t vdenom, v;    point_t Z;    int i;    mpz_t z, temp;    mpz_ptr p = curve->p;    mpz_init(z);    mpz_init(temp);    mpz_set_ui(z, 1);    a = abs(curve->solinasa);    b = abs(curve->solinasb);    point_init(Z);    point_init(bP);    fp2_init(vdenom);    fp2_init(v);    fp2_init(fb);    fp2_init(fbdenom);    fp2_set_1(v);    fp2_set_1(vdenom);    //f_1 = 1    point_set(Z, P);    i = 0;    if (b != 0) {	//work out f_2^b	for(;i<b; i++) {	    fp2_sqr(v, v, p);	    fp2_sqr(vdenom, vdenom, p);	    //g	    //tate_get_tangent(v, Qhat, Z);	    pts_get_tangent(v, Qhat, Z, z, p);	    //h	    //point_add(Z, Z, Z, curve);	    proj_double(Z->x->a, Z->y->a, z, p);	    //tate_get_vertical(vdenom, Qhat, Z);	    pts_get_vertical(vdenom, Qhat, Z, z, p);	}	mpz_invert(z, z, p);	zp_mul(temp, z, z, p);	fp2_mul_mpz(Z->x, Z->x, temp, p);	zp_mul(temp, temp, z, p);	fp2_mul_mpz(Z->y, Z->y, temp, p);	mpz_set_ui(z, 1);	point_set(bP, Z);	if (curve->solinasb < 0) {	    fp2_set(fbdenom, v);	    fp2_set(fb, vdenom);	    tate_get_vertical(fbdenom, Qhat, Z, p);	    fp2_neg(bP->y, bP->y, p);	} else {	    fp2_set(fb, v);	    fp2_set(fbdenom, vdenom);	}    }    //work out f_2^a    for(; i<a; i++) {	fp2_sqr(v, v, p);	fp2_sqr(vdenom, vdenom, p);	//g	pts_get_tangent(v, Qhat, Z, z, p);	//h	proj_double(Z->x->a, Z->y->a, z, p);	pts_get_vertical(vdenom, Qhat, Z, z, p);    }    mpz_invert(z, z, p);    zp_mul(temp, z, z, p);    fp2_mul_mpz(Z->x, Z->x, temp, p);    zp_mul(temp, temp, z, p);    fp2_mul_mpz(Z->y, Z->y, temp, p);    mpz_set_ui(z, 1);    //work out f_(2^a +- 2^b +- 1)