edwards_pairing.cpp

Go to the documentation of this file.

 
#include <cassert>
#include <libff/algebra/curves/edwards/edwards_g1.hpp>
#include <libff/algebra/curves/edwards/edwards_g2.hpp>
#include <libff/algebra/curves/edwards/edwards_init.hpp>
#include <libff/algebra/curves/edwards/edwards_pairing.hpp>
#include <libff/common/profiling.hpp>
 
namespace libff
{
 
bool edwards_Fq_conic_coefficients::operator==(
    const edwards_Fq_conic_coefficients &other) const
{
    return (
        this->c_ZZ == other.c_ZZ && this->c_XY == other.c_XY &&
        this->c_XZ == other.c_XZ);
}
 
std::ostream &operator<<(
    std::ostream &out, const edwards_Fq_conic_coefficients &cc)
{
    out << cc.c_ZZ << OUTPUT_SEPARATOR << cc.c_XY << OUTPUT_SEPARATOR
        << cc.c_XZ;
    return out;
}
 
std::istream &operator>>(std::istream &in, edwards_Fq_conic_coefficients &cc)
{
    in >> cc.c_ZZ;
    consume_OUTPUT_SEPARATOR(in);
    in >> cc.c_XY;
    consume_OUTPUT_SEPARATOR(in);
    in >> cc.c_XZ;
    return in;
}
 
std::ostream &operator<<(
    std::ostream &out, const edwards_tate_G1_precomp &prec_P)
{
    out << prec_P.size() << "\n";
    for (const edwards_Fq_conic_coefficients &cc : prec_P) {
        out << cc << OUTPUT_NEWLINE;
    }
 
    return out;
}
 
std::istream &operator>>(std::istream &in, edwards_tate_G1_precomp &prec_P)
{
    prec_P.clear();
 
    size_t s;
    in >> s;
 
    consume_newline(in);
    prec_P.reserve(s);
 
    for (size_t i = 0; i < s; ++i) {
        edwards_Fq_conic_coefficients cc;
        in >> cc;
        consume_OUTPUT_NEWLINE(in);
        prec_P.emplace_back(cc);
    }
 
    return in;
}
 
bool edwards_tate_G2_precomp::operator==(
    const edwards_tate_G2_precomp &other) const
{
    return (this->y0 == other.y0 && this->eta == other.eta);
}
 
std::ostream &operator<<(
    std::ostream &out, const edwards_tate_G2_precomp &prec_Q)
{
    out << prec_Q.y0 << OUTPUT_SEPARATOR << prec_Q.eta;
    return out;
}
 
std::istream &operator>>(std::istream &in, edwards_tate_G2_precomp &prec_Q)
{
    in >> prec_Q.y0;
    consume_OUTPUT_SEPARATOR(in);
    in >> prec_Q.eta;
    return in;
}
 
bool edwards_Fq3_conic_coefficients::operator==(
    const edwards_Fq3_conic_coefficients &other) const
{
    return (
        this->c_ZZ == other.c_ZZ && this->c_XY == other.c_XY &&
        this->c_XZ == other.c_XZ);
}
 
std::ostream &operator<<(
    std::ostream &out, const edwards_Fq3_conic_coefficients &cc)
{
    out << cc.c_ZZ << OUTPUT_SEPARATOR << cc.c_XY << OUTPUT_SEPARATOR
        << cc.c_XZ;
    return out;
}
 
std::istream &operator>>(std::istream &in, edwards_Fq3_conic_coefficients &cc)
{
    in >> cc.c_ZZ;
    consume_OUTPUT_SEPARATOR(in);
    in >> cc.c_XY;
    consume_OUTPUT_SEPARATOR(in);
    in >> cc.c_XZ;
    return in;
}
 
std::ostream &operator<<(
    std::ostream &out, const edwards_ate_G2_precomp &prec_Q)
{
    out << prec_Q.size() << "\n";
    for (const edwards_Fq3_conic_coefficients &cc : prec_Q) {
        out << cc << OUTPUT_NEWLINE;
    }
 
    return out;
}
 
std::istream &operator>>(std::istream &in, edwards_ate_G2_precomp &prec_Q)
{
    prec_Q.clear();
 
    size_t s;
    in >> s;
 
    consume_newline(in);
 
    prec_Q.reserve(s);
 
    for (size_t i = 0; i < s; ++i) {
        edwards_Fq3_conic_coefficients cc;
        in >> cc;
        consume_OUTPUT_NEWLINE(in);
        prec_Q.emplace_back(cc);
    }
 
    return in;
}
 
bool edwards_ate_G1_precomp::operator==(
    const edwards_ate_G1_precomp &other) const
{
    return (
        this->P_XY == other.P_XY && this->P_XZ == other.P_XZ &&
        this->P_ZZplusYZ == other.P_ZZplusYZ);
}
 
std::ostream &operator<<(
    std::ostream &out, const edwards_ate_G1_precomp &prec_P)
{
    out << prec_P.P_XY << OUTPUT_SEPARATOR << prec_P.P_XZ << OUTPUT_SEPARATOR
        << prec_P.P_ZZplusYZ;
 
    return out;
}
 
std::istream &operator>>(std::istream &in, edwards_ate_G1_precomp &prec_P)
{
    in >> prec_P.P_XY >> prec_P.P_XZ >> prec_P.P_ZZplusYZ;
 
    return in;
}
 
/* final exponentiations */
edwards_Fq6 edwards_final_exponentiation_last_chunk(
    const edwards_Fq6 &elt, const edwards_Fq6 &elt_inv)
{
    enter_block("Call to edwards_final_exponentiation_last_chunk");
    const edwards_Fq6 elt_q = elt.Frobenius_map(1);
    edwards_Fq6 w1_part =
        elt_q.cyclotomic_exp(edwards_final_exponent_last_chunk_w1);
    edwards_Fq6 w0_part;
    if (edwards_final_exponent_last_chunk_is_w0_neg) {
        w0_part =
            elt_inv.cyclotomic_exp(edwards_final_exponent_last_chunk_abs_of_w0);
    } else {
        w0_part =
            elt.cyclotomic_exp(edwards_final_exponent_last_chunk_abs_of_w0);
    }
    edwards_Fq6 result = w1_part * w0_part;
    leave_block("Call to edwards_final_exponentiation_last_chunk");
 
    return result;
}
 
edwards_Fq6 edwards_final_exponentiation_first_chunk(
    const edwards_Fq6 &elt, const edwards_Fq6 &elt_inv)
{
    enter_block("Call to edwards_final_exponentiation_first_chunk");
 
    /* (q^3-1)*(q+1) */
 
    /* elt_q3 = elt^(q^3) */
    const edwards_Fq6 elt_q3 = elt.Frobenius_map(3);
    /* elt_q3_over_elt = elt^(q^3-1) */
    const edwards_Fq6 elt_q3_over_elt = elt_q3 * elt_inv;
    /* alpha = elt^((q^3-1) * q) */
    const edwards_Fq6 alpha = elt_q3_over_elt.Frobenius_map(1);
    /* beta = elt^((q^3-1)*(q+1) */
    const edwards_Fq6 beta = alpha * elt_q3_over_elt;
    leave_block("Call to edwards_final_exponentiation_first_chunk");
    return beta;
}
 
edwards_GT edwards_final_exponentiation(const edwards_Fq6 &elt)
{
    enter_block("Call to edwards_final_exponentiation");
    const edwards_Fq6 elt_inv = elt.inverse();
    const edwards_Fq6 elt_to_first_chunk =
        edwards_final_exponentiation_first_chunk(elt, elt_inv);
    const edwards_Fq6 elt_inv_to_first_chunk =
        edwards_final_exponentiation_first_chunk(elt_inv, elt);
    edwards_GT result = edwards_final_exponentiation_last_chunk(
        elt_to_first_chunk, elt_inv_to_first_chunk);
    leave_block("Call to edwards_final_exponentiation");
 
    return result;
}
 
edwards_tate_G2_precomp edwards_tate_precompute_G2(const edwards_G2 &Q)
{
    enter_block("Call to edwards_tate_precompute_G2");
    edwards_G2 Qcopy = Q;
    Qcopy.to_affine_coordinates();
    edwards_tate_G2_precomp result;
    result.y0 = Qcopy.Y * Qcopy.Z.inverse(); // Y/Z
    result.eta =
        (Qcopy.Z + Qcopy.Y) *
        edwards_Fq6::mul_by_non_residue(Qcopy.X).inverse(); // (Z+Y)/(nqr*X)
    leave_block("Call to edwards_tate_precompute_G2");
 
    return result;
}
 
struct extended_edwards_G1_projective {
    edwards_Fq X;
    edwards_Fq Y;
    edwards_Fq Z;
    edwards_Fq T;
 
    void print() const
    {
        printf("extended edwards_G1 projective X/Y/Z/T:\n");
        X.print();
        Y.print();
        Z.print();
        T.print();
    }
 
    void test_invariant() const { assert(T * Z == X * Y); }
};
 
void doubling_step_for_miller_loop(
    extended_edwards_G1_projective &current, edwards_Fq_conic_coefficients &cc)
{
    const edwards_Fq &X = current.X, &Y = current.Y, &Z = current.Z,
                     &T = current.T;
    const edwards_Fq A = X.squared();       // A    = X1^2
    const edwards_Fq B = Y.squared();       // B    = Y1^2
    const edwards_Fq C = Z.squared();       // C    = Z1^2
    const edwards_Fq D = (X + Y).squared(); // D    = (X1+Y1)^2
    const edwards_Fq E = (Y + Z).squared(); // E    = (Y1+Z1)^2
    const edwards_Fq F = D - (A + B);       // F    = D-(A+B)
    const edwards_Fq G = E - (B + C);       // G    = E-(B+C)
    const edwards_Fq &H = A;                // H    = A (edwards_a=1)
    const edwards_Fq I = H + B;             // I    = H+B
    const edwards_Fq J = C - I;             // J    = C-I
    const edwards_Fq K = J + C;             // K    = J+C
 
    cc.c_ZZ = Y * (T - X); // c_ZZ = 2*Y1*(T1-X1)
    cc.c_ZZ = cc.c_ZZ + cc.c_ZZ;
 
    cc.c_XY = J + J + G; // c_XY = 2*J+G
    cc.c_XZ = X * T - B; // c_XZ = 2*(X1*T1-B) (edwards_a=1)
    cc.c_XZ = cc.c_XZ + cc.c_XZ;
 
    current.X = F * K;       // X3 = F*K
    current.Y = I * (B - H); // Y3 = I*(B-H)
    current.Z = I * K;       // Z3 = I*K
    current.T = F * (B - H); // T3 = F*(B-H)
 
#ifdef DEBUG
    current.test_invariant();
#endif
}
 
void full_addition_step_for_miller_loop(
    const extended_edwards_G1_projective &base,
    extended_edwards_G1_projective &current,
    edwards_Fq_conic_coefficients &cc)
{
    const edwards_Fq &X1 = current.X, &Y1 = current.Y, &Z1 = current.Z,
                     &T1 = current.T;
    const edwards_Fq &X2 = base.X, &Y2 = base.Y, &Z2 = base.Z, &T2 = base.T;
 
    const edwards_Fq A = X1 * X2; // A    = X1*X2
    const edwards_Fq B = Y1 * Y2; // B    = Y1*Y2
    const edwards_Fq C = Z1 * T2; // C    = Z1*T2
    const edwards_Fq D = T1 * Z2; // D    = T1*Z2
    const edwards_Fq E = D + C;   // E    = D+C
    const edwards_Fq F =
        (X1 - Y1) * (X2 + Y2) + B - A; // F    = (X1-Y1)*(X2+Y2)+B-A
    const edwards_Fq G = B + A;        // G    = B + A (edwards_a=1)
    const edwards_Fq H = D - C;        // H    = D-C
    const edwards_Fq I = T1 * T2;      // I    = T1*T2
 
    cc.c_ZZ = (T1 - X1) * (T2 + X2) - I + A;     // c_ZZ = (T1-X1)*(T2+X2)-I+A
    cc.c_XY = X1 * Z2 - X2 * Z1 + F;             // c_XY = X1*Z2-X2*Z1+F
    cc.c_XZ = (Y1 - T1) * (Y2 + T2) - B + I - H; // c_XZ = (Y1-T1)*(Y2+T2)-B+I-H
    current.X = E * F;                           // X3   = E*F
    current.Y = G * H;                           // Y3   = G*H
    current.Z = F * G;                           // Z3   = F*G
    current.T = E * H;                           // T3   = E*H
 
#ifdef DEBUG
    current.test_invariant();
#endif
}
 
void mixed_addition_step_for_miller_loop(
    const extended_edwards_G1_projective &base,
    extended_edwards_G1_projective &current,
    edwards_Fq_conic_coefficients &cc)
{
    const edwards_Fq &X1 = current.X, &Y1 = current.Y, &Z1 = current.Z,
                     &T1 = current.T;
    const edwards_Fq &X2 = base.X, &Y2 = base.Y, &T2 = base.T;
 
    const edwards_Fq A = X1 * X2; // A    = X1*X2
    const edwards_Fq B = Y1 * Y2; // B    = Y1*Y2
    const edwards_Fq C = Z1 * T2; // C    = Z1*T2
    const edwards_Fq D = T1;      // D    = T1*Z2
    const edwards_Fq E = D + C;   // E    = D+C
    const edwards_Fq F =
        (X1 - Y1) * (X2 + Y2) + B - A; // F    = (X1-Y1)*(X2+Y2)+B-A
    const edwards_Fq G = B + A;        // G    = B + A (edwards_a=1)
    const edwards_Fq H = D - C;        // H    = D-C
    const edwards_Fq I = T1 * T2;      // I    = T1*T2
 
    cc.c_ZZ = (T1 - X1) * (T2 + X2) - I + A;     // c_ZZ = (T1-X1)*(T2+X2)-I+A
    cc.c_XY = X1 - X2 * Z1 + F;                  // c_XY = X1*Z2-X2*Z1+F
    cc.c_XZ = (Y1 - T1) * (Y2 + T2) - B + I - H; // c_XZ = (Y1-T1)*(Y2+T2)-B+I-H
    current.X = E * F;                           // X3   = E*F
    current.Y = G * H;                           // Y3   = G*H
    current.Z = F * G;                           // Z3   = F*G
    current.T = E * H;                           // T3   = E*H
 
#ifdef DEBUG
    current.test_invariant();
#endif
}
 
edwards_tate_G1_precomp edwards_tate_precompute_G1(const edwards_G1 &P)
{
    enter_block("Call to edwards_tate_precompute_G1");
    edwards_tate_G1_precomp result;
 
    edwards_G1 Pcopy = P;
    Pcopy.to_affine_coordinates();
 
    extended_edwards_G1_projective P_ext;
    P_ext.X = Pcopy.X;
    P_ext.Y = Pcopy.Y;
    P_ext.Z = Pcopy.Z;
    P_ext.T = Pcopy.X * Pcopy.Y;
 
    extended_edwards_G1_projective R = P_ext;
 
    bool found_one = false;
    for (long i = edwards_modulus_r.max_bits(); i >= 0; --i) {
        const bool bit = edwards_modulus_r.test_bit(i);
        if (!found_one) {
            /* this skips the MSB itself */
            found_one |= bit;
            continue;
        }
 
        /* code below gets executed for all bits (EXCEPT the MSB itself) of
           edwards_modulus_r (skipping leading zeros) in MSB to LSB
           order */
        edwards_Fq_conic_coefficients cc;
        doubling_step_for_miller_loop(R, cc);
        result.push_back(cc);
 
        if (bit) {
            mixed_addition_step_for_miller_loop(P_ext, R, cc);
            result.push_back(cc);
        }
    }
 
    leave_block("Call to edwards_tate_precompute_G1");
    return result;
}
 
edwards_Fq6 edwards_tate_miller_loop(
    const edwards_tate_G1_precomp &prec_P,
    const edwards_tate_G2_precomp &prec_Q)
{
    enter_block("Call to edwards_tate_miller_loop");
 
    edwards_Fq6 f = edwards_Fq6::one();
 
    bool found_one = false;
    size_t idx = 0;
    for (long i = edwards_modulus_r.max_bits() - 1; i >= 0; --i) {
        const bool bit = edwards_modulus_r.test_bit(i);
        if (!found_one) {
            /* this skips the MSB itself */
            found_one |= bit;
            continue;
        }
 
        /* code below gets executed for all bits (EXCEPT the MSB itself) of
           edwards_modulus_r (skipping leading zeros) in MSB to LSB
           order */
        edwards_Fq_conic_coefficients cc = prec_P[idx++];
        edwards_Fq6 g_RR_at_Q = edwards_Fq6(
            edwards_Fq3(cc.c_XZ, edwards_Fq(0l), edwards_Fq(0l)) +
                cc.c_XY * prec_Q.y0,
            cc.c_ZZ * prec_Q.eta);
        f = f.squared() * g_RR_at_Q;
        if (bit) {
            cc = prec_P[idx++];
 
            edwards_Fq6 g_RP_at_Q = edwards_Fq6(
                edwards_Fq3(cc.c_XZ, edwards_Fq(0l), edwards_Fq(0l)) +
                    cc.c_XY * prec_Q.y0,
                cc.c_ZZ * prec_Q.eta);
            f = f * g_RP_at_Q;
        }
    }
    leave_block("Call to edwards_tate_miller_loop");
 
    return f;
}
 
edwards_Fq6 edwards_tate_pairing(const edwards_G1 &P, const edwards_G2 &Q)
{
    enter_block("Call to edwards_tate_pairing");
    edwards_tate_G1_precomp prec_P = edwards_tate_precompute_G1(P);
    edwards_tate_G2_precomp prec_Q = edwards_tate_precompute_G2(Q);
    edwards_Fq6 result = edwards_tate_miller_loop(prec_P, prec_Q);
    leave_block("Call to edwards_tate_pairing");
    return result;
}
 
edwards_GT edwards_tate_reduced_pairing(
    const edwards_G1 &P, const edwards_G2 &Q)
{
    enter_block("Call to edwards_tate_reduced_pairing");
    const edwards_Fq6 f = edwards_tate_pairing(P, Q);
    const edwards_GT result = edwards_final_exponentiation(f);
    leave_block("Call to edwards_tate_reduce_pairing");
    return result;
}
 
struct extended_edwards_G2_projective {
    edwards_Fq3 X;
    edwards_Fq3 Y;
    edwards_Fq3 Z;
    edwards_Fq3 T;
 
    void print() const
    {
        printf("extended edwards_G2 projective X/Y/Z/T:\n");
        X.print();
        Y.print();
        Z.print();
        T.print();
    }
 
    void test_invariant() const { assert(T * Z == X * Y); }
};
 
void doubling_step_for_flipped_miller_loop(
    extended_edwards_G2_projective &current, edwards_Fq3_conic_coefficients &cc)
{
    const edwards_Fq3 &X = current.X, &Y = current.Y, &Z = current.Z,
                      &T = current.T;
    const edwards_Fq3 A = X.squared();       // A    = X1^2
    const edwards_Fq3 B = Y.squared();       // B    = Y1^2
    const edwards_Fq3 C = Z.squared();       // C    = Z1^2
    const edwards_Fq3 D = (X + Y).squared(); // D    = (X1+Y1)^2
    const edwards_Fq3 E = (Y + Z).squared(); // E    = (Y1+Z1)^2
    const edwards_Fq3 F = D - (A + B);       // F    = D-(A+B)
    const edwards_Fq3 G = E - (B + C);       // G    = E-(B+C)
    const edwards_Fq3 H =
        edwards_G2::mul_by_a(A); // edwards_param_twist_coeff_a is 1 * X for us
                                 // H    = twisted_a * A
    const edwards_Fq3 I = H + B; // I    = H+B
    const edwards_Fq3 J = C - I; // J    = C-I
    const edwards_Fq3 K = J + C; // K    = J+C
 
    cc.c_ZZ = Y * (T - X); // c_ZZ = 2*Y1*(T1-X1)
    cc.c_ZZ = cc.c_ZZ + cc.c_ZZ;
 
    // c_XY = 2*(C-edwards_a * A * delta_3-B)+G (edwards_a = 1 for us)
    cc.c_XY = C - edwards_G2::mul_by_a(A) -
              B; // edwards_param_twist_coeff_a is 1 * X for us
    cc.c_XY = cc.c_XY + cc.c_XY + G;
 
    // c_XZ = 2*(edwards_a*X1*T1*delta_3-B) (edwards_a = 1 for us)
    cc.c_XZ = edwards_G2::mul_by_a(X * T) -
              B; // edwards_param_twist_coeff_a is 1 * X for us
    cc.c_XZ = cc.c_XZ + cc.c_XZ;
 
    current.X = F * K;       // X3 = F*K
    current.Y = I * (B - H); // Y3 = I*(B-H)
    current.Z = I * K;       // Z3 = I*K
    current.T = F * (B - H); // T3 = F*(B-H)
#ifdef DEBUG
    current.test_invariant();
#endif
}
 
void full_addition_step_for_flipped_miller_loop(
    const extended_edwards_G2_projective &base,
    extended_edwards_G2_projective &current,
    edwards_Fq3_conic_coefficients &cc)
{
    const edwards_Fq3 &X1 = current.X, &Y1 = current.Y, &Z1 = current.Z,
                      &T1 = current.T;
    const edwards_Fq3 &X2 = base.X, &Y2 = base.Y, &Z2 = base.Z, &T2 = base.T;
 
    const edwards_Fq3 A = X1 * X2; // A    = X1*X2
    const edwards_Fq3 B = Y1 * Y2; // B    = Y1*Y2
    const edwards_Fq3 C = Z1 * T2; // C    = Z1*T2
    const edwards_Fq3 D = T1 * Z2; // D    = T1*Z2
    const edwards_Fq3 E = D + C;   // E    = D+C
    const edwards_Fq3 F =
        (X1 - Y1) * (X2 + Y2) + B - A; // F    = (X1-Y1)*(X2+Y2)+B-A
    // G = B + twisted_edwards_a * A
    const edwards_Fq3 G =
        B +
        edwards_G2::mul_by_a(A);   // edwards_param_twist_coeff_a is 1*X for us
    const edwards_Fq3 H = D - C;   // H    = D-C
    const edwards_Fq3 I = T1 * T2; // I    = T1*T2
 
    // c_ZZ = delta_3* ((T1-X1)*(T2+X2)-I+A)
    cc.c_ZZ = edwards_G2::mul_by_a(
        (T1 - X1) * (T2 + X2) - I +
        A); // edwards_param_twist_coeff_a is 1*X for us
 
    cc.c_XY = X1 * Z2 - X2 * Z1 + F;             // c_XY = X1*Z2-X2*Z1+F
    cc.c_XZ = (Y1 - T1) * (Y2 + T2) - B + I - H; // c_XZ = (Y1-T1)*(Y2+T2)-B+I-H
    current.X = E * F;                           // X3   = E*F
    current.Y = G * H;                           // Y3   = G*H
    current.Z = F * G;                           // Z3   = F*G
    current.T = E * H;                           // T3   = E*H
 
#ifdef DEBUG
    current.test_invariant();
#endif
}
 
void mixed_addition_step_for_flipped_miller_loop(
    const extended_edwards_G2_projective &base,
    extended_edwards_G2_projective &current,
    edwards_Fq3_conic_coefficients &cc)
{
    const edwards_Fq3 &X1 = current.X, &Y1 = current.Y, &Z1 = current.Z,
                      &T1 = current.T;
    const edwards_Fq3 &X2 = base.X, &Y2 = base.Y, &T2 = base.T;
 
    const edwards_Fq3 A = X1 * X2; // A    = X1*X2
    const edwards_Fq3 B = Y1 * Y2; // B    = Y1*Y2
    const edwards_Fq3 C = Z1 * T2; // C    = Z1*T2
    const edwards_Fq3 E = T1 + C;  // E    = T1+C
    const edwards_Fq3 F =
        (X1 - Y1) * (X2 + Y2) + B - A; // F    = (X1-Y1)*(X2+Y2)+B-A
    // G = B + twisted_edwards_a * A
    const edwards_Fq3 G =
        B +
        edwards_G2::mul_by_a(A);   // edwards_param_twist_coeff_a is 1*X for us
    const edwards_Fq3 H = T1 - C;  // H    = T1-C
    const edwards_Fq3 I = T1 * T2; // I    = T1*T2
 
    // c_ZZ = delta_3* ((T1-X1)*(T2+X2)-I+A)
    cc.c_ZZ = edwards_G2::mul_by_a(
        (T1 - X1) * (T2 + X2) - I +
        A); // edwards_param_twist_coeff_a is 1*X for us
 
    cc.c_XY = X1 - X2 * Z1 + F;                  // c_XY = X1*Z2-X2*Z1+F
    cc.c_XZ = (Y1 - T1) * (Y2 + T2) - B + I - H; // c_XZ = (Y1-T1)*(Y2+T2)-B+I-H
    current.X = E * F;                           // X3   = E*F
    current.Y = G * H;                           // Y3   = G*H
    current.Z = F * G;                           // Z3   = F*G
    current.T = E * H;                           // T3   = E*H
 
#ifdef DEBUG
    current.test_invariant();
#endif
}
 
edwards_ate_G1_precomp edwards_ate_precompute_G1(const edwards_G1 &P)
{
    enter_block("Call to edwards_ate_precompute_G1");
    edwards_G1 Pcopy = P;
    Pcopy.to_affine_coordinates();
    edwards_ate_G1_precomp result;
    result.P_XY = Pcopy.X * Pcopy.Y;
    result.P_XZ = Pcopy.X; // P.X * P.Z but P.Z = 1
    result.P_ZZplusYZ =
        (edwards_Fq::one() + Pcopy.Y); // (P.Z + P.Y) * P.Z but P.Z = 1
    leave_block("Call to edwards_ate_precompute_G1");
    return result;
}
 
edwards_ate_G2_precomp edwards_ate_precompute_G2(const edwards_G2 &Q)
{
    enter_block("Call to edwards_ate_precompute_G2");
    const bigint<edwards_Fr::num_limbs> &loop_count = edwards_ate_loop_count;
    edwards_ate_G2_precomp result;
 
    edwards_G2 Qcopy(Q);
    Qcopy.to_affine_coordinates();
 
    extended_edwards_G2_projective Q_ext;
    Q_ext.X = Qcopy.X;
    Q_ext.Y = Qcopy.Y;
    Q_ext.Z = Qcopy.Z;
    Q_ext.T = Qcopy.X * Qcopy.Y;
 
    extended_edwards_G2_projective R = Q_ext;
 
    bool found_one = false;
    for (long i = loop_count.max_bits() - 1; i >= 0; --i) {
        const bool bit = loop_count.test_bit(i);
        if (!found_one) {
            /* this skips the MSB itself */
            found_one |= bit;
            continue;
        }
 
        edwards_Fq3_conic_coefficients cc;
        doubling_step_for_flipped_miller_loop(R, cc);
        result.push_back(cc);
        if (bit) {
            mixed_addition_step_for_flipped_miller_loop(Q_ext, R, cc);
            result.push_back(cc);
        }
    }
 
    leave_block("Call to edwards_ate_precompute_G2");
    return result;
}
 
edwards_Fq6 edwards_ate_miller_loop(
    const edwards_ate_G1_precomp &prec_P, const edwards_ate_G2_precomp &prec_Q)
{
    enter_block("Call to edwards_ate_miller_loop");
    const bigint<edwards_Fr::num_limbs> &loop_count = edwards_ate_loop_count;
 
    edwards_Fq6 f = edwards_Fq6::one();
 
    bool found_one = false;
    size_t idx = 0;
    for (long i = loop_count.max_bits() - 1; i >= 0; --i) {
        const bool bit = loop_count.test_bit(i);
        if (!found_one) {
            /* this skips the MSB itself */
            found_one |= bit;
            continue;
        }
 
        /* code below gets executed for all bits (EXCEPT the MSB itself) of
           edwards_param_p (skipping leading zeros) in MSB to LSB
           order */
        edwards_Fq3_conic_coefficients cc = prec_Q[idx++];
 
        edwards_Fq6 g_RR_at_P = edwards_Fq6(
            prec_P.P_XY * cc.c_XY + prec_P.P_XZ * cc.c_XZ,
            prec_P.P_ZZplusYZ * cc.c_ZZ);
        f = f.squared() * g_RR_at_P;
        if (bit) {
            cc = prec_Q[idx++];
            edwards_Fq6 g_RQ_at_P = edwards_Fq6(
                prec_P.P_ZZplusYZ * cc.c_ZZ,
                prec_P.P_XY * cc.c_XY + prec_P.P_XZ * cc.c_XZ);
            f = f * g_RQ_at_P;
        }
    }
    leave_block("Call to edwards_ate_miller_loop");
 
    return f;
}
 
edwards_Fq6 edwards_ate_double_miller_loop(
    const edwards_ate_G1_precomp &prec_P1,
    const edwards_ate_G2_precomp &prec_Q1,
    const edwards_ate_G1_precomp &prec_P2,
    const edwards_ate_G2_precomp &prec_Q2)
{
    enter_block("Call to edwards_ate_double_miller_loop");
    const bigint<edwards_Fr::num_limbs> &loop_count = edwards_ate_loop_count;
 
    edwards_Fq6 f = edwards_Fq6::one();
 
    bool found_one = false;
    size_t idx = 0;
    for (long i = loop_count.max_bits() - 1; i >= 0; --i) {
        const bool bit = loop_count.test_bit(i);
        if (!found_one) {
            /* this skips the MSB itself */
            found_one |= bit;
            continue;
        }
 
        /* code below gets executed for all bits (EXCEPT the MSB itself) of
           edwards_param_p (skipping leading zeros) in MSB to LSB
           order */
        edwards_Fq3_conic_coefficients cc1 = prec_Q1[idx];
        edwards_Fq3_conic_coefficients cc2 = prec_Q2[idx];
        ++idx;
 
        edwards_Fq6 g_RR_at_P1 = edwards_Fq6(
            prec_P1.P_XY * cc1.c_XY + prec_P1.P_XZ * cc1.c_XZ,
            prec_P1.P_ZZplusYZ * cc1.c_ZZ);
 
        edwards_Fq6 g_RR_at_P2 = edwards_Fq6(
            prec_P2.P_XY * cc2.c_XY + prec_P2.P_XZ * cc2.c_XZ,
            prec_P2.P_ZZplusYZ * cc2.c_ZZ);
        f = f.squared() * g_RR_at_P1 * g_RR_at_P2;
 
        if (bit) {
            cc1 = prec_Q1[idx];
            cc2 = prec_Q2[idx];
            ++idx;
            edwards_Fq6 g_RQ_at_P1 = edwards_Fq6(
                prec_P1.P_ZZplusYZ * cc1.c_ZZ,
                prec_P1.P_XY * cc1.c_XY + prec_P1.P_XZ * cc1.c_XZ);
            edwards_Fq6 g_RQ_at_P2 = edwards_Fq6(
                prec_P2.P_ZZplusYZ * cc2.c_ZZ,
                prec_P2.P_XY * cc2.c_XY + prec_P2.P_XZ * cc2.c_XZ);
            f = f * g_RQ_at_P1 * g_RQ_at_P2;
        }
    }
    leave_block("Call to edwards_ate_double_miller_loop");
 
    return f;
}
 
edwards_Fq6 edwards_ate_pairing(const edwards_G1 &P, const edwards_G2 &Q)
{
    enter_block("Call to edwards_ate_pairing");
    edwards_ate_G1_precomp prec_P = edwards_ate_precompute_G1(P);
    edwards_ate_G2_precomp prec_Q = edwards_ate_precompute_G2(Q);
    edwards_Fq6 result = edwards_ate_miller_loop(prec_P, prec_Q);
    leave_block("Call to edwards_ate_pairing");
    return result;
}
 
edwards_GT edwards_ate_reduced_pairing(const edwards_G1 &P, const edwards_G2 &Q)
{
    enter_block("Call to edwards_ate_reduced_pairing");
    const edwards_Fq6 f = edwards_ate_pairing(P, Q);
    const edwards_GT result = edwards_final_exponentiation(f);
    leave_block("Call to edwards_ate_reduced_pairing");
    return result;
}
 
edwards_G1_precomp edwards_precompute_G1(const edwards_G1 &P)
{
    return edwards_ate_precompute_G1(P);
}
 
edwards_G2_precomp edwards_precompute_G2(const edwards_G2 &Q)
{
    return edwards_ate_precompute_G2(Q);
}
 
edwards_Fq6 edwards_miller_loop(
    const edwards_G1_precomp &prec_P, const edwards_G2_precomp &prec_Q)
{
    return edwards_ate_miller_loop(prec_P, prec_Q);
}
 
edwards_Fq6 edwards_double_miller_loop(
    const edwards_G1_precomp &prec_P1,
    const edwards_G2_precomp &prec_Q1,
    const edwards_G1_precomp &prec_P2,
    const edwards_G2_precomp &prec_Q2)
{
    return edwards_ate_double_miller_loop(prec_P1, prec_Q1, prec_P2, prec_Q2);
}
 
edwards_Fq6 edwards_pairing(const edwards_G1 &P, const edwards_G2 &Q)
{
    return edwards_ate_pairing(P, Q);
}
 
edwards_GT edwards_reduced_pairing(const edwards_G1 &P, const edwards_G2 &Q)
{
    return edwards_ate_reduced_pairing(P, Q);
}
} // namespace libff