From 9570346a59fc16fb0251e752ba08c04d22b282c3 Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Fri, 21 Apr 2023 08:17:11 +0200
Subject: [PATCH] experimenting with 3 times faster EEA

https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Simple_algebraic_field_extensions
---
 galois_field.hh | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/galois_field.hh b/galois_field.hh
index 229178f..8ac3713 100644
--- a/galois_field.hh
+++ b/galois_field.hh
@@ -392,10 +392,32 @@ template <int M, int64_t POLY, typename TYPE>
 GaloisFieldReference<M, POLY, TYPE> rcp(GaloisFieldReference<M, POLY, TYPE> a)
 {
 	assert(a.v);
+#if 1
+	uint32_t r = POLY, newr = a.v;
+	uint32_t t = 0, newt = 1;
+	auto degree = [](uint32_t a) {
+		int d = 0;
+		while (a >>= 1)
+			++d;
+		return d;
+	};
+	while (newr != 1) {
+		int j = degree(newr) - degree(r);
+		if (j < 0) {
+			j = -j;
+			std::swap(newr, r);
+			std::swap(newt, t);
+		}
+		newr ^= r << j;
+		newt ^= t << j;
+	}
+	return GaloisFieldReference<M, POLY, TYPE>(newt);
+#else
 	GaloisFieldReference<M, POLY, TYPE> t(a *= a);
 	for (int i = 0; i < M - 2; ++i)
 		t *= a *= a;
 	return t;
+#endif
 }
 
 template <int M, int64_t POLY, typename TYPE>