start to imple GLV for G1

Makefile

@@ -3,7 +3,7 @@ LIB_DIR=lib
 OBJ_DIR=obj
 EXE_DIR=bin
 SRC_SRC=fp.cpp
-TEST_SRC=fp_test.cpp ec_test.cpp fp_util_test.cpp window_method_test.cpp elgamal_test.cpp fp_tower_test.cpp gmp_test.cpp bn_test.cpp bn256_test.cpp bn384_test.cpp
+TEST_SRC=fp_test.cpp ec_test.cpp fp_util_test.cpp window_method_test.cpp elgamal_test.cpp fp_tower_test.cpp gmp_test.cpp bn_test.cpp bn256_test.cpp bn384_test.cpp glv_test.cpp
 ifeq ($(CPU),x86-64)
   MCL_USE_XBYAK?=1
   TEST_SRC+=mont_fp_test.cpp sq_test.cpp

include/mcl/bn.hpp

@@ -199,6 +199,111 @@ struct MapToT {
 	}
 };
 
+/*
+	Skew Frobenius Map and Efficient Scalar Multiplication for Pairing-Based Cryptography
+	Y. Sakemi, Y. Nogami, K. Okeya, H. Kato, Y. Morikawa
+*/
+template<class Fp>
+struct GLV {
+	typedef mcl::EcT<Fp> G1;
+	Fp w; // (-1 + sqrt(-3)) / 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)
+	{
+		if (!Fp::squareRoot(w, -3)) throw cybozu::Exception("GLV:init");
+		w = (w - 1) / 2;
+		this->r = r;
+		this->z = z;
+		this->isNegative = isNegative;
+		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;
+	}
+	/*
+		s = ap^2 + b mod r
+		assume(s < r);
+	*/
+	void getAB(mpz_class& a, mpz_class& b, const mpz_class& s) 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;
+	}
+	void mul(G1& Q, G1 P, mpz_class x) const
+	{
+		if (x == 0) {
+			Q.clear();
+			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
+		G1 A, B;
+		mulP2(A, P);
+		if (b < 0) {
+			b = -b;
+			G1::neg(P, P);
+		}
+		assert(a >= 0);
+		assert(b >= 0);
+#if 1
+		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[4];
+		tbl[1] = A; tbl[1].normalize();
+		tbl[2] = P; tbl[2].normalize();
+		tbl[3] = A + P; tbl[3].normalize();
+		Q.clear();
+		for (int i = (int)n - 1; i >= 0; i--) {
+			G1::dbl(Q, Q);
+			bool ai = mcl::gmp::testBit(a, i);
+			bool bi = mcl::gmp::testBit(b, i);
+			int c = bi * 2 + ai;
+			if (c > 0) {
+				Q += tbl[c];
+			}
+		}
+#else
+		G1::mul(A, A, a);
+		G1::mul(B, P, b);
+		G1::add(Q, A, B);
+#endif
+	}
+};
+
 template<class Fp>
 struct ParamT {
 	typedef Fp2T<Fp> Fp2;

include/mcl/ec.hpp

@@ -61,6 +61,7 @@ public:
 	*/
 	static bool verifyOrder_;
 	static mpz_class order_;
+	static void (*mulArrayGLV)(EcT& z, const EcT& x, const fp::Unit *y, size_t yn, bool isNegative, bool constTime);
 	/* default constructor is undefined value */
 	EcT() {}
 	EcT(const Fp& _x, const Fp& _y)
@@ -786,6 +787,7 @@ template<class Fp> int EcT<Fp>::specialA_;
 template<class Fp> bool EcT<Fp>::compressedExpression_;
 template<class Fp> bool EcT<Fp>::verifyOrder_;
 template<class Fp> mpz_class EcT<Fp>::order_;
+template<class Fp> static void (*mulArrayGLV)(EcT<Fp>& z, const EcT<Fp>& x, const fp::Unit *y, size_t yn, bool isNegative, bool constTime);
 #ifndef MCL_EC_USE_AFFINE
 template<class Fp> int EcT<Fp>::mode_;
 #endif

readme.md

@@ -215,6 +215,8 @@ This library contains [mie](https://github.com/herumi/mie/) and [Lifted-ElGamal]
   Pairing 2010, ([preprint](http://eprint.iacr.org/2010/354))
 * [_Faster hashing to G2_](http://dx.doi.org/10.1007/978-3-642-28496-0_25),Laura Fuentes-Castañeda,  Edward Knapp,  Francisco Rodríguez-Henríquez,
   SAC 2011, ([preprint](https://eprint.iacr.org/2008/530))
+* [_Skew Frobenius Map and Efficient Scalar Multiplication for Pairing–Based Cryptography_](https://www.researchgate.net/publication/221282560_Skew_Frobenius_Map_and_Efficient_Scalar_Multiplication_for_Pairing-Based_Cryptography),
+Y. Sakemi, Y. Nogami, K. Okeya, Y. Morikawa, CANS 2008.
 
 # Author
 

test/glv_test.cpp

@@ -0,0 +1,44 @@
+#include <cybozu/test.hpp>
+#include <cybozu/xorshift.hpp>
+#include <cybozu/benchmark.hpp>
+
+#if 1
+#include <mcl/bn384.hpp>
+using namespace mcl::bn384;
+#else
+#include <mcl/bn256.hpp>
+using namespace mcl::bn256;
+#endif
+
+#define PUT(x) std::cout << #x "=" << (x) << std::endl;
+
+void testGLV(const mcl::bn::CurveParam& cp)
+{
+	bn384init(cp);
+	G1::setCompressedExpression(false);
+
+	G1 P0, P1, P2;
+	cybozu::XorShift rg;
+	mcl::bn::GLV<Fp> glv;
+	glv.init(BN::param.r, BN::param.z, BN::param.isNegative);
+	for (int i = 1; i < 100; i++) {
+		BN::mapToG1(P0, i);
+		Fr s;
+		s.setRand(rg);
+		mpz_class ss = s.getMpz();
+		G1::mul(P1, P0, ss);
+		glv.mul(P2, P0, ss);
+		CYBOZU_TEST_EQUAL(P1, P2);
+	}
+	Fr s;
+	BN::mapToG1(P0, 123);
+	CYBOZU_BENCH_C("Ec::mul", 100, P1 = P0; s.setRand(rg); G1::mul, P2, P1, s.getMpz());
+	CYBOZU_BENCH_C("Ec::glv", 100, P1 = P0; s.setRand(rg); glv.mul, P2, P1, s.getMpz());
+}
+
+CYBOZU_TEST_AUTO(glv)
+{
+	testGLV(mcl::bn::CurveFp382_1);
+	testGLV(mcl::bn::CurveFp382_2);
+	testGLV(mcl::bn::CurveFp254BNb);
+}