diff --git a/field/eisenstein/doc.go b/field/eisenstein/doc.go
new file mode 100644
index 000000000..3d353e2ab
--- /dev/null
+++ b/field/eisenstein/doc.go
@@ -0,0 +1,7 @@
+// Package Eisenstein provides Eisenstein integer arithmetic.
+//
+// The Eisenstein integers form a commutative ring of algebraic integers in the
+// algebraic number field Q(ω) – the third cyclotomic field.  These are of the
+// form z = a + bω, where a and b are integers and ω is a primitive third root
+// of unity i.e. ω²+ω+1 = 0.
+package eisenstein
diff --git a/field/eisenstein/eisenstein.go b/field/eisenstein/eisenstein.go
new file mode 100644
index 000000000..e80915413
--- /dev/null
+++ b/field/eisenstein/eisenstein.go
@@ -0,0 +1,174 @@
+package eisenstein
+
+import (
+	"math/big"
+)
+
+// A ComplexNumber represents an arbitrary-precision Eisenstein integer.
+type ComplexNumber struct {
+	A0, A1 *big.Int
+}
+
+func (z *ComplexNumber) init() {
+	if z.A0 == nil {
+		z.A0 = new(big.Int)
+
+	}
+	if z.A1 == nil {
+		z.A1 = new(big.Int)
+
+	}
+}
+
+// String implements Stringer interface for fancy printing
+func (z *ComplexNumber) String() string {
+	return z.A0.String() + "+(" + z.A1.String() + "*ω)"
+}
+
+// Equal returns true if z equals x, false otherwise
+func (z *ComplexNumber) Equal(x *ComplexNumber) bool {
+	return z.A0.Cmp(x.A0) == 0 && z.A1.Cmp(x.A1) == 0
+}
+
+// Set sets z to x, and returns z.
+func (z *ComplexNumber) Set(x *ComplexNumber) *ComplexNumber {
+	z.init()
+	z.A0.Set(x.A0)
+	z.A1.Set(x.A1)
+	return z
+}
+
+// SetZero sets z to 0, and returns z.
+func (z *ComplexNumber) SetZero() *ComplexNumber {
+	z.A0 = big.NewInt(0)
+	z.A1 = big.NewInt(0)
+	return z
+}
+
+// SetOne sets z to 1, and returns z.
+func (z *ComplexNumber) SetOne() *ComplexNumber {
+	z.A0 = big.NewInt(1)
+	z.A1 = big.NewInt(0)
+	return z
+}
+
+// Neg sets z to the negative of x, and returns z.
+func (z *ComplexNumber) Neg(x *ComplexNumber) *ComplexNumber {
+	z.init()
+	z.A0.Neg(x.A0)
+	z.A1.Neg(x.A1)
+	return z
+}
+
+// Conjugate sets z to the conjugate of x, and returns z.
+func (z *ComplexNumber) Conjugate(x *ComplexNumber) *ComplexNumber {
+	z.init()
+	z.A0.Sub(x.A0, x.A1)
+	z.A1.Neg(x.A1)
+	return z
+}
+
+// Add sets z to the sum of x and y, and returns z.
+func (z *ComplexNumber) Add(x, y *ComplexNumber) *ComplexNumber {
+	z.init()
+	z.A0.Add(x.A0, y.A0)
+	z.A1.Add(x.A1, y.A1)
+	return z
+}
+
+// Sub sets z to the difference of x and y, and returns z.
+func (z *ComplexNumber) Sub(x, y *ComplexNumber) *ComplexNumber {
+	z.init()
+	z.A0.Sub(x.A0, y.A0)
+	z.A1.Sub(x.A1, y.A1)
+	return z
+}
+
+// Mul sets z to the product of x and y, and returns z.
+//
+// Given that ω²+ω+1=0, the explicit formula is:
+//
+//	(x0+x1ω)(y0+y1ω) = (x0y0-x1y1) + (x0y1+x1y0-x1y1)ω
+func (z *ComplexNumber) Mul(x, y *ComplexNumber) *ComplexNumber {
+	z.init()
+	var t [3]big.Int
+	var z0, z1 big.Int
+	t[0].Mul(x.A0, y.A0)
+	t[1].Mul(x.A1, y.A1)
+	z0.Sub(&t[0], &t[1])
+	t[0].Mul(x.A0, y.A1)
+	t[2].Mul(x.A1, y.A0)
+	t[0].Add(&t[0], &t[2])
+	z1.Sub(&t[0], &t[1])
+	z.A0.Set(&z0)
+	z.A1.Set(&z1)
+	return z
+}
+
+// Norm returns the norm of z.
+//
+// The explicit formula is:
+//
+//	N(x0+x1ω) = x0² + x1² - x0*x1
+func (z *ComplexNumber) Norm() *big.Int {
+	norm := new(big.Int)
+	temp := new(big.Int)
+	norm.Add(
+		norm.Mul(z.A0, z.A0),
+		temp.Mul(z.A1, z.A1),
+	)
+	norm.Sub(
+		norm,
+		temp.Mul(z.A0, z.A1),
+	)
+	return norm
+}
+
+// QuoRem sets z to the quotient of x and y, r to the remainder, and returns z and r.
+func (z *ComplexNumber) QuoRem(x, y, r *ComplexNumber) (*ComplexNumber, *ComplexNumber) {
+	norm := y.Norm()
+	if norm.Cmp(big.NewInt(0)) == 0 {
+		panic("division by zero")
+	}
+	z.Conjugate(y)
+	z.Mul(x, z)
+	z.A0.Div(z.A0, norm)
+	z.A1.Div(z.A1, norm)
+	r.Mul(y, z)
+	r.Sub(x, r)
+
+	return z, r
+}
+
+// HalfGCD returns the rational reconstruction of a, b.
+// This outputs w, v, u s.t. w = a*u + b*v.
+func HalfGCD(a, b *ComplexNumber) [3]*ComplexNumber {
+
+	var aRun, bRun, u, v, u_, v_, quotient, remainder, t, t1, t2 ComplexNumber
+	var sqrt big.Int
+
+	aRun.Set(a)
+	bRun.Set(b)
+	u.SetOne()
+	v.SetZero()
+	u_.SetZero()
+	v_.SetOne()
+
+	// Eisenstein integers form an Euclidean domain for the norm
+	sqrt.Sqrt(a.Norm())
+	for bRun.Norm().Cmp(&sqrt) >= 0 {
+		quotient.QuoRem(&aRun, &bRun, &remainder)
+		t.Mul(&u_, &quotient)
+		t1.Sub(&u, &t)
+		t.Mul(&v_, &quotient)
+		t2.Sub(&v, &t)
+		aRun.Set(&bRun)
+		u.Set(&u_)
+		v.Set(&v_)
+		bRun.Set(&remainder)
+		u_.Set(&t1)
+		v_.Set(&t2)
+	}
+
+	return [3]*ComplexNumber{&bRun, &v_, &u_}
+}
diff --git a/field/eisenstein/eisenstein_test.go b/field/eisenstein/eisenstein_test.go
new file mode 100644
index 000000000..6aff795f9
--- /dev/null
+++ b/field/eisenstein/eisenstein_test.go
@@ -0,0 +1,279 @@
+package eisenstein
+
+import (
+	"crypto/rand"
+	"math/big"
+	"testing"
+
+	"github.com/leanovate/gopter"
+	"github.com/leanovate/gopter/prop"
+)
+
+const (
+	nbFuzzShort = 10
+	nbFuzz      = 50
+	boundSize   = 128
+)
+
+func TestEisensteinReceiverIsOperand(t *testing.T) {
+
+	t.Parallel()
+	parameters := gopter.DefaultTestParameters()
+	if testing.Short() {
+		parameters.MinSuccessfulTests = nbFuzzShort
+	} else {
+		parameters.MinSuccessfulTests = nbFuzz
+	}
+
+	properties := gopter.NewProperties(parameters)
+
+	genE := GenComplexNumber(boundSize)
+
+	properties.Property("Having the receiver as operand (addition) should output the same result", prop.ForAll(
+		func(a, b *ComplexNumber) bool {
+			var c, d ComplexNumber
+			d.Set(a)
+			c.Add(a, b)
+			a.Add(a, b)
+			b.Add(&d, b)
+			return a.Equal(b) && a.Equal(&c) && b.Equal(&c)
+		},
+		genE,
+		genE,
+	))
+
+	properties.Property("Having the receiver as operand (sub) should output the same result", prop.ForAll(
+		func(a, b *ComplexNumber) bool {
+			var c, d ComplexNumber
+			d.Set(a)
+			c.Sub(a, b)
+			a.Sub(a, b)
+			b.Sub(&d, b)
+			return a.Equal(b) && a.Equal(&c) && b.Equal(&c)
+		},
+		genE,
+		genE,
+	))
+
+	properties.Property("Having the receiver as operand (mul) should output the same result", prop.ForAll(
+		func(a, b *ComplexNumber) bool {
+			var c, d ComplexNumber
+			d.Set(a)
+			c.Mul(a, b)
+			a.Mul(a, b)
+			b.Mul(&d, b)
+			return a.Equal(b) && a.Equal(&c) && b.Equal(&c)
+		},
+		genE,
+		genE,
+	))
+
+	properties.Property("Having the receiver as operand (neg) should output the same result", prop.ForAll(
+		func(a *ComplexNumber) bool {
+			var b ComplexNumber
+			b.Neg(a)
+			a.Neg(a)
+			return a.Equal(&b)
+		},
+		genE,
+	))
+
+	properties.Property("Having the receiver as operand (conjugate) should output the same result", prop.ForAll(
+		func(a *ComplexNumber) bool {
+			var b ComplexNumber
+			b.Conjugate(a)
+			a.Conjugate(a)
+			return a.Equal(&b)
+		},
+		genE,
+	))
+
+	properties.TestingRun(t, gopter.ConsoleReporter(false))
+}
+
+func TestEisensteinArithmetic(t *testing.T) {
+
+	t.Parallel()
+	parameters := gopter.DefaultTestParameters()
+	if testing.Short() {
+		parameters.MinSuccessfulTests = nbFuzzShort
+	} else {
+		parameters.MinSuccessfulTests = nbFuzz
+	}
+
+	properties := gopter.NewProperties(parameters)
+
+	genE := GenComplexNumber(boundSize)
+
+	properties.Property("sub & add should leave an element invariant", prop.ForAll(
+		func(a, b *ComplexNumber) bool {
+			var c ComplexNumber
+			c.Set(a)
+			c.Add(&c, b).Sub(&c, b)
+			return c.Equal(a)
+		},
+		genE,
+		genE,
+	))
+
+	properties.Property("neg twice should leave an element invariant", prop.ForAll(
+		func(a *ComplexNumber) bool {
+			var b ComplexNumber
+			b.Neg(a).Neg(&b)
+			return a.Equal(&b)
+		},
+		genE,
+	))
+
+	properties.Property("conj twice should leave an element invariant", prop.ForAll(
+		func(a *ComplexNumber) bool {
+			var b ComplexNumber
+			b.Conjugate(a).Conjugate(&b)
+			return a.Equal(&b)
+		},
+		genE,
+	))
+
+	properties.Property("add zero should leave element invariant", prop.ForAll(
+		func(a *ComplexNumber) bool {
+			var b, zero ComplexNumber
+			zero.SetZero()
+			b.Add(a, &zero)
+			return a.Equal(&b)
+		},
+		genE,
+	))
+
+	properties.Property("mul by one should leave element invariant", prop.ForAll(
+		func(a *ComplexNumber) bool {
+			var b, one ComplexNumber
+			one.SetOne()
+			b.Mul(a, &one)
+			return a.Equal(&b)
+		},
+		genE,
+	))
+
+	properties.Property("add should be commutative", prop.ForAll(
+		func(a, b *ComplexNumber) bool {
+			var c, d ComplexNumber
+			c.Add(a, b)
+			d.Add(b, a)
+			return c.Equal(&d)
+		},
+		genE,
+		genE,
+	))
+
+	properties.Property("add should be assiocative", prop.ForAll(
+		func(a, b, c *ComplexNumber) bool {
+			var d, e ComplexNumber
+			d.Add(a, b).Add(&d, c)
+			e.Add(c, b).Add(&e, a)
+			return e.Equal(&d)
+		},
+		genE,
+		genE,
+		genE,
+	))
+
+	properties.Property("mul should be commutative", prop.ForAll(
+		func(a, b *ComplexNumber) bool {
+			var c, d ComplexNumber
+			c.Mul(a, b)
+			d.Mul(b, a)
+			return c.Equal(&d)
+		},
+		genE,
+		genE,
+	))
+
+	properties.Property("mul should be assiocative", prop.ForAll(
+		func(a, b, c *ComplexNumber) bool {
+			var d, e ComplexNumber
+			d.Mul(a, b).Mul(&d, c)
+			e.Mul(c, b).Mul(&e, a)
+			return e.Equal(&d)
+		},
+		genE,
+		genE,
+		genE,
+	))
+
+	properties.Property("norm should always be positive", prop.ForAll(
+		func(a *ComplexNumber) bool {
+			return a.Norm().Sign() >= 0
+		},
+		genE,
+	))
+
+	properties.TestingRun(t, gopter.ConsoleReporter(false))
+}
+
+func TestEisensteinHalfGCD(t *testing.T) {
+
+	t.Parallel()
+	parameters := gopter.DefaultTestParameters()
+	if testing.Short() {
+		parameters.MinSuccessfulTests = nbFuzzShort
+	} else {
+		parameters.MinSuccessfulTests = nbFuzz
+	}
+
+	properties := gopter.NewProperties(parameters)
+
+	genE := GenComplexNumber(boundSize)
+
+	properties.Property("half-GCD", prop.ForAll(
+		func(a, b *ComplexNumber) bool {
+			res := HalfGCD(a, b)
+			var c, d ComplexNumber
+			c.Mul(b, res[1])
+			d.Mul(a, res[2])
+			d.Add(&c, &d)
+			return d.Equal(res[0])
+		},
+		genE,
+		genE,
+	))
+
+	properties.TestingRun(t, gopter.ConsoleReporter(false))
+}
+
+// GenNumber generates a random integer
+func GenNumber(boundSize int64) gopter.Gen {
+	return func(genParams *gopter.GenParameters) *gopter.GenResult {
+		var bound big.Int
+		bound.Exp(big.NewInt(2), big.NewInt(boundSize), nil)
+		elmt, _ := rand.Int(genParams.Rng, &bound)
+		genResult := gopter.NewGenResult(elmt, gopter.NoShrinker)
+		return genResult
+	}
+}
+
+// GenComplexNumber generates a random integer
+func GenComplexNumber(boundSize int64) gopter.Gen {
+	return gopter.CombineGens(
+		GenNumber(boundSize),
+		GenNumber(boundSize),
+	).Map(func(values []interface{}) *ComplexNumber {
+		return &ComplexNumber{A0: values[0].(*big.Int), A1: values[1].(*big.Int)}
+	})
+}
+
+// bench
+var benchRes [3]*ComplexNumber
+
+func BenchmarkHalfGCD(b *testing.B) {
+	var n, _ = new(big.Int).SetString("100000000000000000000000000000000", 16) // 2^128
+	a0, _ := rand.Int(rand.Reader, n)
+	a1, _ := rand.Int(rand.Reader, n)
+	c0, _ := rand.Int(rand.Reader, n)
+	c1, _ := rand.Int(rand.Reader, n)
+	a := ComplexNumber{A0: a0, A1: a1}
+	c := ComplexNumber{A0: c0, A1: c1}
+	b.ResetTimer()
+	for i := 0; i < b.N; i++ {
+		benchRes = HalfGCD(&a, &c)
+	}
+}