a little optimization of Ec::addJacobi if both points are affine

3 files changed, 90 insertions, 40 deletions

diff --git a/include/mcl/ec.hpp b/include/mcl/ec.hpp
index 8ebf7e7..237b6c3 100644
--- a/include/mcl/ec.hpp
+++ b/include/mcl/ec.hpp
@@ -23,11 +23,11 @@ namespace mcl {
 namespace ec {
 
 enum Mode {
-	Jacobi,
-	Proj
+	Jacobi = 0,
+	Proj = 1
 };
 
-} // mcl::ecl
+} // mcl::ec
 
 /*
 	elliptic curve
@@ -423,27 +423,41 @@ public:
 		dblNoVerifyInf(R, P);
 	}
 #ifndef MCL_EC_USE_AFFINE
-	static inline void addJacobi(EcT& R, const EcT& P, const EcT& Q)
+	static inline void addJacobi(EcT& R, const EcT& P, const EcT& Q, bool isPzOne, bool isQzOne)
 	{
-		const bool isQzOne = Q.z.isOne();
 		Fp r, U1, S1, H, H3;
-		Fp::sqr(r, P.z);
+		if (isPzOne) {
+			// r = 1;
+		} else {
+			Fp::sqr(r, P.z);
+		}
 		if (isQzOne) {
 			U1 = P.x;
-			Fp::mul(H, Q.x, r);
+			if (isPzOne) {
+				H = Q.x;
+			} else {
+				Fp::mul(H, Q.x, r);
+			}
 			H -= U1;
-			r *= P.z;
 			S1 = P.y;
 		} else {
 			Fp::sqr(S1, Q.z);
 			Fp::mul(U1, P.x, S1);
-			Fp::mul(H, Q.x, r);
+			if (isPzOne) {
+				H = Q.x;
+			} else {
+				Fp::mul(H, Q.x, r);
+			}
 			H -= U1;
-			r *= P.z;
 			S1 *= Q.z;
 			S1 *= P.y;
 		}
-		r *= Q.y;
+		if (isPzOne) {
+			r = Q.y;
+		} else {
+			r *= P.z;
+			r *= Q.y;
+		}
 		r -= S1;
 		if (H.isZero()) {
 			if (r.isZero()) {
@@ -453,11 +467,13 @@ public:
 			}
 			return;
 		}
-		if (isQzOne) {
-			Fp::mul(R.z, P.z, H);
+		if (isPzOne) {
+			R.z = H;
 		} else {
-			Fp::mul(R.z, P.z, Q.z);
-			R.z *= H;
+			Fp::mul(R.z, P.z, H);
+		}
+		if (!isQzOne) {
+			R.z *= Q.z;
 		}
 		Fp::sqr(H3, H); // H^2
 		Fp::sqr(R.y, r); // r^2
@@ -549,14 +565,17 @@ public:
 #else
 		const EcT *pP = &P0;
 		const EcT *pQ = &Q0;
+		bool isPzOne = P0.z.isOne();
+		bool isQzOne = Q0.z.isOne();
 		if (pP->z.isOne()) {
 			fp::swap_(pP, pQ);
+			std::swap(isPzOne, isQzOne);
 		}
 		const EcT& P(*pP);
 		const EcT& Q(*pQ);
 		switch (mode_) {
 		case ec::Jacobi:
-			addJacobi(R, P, Q);
+			addJacobi(R, P, Q, isPzOne, isQzOne);
 			break;
 		case ec::Proj:
 			addProj(R, P, Q);
diff --git a/test/bench.hpp b/test/bench.hpp
index e703194..b12947a 100644
--- a/test/bench.hpp
+++ b/test/bench.hpp
@@ -1,5 +1,58 @@
 #include <mcl/lagrange.hpp>
 
+void benchAddDblG1()
+{
+	puts("benchAddDblG1");
+	const int C = 100000;
+	G1 P1, P2, P3;
+	hashAndMapToG1(P1, "a");
+	hashAndMapToG1(P2, "b");
+	P1 += P2;
+	P2 += P1;
+	printf("z.isOne()=%d %d\n", P1.z.isOne(), P2.z.isOne());
+	CYBOZU_BENCH_C("G1::add(1)", C, G1::add, P3, P1, P2);
+	P1.normalize();
+	printf("z.isOne()=%d %d\n", P1.z.isOne(), P2.z.isOne());
+	CYBOZU_BENCH_C("G1::add(2)", C, G1::add, P3, P1, P2);
+	CYBOZU_BENCH_C("G1::add(3)", C, G1::add, P3, P2, P1);
+	P2.normalize();
+	printf("z.isOne()=%d %d\n", P1.z.isOne(), P2.z.isOne());
+	CYBOZU_BENCH_C("G1::add(4)", C, G1::add, P3, P1, P2);
+	P1 = P3;
+	printf("z.isOne()=%d\n", P1.z.isOne());
+	CYBOZU_BENCH_C("G1::dbl(1)", C, G1::dbl, P3, P1);
+	P1.normalize();
+	printf("z.isOne()=%d\n", P1.z.isOne());
+	CYBOZU_BENCH_C("G1::dbl(2)", C, G1::dbl, P3, P1);
+}
+
+void benchAddDblG2()
+{
+	puts("benchAddDblG2");
+	const int C = 100000;
+	G2 P1, P2, P3;
+	hashAndMapToG2(P1, "a");
+	hashAndMapToG2(P2, "b");
+	P1 += P2;
+	P2 += P1;
+	printf("z.isOne()=%d %d\n", P1.z.isOne(), P2.z.isOne());
+	CYBOZU_BENCH_C("G2::add(1)", C, G2::add, P3, P1, P2);
+	P1.normalize();
+	printf("z.isOne()=%d %d\n", P1.z.isOne(), P2.z.isOne());
+	CYBOZU_BENCH_C("G2::add(2)", C, G2::add, P3, P1, P2);
+	CYBOZU_BENCH_C("G2::add(3)", C, G2::add, P3, P2, P1);
+	P2.normalize();
+	printf("z.isOne()=%d %d\n", P1.z.isOne(), P2.z.isOne());
+	CYBOZU_BENCH_C("G2::add(4)", C, G2::add, P3, P1, P2);
+	P1 = P3;
+	printf("z.isOne()=%d\n", P1.z.isOne());
+	CYBOZU_BENCH_C("G2::dbl(1)", C, G2::dbl, P3, P1);
+	P1.normalize();
+	printf("z.isOne()=%d\n", P1.z.isOne());
+	CYBOZU_BENCH_C("G2::dbl(2)", C, G2::dbl, P3, P1);
+}
+
+
 void testBench(const G1& P, const G2& Q)
 {
 	G1 Pa;
@@ -19,34 +72,10 @@ void testBench(const G1& P, const G2& Q)
 	CYBOZU_BENCH_C("G1::mul       ", C, G1::mul, Pa, Pa, a);
 	CYBOZU_BENCH_C("G1::add       ", C, G1::add, Pa, Pa, P);
 	CYBOZU_BENCH_C("G1::dbl       ", C, G1::dbl, Pa, Pa);
-	{
-		G1 P1 = Pa, P2 = Pa + P, P3;
-		CYBOZU_BENCH_C("G1::add(1)    ", C3, G1::add, P3, P1, P2);
-		P1.normalize();
-		CYBOZU_BENCH_C("G1::add(2)    ", C3, G1::add, P3, P1, P2);
-		P2.normalize();
-		CYBOZU_BENCH_C("G1::add(3)    ", C3, G1::add, P3, P1, P2);
-		P1 = P3;
-		CYBOZU_BENCH_C("G1::dbl(1)    ", C3, G1::dbl, P3, P1);
-		P1.normalize();
-		CYBOZU_BENCH_C("G1::dbl(2)    ", C3, G1::dbl, P3, P1);
-	}
 	CYBOZU_BENCH_C("G2::mulCT     ", C, G2::mulCT, Qa, Q, a);
 	CYBOZU_BENCH_C("G2::mul       ", C, G2::mul, Qa, Qa, a);
 	CYBOZU_BENCH_C("G2::add       ", C, G2::add, Qa, Qa, Q);
 	CYBOZU_BENCH_C("G2::dbl       ", C, G2::dbl, Qa, Qa);
-	{
-		G2 Q1 = Qa, Q2 = Qa + Q, Q3;
-		CYBOZU_BENCH_C("G2::add(1)    ", C3, G2::add, Q3, Q1, Q2);
-		Q1.normalize();
-		CYBOZU_BENCH_C("G2::add(2)    ", C3, G2::add, Q3, Q1, Q2);
-		Q2.normalize();
-		CYBOZU_BENCH_C("G2::add(3)    ", C3, G2::add, Q3, Q1, Q2);
-		Q1 = Q3;
-		CYBOZU_BENCH_C("G2::dbl(1)    ", C3, G2::dbl, Q3, Q1);
-		Q1.normalize();
-		CYBOZU_BENCH_C("G2::dbl(2)    ", C3, G2::dbl, Q3, Q1);
-	}
 	CYBOZU_BENCH_C("GT::pow       ", C, GT::pow, e1, e1, a);
 //	CYBOZU_BENCH_C("GT::powGLV    ", C, BN::param.glv2.pow, e1, e1, a);
 	G1 PP;
diff --git a/test/bn_test.cpp b/test/bn_test.cpp
index 929e235..4ee2980 100644
--- a/test/bn_test.cpp
+++ b/test/bn_test.cpp
@@ -375,6 +375,8 @@ CYBOZU_TEST_AUTO(naive)
 		testPrecomputed(P, Q);
 		testMillerLoop2(P, Q);
 		testBench(P, Q);
+		benchAddDblG1();
+		benchAddDblG2();
 	}
 	int count = (int)clk.getCount();
 	if (count) {