From 0ab371bbe9b0efea55a8a5fb47cfbb74b9ced525 Mon Sep 17 00:00:00 2001
From: MITSUNARI Shigeo <herumi@nifty.com>
Date: Wed, 13 Apr 2016 17:59:40 +0900
Subject: [PATCH] add Fp6::inv

---
 include/mcl/fp_tower.hpp | 42 ++++++++++++++++++++++++++++++++++++++++
 test/fp_tower_test.cpp   | 11 +++++++++++
 2 files changed, 53 insertions(+)

diff --git a/include/mcl/fp_tower.hpp b/include/mcl/fp_tower.hpp
index 8d0c887..0ecf7ad 100644
--- a/include/mcl/fp_tower.hpp
+++ b/include/mcl/fp_tower.hpp
@@ -528,6 +528,48 @@ struct Fp6T {
 		t1 -= cf;
 		Fp2::add(z.c, t1, be);
 	}
+	/*
+		x = a + bv + cv^2, v^3 = xi
+		y = 1/x = p/q where
+		p = (a^2 - bc xi) + (c^2 xi - ab)v + (b^2 - ac)v^2
+		q = c^3 xi^2 + b(b^2 - 3ac)xi + a^3
+	*/
+	static inline void inv(Fp6T& y, const Fp6T& x)
+	{
+		const Fp2& a = x.a;
+		const Fp2& b = x.b;
+		const Fp2& c = x.c;
+		Fp2 aa, bb, cc, ab, bc, ac;
+		Fp2::sqr(aa, a);
+		Fp2::sqr(bb, b);
+		Fp2::sqr(cc, c);
+		Fp2::mul(ab, a, b);
+		Fp2::mul(bc, b, c);
+		Fp2::mul(ac, c, a);
+		Fp2 t;
+		Fp2::mul_xi(t, bc);
+		Fp6T p;
+		Fp2::sub(p.a, aa, t);
+		Fp2::mul_xi(t, cc);
+		Fp2::sub(p.b, t, ab);
+		Fp2::sub(p.c, bb, ac);
+		Fp2 q;
+		Fp2::mul(q, aa, a); // a^3
+		Fp2::add(t, ac, ac); // t = 2ac
+		t += ac; // t = 3ac
+		Fp2::sub(t, bb, t); // t = b^2 - 3ac
+		t *= b;
+		Fp2::mul_xi(t, t); // b(b^2 - 3ac)xi
+		q += t;
+		Fp2::mul(t, cc, c);
+		Fp2::mul_xi(t, t);
+		Fp2::mul_xi(t, t); // QQQ : c^3 xi^2
+		q += t;
+		Fp2::inv(q, q);
+		Fp2::mul(y.a, p.a, q);
+		Fp2::mul(y.b, p.b, q);
+		Fp2::mul(y.c, p.c, q);
+	}
 };
 
 /*
diff --git a/test/fp_tower_test.cpp b/test/fp_tower_test.cpp
index 9901b88..6d2c062 100644
--- a/test/fp_tower_test.cpp
+++ b/test/fp_tower_test.cpp
@@ -117,6 +117,17 @@ void testFp6()
 	Fp6::mul(w, x, x);
 	testFp6sqr(a, b, c, z);
 	testFp6sqr(a, b, c, w);
+	for (int i = 0; i < 10; i++) {
+		Fp6::inv(y, x);
+		Fp6::mul(z, y, x);
+		CYBOZU_TEST_EQUAL(z, 1);
+		x += y;
+		y = x;
+		Fp6::inv(y, x);
+		Fp6::mul(z, y, x);
+		y *= x;
+		CYBOZU_TEST_EQUAL(z, 1);
+	}
 }
 
 void testFp12()