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 gf32m_init(element_t e) {
    e->data = pbc_malloc(sizeof(gf32m_s));
    gf32m_ptr p = (gf32m_ptr) e->data;
    field_ptr base = BASE(e);
    element_init(p->_0, base);
    element_init(p->_1, base);
}

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

int main(void)
{
    field_t c;
    field_t Z19;
    element_t P, Q, R;
    mpz_t q, z;
    element_t a;
    int i;

    mpz_init(q);
    mpz_init(z);

    mpz_set_ui(q, 19);

    field_init_fp(Z19, q);
    element_init(a, Z19);

    field_init_curve_singular_with_node(c, Z19);

    element_init(P, c);
    element_init(Q, c);
    element_init(R, c);

    //(3,+/-6) is a generator
    //we have an isomorphism from E_ns to F_19^*
    // (3,6) --> 3
    //(generally (x,y) --> (y+x)/(y-x)

    curve_set_si(R, 3, 6);

    for (i=1; i<=18; i++) {
        mpz_set_si(z, i);
        element_mul_mpz(Q, R, z);
        element_printf("%dR = %B\n", i, Q);
    }

    mpz_set_ui(z, 6);
    element_mul_mpz(P, R, z);
    //P has order 3
    element_printf("P = %B\n", P);

    for (i=1; i<=3; i++) {
        mpz_set_si(z, i);
        element_mul_mpz(Q, R, z);
        tate_3(a, P, Q, R);
        element_printf("e_3(P,%dP) = %B\n", i, a);
    }

    element_double(P, R);
    //P has order 9
    element_printf("P = %B\n", P);
    for (i=1; i<=9; i++) {
        mpz_set_si(z, i);
        element_mul_mpz(Q, P, z);
        tate_9(a, P, Q, R);
        element_printf("e_9(P,%dP) = %B\n", i, a);
    }

    return 0;
}

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

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
}

/* computing $c <- U^M, M=(3^{3m}-1)*(3^m+1)*(3^m+1-\mu*b*3^{(m+1)//2})$
 * This is the algorithm 8 in the paper above. */
static void algorithm8(element_t c, element_t u) {
    field_ptr f6 = FIELD(u), f = FIELD(ITEM(u,0,0));
    params *p = (params *) f->data;
    element_t v, w;
    element_init(v, f6);
    element_init(w, f6);
    algorithm6(v, u);
    algorithm7(v, v);
    element_set(w, v);
    int i;
    for (i = 0; i < (p->m + 1) / 2; i++)
        element_cubic(w, w);
    algorithm7(v, v);
    if (p->m % 12 == 1 || p->m % 12 == 11) { // w <= w^{-\mu*b}
        element_ptr e;
        e = ITEM(w,0,1);
        element_neg(e, e);
        e = ITEM(w,1,1);
        element_neg(e, e);
        e = ITEM(w,2,1);
        element_neg(e, e);
    }
    element_mul(c, v, w);
    element_clear(v);
    element_clear(w);
}

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

//============================================
//   Frobenius Map \phi_p
//============================================
void test_frob(Field f)
{
    int i;
    unsigned long long int t1, t2;
    mpz_t p;
    Element a, b, c;

    mpz_init_set(p, *field_get_char(f));

    element_init(a, f);
    element_init(b, f);
    element_init(c, f);

    for (i = 0; i < 100; i++)
    {
        element_random(a);
        element_pow(b, a, p);

        bn254_fp2_frob_p(c, a);

        assert(element_cmp(b, c) == 0);
    }

    t1 = rdtsc();
    for (i = 0; i < N; i++) { bn254_fp2_frob_p(c, a); }
    t2 = rdtsc();

    printf("element frob: %.2lf [clock]\n", (double)(t2 - t1) / N);

    mpz_clear(p);

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

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

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

void test_list_remove_data (void)
{
  list *l = list_empty();
  assert(list_remove_head(l) == NULL);
  element *tmp = NULL;
  
  element* e1 = element_init("plop");
  element* e2 = element_init("plap");
  element* e3 = element_init("plup");
  
  list_add_head(l,e1);
  list_add_head(l,e2);
  assert(list_remove_data(l,"plup") == NULL);
  list_add_head(l,e3);
  
  tmp = list_remove_data(l,"plup");
  assert( tmp == e3);
  assert(l->head == e2);
  free(tmp);

  tmp = list_remove_data(l,"plop");
  assert( tmp == e1);
  assert(l->head == e2);
  assert(e2->next == NULL);
  
  free(tmp);
  list_free(l);
  
  printf("Les tests de list_remove_data sont réussis!\n");
}

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

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

// start_index从1开始
int BCEChangeDecryptionProduct(byte *global_params_path, int start_index, int length, int *adds, int n_adds, int *rems, int n_rems, byte *decr_prod, byte *decr_prod_out)
{
    global_broadcast_params_t gbs;
    struct single_priv_key_s *priv_key;
    int i, writelen = 0;

    if (global_params_path == NULL)
        return 1;
    if (start_index % NUM_USER_DIVISOR != 1)
        return 2;
    if (adds == NULL && rems == NULL)
        return 4;
    if (decr_prod == NULL)
        return 8;
    if (decr_prod_out == NULL)
        return 9;

    LoadGlobalParams((char *) global_params_path, &gbs);

    if (n_adds > gbs->num_users)
        return 5;
    if (n_rems > gbs->num_users)
        return 7;

    priv_key = (priv_key_t) malloc(length * sizeof(struct single_priv_key_s));

    for (i = 0; i < length; i++)
    {
        // restore index
        priv_key[i].index = start_index + i;

        // restore fake g_i_gamma
        element_init(priv_key[i].g_i_gamma, gbs->pairing->G1);

        // restore fake g_i
        element_init(priv_key[i].g_i, gbs->pairing->G1);

        // restore fake h_i
        element_init(priv_key[i].h_i, gbs->pairing->G2);

        // restore real decr_prod
        element_init(priv_key[i].decr_prod, gbs->pairing->G1);
        decr_prod += element_from_bytes(priv_key[i].decr_prod, decr_prod);
    }

    for (i = 0; i < length; i++)
    {
        Change_decr_prod_indicies(gbs, priv_key[i].index, adds, n_adds, rems, n_rems, &priv_key[i]);
        writelen += element_to_bytes(decr_prod_out + writelen, priv_key[i].decr_prod);
    }

    for (i = 0; i < length; i++)
        FreePK(&priv_key[i]);
    free(priv_key);
    FreeGBP(gbs);
    pbc_free(gbs);

    return 0;
}

static void point_init(element_t e) {
    field_ptr f = BASE(e);
    e->data = pbc_malloc(sizeof(struct point_s));
    point_ptr p = DATA(e);
    element_init(p->x, f);
    element_init(p->y, f);
    p->isinf = 1;
}

static void sn_init(element_ptr e) {
  field_ptr f = e->field->data;
  e->data = pbc_malloc(sizeof(point_t));
  point_ptr p = e->data;
  element_init(p->x, f);
  element_init(p->y, f);
  p->inf_flag = 1;
}

static void curve_init(element_ptr e) {
	curve_data_ptr cdp = (curve_data_ptr)e->field->data;
	point_ptr p; 
	e->data = pbc_malloc(sizeof(*p));
	p = (point_ptr)e->data;
  element_init(p->x, cdp->field);
  element_init(p->y, cdp->field);
  p->inf_flag = 1;
}