#ifndef FP3_TCC_
#define FP3_TCC_

#include <libff/algebra/fields/field_utils.hpp>

namespace libff
{

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::zero()
{
    return Fp3_model<n, modulus>(my_Fp::zero(), my_Fp::zero(), my_Fp::zero());
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::one()
{
    return Fp3_model<n, modulus>(my_Fp::one(), my_Fp::zero(), my_Fp::zero());
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::random_element()
{
    Fp3_model<n, modulus> r;
    r.coeffs[0] = my_Fp::random_element();
    r.coeffs[1] = my_Fp::random_element();
    r.coeffs[2] = my_Fp::random_element();

    return r;
}

template<mp_size_t n, const bigint<n> &modulus>
bool Fp3_model<n, modulus>::operator==(const Fp3_model<n, modulus> &other) const
{
    return (
        this->coeffs[0] == other.coeffs[0] &&
        this->coeffs[1] == other.coeffs[1] &&
        this->coeffs[2] == other.coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
bool Fp3_model<n, modulus>::operator!=(const Fp3_model<n, modulus> &other) const
{
    return !(operator==(other));
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator+(
    const Fp3_model<n, modulus> &other) const
{
    return Fp3_model<n, modulus>(
        this->coeffs[0] + other.coeffs[0],
        this->coeffs[1] + other.coeffs[1],
        this->coeffs[2] + other.coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator-(
    const Fp3_model<n, modulus> &other) const
{
    return Fp3_model<n, modulus>(
        this->coeffs[0] - other.coeffs[0],
        this->coeffs[1] - other.coeffs[1],
        this->coeffs[2] - other.coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> operator*(
    const Fp_model<n, modulus> &lhs, const Fp3_model<n, modulus> &rhs)
{
    return Fp3_model<n, modulus>(
        lhs * rhs.coeffs[0], lhs * rhs.coeffs[1], lhs * rhs.coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator*(
    const Fp3_model<n, modulus> &other) const
{

    const my_Fp &A = other.coeffs[0], &B = other.coeffs[1],
                &C = other.coeffs[2], &a = this->coeffs[0],
                &b = this->coeffs[1], &c = this->coeffs[2];
    const my_Fp aA = a * A;
    const my_Fp bB = b * B;
    const my_Fp cC = c * C;

    return Fp3_model<n, modulus>(
        aA + non_residue * ((b + c) * (B + C) - bB - cC),
        (a + b) * (A + B) - aA - bB + non_residue * cC,
        (a + c) * (A + C) - aA + bB - cC);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator-() const
{
    return Fp3_model<n, modulus>(
        -this->coeffs[0], -this->coeffs[1], -this->coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::squared() const
{

    const my_Fp &a = this->coeffs[0], &b = this->coeffs[1],
                &c = this->coeffs[2];
    const my_Fp s0 = a.squared();
    const my_Fp ab = a * b;
    const my_Fp s1 = ab + ab;
    const my_Fp s2 = (a - b + c).squared();
    const my_Fp bc = b * c;
    const my_Fp s3 = bc + bc;
    const my_Fp s4 = c.squared();

    return Fp3_model<n, modulus>(
        s0 + non_residue * s3, s1 + non_residue * s4, s1 + s2 + s3 - s0 - s4);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::inverse() const
{
    const my_Fp &a = this->coeffs[0], &b = this->coeffs[1],
                &c = this->coeffs[2];

    const my_Fp t0 = a.squared();
    const my_Fp t1 = b.squared();
    const my_Fp t2 = c.squared();
    const my_Fp t3 = a * b;
    const my_Fp t4 = a * c;
    const my_Fp t5 = b * c;
    const my_Fp c0 = t0 - non_residue * t5;
    const my_Fp c1 = non_residue * t2 - t3;
    const my_Fp c2 = t1 - t4; 

    const my_Fp t6 = (a * c0 + non_residue * (c * c1 + b * c2)).inverse();
    return Fp3_model<n, modulus>(t6 * c0, t6 * c1, t6 * c2);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::Frobenius_map(
    unsigned long power) const
{
    return Fp3_model<n, modulus>(
        coeffs[0],
        Frobenius_coeffs_c1[power % 3] * coeffs[1],
        Frobenius_coeffs_c2[power % 3] * coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::sqrt() const
{
    Fp3_model<n, modulus> one = Fp3_model<n, modulus>::one();

    size_t v = Fp3_model<n, modulus>::s;
    Fp3_model<n, modulus> z = Fp3_model<n, modulus>::nqr_to_t;
    Fp3_model<n, modulus> w = (*this) ^ Fp3_model<n, modulus>::t_minus_1_over_2;
    Fp3_model<n, modulus> x = (*this) * w;
    Fp3_model<n, modulus> b = x * w; 

#if DEBUG

    Fp3_model<n, modulus> check = b;
    for (size_t i = 0; i < v - 1; ++i) {
        check = check.squared();
    }
    if (check != one) {
        assert(0);
    }
#endif

    while (b != one) {
        size_t m = 0;
        Fp3_model<n, modulus> b2m = b;
        while (b2m != one) {

            b2m = b2m.squared();
            m += 1;
        }

        int j = v - m - 1;
        w = z;
        while (j > 0) {
            w = w.squared();
            --j;
        } 

        z = w.squared();
        b = b * z;
        x = x * w;
        v = m;
    }

    return x;
}

template<mp_size_t n, const bigint<n> &modulus>
template<mp_size_t m>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator^(
    const bigint<m> &pow) const
{
    return power<Fp3_model<n, modulus>>(*this, pow);
}

template<mp_size_t n, const bigint<n> &modulus>
std::ostream &operator<<(std::ostream &out, const Fp3_model<n, modulus> &el)
{
    out << el.coeffs[0] << OUTPUT_SEPARATOR << el.coeffs[1] << OUTPUT_SEPARATOR
        << el.coeffs[2];
    return out;
}

template<mp_size_t n, const bigint<n> &modulus>
std::istream &operator>>(std::istream &in, Fp3_model<n, modulus> &el)
{
    in >> el.coeffs[0] >> el.coeffs[1] >> el.coeffs[2];
    return in;
}

template<mp_size_t n, const bigint<n> &modulus>
std::ostream &operator<<(
    std::ostream &out, const std::vector<Fp3_model<n, modulus>> &v)
{
    out << v.size() << "\n";
    for (const Fp3_model<n, modulus> &t : v) {
        out << t << OUTPUT_NEWLINE;
    }

    return out;
}

template<mp_size_t n, const bigint<n> &modulus>
std::istream &operator>>(
    std::istream &in, std::vector<Fp3_model<n, modulus>> &v)
{
    v.clear();

    size_t s;
    in >> s;

    char b;
    in.read(&b, 1);

    v.reserve(s);

    for (size_t i = 0; i < s; ++i) {
        Fp3_model<n, modulus> el;
        in >> el;
        v.emplace_back(el);
    }

    return in;
}

} 
#endif
/** @file
 *****************************************************************************
 Implementation of arithmetic in the finite field F[p^3].
 *****************************************************************************
 * @author     This file is part of libff, developed by SCIPR Lab
 *             and contributors (see AUTHORS).
 * @copyright  MIT license (see LICENSE file)
 *****************************************************************************/

#ifndef FP3_TCC_
#define FP3_TCC_

#include <libff/algebra/fields/field_utils.hpp>

namespace libff
{

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::zero()
{
    return Fp3_model<n, modulus>(my_Fp::zero(), my_Fp::zero(), my_Fp::zero());
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::one()
{
    return Fp3_model<n, modulus>(my_Fp::one(), my_Fp::zero(), my_Fp::zero());
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::random_element()
{
    Fp3_model<n, modulus> r;
    r.coeffs[0] = my_Fp::random_element();
    r.coeffs[1] = my_Fp::random_element();
    r.coeffs[2] = my_Fp::random_element();

    return r;
}

template<mp_size_t n, const bigint<n> &modulus>
bool Fp3_model<n, modulus>::operator==(const Fp3_model<n, modulus> &other) const
{
    return (
        this->coeffs[0] == other.coeffs[0] &&
        this->coeffs[1] == other.coeffs[1] &&
        this->coeffs[2] == other.coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
bool Fp3_model<n, modulus>::operator!=(const Fp3_model<n, modulus> &other) const
{
    return !(operator==(other));
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator+(
    const Fp3_model<n, modulus> &other) const
{
    return Fp3_model<n, modulus>(
        this->coeffs[0] + other.coeffs[0],
        this->coeffs[1] + other.coeffs[1],
        this->coeffs[2] + other.coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator-(
    const Fp3_model<n, modulus> &other) const
{
    return Fp3_model<n, modulus>(
        this->coeffs[0] - other.coeffs[0],
        this->coeffs[1] - other.coeffs[1],
        this->coeffs[2] - other.coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> operator*(
    const Fp_model<n, modulus> &lhs, const Fp3_model<n, modulus> &rhs)
{
    return Fp3_model<n, modulus>(
        lhs * rhs.coeffs[0], lhs * rhs.coeffs[1], lhs * rhs.coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator*(
    const Fp3_model<n, modulus> &other) const
{
    /* Devegili OhEig Scott Dahab --- Multiplication and Squaring on
     * Pairing-Friendly Fields.pdf; Section 4 (Karatsuba) */
    const my_Fp &A = other.coeffs[0], &B = other.coeffs[1],
                &C = other.coeffs[2], &a = this->coeffs[0],
                &b = this->coeffs[1], &c = this->coeffs[2];
    const my_Fp aA = a * A;
    const my_Fp bB = b * B;
    const my_Fp cC = c * C;

    return Fp3_model<n, modulus>(
        aA + non_residue * ((b + c) * (B + C) - bB - cC),
        (a + b) * (A + B) - aA - bB + non_residue * cC,
        (a + c) * (A + C) - aA + bB - cC);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator-() const
{
    return Fp3_model<n, modulus>(
        -this->coeffs[0], -this->coeffs[1], -this->coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::squared() const
{
    /* Devegili OhEig Scott Dahab --- Multiplication and Squaring on
     * Pairing-Friendly Fields.pdf; Section 4 (CH-SQR2) */
    const my_Fp &a = this->coeffs[0], &b = this->coeffs[1],
                &c = this->coeffs[2];
    const my_Fp s0 = a.squared();
    const my_Fp ab = a * b;
    const my_Fp s1 = ab + ab;
    const my_Fp s2 = (a - b + c).squared();
    const my_Fp bc = b * c;
    const my_Fp s3 = bc + bc;
    const my_Fp s4 = c.squared();

    return Fp3_model<n, modulus>(
        s0 + non_residue * s3, s1 + non_residue * s4, s1 + s2 + s3 - s0 - s4);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::inverse() const
{
    const my_Fp &a = this->coeffs[0], &b = this->coeffs[1],
                &c = this->coeffs[2];

    /* From "High-Speed Software Implementation of the Optimal Ate Pairing over
     * Barreto-Naehrig Curves"; Algorithm 17 */
    const my_Fp t0 = a.squared();
    const my_Fp t1 = b.squared();
    const my_Fp t2 = c.squared();
    const my_Fp t3 = a * b;
    const my_Fp t4 = a * c;
    const my_Fp t5 = b * c;
    const my_Fp c0 = t0 - non_residue * t5;
    const my_Fp c1 = non_residue * t2 - t3;
    const my_Fp c2 = t1 - t4; // typo in paper referenced above. should be "-"
                              // as per Scott, but is "*"
    const my_Fp t6 = (a * c0 + non_residue * (c * c1 + b * c2)).inverse();
    return Fp3_model<n, modulus>(t6 * c0, t6 * c1, t6 * c2);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::Frobenius_map(
    unsigned long power) const
{
    return Fp3_model<n, modulus>(
        coeffs[0],
        Frobenius_coeffs_c1[power % 3] * coeffs[1],
        Frobenius_coeffs_c2[power % 3] * coeffs[2]);
}

template<mp_size_t n, const bigint<n> &modulus>
Fp3_model<n, modulus> Fp3_model<n, modulus>::sqrt() const
{
    Fp3_model<n, modulus> one = Fp3_model<n, modulus>::one();

    size_t v = Fp3_model<n, modulus>::s;
    Fp3_model<n, modulus> z = Fp3_model<n, modulus>::nqr_to_t;
    Fp3_model<n, modulus> w = (*this) ^ Fp3_model<n, modulus>::t_minus_1_over_2;
    Fp3_model<n, modulus> x = (*this) * w;
    Fp3_model<n, modulus> b = x * w; // b = (*this)^t

#if DEBUG
    // check if square with euler's criterion
    Fp3_model<n, modulus> check = b;
    for (size_t i = 0; i < v - 1; ++i) {
        check = check.squared();
    }
    if (check != one) {
        assert(0);
    }
#endif

    // compute square root with Tonelli--Shanks
    // (does not terminate if not a square!)

    while (b != one) {
        size_t m = 0;
        Fp3_model<n, modulus> b2m = b;
        while (b2m != one) {
            /* invariant: b2m = b^(2^m) after entering this loop */
            b2m = b2m.squared();
            m += 1;
        }

        int j = v - m - 1;
        w = z;
        while (j > 0) {
            w = w.squared();
            --j;
        } // w = z^2^(v-m-1)

        z = w.squared();
        b = b * z;
        x = x * w;
        v = m;
    }

    return x;
}

template<mp_size_t n, const bigint<n> &modulus>
template<mp_size_t m>
Fp3_model<n, modulus> Fp3_model<n, modulus>::operator^(
    const bigint<m> &pow) const
{
    return power<Fp3_model<n, modulus>>(*this, pow);
}

template<mp_size_t n, const bigint<n> &modulus>
std::ostream &operator<<(std::ostream &out, const Fp3_model<n, modulus> &el)
{
    out << el.coeffs[0] << OUTPUT_SEPARATOR << el.coeffs[1] << OUTPUT_SEPARATOR
        << el.coeffs[2];
    return out;
}

template<mp_size_t n, const bigint<n> &modulus>
std::istream &operator>>(std::istream &in, Fp3_model<n, modulus> &el)
{
    in >> el.coeffs[0] >> el.coeffs[1] >> el.coeffs[2];
    return in;
}

template<mp_size_t n, const bigint<n> &modulus>
std::ostream &operator<<(
    std::ostream &out, const std::vector<Fp3_model<n, modulus>> &v)
{
    out << v.size() << "\n";
    for (const Fp3_model<n, modulus> &t : v) {
        out << t << OUTPUT_NEWLINE;
    }

    return out;
}

template<mp_size_t n, const bigint<n> &modulus>
std::istream &operator>>(
    std::istream &in, std::vector<Fp3_model<n, modulus>> &v)
{
    v.clear();

    size_t s;
    in >> s;

    char b;
    in.read(&b, 1);

    v.reserve(s);

    for (size_t i = 0; i < s; ++i) {
        Fp3_model<n, modulus> el;
        in >> el;
        v.emplace_back(el);
    }

    return in;
}

} // namespace libff
#endif // FP3_TCC_
