diff --git a/include/mcl/bn.hpp b/include/mcl/bn.hpp
index 58bfc27..1a51c03 100644
--- a/include/mcl/bn.hpp
+++ b/include/mcl/bn.hpp
@@ -541,6 +541,13 @@ struct BNT {
 		addLineWithoutP(l, R, Q);
 		updateLine(l, P);
 	}
+	static void Fp6_cb_mul_G1_xy(Fp6& l, const G1& P)
+	{
+		assert(P.isNormalized());
+		Fp2::mulFp(l.c, l.c, P.x);
+		Fp2::mulFp(l.b, l.b, P.y);
+	}
+
 	static void convertFp6toFp12(Fp12& y, const Fp6& x)
 	{
 		y.clear();
@@ -1001,21 +1008,22 @@ struct BNT {
 #endif
 		exp_d1(y, y);
 	}
-	static void pairing(Fp12& f, const G2& Q, const G1& P)
+	static void millerLoop(Fp12& f, const G2& Q, const G1& P)
 	{
 		P.normalize();
 		Q.normalize();
-		Fp6 l;
 		G2 T = Q;
-		f = 1;
 		G2 negQ;
-		G2::neg(negQ, Q);
+		if (param.useNAF) {
+			G2::neg(negQ, Q);
+		}
 		Fp6 d;
 		dblLine(d, T, P);
 		Fp6 e;
 		assert(param.siTbl[1] == 1);
 		addLine(e, T, Q, P);
 		mul_024_024(f, d, e);
+		Fp6 l;
 		for (size_t i = 2; i < param.siTbl.size(); i++) {
 			dblLine(l, T, P);
 			Fp12::sqr(f, f);
@@ -1042,8 +1050,86 @@ struct BNT {
 		Fp12 ft;
 		mul_024_024(ft, d, e);
 		f *= ft;
+	}
+	static void pairing(Fp12& f, const G2& Q, const G1& P)
+	{
+		millerLoop(f, Q, P);
 		finalExp(f, f);
 	}
+	static void precomputeG2(std::vector<Fp6>& Qcoeff, const G2& Q)
+	{
+		Qcoeff.clear();
+		Q.normalize();
+		G2 T = Q;
+		G2 negQ;
+		if (param.useNAF) {
+			G2::neg(negQ, Q);
+		}
+		Fp6 d;
+		dblLineWithoutP(d, T);
+		Qcoeff.push_back(d);
+		Fp6 e;
+		assert(param.siTbl[1] == 1);
+		addLineWithoutP(e, T, Q);
+		Qcoeff.push_back(e);
+		Fp6 l;
+		for (size_t i = 2; i < param.siTbl.size(); i++) {
+			dblLineWithoutP(l, T);
+			Qcoeff.push_back(l);
+			if (param.siTbl[i]) {
+				if (param.siTbl[i] > 0) {
+					addLineWithoutP(l, T, Q);
+				} else {
+					addLineWithoutP(l, T, negQ);
+				}
+				Qcoeff.push_back(l);
+			}
+		}
+		G2 Q1, Q2;
+		FrobeniusOnTwist(Q1, Q);
+		FrobeniusOnTwist(Q2, Q1);
+		G2::neg(Q2, Q2);
+		if (param.z < 0) {
+			G2::neg(T, T);
+		}
+		addLineWithoutP(d, T, Q1);
+		Qcoeff.push_back(d);
+		addLineWithoutP(e, T, Q2);
+		Qcoeff.push_back(e);
+	}
+	static void precomputedMillerLoop(Fp12& f, const std::vector<Fp6>& Qcoeff, const G1& P)
+	{
+		P.normalize();
+		size_t idx = 0;
+		Fp6 d = Qcoeff[idx++];
+		Fp6_cb_mul_G1_xy(d, P);
+
+		Fp6 e = Qcoeff[idx++];
+		Fp6_cb_mul_G1_xy(e, P);
+		mul_024_024(f, d, e);
+		Fp6 l;
+		for (size_t i = 2; i < param.siTbl.size(); i++) {
+			l = Qcoeff[idx++];
+			Fp6_cb_mul_G1_xy(l, P);
+			Fp12::sqr(f, f);
+			mul_024(f, f, l);
+			if (param.siTbl[i]) {
+				l = Qcoeff[idx++];
+				Fp6_cb_mul_G1_xy(l, P);
+				mul_024(f, f, l);
+			}
+		}
+		if (param.z < 0) {
+			Fp6::neg(f.b, f.b);
+		}
+		d = Qcoeff[idx++];
+		Fp6_cb_mul_G1_xy(d, P);
+		e = Qcoeff[idx++];
+		Fp6_cb_mul_G1_xy(e, P);
+		Fp12 ft;
+		mul_024_024(ft, d, e);
+		f *= ft;
+	}
 };
 
 template<class Fp>
diff --git a/test/bn_test.cpp b/test/bn_test.cpp
index 01f4bcf..3c43aa6 100644
--- a/test/bn_test.cpp
+++ b/test/bn_test.cpp
@@ -161,7 +161,18 @@ void testCompress()
 	CYBOZU_TEST_EQUAL(b, c);
 }
 
-void test(const TestSet& ts)
+void testPrecomputed(const G2& Q, const G1& P)
+{
+	Fp12 e1, e2;
+	BN::pairing(e1, Q, P);
+	std::vector<Fp6> Qcoeff;
+	BN::precomputeG2(Qcoeff, Q);
+	BN::precomputedMillerLoop(e2, Qcoeff, P);
+	BN::finalExp(e2, e2);
+	CYBOZU_TEST_EQUAL(e1, e2);
+}
+
+void testPairing(const TestSet& ts)
 {
 	G1 P(ts.g1.a, ts.g1.b);
 	G2 Q(Fp2(ts.g2.aa, ts.g2.ab), Fp2(ts.g2.ba, ts.g2.bb));
@@ -171,16 +182,8 @@ void test(const TestSet& ts)
 	{
 		std::stringstream ss(ts.e);
 		ss >> e2;
-//		mpz_class x = BN::param.z;
-//		x = 2 * x * (6 * x * x + 3 * x + 1);
-//		Fp12::pow(e1, e1, x);
 	}
 	CYBOZU_TEST_EQUAL(e1, e2);
-	/*
-		ate-pairing on Haswell
-		miller loop : 700Kclk
-		final exp   : 460Kclk
-	*/
 #if 0
 	for (int i = 0; i < 1000; i++) BN::pairing(e1, Q, P);
 //	CYBOZU_BENCH_C("pairing", 1000, BN::pairing, e1, Q, P); // 2.4Mclk
@@ -209,6 +212,7 @@ void test(const TestSet& ts)
 	CYBOZU_BENCH("pairing", BN::pairing, e1, Q, P); // 2.4Mclk
 	CYBOZU_BENCH("finalExp", BN::finalExp, e1, e1); // 1.3Mclk
 #endif
+	testPrecomputed(Q, P);
 }
 
 CYBOZU_TEST_AUTO(naive)
@@ -222,7 +226,7 @@ CYBOZU_TEST_AUTO(naive)
 		testMapToG2();
 		testCyclotomic();
 		testCompress();
-		test(ts);
+		testPairing(ts);
 //break;
 	}
 	int count = (int)clk.getCount();