59.c

这是一个C的源代码

    element_mul(t1, pm2, s0);
    element_mul(t0, p0, sm2);
    element_sub(v0m2, t0, t1);
    element_mul(v0m2, v0m2, W20inv);

    element_mul(t0, p0, sm1);
    element_mul(t1, pm1, s0);
    element_sub(v0m1, t0, t1);

    element_mul(t1, pm1, s1);
    element_mul(t0, p1, sm1);
    element_sub(v00, t0, t1);
    element_mul(v00, v00, W20inv);

    element_mul(t0, p1, s0);
    element_mul(t1, p0, s1);
    element_sub(v01, t0, t1);

    element_mul(t1, p0, s2);
    element_mul(t0, p2, s0);
    element_sub(v02, t0, t1);
    element_mul(v02, v02, W20inv);

    element_mul(t0, p2, s1);
    element_mul(t1, p1, s2);
    element_sub(v03, t0, t1);

    element_mul(t1, p1, s3);
    element_mul(t0, p3, s1);
    element_sub(v04, t0, t1);
    element_mul(v04, v04, W20inv);

    element_square(t0, v10);
    element_mul(t1, v1m1, v11);

    element_mul(t2, pm1, t0);
    element_mul(v1m1, t1, sm1);
    element_sub(v1m1, v1m1, t2);

    element_mul(t2, p0, t0);
    element_mul(v10, t1, s0);
    element_sub(v10, v10, t2);

    element_mul(t2, p1, t0);
    element_mul(v11, t1, s1);
    element_sub(v11, v11, t2);
    element_mul(v11, v11, Wm11inv);

    element_printf("VEC1: %B %B %B\n", v1m1, v10, v11);
    element_printf("VEC0: %B %B %B %B %B %B %B %B\n",
	    v0m3, v0m2, v0m1, v00, v01, v02, v03, v04);

    //DoubleAdd
    element_square(sm2, v0m2);
    element_square(sm1, v0m1);
    element_square(s0, v00);
    element_square(s1, v01);
    element_square(s2, v02);
    element_square(s3, v03);

    element_mul(pm2, v0m3, v0m1);
    element_mul(pm1, v0m2, v00);
    element_mul(p0, v0m1, v01);
    element_mul(p1, v00, v02);
    element_mul(p2, v01, v03);
    element_mul(p3, v02, v04);

    element_mul(t1, pm2, s0);
    element_mul(t0, p0, sm2);
    element_sub(v0m3, t0, t1);
    element_mul(v0m3, v0m3, W20inv);

    element_mul(t0, p0, sm1);
    element_mul(t1, pm1, s0);
    element_sub(v0m2, t0, t1);

    element_mul(t1, pm1, s1);
    element_mul(t0, p1, sm1);
    element_sub(v0m1, t0, t1);
    element_mul(v0m1, v0m1, W20inv);

    element_mul(t0, p1, s0);
    element_mul(t1, p0, s1);
    element_sub(v00, t0, t1);

    element_mul(t1, p0, s2);
    element_mul(t0, p2, s0);
    element_sub(v01, t0, t1);
    element_mul(v01, v01, W20inv);

    element_mul(t0, p2, s1);
    element_mul(t1, p1, s2);
    element_sub(v02, t0, t1);

    element_mul(t1, p1, s3);
    element_mul(t0, p3, s1);
    element_sub(v03, t0, t1);
    element_mul(v03, v03, W20inv);

    element_mul(t0, p3, s2);
    element_mul(t1, p2, s3);
    element_sub(v04, t0, t1);

    element_square(t0, v10);
    element_mul(t1, v1m1, v11);

    element_mul(t2, p0, t0);
    element_mul(v1m1, t1, s0);
    element_sub(v1m1, v1m1, t2);

    element_mul(t2, p1, t0);
    element_mul(v10, t1, s1);
    element_sub(v10, v10, t2);
    element_mul(v10, v10, Wm11inv);

    element_mul(t2, t1, s2);
    element_mul(v11, p2, t0);
    element_sub(v11, v11, t2);
    element_mul(v11, v11, W2m1inv);

    element_printf("VEC1: %B %B %B\n", v1m1, v10, v11);
    element_printf("VEC0: %B %B %B %B %B %B %B %B\n",
	    v0m3, v0m2, v0m1, v00, v01, v02, v03, v04);
    element_div(z, v11, v01);
    element_printf("prepow: %B\n", z);

    element_pow_mpz(z, z, q1r);

    mpz_clear(q1r);
}

void miller(element_t z, element_t PR, element_t R, element_t P, element_t Q)
{
    int m = mpz_sizeinbase(order, 2) - 2;

    element_t Z;
    element_t z1;
    element_t x1;
    element_init_same_as(Z, PR);

    element_set(Z, P);
    element_set1(z);
    element_init_same_as(z1, z);
    element_init_same_as(x1, z);

    do_vert(x1, PR, Q);
    element_printf("vert(P+R) %B\n", x1);
    do_line(z1, P, R, Q);
    element_printf("line(P,R) %B\n", z1);
    element_div(x1, x1, z1);
    element_printf("x1 %B\n", x1);
    element_set(z, x1);

    for (;;) {
	printf("iteration %d: %d\n", m, mpz_tstbit(order,m));
	element_square(z, z);
	element_printf("squared: %B\n", z);
	do_tangent(z1, Z, Q);
	element_mul(z, z, z1);

	element_double(Z, Z);
	do_vert(z1, Z, Q);
	element_div(z, z, z1);
	element_printf("pre-if: %B\n", z);

	if (mpz_tstbit(order, m)) {
	    element_mul(z, z, x1);
	    do_vert(z1, P, Q);
	    element_mul(z, z, z1);
	    element_printf("done %B\n", z);
	    /*
	    do_line(z1, Z, P, Q);
	    element_mul(z, z, z1);
	    element_add(Z, Z, P);
	    do_vert(z1, Z, Q);
	    element_div(z, z, z1);
	    */
	}
	if (!m) break;
	m--;
    }

    element_clear(x1);
    element_clear(z1);
}
/*
*/

void weil(element_t w, element_t g, element_t h)
{
    element_t gr;
    element_t hs;
    element_t r;
    element_t s;
    element_t z, z0, z1;

    element_init(z, Fq2);
    element_init(z0, Fq2);
    element_init(z1, Fq2);

    element_init_same_as(gr, g);
    element_init_same_as(hs, h);
    element_init_same_as(r, g);
    element_init_same_as(s, h);

    element_random(r);
    element_random(s);
    //point_random always takes the same square root
    //why not take the other one for once?
    element_neg(r, r);
    element_set_str(r, "[[40,0],[54,0]]", 0);
    element_set_str(s, "[[48,55],[28,51]]", 0);

    element_printf("chose R = %B\n", r);
    element_printf("chose S = %B\n", s);
    element_add(gr, g, r);
    element_add(hs, h, s);

    element_printf("P+R = %B\n", gr);
    element_printf("Q+S = %B\n", hs);
    miller(z, gr, r, g, hs);
    miller(z0, gr, r, g, s);
    element_div(z1, z, z0);
    element_printf("num: %B\n", z1);

    miller(z, hs, s, h, gr);
    miller(z0, hs, s, h, r);
    element_div(w, z, z0);
    element_printf("denom: %B\n", w);

    element_div(w, z1, w);
}

void fasterweil(element_t w, element_t g, element_t h)
{
    element_t hs;
    element_t s;
    element_t z, z0, z1;

    element_init(z, Fq2);
    element_init(z0, Fq2);
    element_init(z1, Fq2);

    element_init_same_as(hs, h);
    element_init_same_as(s, h);

    element_random(s);
    //point_random always takes the same square root
    //why not take the other one for once?
    element_set_str(s, "[[48,55],[28,51]]", 0);

    element_printf("chose S = %B\n", s);
    element_add(hs, h, s);

    element_printf("Q+S = %B\n", hs);

    millertate(z, g, hs);
    millertate(z0, g, s);
    element_div(z1, z, z0);
    element_printf("num: %B\n", z1);

    miller(w, hs, s, h, g);
    element_printf("denom: %B\n", w);

    element_div(w, z1, w);
}

int main(void)
{
    int i;
    element_t g, h;
    element_t w0, w1;
    element_t a, b;
    mpz_t prime, cofac;

    mpz_init(prime);
    mpz_init(order);
    mpz_init(cofac);
    mpz_set_ui(prime, 59);

    field_init_fp(Fq, prime);

    element_init(a, Fq);
    element_init(b, Fq);

    field_init_fi(Fq2, Fq);

    element_set1(a);
    element_set0(b);
    mpz_set_ui(order, 5);
    mpz_set_ui(cofac, 12);

    field_init_curve_ab(E, a, b, order, cofac);

    element_clear(a);
    element_clear(b);
    element_init(a, Fq2);
    element_init(b, Fq2);
    element_set1(a);
    element_set0(b);

    mpz_mul(cofac, cofac, cofac);
    field_init_curve_ab(E2, a, b, order, NULL);

    element_init(g, E2);
    element_init(h, E2);

    element_init(w0, Fq2);
    element_init(w1, Fq2);

    /*
    do {
	element_random(g);
    } while (element_is1(g));
    for (i=1; i<5; i++) {
	element_mul(h, h, g);
	element_printf("%d: %B\n", i, h);
	element_printf("tangent = ");
	do_tangent(h);
    }
    */
    element_set_str(g, "[[25,0],[30,0]", 0);
    element_set_str(h, "[[34,0],[0,30]", 0);
    weil(w0, g, h);
    element_printf("weil: %B\n", w0);

    element_set1(w1);
    for (i=1; i<6; i++) {
	element_mul(w1, w1, w0);
	element_printf("%d: %B\n", i, w1);
    }

    fasterweil(w0, g, h);
    element_printf("fasterweil: %B\n", w0);

    element_set1(w1);
    for (i=1; i<6; i++) {
	element_mul(w1, w1, w0);
	element_printf("%d: %B\n", i, w1);
    }

    tate(w0, g, h);
    element_printf("tate: %B\n", w0);

    element_set1(w1);
    for (i=1; i<6; i++) {
	element_mul(w1, w1, w0);
	element_printf("%d: %B\n", i, w1);
    }

    shipseystange(w0, g, h);
    element_printf("ss-tate: %B\n", w0);

    element_set1(w1);
    for (i=1; i<6; i++) {
	element_mul(w1, w1, w0);
	element_printf("%d: %B\n", i, w1);
    }
    return 0;
}