diff --git a/arithmetic.go b/arithmetic.go
index 2b03711..7dead6c 100644
--- a/arithmetic.go
+++ b/arithmetic.go
@@ -69,6 +69,12 @@ func compareBigIntArray(a, b []*big.Int) bool {
 	return true
 }
 
+//nolint:deadcode,unused
+func checkArrayOfZeroes(a []*big.Int) bool {
+	z := arrayOfZeroes(len(a))
+	return compareBigIntArray(a, z)
+}
+
 func fAdd(a, b *big.Int) *big.Int {
 	ab := new(big.Int).Add(a, b)
 	return ab.Mod(ab, R)
@@ -89,7 +95,6 @@ func fDiv(a, b *big.Int) *big.Int {
 	return new(big.Int).Mod(ab, R)
 }
 
-//nolint:unused,deadcode // TODO check
 func fNeg(a *big.Int) *big.Int {
 	return new(big.Int).Mod(new(big.Int).Neg(a), R)
 }
@@ -170,6 +175,19 @@ func polynomialDiv(a, b []*big.Int) ([]*big.Int, []*big.Int) {
 	return r, rem
 }
 
+func polynomialMulByConstant(a []*big.Int, c *big.Int) []*big.Int {
+	for i := 0; i < len(a); i++ {
+		a[i] = fMul(a[i], c)
+	}
+	return a
+}
+func polynomialDivByConstant(a []*big.Int, c *big.Int) []*big.Int {
+	for i := 0; i < len(a); i++ {
+		a[i] = fDiv(a[i], c)
+	}
+	return a
+}
+
 // polynomialEval evaluates the polinomial over the Finite Field at the given value x
 func polynomialEval(p []*big.Int, x *big.Int) *big.Int {
 	r := big.NewInt(int64(0))
@@ -202,6 +220,16 @@ func newPolZeroAt(pointPos, totalPoints int, height *big.Int) []*big.Int {
 	return r
 }
 
+// zeroPolynomial returns the zero polynomial:
+// z(x) = (x - z_0) (x - z_1) ... (x - z_{k-1})
+func zeroPolynomial(zs []*big.Int) []*big.Int {
+	z := []*big.Int{fNeg(zs[0]), big.NewInt(1)}
+	for i := 1; i < len(zs); i++ {
+		z = polynomialMul(z, []*big.Int{fNeg(zs[i]), big.NewInt(1)})
+	}
+	return z
+}
+
 var sNums = map[string]string{
 	"0": "⁰",
 	"1": "¹",
@@ -237,4 +265,40 @@ func PolynomialToString(p []*big.Int) string {
 	return s
 }
 
+// LagrangeInterpolation implements the Lagrange interpolation:
+// https://en.wikipedia.org/wiki/Lagrange_polynomial
+func LagrangeInterpolation(x, y []*big.Int) ([]*big.Int, error) {
+	// p(x) will be the interpoled polynomial
+	// var p []*big.Int
+	if len(x) != len(y) {
+		return nil, fmt.Errorf("len(x)!=len(y): %d, %d", len(x), len(y))
+	}
+	p := arrayOfZeroes(len(x))
+	k := len(x)
+
+	for j := 0; j < k; j++ {
+		// jPol is the Lagrange basis polynomial for each point
+		var jPol []*big.Int
+		for m := 0; m < k; m++ {
+			// if x[m] == x[j] {
+			if m == j {
+				continue
+			}
+			// numerator & denominator of the current iteration
+			num := []*big.Int{fNeg(x[m]), big.NewInt(1)} // (x^1 - x_m)
+			den := fSub(x[j], x[m])                      // x_j-x_m
+			mPol := polynomialDivByConstant(num, den)
+			if len(jPol) == 0 {
+				// first j iteration
+				jPol = mPol
+				continue
+			}
+			jPol = polynomialMul(jPol, mPol)
+		}
+		p = polynomialAdd(p, polynomialMulByConstant(jPol, y[j]))
+	}
+
+	return p, nil
+}
+
 // TODO add method to 'clean' the polynomial, to remove right-zeroes
diff --git a/arithmetic_test.go b/arithmetic_test.go
index 51a627e..2d21080 100644
--- a/arithmetic_test.go
+++ b/arithmetic_test.go
@@ -80,6 +80,21 @@ func TestPolynomial(t *testing.T) {
 	assert.Equal(t, polynomialEval(o, big.NewInt(3)), b4)
 	o = newPolZeroAt(2, 4, b3)
 	assert.Equal(t, polynomialEval(o, big.NewInt(2)), b3)
+
+	// polynomialEval
+	// p(x) = x^3 + x + 5
+	p := []*big.Int{
+		big.NewInt(5),
+		big.NewInt(1), // x^1
+		big.NewInt(0), // x^2
+		big.NewInt(1), // x^3
+	}
+	assert.Equal(t, "1x³ + 1x¹ + 5", PolynomialToString(p))
+	assert.Equal(t, "35", polynomialEval(p, big.NewInt(3)).String())
+	assert.Equal(t, "1015", polynomialEval(p, big.NewInt(10)).String())
+	assert.Equal(t, "16777477", polynomialEval(p, big.NewInt(256)).String())
+	assert.Equal(t, "125055", polynomialEval(p, big.NewInt(50)).String())
+	assert.Equal(t, "7", polynomialEval(p, big.NewInt(1)).String())
 }
 
 func BenchmarkArithmetic(b *testing.B) {
@@ -109,3 +124,42 @@ func BenchmarkArithmetic(b *testing.B) {
 		}
 	})
 }
+
+func TestLagrangeInterpolation(t *testing.T) {
+	x0 := big.NewInt(3)
+	y0 := big.NewInt(35)
+	x1 := big.NewInt(10)
+	y1 := big.NewInt(1015)
+	x2 := big.NewInt(256)
+	y2 := big.NewInt(16777477)
+	x3 := big.NewInt(50)
+	y3 := big.NewInt(125055)
+
+	xs := []*big.Int{x0, x1, x2, x3}
+	ys := []*big.Int{y0, y1, y2, y3}
+
+	p, err := LagrangeInterpolation(xs, ys)
+	assert.Nil(t, err)
+	assert.Equal(t, "1x³ + 1x¹ + 5", PolynomialToString(p))
+
+	assert.Equal(t, y0, polynomialEval(p, x0))
+	assert.Equal(t, y1, polynomialEval(p, x1))
+	assert.Equal(t, y2, polynomialEval(p, x2))
+}
+
+func TestZeroPolynomial(t *testing.T) {
+	x0 := big.NewInt(1)
+	x1 := big.NewInt(40)
+	x2 := big.NewInt(512)
+	xs := []*big.Int{x0, x1, x2}
+
+	z := zeroPolynomial(xs)
+	assert.Equal(t, "1x³ "+
+		"+ 21888242871839275222246405745257275088548364400416034343698204186575808495064x² "+
+		"+ 21032x¹ + 21888242871839275222246405745257275088548364400416034343698204186575808475137",
+		PolynomialToString(z))
+
+	assert.Equal(t, "0", polynomialEval(z, x0).String())
+	assert.Equal(t, "0", polynomialEval(z, x1).String())
+	assert.Equal(t, "0", polynomialEval(z, x2).String())
+}