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

void cc_pairings_affine(element_ptr out, element_t in1[], element_t in2[],
        int n_prod, pairing_t pairing) {
  element_ptr Qbase;
  element_t* Qx = pbc_malloc(sizeof(element_t)*n_prod);
  element_t* Qy = pbc_malloc(sizeof(element_t)*n_prod);
  pptr p = pairing->data;
  int i;
  for(i=0; i<n_prod; i++){
          element_init(Qx[i], p->Fqd);
          element_init(Qy[i], p->Fqd);
        Qbase = in2[i];
          // Twist: (x, y) --> (v^-1 x, v^-(3/2) y)
          // where v is the quadratic nonresidue used to construct the twist.
          element_mul(Qx[i], curve_x_coord(Qbase), p->nqrinv);
          // v^-3/2 = v^-2 * v^1/2
          element_mul(Qy[i], curve_y_coord(Qbase), p->nqrinv2);
  }
  cc_millers_no_denom_affine(out, pairing->r, in1, Qx, Qy, n_prod);
  cc_tatepower(out, out, pairing);

  for(i=0; i<n_prod; i++){
          element_clear(Qx[i]);
                element_clear(Qy[i]);
  }
  pbc_free(Qx);
  pbc_free(Qy);
}

static void f_tateexp(element_t out) {
  element_t x, y, epow;
  f_pairing_data_ptr p = out->field->pairing->data;
  element_init(x, p->Fq12);
  element_init(y, p->Fq12);
  element_init(epow, p->Fq2);

  #define qpower(e1, e) {                                         \
    element_set(element_item(e1, 0), element_item(out, 0));       \
    element_mul(element_item(e1, 1), element_item(out, 1), e);    \
    element_square(epow, e);                                      \
    element_mul(element_item(e1, 2), element_item(out, 2), epow); \
    element_mul(epow, epow, e);                                   \
    element_mul(element_item(e1, 3), element_item(out, 3), epow); \
    element_mul(epow, epow, e);                                   \
    element_mul(element_item(e1, 4), element_item(out, 4), epow); \
    element_mul(epow, epow, e);                                   \
    element_mul(element_item(e1, 5), element_item(out, 5), epow); \
  }

  qpower(y, p->xpowq8);
  qpower(x, p->xpowq6);
  element_mul(y, y, x);
  qpower(x, p->xpowq2);
  element_mul(x, x, out);
  element_invert(x, x);
  element_mul(out, y, x);

  element_clear(epow);
  element_clear(x);
  element_clear(y);
  element_pow_mpz(out, out, p->tateexp);
  #undef qpower
}

static void test_gf3m_cubic(void) {
    element_random(a);
    element_mul(b, a, a);
    element_mul(b, a, b);
    element_cubic(a, a);
    EXPECT(!element_cmp(a, b));
}

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++;
  }

static void test_gf32m_cubic(void) {
    element_random(a2);
    element_mul(b2, a2, a2);
    element_mul(b2, b2, a2);
    element_cubic(a2, a2);
    EXPECT(!element_cmp(a2, b2));
}

static void test_gf33m_cubic(void) {
    element_random(a3);
    element_mul(b3, a3, a3);
    element_mul(b3, b3, a3);
    element_cubic(a3, a3);
    EXPECT(!element_cmp(a3, b3));
}

// The final powering, where we standardize the coset representative.
static void cc_tatepower(element_ptr out, element_ptr in, pairing_t pairing) {
  pptr p = pairing->data;
  #define qpower(sign) {                         \
    polymod_const_mul(e2, inre[1], p->xpowq);    \
    element_set(e0re, e2);                       \
    polymod_const_mul(e2, inre[2], p->xpowq2);   \
    element_add(e0re, e0re, e2);                 \
    element_add(e0re0, e0re0, inre[0]);          \
                                                 \
    if (sign > 0) {                              \
      polymod_const_mul(e2, inim[1], p->xpowq);  \
      element_set(e0im, e2);                     \
      polymod_const_mul(e2, inim[2], p->xpowq2); \
      element_add(e0im, e0im, e2);               \
      element_add(e0im0, e0im0, inim[0]);        \
    } else {                                     \
      polymod_const_mul(e2, inim[1], p->xpowq);  \
      element_neg(e0im, e2);                     \
      polymod_const_mul(e2, inim[2], p->xpowq2); \
      element_sub(e0im, e0im, e2);               \
      element_sub(e0im0, e0im0, inim[0]);        \
    }                                            \
  }
  if (p->k == 6) {
    // See thesis, section 6.9, "The Final Powering", which gives a formula
    // for the first step of the final powering when Fq6 has been implemented
    // as a quadratic extension on top of a cubic extension.
    element_t e0, e2, e3;
    element_init(e0, p->Fqk);
    element_init(e2, p->Fqd);
    element_init(e3, p->Fqk);
    element_ptr e0re = element_x(e0);
    element_ptr e0im = element_y(e0);
    element_ptr e0re0 = ((element_t *) e0re->data)[0];
    element_ptr e0im0 = ((element_t *) e0im->data)[0];
    element_t *inre = element_x(in)->data;
    element_t *inim = element_y(in)->data;
    // Expressions in the formula are similar, hence the following function.
    qpower(1);
    element_set(e3, e0);
    element_set(e0re, element_x(in));
    element_neg(e0im, element_y(in));
    element_mul(e3, e3, e0);
    qpower(-1);
    element_mul(e0, e0, in);
    element_invert(e0, e0);
    element_mul(in, e3, e0);

    element_set(e0, in);
    // We use Lucas sequences to complete the final powering.
    lucas_even(out, e0, pairing->phikonr);

    element_clear(e0);
    element_clear(e2);
    element_clear(e3);
  } else {
    element_pow_mpz(out, in, p->tateexp);
  }
  #undef qpower
}

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

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

/* $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);
}

static void e_miller_proj(element_t res, element_t P,
    element_ptr QR, element_ptr R,
    e_pairing_data_ptr p) {
  //collate divisions
  int n;
  element_t v, vd;
  element_t v1, vd1;
  element_t Z, Z1;
  element_t a, b, c;
  const element_ptr cca = curve_a_coeff(P);
  element_t e0, e1;
  const element_ptr e2 = a, e3 = b;
  element_t z, z2;
  int i;
  element_ptr Zx, Zy;
  const element_ptr Px = curve_x_coord(P);
  const element_ptr numx = curve_x_coord(QR);
  const element_ptr numy = curve_y_coord(QR);
  const element_ptr denomx = curve_x_coord(R);
  const element_ptr denomy = curve_y_coord(R);

  //convert Z from weighted projective (Jacobian) to affine
  //i.e. (X, Y, Z) --> (X/Z^2, Y/Z^3)
  //also sets z to 1
  void to_affine(void)
  {
    element_invert(z, z);
    element_square(e0, z);
    element_mul(Zx, Zx, e0);
    element_mul(e0, e0, z);
    element_mul(Zy, Zy, e0);
    element_set1(z);
    element_set1(z2);
  }

static void test_gf36m_cubic(void) {
    element_random(a6);
    element_mul(b6, a6, a6);
    element_mul(b6, b6, a6);
    element_cubic(a6, a6);
    EXPECT(!element_cmp(a6, b6));
}

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_iとk_iをペアリングする関数
//¬記号で別々の処理する
//(v_i - x_t)も必要→とりあえず置いておこう……→一応できた？
Element *pairing_c_k(EC_PAIRING p, rho_i *rho, EC_POINT *c, EC_POINT *k, mpz_t *alpha_i) {
	int i;
	Element *result;
	result = (Element*)malloc(sizeof(Element));
	Element egg, tempegg1, tempegg2;
	element_init(egg, p->g3);
	element_init(tempegg1, p->g3);
	element_init(tempegg2, p->g3);
	element_init(*result, p->g3);
	mpz_t temp1;
	mpz_init(temp1);
	mpz_t temp2;
	mpz_init(temp2);
	mpz_t order;
	mpz_init(order);
	mpz_set(order, *pairing_get_order(p));
	element_set_one(*result);

	if (alpha_i == NULL && rho == NULL) { //e(c_0, k_0)
			for (i = 0; i < 5; i++) {
				pairing_map(tempegg1, c[i], k[i], p);
				element_mul(tempegg2, tempegg1, *result);
				element_set(*result, tempegg2);
			}
	}
    else if (mpz_cmp_ui(*alpha_i, 0) == 0) {//return 1
    }
	else if (rho->is_negated == FALSE) {
			for (i = 0; i < 7; i++) {
				pairing_map(tempegg1, c[i], k[i], p);
				element_mul(tempegg2, tempegg1, *result);
				element_set(*result, tempegg2);
			}
			element_pow(tempegg1, *result, *alpha_i);
			element_set(*result, tempegg1);
	}
	else { //is_negated == TRUE
			for (i = 0; i < 7; i++) {
				pairing_map(tempegg1, c[i], k[i], p);
				element_mul(tempegg2, tempegg1, *result);
				element_set(*result, tempegg2);
			}
		mpz_set_ui(temp1, rho->v_t[0]); //v_i - x_t
		mpz_invert(temp2, temp1, order);
		mpz_mul(temp1, temp2, *alpha_i); // alpha_i / (v_i - x_t)
		mpz_mod(*alpha_i, temp1, order);
		element_pow(tempegg1, *result, *alpha_i);
		element_set(*result, tempegg1);
	}

    mpz_clear(order);
    mpz_clear(temp2);
    mpz_clear(temp1);
	element_clear(egg);
	element_clear(tempegg1);
	element_clear(tempegg2);

	return result;
}

static void sn_add(element_t c, element_t a, element_t b) {
  point_ptr r = c->data;
  point_ptr p = a->data;
  point_ptr q = b->data;
  if (p->inf_flag) {
    sn_set(c, b);
    return;
  }
  if (q->inf_flag) {
    sn_set(c, a);
    return;
  }
  if (!element_cmp(p->x, q->x)) {
    if (!element_cmp(p->y, q->y)) {
      if (element_is0(p->y)) {
        r->inf_flag = 1;
        return;
      } else {
        sn_double_no_check(r, p);
        return;
      }
    }
    //points are inverses of each other
    r->inf_flag = 1;
    return;
  } else {
    element_t lambda, e0, e1;

    element_init(lambda, p->x->field);
    element_init(e0, p->x->field);
    element_init(e1, p->x->field);

    //lambda = (y2-y1)/(x2-x1)
    element_sub(e0, q->x, p->x);
    element_invert(e0, e0);
    element_sub(lambda, q->y, p->y);
    element_mul(lambda, lambda, e0);
    //x3 = lambda^2 - x1 - x2 - 1
    element_square(e0, lambda);
    element_sub(e0, e0, p->x);
    element_sub(e0, e0, q->x);
    element_set1(e1);
    element_sub(e0, e0, e1);
    //y3 = (x1-x3)lambda - y1
    element_sub(e1, p->x, e0);
    element_mul(e1, e1, lambda);
    element_sub(e1, e1, p->y);

    element_set(r->x, e0);
    element_set(r->y, e1);
    r->inf_flag = 0;

    element_clear(lambda);
    element_clear(e0);
    element_clear(e1);
  }
}

//--------------------------------------------------
//  square root in extended Fp
//--------------------------------------------------
int bn254_fp2_sqrt(Element z, const Element x)
{
    mpz_t _v;
    int m, r, i;

    field_precomp_sqrt_p ps;

    Element *t = field(z)->tmp;

    if (!element_is_sqr(x)) {
        return FALSE;
    }

    ps = ((field_precomp_p)(field(x)->precomp))->ps;

    element_set(t[0], ps->n_v);  // t0 = n^v

    r = ps->e; // r = e

    mpz_init_set(_v, ps->v);
    mpz_sub_ui(_v, _v, 1);
    mpz_tdiv_q_2exp(_v, _v, 1);

    element_pow(t[1], x, _v);   // t1 = x^{(v-1)/2}
    element_sqr(t[2], t[1]);
    element_mul(t[2], t[2], x); // t2 = x*t1^2
    element_mul(t[1], x, t[1]); // t1 = x*t1

    mpz_clear(_v);

    while (!element_is_one(t[2]))
    {
        m = 0;
        element_set(t[3], t[2]);

        do {
            element_sqr(t[3], t[3]);
            m++;
        }
        while (!element_is_one(t[3]) && m < r);

        r = r - m - 1;
        element_set(t[3], t[0]);
        for (i = 1; i <= r; i++) {
            element_sqr(t[3], t[3]);
        } // t3 = t2^{r-m-1}
        element_sqr(t[0], t[3]);      // t0 = t3^2
        r = m;
        element_mul(t[1], t[1], t[3]);// t1 = t1*t3
        element_mul(t[2], t[2], t[0]);// t2 = t2*t0
    }

    element_set(z, t[1]);

    return TRUE;
}

  void do_vertical(element_t e, element_t edenom, element_ptr Ax)
  {
    element_mul(e0, numx, z2);
    element_sub(e0, e0, Ax);
    element_mul(e, e, e0);

    element_mul(e0, denomx, z2);
    element_sub(e0, e0, Ax);
    element_mul(edenom, edenom, e0);
  }

static void curve_mul(element_ptr c, element_ptr a, element_ptr b) {
	curve_data_ptr cdp = (curve_data_ptr)a->field->data;
  point_ptr r = (point_ptr)c->data, p = (point_ptr)a->data, q = (point_ptr)b->data;

  if (p->inf_flag) {
    curve_set(c, b);
    return;
  }
  if (q->inf_flag) {
    curve_set(c, a);
    return;
  }
  if (!element_cmp(p->x, q->x)) {
    if (!element_cmp(p->y, q->y)) {
      if (element_is0(p->y)) {
        r->inf_flag = 1;
        return;
      } else {
        double_no_check(r, p, cdp->a);
        return;
      }
    }
    //points are inverses of each other
    r->inf_flag = 1;
    return;
  } else {
    element_t lambda, e0, e1;

    element_init(lambda, cdp->field);
    element_init(e0, cdp->field);
    element_init(e1, cdp->field);

    //lambda = (y2-y1)/(x2-x1)
    element_sub(e0, q->x, p->x);
    element_invert(e0, e0);
    element_sub(lambda, q->y, p->y);
    element_mul(lambda, lambda, e0);
    //x3 = lambda^2 - x1 - x2
    element_square(e0, lambda);
    element_sub(e0, e0, p->x);
    element_sub(e0, e0, q->x);
    //y3 = (x1-x3)lambda - y1
    element_sub(e1, p->x, e0);
    element_mul(e1, e1, lambda);
    element_sub(e1, e1, p->y);

    element_set(r->x, e0);
    element_set(r->y, e1);
    r->inf_flag = 0;

    element_clear(lambda);
    element_clear(e0);
    element_clear(e1);
  }
}