bn128_g1.cpp

Go to the documentation of this file.

 
#include <libff/algebra/curves/bn128/bn128_g1.hpp>
#include <libff/algebra/curves/bn128/bn_utils.hpp>
 
namespace libff
{
 
#ifdef PROFILE_OP_COUNTS
long long bn128_G1::add_cnt = 0;
long long bn128_G1::dbl_cnt = 0;
#endif
 
std::vector<size_t> bn128_G1::wnaf_window_table;
std::vector<size_t> bn128_G1::fixed_base_exp_window_table;
bn128_G1 bn128_G1::G1_zero;
bn128_G1 bn128_G1::G1_one;
bigint<bn128_G1::h_limbs> bn128_G1::h;
 
bn::Fp bn128_G1::sqrt(const bn::Fp &el)
{
    size_t v = bn128_Fq_s;
    bn::Fp z = bn128_Fq_nqr_to_t;
    bn::Fp w = mie::power(el, bn128_Fq_t_minus_1_over_2);
    bn::Fp x = el * w;
    bn::Fp b = x * w;
 
#if DEBUG
    // check if square with Euler's criterion
    bn::Fp check = b;
    for (size_t i = 0; i < v - 1; ++i) {
        bn::Fp::square(check, check);
    }
 
    assert(check == bn::Fp(1));
#endif
 
    // compute square root with Tonelli--Shanks
    // (does not terminate if not a square!)
 
    while (b != bn::Fp(1)) {
        size_t m = 0;
        bn::Fp b2m = b;
        while (b2m != bn::Fp(1)) {
            // invariant: b2m = b^(2^m) after entering this loop
            bn::Fp::square(b2m, b2m);
            m += 1;
        }
 
        // w = z^2^(v-m-1)
        int j = v - m - 1;
        w = z;
        while (j > 0) {
            bn::Fp::square(w, w);
            --j;
        }
 
        z = w * w;
        b = b * z;
        x = x * w;
        v = m;
    }
 
    return x;
}
 
bn128_G1::bn128_G1()
{
    this->X = G1_zero.X;
    this->Y = G1_zero.Y;
    this->Z = G1_zero.Z;
}
 
void bn128_G1::print() const
{
    if (this->is_zero()) {
        printf("O\n");
    } else {
        bn128_G1 copy(*this);
        copy.to_affine_coordinates();
        std::cout << "(" << copy.X.toString(10) << " : " << copy.Y.toString(10)
                  << " : " << copy.Z.toString(10) << ")\n";
    }
}
 
void bn128_G1::print_coordinates() const
{
    if (this->is_zero()) {
        printf("O\n");
    } else {
        std::cout << "(" << X.toString(10) << " : " << Y.toString(10) << " : "
                  << Z.toString(10) << ")\n";
    }
}
 
void bn128_G1::to_affine_coordinates()
{
    if (this->is_zero()) {
        X = 0;
        Y = 1;
        Z = 0;
    } else {
        bn::Fp r;
        r = Z;
        r.inverse();
        bn::Fp::square(Z, r);
        X *= Z;
        r *= Z;
        Y *= r;
        Z = 1;
    }
}
 
void bn128_G1::to_special() { this->to_affine_coordinates(); }
 
bool bn128_G1::is_special() const { return (this->is_zero() || this->Z == 1); }
 
bool bn128_G1::is_zero() const { return Z.isZero(); }
 
bool bn128_G1::operator==(const bn128_G1 &other) const
{
    if (this->is_zero()) {
        return other.is_zero();
    }
 
    if (other.is_zero()) {
        return false;
    }
 
    /* now neither is O */
 
    bn::Fp Z1sq, Z2sq, lhs, rhs;
    bn::Fp::square(Z1sq, this->Z);
    bn::Fp::square(Z2sq, other.Z);
    bn::Fp::mul(lhs, Z2sq, this->X);
    bn::Fp::mul(rhs, Z1sq, other.X);
 
    if (lhs != rhs) {
        return false;
    }
 
    bn::Fp Z1cubed, Z2cubed;
    bn::Fp::mul(Z1cubed, Z1sq, this->Z);
    bn::Fp::mul(Z2cubed, Z2sq, other.Z);
    bn::Fp::mul(lhs, Z2cubed, this->Y);
    bn::Fp::mul(rhs, Z1cubed, other.Y);
 
    return (lhs == rhs);
}
 
bool bn128_G1::operator!=(const bn128_G1 &other) const
{
    return !(operator==(other));
}
 
bn128_G1 bn128_G1::operator+(const bn128_G1 &other) const
{
    // handle special cases having to do with O
    if (this->is_zero()) {
        return other;
    }
 
    if (other.is_zero()) {
        return *this;
    }
 
    // no need to handle points of order 2,4
    // (they cannot exist in a prime-order subgroup)
 
    // handle double case, and then all the rest
    if (this->operator==(other)) {
        return this->dbl();
    } else {
        return this->add(other);
    }
}
 
bn128_G1 bn128_G1::operator-() const
{
    bn128_G1 result(*this);
    bn::Fp::neg(result.Y, result.Y);
    return result;
}
 
bn128_G1 bn128_G1::operator-(const bn128_G1 &other) const
{
    return (*this) + (-other);
}
 
bn128_G1 bn128_G1::add(const bn128_G1 &other) const
{
#ifdef PROFILE_OP_COUNTS
    this->add_cnt++;
#endif
 
    bn::Fp this_coord[3], other_coord[3], result_coord[3];
    this->fill_coord(this_coord);
    other.fill_coord(other_coord);
    bn::ecop::ECAdd(result_coord, this_coord, other_coord);
 
    bn128_G1 result(result_coord);
    return result;
}
 
bn128_G1 bn128_G1::mixed_add(const bn128_G1 &other) const
{
    if (this->is_zero()) {
        return other;
    }
 
    if (other.is_zero()) {
        return *this;
    }
 
    // no need to handle points of order 2,4
    // (they cannot exist in a prime-order subgroup)
 
#ifdef DEBUG
    assert(other.is_special());
#endif
 
    // check for doubling case
 
    // using Jacobian coordinates so:
    // (X1:Y1:Z1) = (X2:Y2:Z2)
    // iff
    // X1/Z1^2 == X2/Z2^2 and Y1/Z1^3 == Y2/Z2^3
    // iff
    // X1 * Z2^2 == X2 * Z1^2 and Y1 * Z2^3 == Y2 * Z1^3
 
    // we know that Z2 = 1
 
    bn::Fp Z1Z1;
    bn::Fp::square(Z1Z1, this->Z);
    const bn::Fp &U1 = this->X;
    bn::Fp U2;
    bn::Fp::mul(U2, other.X, Z1Z1);
    bn::Fp Z1_cubed;
    bn::Fp::mul(Z1_cubed, this->Z, Z1Z1);
 
    const bn::Fp &S1 = this->Y;
    bn::Fp S2;
    bn::Fp::mul(S2, other.Y, Z1_cubed); // S2 = Y2*Z1*Z1Z1
 
    if (U1 == U2 && S1 == S2) {
        // dbl case; nothing of above can be reused
        return this->dbl();
    }
 
#ifdef PROFILE_OP_COUNTS
    this->add_cnt++;
#endif
 
    bn128_G1 result;
    bn::Fp H, HH, I, J, r, V, tmp;
    // H = U2-X1
    bn::Fp::sub(H, U2, this->X);
    // HH = H^2
    bn::Fp::square(HH, H);
    // I = 4*HH
    bn::Fp::add(tmp, HH, HH);
    bn::Fp::add(I, tmp, tmp);
    // J = H*I
    bn::Fp::mul(J, H, I);
    // r = 2*(S2-Y1)
    bn::Fp::sub(tmp, S2, this->Y);
    bn::Fp::add(r, tmp, tmp);
    // V = X1*I
    bn::Fp::mul(V, this->X, I);
    // X3 = r^2-J-2*V
    bn::Fp::square(result.X, r);
    bn::Fp::sub(result.X, result.X, J);
    bn::Fp::sub(result.X, result.X, V);
    bn::Fp::sub(result.X, result.X, V);
    // Y3 = r*(V-X3)-2*Y1*J
    bn::Fp::sub(tmp, V, result.X);
    bn::Fp::mul(result.Y, r, tmp);
    bn::Fp::mul(tmp, this->Y, J);
    bn::Fp::sub(result.Y, result.Y, tmp);
    bn::Fp::sub(result.Y, result.Y, tmp);
    // Z3 = (Z1+H)^2-Z1Z1-HH
    bn::Fp::add(tmp, this->Z, H);
    bn::Fp::square(result.Z, tmp);
    bn::Fp::sub(result.Z, result.Z, Z1Z1);
    bn::Fp::sub(result.Z, result.Z, HH);
    return result;
}
 
bn128_G1 bn128_G1::dbl() const
{
#ifdef PROFILE_OP_COUNTS
    this->dbl_cnt++;
#endif
 
    bn::Fp this_coord[3], result_coord[3];
    this->fill_coord(this_coord);
    bn::ecop::ECDouble(result_coord, this_coord);
 
    bn128_G1 result(result_coord);
    return result;
}
 
bn128_G1 bn128_G1::mul_by_cofactor() const
{
    // Cofactor = 1
    return (*this);
}
 
const bn128_G1 &bn128_G1::zero() { return G1_zero; }
 
const bn128_G1 &bn128_G1::one() { return G1_one; }
 
bn128_G1 bn128_G1::random_element()
{
    return bn128_Fr::random_element().as_bigint() * G1_one;
}
 
bool bn128_G1::is_well_formed() const
{
    if (this->is_zero()) {
        return true;
    } else {
        /*
          y^2 = x^3 + b
 
          We are using Jacobian coordinates, so equation we need to check is
          actually
 
          (y/z^3)^2 = (x/z^2)^3 + b
          y^2 / z^6 = x^3 / z^6 + b
          y^2 = x^3 + b z^6
        */
        bn::Fp X2, Y2, Z2;
        bn::Fp::square(X2, this->X);
        bn::Fp::square(Y2, this->Y);
        bn::Fp::square(Z2, this->Z);
 
        bn::Fp X3, Z3, Z6;
        bn::Fp::mul(X3, X2, this->X);
        bn::Fp::mul(Z3, Z2, this->Z);
        bn::Fp::square(Z6, Z3);
 
        return (Y2 == X3 + bn128_coeff_b * Z6);
    }
}
 
bool bn128_G1::is_in_safe_subgroup() const { return true; }
 
void bn128_G1::write_uncompressed(std::ostream &out) const
{
    bn128_G1 gcopy(*this);
    gcopy.to_affine_coordinates();
 
    out << (gcopy.is_zero() ? '1' : '0') << OUTPUT_SEPARATOR;
 
#ifndef BINARY_OUTPUT
    out << gcopy.X << OUTPUT_SEPARATOR << gcopy.Y;
#else
    out.write((char *)&gcopy.X, sizeof(gcopy.X));
    out.write((char *)&gcopy.Y, sizeof(gcopy.Y));
#endif
}
 
void bn128_G1::write_compressed(std::ostream &out) const
{
    bn128_G1 gcopy(*this);
    gcopy.to_affine_coordinates();
 
    out << (gcopy.is_zero() ? '1' : '0') << OUTPUT_SEPARATOR;
 
#ifndef BINARY_OUTPUT
    out << gcopy.X;
#else
    out.write((char *)&gcopy.X, sizeof(gcopy.X));
#endif
    out << OUTPUT_SEPARATOR << (((unsigned char *)&gcopy.Y)[0] & 1 ? '1' : '0');
}
 
void bn128_G1::read_uncompressed(std::istream &in, bn128_G1 &g)
{
    char is_zero;
    in.read((char *)&is_zero, 1); // this reads is_zero;
    is_zero -= '0';
    consume_OUTPUT_SEPARATOR(in);
 
#ifndef BINARY_OUTPUT
    in >> g.X;
    consume_OUTPUT_SEPARATOR(in);
    in >> g.Y;
#else
    in.read((char *)&g.X, sizeof(g.X));
    in.read((char *)&g.Y, sizeof(g.Y));
#endif
 
    /* finalize */
    if (!is_zero) {
        g.Z = bn::Fp(1);
    } else {
        g = bn128_G1::zero();
    }
}
 
void bn128_G1::read_compressed(std::istream &in, bn128_G1 &g)
{
    char is_zero;
    // this reads is_zero;
    in.read((char *)&is_zero, 1);
    is_zero -= '0';
    consume_OUTPUT_SEPARATOR(in);
 
    bn::Fp tX;
#ifndef BINARY_OUTPUT
    in >> tX;
#else
    in.read((char *)&tX, sizeof(tX));
#endif
    consume_OUTPUT_SEPARATOR(in);
    unsigned char Y_lsb;
    in.read((char *)&Y_lsb, 1);
    Y_lsb -= '0';
 
    // y = +/- sqrt(x^3 + b)
    if (!is_zero) {
        g.X = tX;
        bn::Fp tX2, tY2;
        bn::Fp::square(tX2, tX);
        bn::Fp::mul(tY2, tX2, tX);
        bn::Fp::add(tY2, tY2, bn128_coeff_b);
 
        g.Y = bn128_G1::sqrt(tY2);
        if ((((unsigned char *)&g.Y)[0] & 1) != Y_lsb) {
            bn::Fp::neg(g.Y, g.Y);
        }
    }
 
    /* finalize */
    if (!is_zero) {
        g.Z = bn::Fp(1);
    } else {
        g = bn128_G1::zero();
    }
}
 
std::ostream &operator<<(std::ostream &out, const bn128_G1 &g)
{
#ifdef NO_PT_COMPRESSION
    g.write_uncompressed(out);
#else
    g.write_compressed(out);
#endif
    return out;
}
 
std::istream &operator>>(std::istream &in, bn128_G1 &g)
{
#ifdef NO_PT_COMPRESSION
    bn128_G1::read_uncompressed(in, g);
#else
    bn128_G1::read_compressed(in, g);
#endif
    return in;
}
 
void bn128_G1::batch_to_special_all_non_zeros(std::vector<bn128_G1> &vec)
{
    std::vector<bn::Fp> Z_vec;
    Z_vec.reserve(vec.size());
 
    for (auto &el : vec) {
        Z_vec.emplace_back(el.Z);
    }
    bn_batch_invert<bn::Fp>(Z_vec);
 
    const bn::Fp one = 1;
 
    for (size_t i = 0; i < vec.size(); ++i) {
        bn::Fp Z2, Z3;
        bn::Fp::square(Z2, Z_vec[i]);
        bn::Fp::mul(Z3, Z2, Z_vec[i]);
 
        bn::Fp::mul(vec[i].X, vec[i].X, Z2);
        bn::Fp::mul(vec[i].Y, vec[i].Y, Z3);
        vec[i].Z = one;
    }
}
 
} // namespace libff