diff --git a/include/mcl/bn.hpp b/include/mcl/bn.hpp index 4dae1df..fa2b0fd 100644 --- a/include/mcl/bn.hpp +++ b/include/mcl/bn.hpp @@ -206,54 +206,46 @@ struct MapToT { template struct GLV { typedef mcl::EcT G1; - Fp w; // (-1 + sqrt(-3)) / 2 + Fp rw; // rw = 1 / w = (-1 - sqrt(-3)) / 2 + int m; + mpz_class v0, v1; + mpz_class B[2][2]; mpz_class r; - mpz_class z; - bool isNegative; - mpz_class v; // 6z^2 + 4z + 1 > 0 - mpz_class c; // 2z + 1 - void init(const mpz_class& r, const mpz_class& z, bool isNegative) + void init(const mpz_class& r, const mpz_class& z) { - if (!Fp::squareRoot(w, -3)) throw cybozu::Exception("GLV:init"); - w = (w - 1) / 2; + if (!Fp::squareRoot(rw, -3)) throw cybozu::Exception("GLV:init"); + rw = (-1 - rw) / 2; this->r = r; - this->z = z; - this->isNegative = isNegative; - v = 1 + z * (4 + z * 6); - c = 2 * z + 1; + m = gmp::getBitSize(r); + m = (m + fp::UnitBitSize - 1) & ~(fp::UnitBitSize - 1);// a little better size + v0 = ((6 * z * z + 4 * z + 1) << m) / r; + v1 = ((-2 * z - 1) << m) / r; + B[0][0] = 6 * z * z + 2 * z; + B[0][1] = -2 * z - 1; + B[1][0] = -2 * z - 1; + B[1][1] = -6 * z * z - 4 * z - 1; } /* - (p^2 mod r) (x, y) = (wx, -y) + lambda = 36z^4 - 1 + lambda (x, y) = (rw x, y) */ - void mulP2(G1& Q, const G1& P) const + void mulLambda(G1& Q, const G1& P) const { - Fp::mul(Q.x, P.x, w); - Fp::neg(Q.y, P.y); + Fp::mul(Q.x, P.x, rw); + Q.y = P.y; Q.z = P.z; } /* - s = ap^2 + b mod r - assume(s < r); + lambda = 36 z^4 - 1 + x = a + b * lambda mod r */ - void getAB(mpz_class& a, mpz_class& b, const mpz_class& s) const + void split(mpz_class& a, mpz_class& b, const mpz_class& x) const { - assert(0 < s && s < r); - /* - s = s1 * v + s2 // s1 = s / v, s2 = s % v - = s1 * c * p^2 + s2 // vP = cp^2 P - = (s3 * v + s4) * p^2 + s2 // s3 = (s1 * c) / v, s4 = (s1 * c) % v - = (s3 * c * p^2 + s4) * p^2 + s2 - = (s3 * c) * p^4 + s4 * p^2 + s2 // s5 = s3 * c, p^4 = p^2 - 1 - = s5 * (p^2 - 1) + s4 * p^2 + s2 - = (s4 + s5) * p^2 + (s2 - s5) - */ mpz_class t; - mcl::gmp::divmod(a, t, s, v); // a = t / v, t = t % v - a *= c; - mcl::gmp::divmod(b, a, a, v); // b = a / v, a = a % v - b *= c; - a += b; - b = t - b; + t = (x * v0) >> m; + b = (x * v1) >> m; + a = x - (t * B[0][0] + b * B[1][0]); + b = - (t * B[0][1] + b * B[1][1]); } void mul(G1& Q, G1 P, mpz_class x, bool constTime = false) const { @@ -263,72 +255,32 @@ struct GLV { return; } if (x < 0) { -// G1::neg(P, P); -// x = -x; x += r; } mpz_class a, b; - getAB(a, b, x); - // Q = (ap^2 + b)P + split(a, b, x); G1 A; - mulP2(A, P); + if (a < 0) { + G1::neg(A, P); + a = -a; + } else { + A = P; + } if (b < 0) { + G1::neg(Q, P); b = -b; - G1::neg(P, P); - } - assert(a >= 0); - assert(b >= 0); -#if 0 - G1::mul(A, A, a); - G1 B; - G1::mul(B, P, b); - G1::add(Q, A, B); - return; -#endif -#if 0 // slow - size_t nA = mcl::gmp::getBitSize(a); - size_t nB = mcl::gmp::getBitSize(b); - size_t n = std::max(nA, nB); - assert(n > 0); - G1 tbl[16]; - tbl[0].clear(); - tbl[1] = A; - G1::dbl(tbl[2], tbl[1]); - G1::add(tbl[3], tbl[2], tbl[1]); - for (int i = 1; i < 4; i++) { - for (int j = 0; j < 4; j++) { - G1::add(tbl[i * 4 + j], tbl[(i - 1) * 4 + j], P); - } - } - for (int i = 1; i < 16; i++) { - tbl[i].normalize(); - } - if (n & 1) { - n--; - bool ai = mcl::gmp::testBit(a, n); - bool bi = mcl::gmp::testBit(b, n); - unsigned int idx = bi * 4 + ai; - Q = tbl[idx]; - if (n == 0) return; } else { - Q.clear(); - } - for (int i = (int)n - 2; i >= 0; i -= 2) { - G1::dbl(Q, Q); - G1::dbl(Q, Q); - bool a0 = mcl::gmp::testBit(a, i + 0); - bool a1 = mcl::gmp::testBit(a, i + 1); - bool b0 = mcl::gmp::testBit(b, i + 0); - bool b1 = mcl::gmp::testBit(b, i + 1); - unsigned int c = b1 * 8 + b0 * 4 + a1 * 2 + a0; - if (c > 0) { - Q += tbl[c]; - } + Q = P; } + mulLambda(Q, Q); +#if 0 + G1::mulBase(A, A, a); + G1::mulBase(Q, Q, b); + Q += A; #else A.normalize(); - P.normalize(); - G1 tbl[4] = { A, A, P, A + P }; // tbl[0] : dummy + Q.normalize(); + G1 tbl[4] = { A, A, Q, A + Q }; // tbl[0] : dummy tbl[3].normalize(); typedef mcl::fp::Unit Unit; const int aN = (int)mcl::gmp::getUnitSize(a); @@ -442,7 +394,7 @@ struct ParamT { G2::init(0, b_div_xi, mcl::ec::Proj); G2::setOrder(r); mapTo.init(2 * p - r); - glv.init(r, z, isNegative); + glv.init(r, z); Fp2::pow(g[0], xi, (p - 1) / 6); // g = xi^((p-1)/6) for (size_t i = 1; i < gN; i++) { diff --git a/test/glv_test.cpp b/test/glv_test.cpp index 732b1c3..b2df070 100644 --- a/test/glv_test.cpp +++ b/test/glv_test.cpp @@ -13,6 +13,103 @@ using namespace mcl::bn256; #define PUT(x) std::cout << #x "=" << (x) << std::endl; +/* + Skew Frobenius Map and Efficient Scalar Multiplication for Pairing-Based Cryptography + Y. Sakemi, Y. Nogami, K. Okeya, H. Kato, Y. Morikawa +*/ +struct oldGLV { + Fp w; // (-1 + sqrt(-3)) / 2 + mpz_class r; + mpz_class v; // 6z^2 + 4z + 1 > 0 + mpz_class c; // 2z + 1 + void init(const mpz_class& r, const mpz_class& z) + { + if (!Fp::squareRoot(w, -3)) throw cybozu::Exception("GLV:init"); + w = (w - 1) / 2; + this->r = r; + v = 1 + z * (4 + z * 6); + c = 2 * z + 1; + } + /* + (p^2 mod r) (x, y) = (wx, -y) + */ + void mulP2(G1& Q, const G1& P) const + { + Fp::mul(Q.x, P.x, w); + Fp::neg(Q.y, P.y); + Q.z = P.z; + } + /* + x = ap^2 + b mod r + assume(x < r); + */ + void split(mpz_class& a, mpz_class& b, const mpz_class& x) const + { + assert(0 < x && x < r); + /* + x = s1 * v + s2 // s1 = x / v, s2 = x % v + = s1 * c * p^2 + s2 // vP = cp^2 P + = (s3 * v + s4) * p^2 + s2 // s3 = (s1 * c) / v, s4 = (s1 * c) % v + = (s3 * c * p^2 + s4) * p^2 + s2 + = (s3 * c) * p^4 + s4 * p^2 + s2 // s5 = s3 * c, p^4 = p^2 - 1 + = s5 * (p^2 - 1) + s4 * p^2 + s2 + = (s4 + s5) * p^2 + (s2 - s5) + */ + mpz_class t; + mcl::gmp::divmod(a, t, x, v); // a = t / v, t = t % v + a *= c; + mcl::gmp::divmod(b, a, a, v); // b = a / v, a = a % v + b *= c; + a += b; + b = t - b; + } + template + void mul(G1& Q, const G1& P, const mpz_class& x) const + { + G1 A, B; + mpz_class a, b; + split(a, b, x); + mulP2(A, P); + G1::mul(A, A, a); + G1::mul(B, P, b); + G1::add(Q, A, B); + } +}; + +template +void compareLength(const GLV1& rhs, const GLV2& lhs) +{ + cybozu::XorShift rg; + int Rc = 0; + int Lc = 0; + int eq = 0; + mpz_class R0, R1, L0, L1, x; + Fr r; + for (int i = 1; i < 1000; i++) { + r.setRand(rg); + x = r.getMpz(); + rhs.split(R0, R1, x); + lhs.split(L0, L1, x); + + size_t R0n = mcl::gmp::getBitSize(R0); + size_t R1n = mcl::gmp::getBitSize(R1); + size_t L0n = mcl::gmp::getBitSize(L0); + size_t L1n = mcl::gmp::getBitSize(L1); + size_t Rn = std::max(R0n, R1n); + size_t Ln = std::max(L0n, L1n); + if (Rn == Ln) { + eq++; + } + if (Rn > Ln) { + Rc++; + } + if (Rn < Ln) { + Lc++; + } + } + printf("eq=%d small is better rhs=%d, lhs=%d\n", eq, Rc, Lc); +} + void testGLV(const mcl::bn::CurveParam& cp) { bn384init(cp); @@ -21,8 +118,14 @@ void testGLV(const mcl::bn::CurveParam& cp) G1 P0, P1, P2; BN::mapToG1(P0, 1); cybozu::XorShift rg; + + oldGLV oldGlv; + oldGlv.init(BN::param.r, BN::param.z); + mcl::bn::GLV glv; - glv.init(BN::param.r, BN::param.z, BN::param.isNegative); + glv.init(BN::param.r, BN::param.z); + compareLength(glv, oldGlv); + for (int i = 1; i < 100; i++) { BN::mapToG1(P0, i); Fr s; @@ -33,25 +136,11 @@ void testGLV(const mcl::bn::CurveParam& cp) CYBOZU_TEST_EQUAL(P1, P2); glv.mul(P2, P0, ss, true); CYBOZU_TEST_EQUAL(P1, P2); - } - for (int i = -100; i < 100; i++) { - mpz_class ss = i; - G1::mulBase(P1, P0, ss); - glv.mul(P2, P0, ss); - CYBOZU_TEST_EQUAL(P1, P2); - glv.mul(P2, P0, ss, true); + oldGlv.mul(P2, P0, ss); CYBOZU_TEST_EQUAL(P1, P2); } for (int i = -100; i < 100; i++) { - mpz_class ss = i * glv.v; - G1::mulBase(P1, P0, ss); - glv.mul(P2, P0, ss); - CYBOZU_TEST_EQUAL(P1, P2); - glv.mul(P2, P0, ss, true); - CYBOZU_TEST_EQUAL(P1, P2); - } - for (int i = -100; i < 100; i++) { - mpz_class ss = -i * glv.v + -i * glv.c; + mpz_class ss = i; G1::mulBase(P1, P0, ss); glv.mul(P2, P0, ss); CYBOZU_TEST_EQUAL(P1, P2); @@ -66,7 +155,7 @@ void testGLV(const mcl::bn::CurveParam& cp) CYBOZU_TEST_AUTO(glv) { + testGLV(mcl::bn::CurveFp254BNb); testGLV(mcl::bn::CurveFp382_1); testGLV(mcl::bn::CurveFp382_2); - testGLV(mcl::bn::CurveFp254BNb); }