From 01438be6e165c1a38727721b8eb134220de25b4e Mon Sep 17 00:00:00 2001
From: Ahmet Inan <xdsopl@gmail.com>
Date: Mon, 27 Apr 2026 04:43:35 +0200
Subject: [PATCH] unroll the pow loop in rcp

---
 mersenne_31.hh    |  6 +++++-
 tests/m31_test.cc | 12 ++++++++++--
 2 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/mersenne_31.hh b/mersenne_31.hh
index 8381e85..7f0d323 100644
--- a/mersenne_31.hh
+++ b/mersenne_31.hh
@@ -114,7 +114,11 @@ Mersenne31 pow(Mersenne31 a, uint32_t m)
 
 Mersenne31 rcp(Mersenne31 a)
 {
-	return pow(a, a.P - 2);
+	Mersenne31 t = a;
+	a *= a;
+	for (int i = 0; i < 29; ++i)
+		t *= a *= a;
+	return t;
 }
 
 Mersenne31 div(Mersenne31 a, Mersenne31 b)
diff --git a/tests/m31_test.cc b/tests/m31_test.cc
index b46c586..865a651 100644
--- a/tests/m31_test.cc
+++ b/tests/m31_test.cc
@@ -12,9 +12,10 @@ Copyright 2026 Ahmet Inan <inan@aicodix.de>
 #include <functional>
 #include "mersenne_31.hh"
 
+typedef CODE::Mersenne31 M31;
+
 void random_test(uint32_t count)
 {
-	typedef CODE::Mersenne31 M31;
 	std::random_device rd;
 	typedef std::default_random_engine generator;
 	typedef std::uniform_int_distribution<int> distribution;
@@ -27,7 +28,7 @@ void random_test(uint32_t count)
 			assert((M31(a) * M31(b))() == uint32_t((uint64_t(a) * b) % M31::P));
 		}
 	}
-	for (uint32_t i = 0; i < count; ++i) {
+	for (uint32_t i = 0; i < count * count; ++i) {
 		uint32_t a = rand1();
 		assert(rcp(M31(a)) * M31(a) == M31(1));
 	}
@@ -47,8 +48,15 @@ void random_test(uint32_t count)
 	}
 }
 
+void exhaustive_test()
+{
+	for (uint32_t a = 1; a < M31::P; ++a)
+		assert(rcp(M31(a)) * M31(a) == M31(1));
+}
+
 int main()
 {
+	//exhaustive_test();
 	random_test(10000);
 	std::cerr << "Mersenne 2^31-1 prime field arithmetic test passed!" << std::endl;
 	return 0;