diff --git a/prime_field.hh b/prime_field.hh
new file mode 100644
index 0000000..88996bd
--- /dev/null
+++ b/prime_field.hh
@@ -0,0 +1,145 @@
+/*
+Prime field arithmetic
+
+Copyright 2024 Ahmet Inan <inan@aicodix.de>
+*/
+
+#pragma once
+
+#include <cassert>
+
+namespace CODE {
+
+template <typename TYPE, TYPE PRIME>
+struct PrimeField
+{
+	static_assert(std::is_unsigned<TYPE>::value, "TYPE must be unsigned");
+	static_assert(std::numeric_limits<TYPE>::max() / (PRIME-1) >= (PRIME-1), "Type not wide enough");
+	static constexpr TYPE P = PRIME;
+	TYPE v;
+	PrimeField() = default;
+	explicit PrimeField(TYPE v) : v(v)
+	{
+	}
+	PrimeField<TYPE, PRIME> operator *= (PrimeField<TYPE, PRIME> a)
+	{
+		return *this = *this * a;
+	}
+	PrimeField<TYPE, PRIME> operator /= (PrimeField<TYPE, PRIME> a)
+	{
+		return *this = *this / a;
+	}
+	PrimeField<TYPE, PRIME> operator += (PrimeField<TYPE, PRIME> a)
+	{
+		return *this = *this + a;
+	}
+	PrimeField<TYPE, PRIME> operator -= (PrimeField<TYPE, PRIME> a)
+	{
+		return *this = *this - a;
+	}
+	TYPE operator () ()
+	{
+		return v;
+	}
+};
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> reduce(PrimeField<TYPE, PRIME> a)
+{
+	return PrimeField<TYPE, PRIME>(a.v % a.P);
+}
+
+template <typename TYPE, TYPE PRIME>
+bool operator == (PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return a.v == b.v;
+}
+
+template <typename TYPE, TYPE PRIME>
+bool operator != (PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return a.v != b.v;
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> add(PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return PrimeField<TYPE, PRIME>(a.v + b.v);
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> operator + (PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return reduce(add(a, b));
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> sub(PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return PrimeField<TYPE, PRIME>(a.v - b.v + a.P);
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> operator - (PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return reduce(sub(a, b));
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> neg(PrimeField<TYPE, PRIME> a)
+{
+	return PrimeField<TYPE, PRIME>(a.P - a.v);
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> operator - (PrimeField<TYPE, PRIME> a)
+{
+	return reduce(neg(a));
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> mul(PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return PrimeField<TYPE, PRIME>(a.v * b.v);
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> operator * (PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return reduce(mul(a, b));
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> rcp(PrimeField<TYPE, PRIME> a)
+{
+	assert(a.v);
+	if (a.v == 1)
+		return a;
+	TYPE t = 0, newt = 1;
+	TYPE r = a.P, newr = a.v;
+	while (newr) {
+		TYPE quotient = r / newr;
+		t -= quotient * newt;
+		r -= quotient * newr;
+		std::swap(newt, t);
+		std::swap(newr, r);
+	}
+	assert(r == 1);
+	if (t >= a.P)
+		t += a.P;
+	return PrimeField<TYPE, PRIME>(t);
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> div(PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return mul(a, rcp(b));
+}
+
+template <typename TYPE, TYPE PRIME>
+PrimeField<TYPE, PRIME> operator / (PrimeField<TYPE, PRIME> a, PrimeField<TYPE, PRIME> b)
+{
+	return reduce(div(a, b));
+}
+
+}
diff --git a/tests/pf_test.cc b/tests/pf_test.cc
new file mode 100644
index 0000000..ee38696
--- /dev/null
+++ b/tests/pf_test.cc
@@ -0,0 +1,36 @@
+/*
+Test for the prime field arithmetic
+
+Copyright 2024 Ahmet Inan <inan@aicodix.de>
+*/
+
+#include <cstdint>
+#include <cassert>
+#include <iostream>
+#include "prime_field.hh"
+
+template <typename TYPE, TYPE PRIME>
+void exhaustive_test()
+{
+	typedef CODE::PrimeField<TYPE, PRIME> PF;
+	for (TYPE a = 0; a < PRIME; ++a)
+		for (TYPE b = 0; b < PRIME; ++b)
+			assert((PF(a) * PF(b))() == (a * b) % PRIME);
+	for (TYPE a = 1; a < PRIME; ++a)
+		assert(rcp(PF(a)) * PF(a) == PF(1));
+	for (TYPE a = 0; a < PRIME; ++a)
+		for (TYPE b = 0; b < PRIME; ++b)
+			assert((PF(a) + PF(b))() == (a + b) % PRIME);
+	for (TYPE a = 0; a < PRIME; ++a)
+		for (TYPE b = 0; b < PRIME; ++b)
+			assert((PF(a) - PF(b))() == (a - b + PRIME) % PRIME);
+}
+
+int main()
+{
+	exhaustive_test<uint32_t, 257>();
+	exhaustive_test<uint64_t, 65537>();
+	std::cerr << "Prime field arithmetic test passed!" << std::endl;
+	return 0;
+}
+