From f7ed8e13d546e2d5acea3ac82d7745112b1bbcf1 Mon Sep 17 00:00:00 2001
From: Ahmet Inan <xdsopl@gmail.com>
Date: Sun, 26 Apr 2026 12:38:25 +0200
Subject: [PATCH] added Mersenne 31 prime field

---
 mersenne_31.hh    | 130 ++++++++++++++++++++++++++++++++++++++++++++++
 tests/m31_test.cc |  56 ++++++++++++++++++++
 2 files changed, 186 insertions(+)
 create mode 100644 mersenne_31.hh
 create mode 100644 tests/m31_test.cc

diff --git a/mersenne_31.hh b/mersenne_31.hh
new file mode 100644
index 0000000..6877b5e
--- /dev/null
+++ b/mersenne_31.hh
@@ -0,0 +1,130 @@
+/*
+Mersenne 2^31-1 prime field arithmetic
+
+Copyright 2026 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+namespace CODE {
+
+struct Mersenne31;
+
+Mersenne31 operator + (Mersenne31, Mersenne31);
+Mersenne31 operator - (Mersenne31, Mersenne31);
+Mersenne31 operator * (Mersenne31, Mersenne31);
+Mersenne31 operator / (Mersenne31, Mersenne31);
+
+struct Mersenne31
+{
+	static constexpr uint32_t P = 0x7FFFFFFF;
+	uint32_t v;
+	Mersenne31() = default;
+	explicit Mersenne31(uint32_t v) : v(v)
+	{
+	}
+	Mersenne31 operator *= (Mersenne31 a)
+	{
+		return *this = *this * a;
+	}
+	Mersenne31 operator /= (Mersenne31 a)
+	{
+		return *this = *this / a;
+	}
+	Mersenne31 operator += (Mersenne31 a)
+	{
+		return *this = *this + a;
+	}
+	Mersenne31 operator -= (Mersenne31 a)
+	{
+		return *this = *this - a;
+	}
+	uint32_t operator () ()
+	{
+		return v == P ? 0 : v;
+	}
+};
+
+Mersenne31 reduce(Mersenne31 a)
+{
+	return Mersenne31((a.v & 0x7FFFFFFF) + (a.v >> 31));
+}
+
+bool operator == (Mersenne31 a, Mersenne31 b)
+{
+	return a.v == b.v;
+}
+
+bool operator != (Mersenne31 a, Mersenne31 b)
+{
+	return a.v != b.v;
+}
+
+Mersenne31 add(Mersenne31 a, Mersenne31 b)
+{
+	return Mersenne31(a.v + b.v);
+}
+
+Mersenne31 operator + (Mersenne31 a, Mersenne31 b)
+{
+	return reduce(add(a, b));
+}
+
+Mersenne31 neg(Mersenne31 a)
+{
+	return Mersenne31(a.P - a.v);
+}
+
+Mersenne31 operator - (Mersenne31 a)
+{
+	return reduce(neg(a));
+}
+
+Mersenne31 sub(Mersenne31 a, Mersenne31 b)
+{
+	return add(a, -b);
+}
+
+Mersenne31 operator - (Mersenne31 a, Mersenne31 b)
+{
+	return reduce(sub(a, b));
+}
+
+Mersenne31 mul(Mersenne31 a, Mersenne31 b)
+{
+	uint64_t v = uint64_t(a.v) * uint64_t(b.v);
+	uint32_t l = v & 0x7FFFFFFF;
+	uint32_t h = v >> 31;
+	return Mersenne31(l + h);
+}
+
+Mersenne31 operator * (Mersenne31 a, Mersenne31 b)
+{
+	return reduce(mul(a, b));
+}
+
+Mersenne31 pow(Mersenne31 a, uint32_t m)
+{
+	Mersenne31 t(1);
+	for (;m; m >>= 1, a *= a)
+		if (m & 1)
+			t *= a;
+	return t;
+}
+
+Mersenne31 rcp(Mersenne31 a)
+{
+	return pow(a, a.P - 2);
+}
+
+Mersenne31 div(Mersenne31 a, Mersenne31 b)
+{
+	return mul(a, rcp(b));
+}
+
+Mersenne31 operator / (Mersenne31 a, Mersenne31 b)
+{
+	return reduce(div(a, b));
+}
+
+}
diff --git a/tests/m31_test.cc b/tests/m31_test.cc
new file mode 100644
index 0000000..b46c586
--- /dev/null
+++ b/tests/m31_test.cc
@@ -0,0 +1,56 @@
+/*
+Test for the Mersenne 2^31-1 prime field arithmetic
+
+Copyright 2026 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <random>
+#include <chrono>
+#include <cstdint>
+#include <cassert>
+#include <iostream>
+#include <functional>
+#include "mersenne_31.hh"
+
+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;
+	auto rand0 = std::bind(distribution(0, M31::P-1), generator(rd()));
+	auto rand1 = std::bind(distribution(1, M31::P-1), generator(rd()));
+	for (uint32_t i = 0; i < count; ++i) {
+		uint32_t a = rand0();
+		for (uint32_t j = 0; j < count; ++j) {
+			uint32_t b = rand0();
+			assert((M31(a) * M31(b))() == uint32_t((uint64_t(a) * b) % M31::P));
+		}
+	}
+	for (uint32_t i = 0; i < count; ++i) {
+		uint32_t a = rand1();
+		assert(rcp(M31(a)) * M31(a) == M31(1));
+	}
+	for (uint32_t i = 0; i < count; ++i) {
+		uint32_t a = rand0();
+		for (uint32_t j = 0; j < count; ++j) {
+			uint32_t b = rand0();
+			assert((M31(a) + M31(b))() == (a + b) % M31::P);
+		}
+	}
+	for (uint32_t i = 0; i < count; ++i) {
+		uint32_t a = rand0();
+		for (uint32_t j = 0; j < count; ++j) {
+			uint32_t b = rand0();
+			assert((M31(a) - M31(b))() == (a - b + M31::P) % M31::P);
+		}
+	}
+}
+
+int main()
+{
+	random_test(10000);
+	std::cerr << "Mersenne 2^31-1 prime field arithmetic test passed!" << std::endl;
+	return 0;
+}
+