added Cauchy Mersenne 31 prime field erasure encoding

cauchy_mersenne_erasure_coding.hh

@@ -0,0 +1,89 @@
+/*
+Cauchy Mersenne 2^31-1 Prime Field Erasure Coding
+
+Copyright 2026 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include "mersenne_31.hh"
+
+namespace CODE {
+
+struct CauchyMersenneErasureCoding
+{
+	typedef Mersenne31 M31;
+	M31 row_num, row_den;
+	// $a_{ij} = \frac{1}{x_i + y_j}$
+	M31 cauchy_matrix(int i, int j)
+	{
+		M31 row(i), col(j);
+		return rcp(row + col);
+	}
+	// $b_{ij} = \frac{\prod_{k=1}^{n}{(x_j + y_k)(x_k + y_i)}}{(x_j + y_i)\prod_{k \ne j}^{n}{(x_j - x_k)}\prod_{k \ne i}^{n}{(y_i - y_k)}}$
+	M31 inverse_cauchy_matrix(const int *rows, int i, int j, int n)
+	{
+#if 0
+		M31 row_j(rows[j]), col_i(i);
+		M31 prod_xy(1), prod_x(1), prod_y(1);
+		for (int k = 0; k < n; k++) {
+			M31 row_k(rows[k]), col_k(k);
+			prod_xy *= (row_j + col_k) * (row_k + col_i);
+			if (k != j)
+				prod_x *= (row_j - row_k);
+			if (k != i)
+				prod_y *= (col_i - col_k);
+		}
+		return prod_xy / ((row_j + col_i) * prod_x * prod_y);
+#else
+		M31 row_j(rows[j]), col_i(i);
+		if (j == 0) {
+			M31 num(1), den(1);
+			for (int k = 0; k < n; k++) {
+				M31 row_k(rows[k]), col_k(k);
+				num *= row_k + col_i;
+				if (k != i)
+					den *= col_i - col_k;
+			}
+			row_num = num;
+			row_den = den;
+		}
+		M31 num(row_num), den(row_den);
+		for (int k = 0; k < n; k++) {
+			M31 row_k(rows[k]), col_k(k);
+			num *= row_j + col_k;
+			if (k != j)
+				den *= row_j - row_k;
+		}
+		return num / ((row_j + col_i) * den);
+#endif
+	}
+	void mac(M31 *c, const M31 *a, M31 b, int len, bool init)
+	{
+		if (init) {
+			for (int i = 0; i < len; i++)
+				c[i] = b * a[i];
+		} else {
+			for (int i = 0; i < len; i++)
+				c[i] += b * a[i];
+		}
+	}
+	void encode(const M31 *data, M31 *block, int block_id, int block_len, int block_cnt)
+	{
+		assert(block_id >= block_cnt && block_id < int(M31::P) / 2);
+		for (int k = 0; k < block_cnt; k++) {
+			M31 a_ik = cauchy_matrix(block_id, k);
+			mac(block, data + block_len * k, a_ik, block_len, !k);
+		}
+	}
+	void decode(M31 *data, const M31 *blocks, const int *block_ids, int block_idx, int block_len, int block_cnt)
+	{
+		for (int k = 0; k < block_cnt; k++) {
+			M31 b_ik = inverse_cauchy_matrix(block_ids, block_idx, k, block_cnt);
+			mac(data, blocks + block_len * k, b_ik, block_len, !k);
+		}
+	}
+};
+
+}
+

cauchy_prime_field_erasure_coding.hh

@@ -15,7 +15,7 @@ struct CauchyPrimeFieldErasureCoding
 	PF temp[MAX_LEN];
 	typedef unsigned used_word;
 	static constexpr int used_width = 8 * sizeof(used_word);
-	static constexpr int used_length = (PF::P + used_width - 1) / used_width;
+	static constexpr int used_length = (MAX_LEN + used_width - 1) / used_width;
 	used_word used_values[used_length];
 	PF row_num, row_den;
 	// $a_{ij} = \frac{1}{x_i + y_j}$
@@ -112,9 +112,10 @@ struct CauchyPrimeFieldErasureCoding
 		for (int i = 0; i < used_length; ++i)
 			used_values[i] = 0;
 		for (int i = 0; i < block_len; ++i)
-			used_values[temp[i]()/used_width] |= 1 << temp[i]()%used_width;
+			if (temp[i]()/used_width < used_length)
+				used_values[temp[i]()/used_width] |= 1 << temp[i]()%used_width;
 		int s = 0;
-		while (used_values[s/used_width] & 1 << s%used_width)
+		while (s/used_width < used_length && used_values[s/used_width] & 1 << s%used_width)
 			++s;
 		return s;
 	}

mersenne_31.hh

@@ -0,0 +1,134 @@
+/*
+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() == b();
+}
+
+bool operator != (Mersenne31 a, Mersenne31 b)
+{
+	return a() != b();
+}
+
+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 neg(a);
+}
+
+Mersenne31 sub(Mersenne31 a, Mersenne31 b)
+{
+	return add(a, neg(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)
+{
+	Mersenne31 t = a;
+	a *= a;
+	for (int i = 0; i < 29; ++i)
+		t *= a *= a;
+	return t;
+}
+
+Mersenne31 div(Mersenne31 a, Mersenne31 b)
+{
+	return mul(a, rcp(b));
+}
+
+Mersenne31 operator / (Mersenne31 a, Mersenne31 b)
+{
+	return reduce(div(a, b));
+}
+
+}

pac_encoder.hh

@@ -22,7 +22,7 @@ class PACEncoder
 		bool b4 = (*state >> 4) & 1;
 		bool b6 = (*state >> 6) & 1;
 		bool output = input ^ b1 ^ b3 ^ b4 ^ b6;
-		*state = ((*state & 126) << 1) | (input ? 2 : 0) | (output ? 1 : 0);
+		*state = ((*state & 62) << 1) | (input ? 2 : 0) | (output ? 1 : 0);
 		return output;
 	}
 public:

pac_list_decoder.hh

@@ -46,7 +46,7 @@ struct PACListTree<TYPE, 1>
 		bool b4 = (*state >> 4) & 1;
 		bool b6 = (*state >> 6) & 1;
 		bool output = input ^ b1 ^ b3 ^ b4 ^ b6;
-		*state = ((*state & 126) << 1) | (input ? 2 : 0) | (output ? 1 : 0);
+		*state = ((*state & 62) << 1) | (input ? 2 : 0) | (output ? 1 : 0);
 		return output;
 	}
 	static MAP rate0(PATH *metric, int *state, TYPE *hard, TYPE *soft)

tests/cme_regression_test.cc

@@ -0,0 +1,77 @@
+/*
+Regression Test for the Cauchy Mersenne 2^31-1 Prime Field Encoder and Decoder
+
+Copyright 2026 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <cstdlib>
+#include <cassert>
+#include <chrono>
+#include <random>
+#include <iostream>
+#include <functional>
+#include "mersenne_31.hh"
+#include "cauchy_mersenne_erasure_coding.hh"
+
+void cme_test(int trials)
+{
+	typedef CODE::Mersenne31 M31;
+	int value_bytes = sizeof(M31);
+	//int value_bits = value_bytes * 8;
+	const int MAX_LEN = 1024;
+	CODE::CauchyMersenneErasureCoding cme;
+	std::random_device rd;
+	std::default_random_engine generator(rd());
+	typedef std::uniform_int_distribution<int> distribution;
+	auto rnd_cnt = std::bind(distribution(1, 256), generator);
+	auto rnd_len = std::bind(distribution(1, MAX_LEN), generator);
+	auto rnd_dat = std::bind(distribution(0, M31::P-1), generator);
+	//auto rnd_dat = std::bind(distribution(0, (1 << value_bits) - 1), generator);
+	while (--trials) {
+		int block_count = rnd_cnt();
+		int idents_total = 1000000; // M31::P / 2 - block_count;
+		int block_values = rnd_len();
+		int block_bytes = block_values * value_bytes;
+		int data_values = block_count * block_values;
+		int data_bytes = data_values * value_bytes;
+		M31 *orig = new M31[data_values];
+		M31 *data = new M31[data_values];
+		M31 *blocks = new M31[data_values];
+		int *idents = new int[idents_total];
+		for (int i = 0; i < data_values; ++i)
+			orig[i] = M31(rnd_dat());
+		for (int i = 0; i < idents_total; ++i)
+			idents[i] = block_count + i;
+		for (int i = 0; i < block_count; i++) {
+			std::uniform_int_distribution<int> hat(i, idents_total - 1);
+			std::swap(idents[i], idents[hat(generator)]);
+		}
+		auto enc_start = std::chrono::system_clock::now();
+		for (int i = 0; i < block_count; ++i)
+			cme.encode(orig, blocks + block_values * i, idents[i], block_values, block_count);
+		auto enc_end = std::chrono::system_clock::now();
+		auto enc_usec = std::chrono::duration_cast<std::chrono::microseconds>(enc_end - enc_start);
+		double enc_mbs = double(data_bytes) / enc_usec.count();
+		auto dec_start = std::chrono::system_clock::now();
+		for (int i = 0; i < block_count; ++i)
+			cme.decode(data + block_values * i, blocks, idents, i, block_values, block_count);
+		auto dec_end = std::chrono::system_clock::now();
+		auto dec_usec = std::chrono::duration_cast<std::chrono::microseconds>(dec_end - dec_start);
+		double dec_mbs = double(data_bytes) / dec_usec.count();
+		std::cout << "block count = " << block_count << ", block size = " << block_bytes << " bytes, encoding speed = " << enc_mbs << " megabyte per second, decoding speed = " << dec_mbs << " megabyte per second" << std::endl;
+		for (int i = 0; i < data_values; ++i)
+			assert(data[i] == orig[i]);
+		delete[] idents;
+		delete[] blocks;
+		delete[] orig;
+		delete[] data;
+	}
+}
+
+int main()
+{
+	cme_test(100);
+	std::cerr << "Cauchy Mersenne 2^31-1 prime field regression test passed!" << std::endl;
+	return 0;
+}
+

tests/m31_test.cc

@@ -0,0 +1,64 @@
+/*
+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"
+
+typedef CODE::Mersenne31 M31;
+
+void random_test(uint32_t count)
+{
+	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 * 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);
+		}
+	}
+}
+
+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;
+}
+