static void do_vert(element_ptr z, element_ptr V, element_ptr Q)
{
    element_ptr Vx = curve_x_coord(V);
    element_ptr Qx = curve_x_coord(Q);
    element_ptr Qy = curve_y_coord(Q);

    element_t a, b, c;
    element_init_same_as(a, Vx);
    element_init_same_as(b, Vx);
    element_init_same_as(c, Vx);

    //a = 1
    //b = 0;
    //c = -Vx
    element_set1(a);
    element_set0(b);
    element_neg(c, Vx);

    element_printf("vert at %B: %B %B %B\n", Vx, a, b, c);
    element_mul(a, a, Qx);
    element_mul(b, b, Qy);
    element_add(c, c, a);
    element_add(z, c, b);
    element_printf("vert eval = %B\n", z);
    element_clear(a);
    element_clear(b);
    element_clear(c);
}

void miller(element_t res, mpz_t q, element_t P, element_ptr Qx, element_ptr Qy) {
  int m;
  element_t v;
  element_t Z;
  element_t a, b, c;
  element_t t0;
  element_t e0;
  const element_ptr cca = curve_a_coeff(P);
  const element_ptr Px = curve_x_coord(P);
  const element_ptr Py = curve_y_coord(P);
  element_ptr Zx, Zy;

  void do_tangent(void) {
    // a = -(3 Zx^2 + cc->a)
    // b = 2 * Zy
    // c = -(2 Zy^2 + a Zx);

    element_square(a, Zx); mult1++;
    element_mul_si(a, a, 3); add1++; add1++; add1++;
    element_add(a, a, cca); add1++;
    element_neg(a, a);

    element_add(b, Zy, Zy); add1++;

    element_mul(t0, b, Zy); mult1++;
    element_mul(c, a, Zx); mult1++;
    element_add(c, c, t0); add1++;
    element_neg(c, c);

    d_miller_evalfn(e0, a, b, c, Qx, Qy);
    element_mul(v, v, e0); multk++;
  }

void bb_sign(unsigned char *sig, unsigned int hashlen, unsigned char *hash, bb_public_key_t pk, bb_private_key_t sk)
{
	int len;
	element_t sigma;
	element_t r, z, m;
	bb_sys_param_ptr param = pk->param;
	pairing_ptr pairing = param->pairing;

	element_init(r, pairing->Zr);
	element_init(z, pairing->Zr);
	element_init(m, pairing->Zr);

	element_random(r);
	element_from_hash(m, hash, hashlen);
	element_mul(z, sk->y, r);
	element_add(z, z, sk->x);
	element_add(z, z, m);
	element_invert(z, z);
	element_init(sigma, pairing->G1);
	element_pow_zn(sigma, pk->g1, z);

	len = element_to_bytes_x_only(sig, sigma);
	element_to_bytes(&sig[len], r);

	element_clear(sigma);
	element_clear(r);
	element_clear(z);
	element_clear(m);
}

static void curve_from_hash(element_t a, void *data, int len) {
  element_t t, t1;
  point_ptr p = (point_ptr)a->data;
  curve_data_ptr cdp = (curve_data_ptr)a->field->data;

  element_init(t, cdp->field);
  element_init(t1, cdp->field);
  p->inf_flag = 0;
  element_from_hash(p->x, data, len);
  for(;;) {
    element_square(t, p->x);
    element_add(t, t, cdp->a);
    element_mul(t, t, p->x);
    element_add(t, t, cdp->b);
    if (element_is_sqr(t)) break;
    // Compute x <- x^2 + 1 and try again.
    element_square(p->x, p->x);
    element_set1(t);
    element_add(p->x, p->x, t);
  }
  element_sqrt(p->y, t);
  if (element_sgn(p->y) < 0) element_neg(p->y, p->y);

  if (cdp->cofac) element_mul_mpz(a, a, cdp->cofac);

  element_clear(t);
  element_clear(t1);
}

/* $e<- a*b$ */
static void gf32m_mult(element_t e, element_t a, element_t b) {
    element_ptr a0 = GF32M(a)->_0, a1 = GF32M(a)->_1, b0 = GF32M(b)->_0, b1 =
            GF32M(b)->_1, e0 = GF32M(e)->_0, e1 = GF32M(e)->_1;
    field_ptr base = BASE(a);
    element_t a0b0, a1b1, t0, t1, c1;
    element_init(a0b0, base);
    element_init(a1b1, base);
    element_init(t0, base);
    element_init(t1, base);
    element_init(c1, base);
    element_mul(a0b0, a0, b0);
    element_mul(a1b1, a1, b1);
    element_add(t0, a1, a0);
    element_add(t1, b1, b0);
    element_mul(c1, t0, t1); // c1 == (a1+a0)*(b1+b0)
    element_sub(c1, c1, a1b1);
    element_sub(c1, c1, a0b0);
    element_ptr c0 = a0b0;
    element_sub(c0, c0, a1b1); // c0 == a0*b0 - a1*b1
    element_set(e0, c0);
    element_set(e1, c1);
    element_clear(a0b0);
    element_clear(a1b1);
    element_clear(t0);
    element_clear(t1);
    element_clear(c1);
}

  void do_line(element_ptr e, element_ptr edenom, element_ptr A, element_ptr B)
  {
    element_ptr Ax = curve_x_coord(A);
    element_ptr Ay = curve_y_coord(A);
    element_ptr Bx = curve_x_coord(B);
    element_ptr By = curve_y_coord(B);

    element_sub(b, Bx, Ax);
    element_sub(a, Ay, By);
    element_mul(c, Ax, By);
    element_mul(e0, Ay, Bx);
    element_sub(c, c, e0);

    element_mul(e0, a, numx);
    element_mul(e1, b, numy);
    element_add(e0, e0, e1);
    element_add(e0, e0, c);
    element_mul(e, e, e0);

    element_mul(e0, a, denomx);
    element_mul(e1, b, denomy);
    element_add(e0, e0, e1);
    element_add(e0, e0, c);
    element_mul(edenom, edenom, e0);
  }

static void fq_add(element_ptr n, element_ptr a, element_ptr b) {
  eptr p = a->data;
  eptr q = b->data;
  eptr r = n->data;
  element_add(r->x, p->x, q->x);
  element_add(r->y, p->y, q->y);
}

static void fi_mul(element_ptr n, element_ptr a, element_ptr b) {
  eptr p = a->data;
  eptr q = b->data;
  eptr r = n->data;
  element_t e0, e1, e2;

  element_init(e0, p->x->field);
  element_init(e1, e0->field);
  element_init(e2, e0->field);
  /* Naive method:
  element_mul(e0, p->x, q->x);
  element_mul(e1, p->y, q->y);
  element_sub(e0, e0, e1);
  element_mul(e1, p->x, q->y);
  element_mul(e2, p->y, q->x);
  element_add(e1, e1, e2);
  element_set(r->x, e0);
  element_set(r->y, e1);
  */
  // Karatsuba multiplicaiton:
  element_add(e0, p->x, p->y);
  element_add(e1, q->x, q->y);
  element_mul(e2, e0, e1);
  element_mul(e0, p->x, q->x);
  element_sub(e2, e2, e0);
  element_mul(e1, p->y, q->y);
  element_sub(r->x, e0, e1);
  element_sub(r->y, e2, e1);

  element_clear(e0);
  element_clear(e1);
  element_clear(e2);
}

static void fi_sqrt(element_ptr n, element_ptr e) {
  eptr p = e->data;
  eptr r = n->data;
  element_t e0, e1, e2;

  // If (a+bi)^2 = x+yi then 2a^2 = x +- sqrt(x^2 + y^2)
  // where we choose the sign so that a exists, and 2ab = y.
  // Thus 2b^2 = - (x -+ sqrt(x^2 + y^2)).
  element_init(e0, p->x->field);
  element_init(e1, e0->field);
  element_init(e2, e0->field);
  element_square(e0, p->x);
  element_square(e1, p->y);
  element_add(e0, e0, e1);
  element_sqrt(e0, e0);
  // e0 = sqrt(x^2 + y^2)
  element_add(e1, p->x, e0);
  element_set_si(e2, 2);
  element_invert(e2, e2);
  element_mul(e1, e1, e2);
  // e1 = (x + sqrt(x^2 + y^2))/2
  if (!element_is_sqr(e1)) {
    element_sub(e1, e1, e0);
    // e1 should be a square.
  }
  element_sqrt(e0, e1);
  element_add(e1, e0, e0);
  element_invert(e1, e1);
  element_mul(r->y, p->y, e1);
  element_set(r->x, e0);
  element_clear(e0);
  element_clear(e1);
  element_clear(e2);
}

static void point_add(element_t c, element_t a, element_t b) {
    point_ptr p1 = DATA(a), p2 = DATA(b), p3 = DATA(c);
    int inf1 = p1->isinf, inf2 = p2->isinf;
    element_ptr x1 = p1->x, y1 = p1->y, x2 = p2->x, y2 = p2->y;
    field_ptr f = FIELD(x1);
    if (inf1) {
        point_set(c, b);
        return;
    }
    if (inf2) {
        point_set(c, a);
        return;
    }
    element_t v0, v1, v2, v3, v4, ny2;
    element_init(v0, f);
    element_init(v1, f);
    element_init(v2, f);
    element_init(v3, f);
    element_init(v4, f);
    element_init(ny2, f);
    if (!element_cmp(x1, x2)) { // x1 == x2
        element_neg(ny2, y2); // ny2 == -y2
        if (!element_cmp(y1, ny2)) {
            p3->isinf = 1;
            goto end;
        }
        if (!element_cmp(y1, y2)) { // y1 == y2
            element_invert(v0, y1); // v0 == y1^{-1}
            element_mul(v1, v0, v0); // v1 == [y1^{-1}]^2
            element_add(p3->x, v1, x1); // v1 == [y1^{-1}]^2 + x1
            element_cubic(v2, v0); // v2 == [y1^{-1}]^3
            element_add(v2, v2, y1); // v2 == [y1^{-1}]^3 + y1
            element_neg(p3->y, v2); // p3 == -([y1^{-1}]^3 + y1)
            p3->isinf = 0;
            goto end;
        }
    }
    // $P1 \ne \pm P2$
    element_sub(v0, x2, x1); // v0 == x2-x1
    element_invert(v1, v0); // v1 == (x2-x1)^{-1}
    element_sub(v0, y2, y1); // v0 == y2-y1
    element_mul(v2, v0, v1); // v2 == (y2-y1)/(x2-x1)
    element_mul(v3, v2, v2); // v3 == [(y2-y1)/(x2-x1)]^2
    element_cubic(v4, v2); // v4 == [(y2-y1)/(x2-x1)]^3
    element_add(v0, x1, x2); // v0 == x1+x2
    element_sub(v3, v3, v0); // v3 == [(y2-y1)/(x2-x1)]^2 - (x1+x2)
    element_add(v0, y1, y2); // v0 == y1+y2
    element_sub(v4, v0, v4); // v4 == (y1+y2) - [(y2-y1)/(x2-x1)]^3
    p3->isinf = 0;
    element_set(p3->x, v3);
    element_set(p3->y, v4);
end:
    element_clear(v0);
    element_clear(v1);
    element_clear(v2);
    element_clear(v3);
    element_clear(v4);
    element_clear(ny2);
}

/* $c <- a+b$ */
static void gf33m_add(element_t c, element_t a, element_t b) {
    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, b0 =
            GF33M(b)->_0, b1 = GF33M(b)->_1, b2 = GF33M(b)->_2, c0 =
            GF33M(c)->_0, c1 = GF33M(c)->_1, c2 = GF33M(c)->_2;
    element_add(c0, a0, b0);
    element_add(c1, a1, b1);
    element_add(c2, a2, b2);
}

/* $e <- a^{-1}$ */
static void gf33m_invert(element_t e, element_t a) {
    element_ptr a0 = GF33M(a)->_0, a1 = GF33M(a)->_1, a2 = GF33M(a)->_2, e0 =
            GF33M(e)->_0, e1 = GF33M(e)->_1, e2 = GF33M(e)->_2;
    field_ptr base = BASE(e);
    element_t a02, a12, a22;
    element_init(a02, base);
    element_init(a12, base);
    element_init(a22, base);
    element_mul(a02, a0, a0);
    element_mul(a12, a1, a1);
    element_mul(a22, a2, a2);
    element_t v0;
    element_init(v0, base);
    element_sub(v0, a0, a2); // v0 == a0-a2
    element_t delta;
    element_init(delta, base);
    element_mul(delta, v0, a02); // delta = (a0-a2)*(a0^2), free
    element_sub(v0, a1, a0); // v0 == a1-a0
    element_t c0;
    element_init(c0, base);
    element_mul(c0, v0, a12); // c0 == (a1-a0)*(a1^2)
    element_add(delta, delta, c0); // delta = (a0-a2)*(a0^2) + (a1-a0)*(a1^2)
    element_sub(v0, a2, v0); // v0 == a2-(a1-a0) = a0-a1+a2
    element_t c1;
    element_init(c1, base);
    element_mul(c1, v0, a22); // c1 == (a0-a1+a2)*(a2^2)
    element_add(delta, delta, c1); // delta = (a0-a2)*(a0^2) + (a1-a0)*(a1^2) + (a0-a1+a2)*(a2^2)
    element_invert(delta, delta); // delta = [(a0-a2)*(a0^2) + (a1-a0)*(a1^2) + (a0-a1+a2)*(a2^2)] ^ {-1}
    element_add(v0, a02, a22); // v0 == a0^2+a2^2
    element_t c2;
    element_init(c2, base);
    element_mul(c2, a0, a2); // c2 == a0*a2
    element_sub(c0, v0, c2); // c0 == a0^2+a2^2-a0*a2
    element_add(v0, a1, a2); // v0 == a1+a2
    element_t c3;
    element_init(c3, base);
    element_mul(c3, a1, v0); // c3 == a1*(a1+a2)
    element_sub(c0, c0, c3); // c0 == a0^2+a2^2-a0*a2-a1*(a1+a2)
    element_mul(c0, c0, delta); // c0 *= delta
    element_mul(c1, a0, a1); // c1 == a0*a1
    element_sub(c1, a22, c1); // c1 == a2^2-a0*a1
    element_mul(c1, c1, delta); // c1 *= delta
    element_sub(c2, a12, c2); // c2 == a1^2-a0*a2
    element_sub(c2, c2, a22); // c2 == a1^2-a0*a2-a2^2
    element_mul(c2, c2, delta); // c2 *= delta
    element_set(e0, c0);
    element_set(e1, c1);
    element_set(e2, c2);
    element_clear(a02);
    element_clear(a12);
    element_clear(a22);
    element_clear(v0);
    element_clear(delta);
    element_clear(c0);
    element_clear(c1);
    element_clear(c2);
    element_clear(c3);
}

static void do_line(element_ptr z, element_ptr V, element_ptr P, element_ptr Q)
{
    element_ptr Vx = curve_x_coord(V);
    element_ptr Vy = curve_y_coord(V);
    element_ptr Px = curve_x_coord(P);
    element_ptr Py = curve_y_coord(P);
    element_ptr Qx = curve_x_coord(Q);
    element_ptr Qy = curve_y_coord(Q);

    element_t a, b, c, e0;
    element_init_same_as(a, Vx);
    element_init_same_as(b, Vx);
    element_init_same_as(c, Vx);
    element_init_same_as(e0, Vx);

    //a = -(B.y - A.y) / (B.x - A.x);
    //b = 1;
    //c = -(A.y + a * A.x);

    element_sub(a, Py, Vy);
    element_sub(b, Vx, Px);
    element_div(a, a, b);
    element_set1(b);
    element_mul(c, a, Vx);
    element_add(c, c, Vy);
    element_neg(c, c);

    /*
    //but we could multiply by B.x - A.x to avoid division, so
    //a = -(By - Ay)
    //b = Bx - Ax
    //c = -(Ay b + a Ax);
    element_sub(a, Vy, Py);
    element_sub(b, Px, Vx);
    element_mul(c, Vx, Py);
    element_mul(e0, Vy, Px);
    element_sub(c, c, e0);
    //
    //actually no, since fasterweil won't work if we do this
    */

    element_printf("line at %B: %B %B %B\n", V, a, b, c);
    element_mul(a, a, Qx);
    element_mul(b, b, Qy);
    element_add(c, c, a);
    element_add(z, c, b);
    element_printf(" = %B\n", z);

    element_clear(a);
    element_clear(b);
    element_clear(c);
    element_clear(e0);
}

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

    element_clear(gr);
    element_clear(r);
    element_clear(hs);
    element_clear(s);
    element_clear(z);
    element_clear(z0);
    element_clear(z1);
}

static void test_gf3m_add(void) {
    element_random(a);
    element_add(b, a, a);
    element_add(b, b, b);
    element_sub(b, b, a);
    element_sub(b, b, a);
    element_sub(b, b, a);
    EXPECT(!element_cmp(a, b));

    element_add(b, params(a)->p, a);
    element_sub(b, b, params(a)->p);
    EXPECT(!element_cmp(a, b));
}

// Computes a point on the elliptic curve Y^2 = X^3 + a X + b given its
// x-coordinate.
// Requires a solution to exist.
static void point_from_x(point_ptr p, element_t x, element_t a, element_t b) {
  element_t t;

  element_init(t, x->field);
  p->inf_flag = 0;
  element_square(t, x);
  element_add(t, t, a);
  element_mul(t, t, x);
  element_add(t, t, b);
  element_sqrt(p->y, t);
  element_set(p->x, x);

  element_clear(t);
}

static void point_random(element_t a) {
    point_ptr p = DATA(a);
    element_ptr x = p->x, y = p->y;
    field_ptr f = x->field;
    p->isinf = 0;
    element_t t, t2, e1;
    element_init(t, f);
    element_init(e1, f);
    element_set1(e1);
    element_init(t2, f);
    do {
        element_random(x);
        if (element_is0(x))
            continue;
        element_cubic(t, x); // t == x^3
        element_sub(t, t, x); // t == x^3 - x
        element_add(t, t, e1); // t == x^3 - x + 1
        element_sqrt(y, t);  // y == sqrt(x^3 - x + 1)
        element_mul(t2, y, y); // t2 == x^3 - x + 1
    } while (element_cmp(t2, t)); // t2 != t

    // make sure order of $a$ is order of $G_1$
    pairing_ptr pairing = FIELD(a)->pairing;
    pairing_data_ptr dp = pairing->data;
    element_pow_mpz(a, a, dp->n2);

    element_clear(t);
    element_clear(t2);
    element_clear(e1);
}

void tate(element_t z, element_t P, element_t Q)
{
    mpz_t q1r;

    mpz_init(q1r);
    mpz_set_ui(q1r, 696);

    /*
    millertate(z, P, Q);
    element_printf("prepow: z = %B\n", z);
    element_pow_mpz(z, z, q1r);
    */
    {
	element_t R, QR;
	element_t z0;

	element_init_same_as(R, P);
	element_init_same_as(QR, P);
	element_init_same_as(z0, z);

	element_random(R);
	element_add(QR, Q, R);

	millertate(z, P, QR);
	millertate(z0, P, R);
	element_div(z, z, z0);
	element_pow_mpz(z, z, q1r);
	element_clear(R);
	element_clear(QR);
    }

    mpz_clear(q1r);
}

void CipherText::langrange(element_t* ys, int index, int k, int num){
  element_t delta;
  element_t numerator;
  element_t denominator;
  element_t temp;
  element_init_Zr(delta, *(this->p));
  element_init_Zr(numerator, *(this->p));
  element_init_Zr(denominator, *(this->p));
  element_init_Zr(temp, *(this->p));
  element_init_Zr(ys[index], *(this->p));
  element_set0(ys[index]);
  int i, j;
  for(i = 0; i < k; i++){
    //compute the langrange coefficent l
    element_set1(delta);
    for(j = 0; j < k; j++){
      if( j != i){
        element_set_si(numerator, index - j);
        element_set_si(denominator, i - j);
        element_div(numerator, numerator, denominator);
        element_mul(delta, delta, numerator);
      }
    }

    element_mul(temp, ys[i], delta);
    element_add(ys[index], ys[index], temp);
  }
}

static int fi_is_sqr(element_ptr e) {
  // x + yi is a square <=> x^2 + y^2 is (in the base field).

  // Proof: (=>) if x+yi = (a+bi)^2, then a^2 - b^2 = x, 2ab = y,
  // thus (a^2 + b^2)^2 = (a^2 - b^2)^2 + (2ab)^2 =  x^2 + y^2

  // (<=) Suppose A^2 = x^2 + y^2. If there exist a, b satisfying:
  //   a^2 = (+-A + x)/2, b^2 = (+-A - x)/2
  // then (a + bi)^2 = x + yi.
  //
  // We show that exactly one of (A + x)/2, (-A + x)/2 is a quadratic residue
  // (thus a, b do exist). Suppose not. Then the product (x^2 - A^2) / 4 is
  // some quadratic residue, a contradiction since this would imply x^2 - A^2 =
  // -y^2 is also a quadratic residue, but we know -1 is not a quadratic
  // residue. QED.
  eptr p = e->data;
  element_t e0, e1;
  int result;
  element_init(e0, p->x->field);
  element_init(e1, e0->field);
  element_square(e0, p->x);
  element_square(e1, p->y);
  element_add(e0, e0, e1);
  result = element_is_sqr(e0);
  element_clear(e0);
  element_clear(e1);
  return result;
}