From 9774948cd3756eb1706b5f2cb60dee70542e477b Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Wed, 19 Apr 2023 13:44:01 +0200
Subject: [PATCH] added reference implementation and exhaustive testing

---
 galois_field.hh  | 87 ++++++++++++++++++++++++++++++++++++++++++++++++
 tests/gf_test.cc | 24 +++++++++++++
 2 files changed, 111 insertions(+)

diff --git a/galois_field.hh b/galois_field.hh
index d6e3086..229178f 100644
--- a/galois_field.hh
+++ b/galois_field.hh
@@ -318,4 +318,91 @@ private:
 	GF::Tables<M, POLY, TYPE> Tables;
 };
 
+template <int M, int64_t POLY, typename TYPE>
+struct GaloisFieldReference
+{
+	static const int64_t Q = int64_t(1) << M;
+	static const TYPE N = Q - 1;
+	static_assert(M <= 8 * sizeof(TYPE), "TYPE not wide enough");
+	static_assert(Q == (POLY & ~(Q - 1)), "POLY not of degree Q");
+	static const TYPE P = POLY;
+	TYPE v;
+	GaloisFieldReference() = default;
+	explicit GaloisFieldReference(TYPE v) : v(v)
+	{
+		assert(v <= N);
+	}
+	GaloisFieldReference<M, POLY, TYPE> operator *= (GaloisFieldReference<M, POLY, TYPE> a)
+	{
+		return *this = *this * a;
+	}
+	GaloisFieldReference<M, POLY, TYPE> operator += (GaloisFieldReference<M, POLY, TYPE> a)
+	{
+		return *this = *this + a;
+	}
+};
+
+template <int M, int64_t POLY, typename TYPE>
+bool operator == (GaloisFieldReference<M, POLY, TYPE> a, GaloisFieldReference<M, POLY, TYPE> b)
+{
+	return a.v == b.v;
+}
+
+template <int M, int64_t POLY, typename TYPE>
+bool operator != (GaloisFieldReference<M, POLY, TYPE> a, GaloisFieldReference<M, POLY, TYPE> b)
+{
+	return a.v != b.v;
+}
+
+template <int M, int64_t POLY, typename TYPE>
+GaloisFieldReference<M, POLY, TYPE> operator + (GaloisFieldReference<M, POLY, TYPE> a, GaloisFieldReference<M, POLY, TYPE> b)
+{
+	return GaloisFieldReference<M, POLY, TYPE>(a.v ^ b.v);
+}
+
+template <int M, int64_t POLY, typename TYPE>
+GaloisFieldReference<M, POLY, TYPE> operator - (GaloisFieldReference<M, POLY, TYPE> a, GaloisFieldReference<M, POLY, TYPE> b)
+{
+	return a + b;
+}
+
+template <int M, int64_t POLY, typename TYPE>
+GaloisFieldReference<M, POLY, TYPE> operator - (GaloisFieldReference<M, POLY, TYPE> a)
+{
+	return a;
+}
+
+template <int M, int64_t POLY, typename TYPE>
+GaloisFieldReference<M, POLY, TYPE> operator * (GaloisFieldReference<M, POLY, TYPE> a, GaloisFieldReference<M, POLY, TYPE> b)
+{
+	GaloisFieldReference<M, POLY, TYPE> p(0);
+	for (int i = 0; i < M; ++i) {
+		if (b.v & 1)
+			p.v ^= a.v;
+		if (a.v & (TYPE(1) << (M - 1)))
+			a.v = (a.v << 1) ^ a.P;
+		else
+			a.v <<= 1;
+		b.v >>= 1;
+	}
+	return p;
+}
+
+template <int M, int64_t POLY, typename TYPE>
+GaloisFieldReference<M, POLY, TYPE> rcp(GaloisFieldReference<M, POLY, TYPE> a)
+{
+	assert(a.v);
+	GaloisFieldReference<M, POLY, TYPE> t(a *= a);
+	for (int i = 0; i < M - 2; ++i)
+		t *= a *= a;
+	return t;
+}
+
+template <int M, int64_t POLY, typename TYPE>
+GaloisFieldReference<M, POLY, TYPE> operator / (GaloisFieldReference<M, POLY, TYPE> a, GaloisFieldReference<M, POLY, TYPE> b)
+{
+	assert(b.v);
+	return a * rcp(b);
+}
+
 }
diff --git a/tests/gf_test.cc b/tests/gf_test.cc
index 1ca0e5b..4c9c214 100644
--- a/tests/gf_test.cc
+++ b/tests/gf_test.cc
@@ -8,36 +8,60 @@ Copyright 2018 Ahmet Inan <inan@aicodix.de>
 #include <iostream>
 #include "galois_field.hh"
 
+template <typename GF, typename GFR>
+void exhaustive_test()
+{
+	for (int j = 0; j < GF::Q; ++j) {
+		typename GF::ValueType a(j);
+		for (int i = 0; i < GF::Q; ++i) {
+			typename GF::ValueType b(i);
+			assert((a * b).v == (GFR(j) * GFR(i)).v);
+		}
+	}
+	for (int i = 1; i < GF::Q; ++i) {
+		typename GF::ValueType a(i);
+		assert(rcp(a).v == rcp(GFR(i)).v);
+	}
+}
+
 int main()
 {
 	if (1) {
 		// BBC WHP031 RS(15, 11) T=2
 		typedef CODE::GaloisField<4, 0b10011, uint8_t> GF;
+		typedef CODE::GaloisFieldReference<4, 0b10011, uint8_t> GFR;
 		GF instance;
 		GF::ValueType a(3), b(7), c(15), d(6);
 		assert(a * b + c == d);
+		exhaustive_test<GF, GFR>();
 	}
 	if (1) {
 		// DVB-T RS(255, 239) T=8
 		typedef CODE::GaloisField<8, 0b100011101, uint8_t> GF;
+		typedef CODE::GaloisFieldReference<8, 0b100011101, uint8_t> GFR;
 		GF instance;
 		GF::ValueType a(154), b(83), c(144), d(63);
 		assert(fma(a, b, c) == d);
+		exhaustive_test<GF, GFR>();
 	}
 	if (1) {
 		// FUN RS(65535, 65471) T=32
 		typedef CODE::GaloisField<16, 0b10001000000001011, uint16_t> GF;
+		typedef CODE::GaloisFieldReference<16, 0b10001000000001011, uint16_t> GFR;
 		GF *instance = new GF();
 		GF::ValueType a(42145), b(13346), c(40958), d(35941);
 		assert(a / b + c == d);
+		exhaustive_test<GF, GFR>();
 		delete instance;
 	}
 	if (1) {
 		// DVB-S2 FULL BCH(65535, 65343) T=12
 		typedef CODE::GaloisField<16, 0b10000000000101101, uint16_t> GF;
+		typedef CODE::GaloisFieldReference<16, 0b10000000000101101, uint16_t> GFR;
 		GF *instance = new GF();
 		GF::ValueType a(38532), b(7932), c(34283), d(22281);
 		assert(a / (rcp(b) + c) == d);
+		exhaustive_test<GF, GFR>();
 		delete instance;
 	}
 	std::cerr << "Galois field arithmetic test passed!" << std::endl;