From 2a7fa2b7b7dbdd41868afa3c99459a84126754d5 Mon Sep 17 00:00:00 2001 From: Ahmet Inan Date: Thu, 20 Sep 2018 14:06:57 +0200 Subject: [PATCH] added Galois field arithmetic --- README.md | 4 + galois_field.hh | 302 +++++++++++++++++++++++++++++++++++++++++++++++ tests/gf_test.cc | 46 ++++++++ 3 files changed, 352 insertions(+) create mode 100644 galois_field.hh create mode 100644 tests/gf_test.cc diff --git a/README.md b/README.md index 9c62fdc..1ae163c 100644 --- a/README.md +++ b/README.md @@ -31,3 +31,7 @@ Here a [Pseudorandom number generator](https://en.wikipedia.org/wiki/Pseudorando When dealing with unaligned and arbitrary-bit-sized elements in a data stream, the bitwise stream container might help avoiding some headaches. +### [galois_field.hh](galois_field.hh) + +We have to thank Évariste Galois' contribution of the [Finite field](https://en.wikipedia.org/wiki/Finite_field) to mathematics, which laid the cornerstone for a variety of applications that we take for granted today. + diff --git a/galois_field.hh b/galois_field.hh new file mode 100644 index 0000000..81fbadf --- /dev/null +++ b/galois_field.hh @@ -0,0 +1,302 @@ +/* +Galois field arithmetic + +Copyright 2018 Ahmet Inan +*/ + +#ifndef GALOIS_FIELD_HH +#define GALOIS_FIELD_HH + +#include + +namespace CODE { +namespace GF { + +template +struct Index; + +template +struct Value +{ + static const int Q = 1 << M, N = Q - 1; + static_assert(M <= 8 * sizeof(TYPE), "TYPE not wide enough"); + static_assert(Q == (POLY & ~N), "POLY not of degree Q"); + TYPE v; + Value() {} + explicit Value(TYPE v) : v(v) + { + assert(v <= N); + } + explicit operator bool () const { return v; } + explicit operator int () const { return v; } + Value operator *= (Index a) + { + assert(a.i < a.modulus()); + return *this = *this * a; + } + Value operator *= (Value a) + { + assert(a.v <= a.N); + return *this = *this * a; + } + Value operator += (Value a) + { + assert(a.v <= a.N); + return *this = *this + a; + } + static const Value zero() + { + return Value(0); + } +}; + +template +struct Index +{ + static const int Q = 1 << M, N = Q - 1; + static_assert(M <= 8 * sizeof(TYPE), "TYPE not wide enough"); + static_assert(Q == (POLY & ~N), "POLY not of degree Q"); + TYPE i; + Index() {} + explicit Index(TYPE i) : i(i) + { + assert(i < modulus()); + } + explicit operator int () const { return i; } + Index operator *= (Index a) + { + assert(a.i < a.modulus()); + assert(i < modulus()); + return *this = *this * a; + } + Index operator /= (Index a) + { + assert(a.i < a.modulus()); + assert(i < modulus()); + return *this = *this / a; + } + static const TYPE modulus() + { + return N; + } +}; + +template +struct Tables +{ + static const int Q = 1 << M, N = Q - 1; + static_assert(M <= 8 * sizeof(TYPE), "TYPE not wide enough"); + static_assert(Q == (POLY & ~N), "POLY not of degree Q"); + static TYPE *LOG, *EXP; + TYPE log_[Q], exp_[Q]; + static TYPE next(TYPE a) + { + return a & (TYPE)(Q >> 1) ? (a << 1) ^ (TYPE)POLY : a << 1; + } + static TYPE log(TYPE a) + { + assert(LOG != nullptr); + assert(a <= N); + return LOG[a]; + } + static TYPE exp(TYPE a) + { + assert(EXP != nullptr); + assert(a <= N); + return EXP[a]; + } + Tables() + { + assert(LOG == nullptr); + LOG = log_; + assert(EXP == nullptr); + EXP = exp_; + log_[exp_[N] = 0] = N; + TYPE a = 1; + for (int i = 0; i < N; ++i, a = next(a)) + log_[exp_[i] = a] = i; + assert(1 == a); + } + ~Tables() + { + assert(LOG != nullptr); + LOG = nullptr; + assert(EXP != nullptr); + EXP = nullptr; + } +}; + +template +TYPE *Tables::LOG; + +template +TYPE *Tables::EXP; + +template +Index index(Value a) +{ + assert(a.v <= a.N); + assert(a.v); + return Index(Tables::log(a.v)); +} + +template +Value value(Index a) { + assert(a.i < a.modulus()); + return Value(Tables::exp(a.i)); +} + +template +bool operator == (Value a, Value b) +{ + assert(a.v <= a.N); + assert(b.v <= b.N); + return a.v == b.v; +} + +template +bool operator != (Value a, Value b) +{ + assert(a.v <= a.N); + assert(b.v <= b.N); + return a.v != b.v; +} + +template +Value operator + (Value a, Value b) +{ + assert(a.v <= a.N); + assert(b.v <= b.N); + return Value(a.v ^ b.v); +} + +template +Index operator * (Index a, Index b) +{ + assert(a.i < a.modulus()); + assert(b.i < b.modulus()); + TYPE tmp = a.i + b.i; + return Index(a.modulus() - a.i <= b.i ? tmp - a.modulus() : tmp); +} + +template +Value operator * (Value a, Value b) +{ + assert(a.v <= a.N); + assert(b.v <= b.N); + return (!a.v || !b.v) ? a.zero() : value(index(a) * index(b)); +} + +template +Value rcp(Value a) +{ + assert(a.v <= a.N); + assert(a.v); + return value(Index(index(a).modulus() - index(a).i)); +} + +template +Index operator / (Index a, Index b) +{ + assert(a.i < a.modulus()); + assert(b.i < b.modulus()); + TYPE tmp = a.i - b.i; + return Index(a.i < b.i ? tmp + a.modulus() : tmp); +} + +template +Value operator / (Value a, Value b) +{ + assert(a.v <= a.N); + assert(b.v <= b.N); + assert(b.v); + return !a.v ? a.zero() : value(index(a) / index(b)); +} + +template +Value operator / (Index a, Value b) +{ + assert(a.i < a.modulus()); + assert(b.v <= b.N); + assert(b.v); + return value(a / index(b)); +} + +template +Value operator / (Value a, Index b) +{ + assert(a.v <= a.N); + assert(b.i < b.modulus()); + return !a.v ? a.zero() : value(index(a) / b); +} + +template +Value operator * (Index a, Value b) +{ + assert(a.i < a.modulus()); + assert(b.v <= b.N); + return !b.v ? b.zero() : value(a * index(b)); +} + +template +Value operator * (Value a, Index b) +{ + assert(a.v <= a.N); + assert(b.i < b.modulus()); + return !a.v ? a.zero() : value(index(a) * b); +} + +template +Value fma(Index a, Index b, Value c) +{ + assert(a.i < a.modulus()); + assert(b.i < b.modulus()); + assert(c.v <= c.N); + return value(a * b) + c; +} + +template +Value fma(Index a, Value b, Value c) +{ + assert(a.i < a.modulus()); + assert(b.v <= b.N); + assert(c.v <= c.N); + return !b.v ? c : (value(a * index(b)) + c); +} + +template +Value fma(Value a, Index b, Value c) +{ + assert(a.v <= a.N); + assert(b.i < b.modulus()); + assert(c.v <= c.N); + return !a.v ? c : (value(index(a) * b) + c); +} + +template +Value fma(Value a, Value b, Value c) +{ + assert(a.v <= a.N); + assert(b.v <= b.N); + assert(c.v <= c.N); + return (!a.v || !b.v) ? c : (value(index(a) * index(b)) + c); +} + +} + +template +class GaloisField +{ +public: + static const int M = WIDTH, Q = 1 << M, N = Q - 1; + static_assert(M <= 8 * sizeof(TYPE), "TYPE not wide enough"); + static_assert(Q == (POLY & ~N), "POLY not of degree Q"); + typedef TYPE value_type; + typedef GF::Value ValueType; + typedef GF::Index IndexType; +private: + GF::Tables Tables; +}; + +} +#endif diff --git a/tests/gf_test.cc b/tests/gf_test.cc new file mode 100644 index 0000000..1ca0e5b --- /dev/null +++ b/tests/gf_test.cc @@ -0,0 +1,46 @@ +/* +Test for the Galois field arithmetic + +Copyright 2018 Ahmet Inan +*/ + +#include +#include +#include "galois_field.hh" + +int main() +{ + if (1) { + // BBC WHP031 RS(15, 11) T=2 + typedef CODE::GaloisField<4, 0b10011, uint8_t> GF; + GF instance; + GF::ValueType a(3), b(7), c(15), d(6); + assert(a * b + c == d); + } + if (1) { + // DVB-T RS(255, 239) T=8 + typedef CODE::GaloisField<8, 0b100011101, uint8_t> GF; + GF instance; + GF::ValueType a(154), b(83), c(144), d(63); + assert(fma(a, b, c) == d); + } + if (1) { + // FUN RS(65535, 65471) T=32 + typedef CODE::GaloisField<16, 0b10001000000001011, uint16_t> GF; + GF *instance = new GF(); + GF::ValueType a(42145), b(13346), c(40958), d(35941); + assert(a / b + c == d); + delete instance; + } + if (1) { + // DVB-S2 FULL BCH(65535, 65343) T=12 + typedef CODE::GaloisField<16, 0b10000000000101101, uint16_t> GF; + GF *instance = new GF(); + GF::ValueType a(38532), b(7932), c(34283), d(22281); + assert(a / (rcp(b) + c) == d); + delete instance; + } + std::cerr << "Galois field arithmetic test passed!" << std::endl; + return 0; +} +