Adds polynomials and Lagrange polynomials.

math/polynomial/polynomial.go

@@ -0,0 +1,123 @@
+// Package polynomial provides representations of polynomials over the scalars
+// of a group.
+package polynomial
+
+import (
+	"fmt"
+
+	"github.com/cloudflare/circl/group"
+)
+
+// Polynomial stores a polynomial over the set of scalars of a group.
+type Polynomial struct {
+	// Internal representation is in polynomial basis:
+	// Thus,
+	//     p(x) = \sum_i^k c[i] x^i,
+	// where k = len(c)-1 is the degree of the polynomial.
+	c []group.Scalar
+}
+
+// New creates a new polynomial given its coefficients in ascending order.
+// Thus,
+//     p(x) = \sum_i^k c[i] x^i,
+// where k = len(c)-1 is the degree of the polynomial.
+func New(coeffs []group.Scalar) Polynomial {
+	if len(coeffs) == 0 {
+		panic("polynomial: invalid length")
+	}
+	return Polynomial{coeffs}
+}
+
+func (p Polynomial) Degree() uint { return uint(len(p.c) - 1) }
+
+func (p Polynomial) Evaluate(x group.Scalar) group.Scalar {
+	l := len(p.c)
+	px := p.c[l-1].Copy()
+	for i := l - 2; i >= 0; i-- {
+		px.Mul(px, x)
+		px.Add(px, p.c[i])
+	}
+	return px
+}
+
+// LagrangePolynomial stores a Lagrange polynomial over the set of scalars of a group.
+type LagrangePolynomial struct {
+	// Internal representation is in Lagrange basis:
+	// Thus,
+	//     p(x) = \sum_i^k y[i] L_j(x), where k is the degree of the polynomial,
+	//     L_j(x) = \prod_i^k (x-x[i])/(x[j]-x[i]),
+	//     y[i] = p(x[i]), and
+	//     all x[i] are different.
+	g    group.Group
+	x, y []group.Scalar
+}
+
+// NewLagrangePolynomial creates a polynomial in Lagrange basis given a list
+// of nodes (x) and values (y), such that:
+//     p(x) = \sum_i^k y[i] L_j(x), where k is the degree of the polynomial,
+//     L_j(x) = \prod_i^k (x-x[i])/(x[j]-x[i]),
+//     y[i] = p(x[i]), and
+//     all x[i] are different.
+// It panics if one of these conditions does not hold.
+func NewLagrangePolynomial(g group.Group, x, y []group.Scalar) LagrangePolynomial {
+	if len(x) != len(y) {
+		panic("lagrange: invalid length")
+	}
+	if !isAllDifferent(x) {
+		panic("lagrange: nodes must be different")
+	}
+
+	return LagrangePolynomial{g, x, y}
+}
+
+func (l LagrangePolynomial) Degree() uint { return uint(len(l.x) - 1) }
+
+func (l LagrangePolynomial) Evaluate(x group.Scalar) group.Scalar {
+	px := l.g.NewScalar()
+	tmp := l.g.NewScalar()
+	for i := range l.x {
+		LjX := baseRatio(uint(i), l.x, x)
+		tmp.Mul(l.y[i], LjX)
+		px.Add(px, tmp)
+	}
+
+	return px
+}
+
+// LagrangeBase returns the j-th Lagrange polynomial base evaluated at x.
+// Thus, L_j(x) = \prod (x - x[i]) / (x[j] - x[i]) for 0 <= i < k, and i != j.
+func LagrangeBase(jth uint, xi []group.Scalar, x group.Scalar) group.Scalar {
+	if jth >= uint(len(xi)) {
+		panic("lagrange: invalid index")
+	}
+	return baseRatio(jth, xi, x)
+}
+
+func baseRatio(jth uint, xi []group.Scalar, x group.Scalar) group.Scalar {
+	num := x.Copy()
+	num.SetUint64(1)
+	den := x.Copy()
+	den.SetUint64(1)
+
+	tmp := x.Copy()
+	for i := range xi {
+		if uint(i) != jth {
+			num.Mul(num, tmp.Sub(x, xi[i]))
+			den.Mul(den, tmp.Sub(xi[jth], xi[i]))
+		}
+	}
+
+	return num.Mul(num, den.Inv(den))
+}
+
+func isAllDifferent(x []group.Scalar) bool {
+	m := make(map[string]struct{})
+	for i := range x {
+		k := x[i].(fmt.Stringer).String()
+		if _, exists := m[k]; exists {
+			return false
+		}
+		m[k] = struct{}{}
+	}
+	return true
+}

math/polynomial/polynomial_test.go

@@ -0,0 +1,101 @@
+package polynomial_test
+
+import (
+	"testing"
+
+	"github.com/cloudflare/circl/group"
+	"github.com/cloudflare/circl/internal/test"
+	"github.com/cloudflare/circl/math/polynomial"
+)
+
+func TestPolyEval(t *testing.T) {
+	g := group.P256
+	c := []group.Scalar{
+		g.NewScalar(),
+		g.NewScalar(),
+		g.NewScalar(),
+	}
+	c[0].SetUint64(5)
+	c[1].SetUint64(5)
+	c[2].SetUint64(2)
+	p := polynomial.New(c)
+
+	x := g.NewScalar()
+	x.SetUint64(10)
+
+	got := p.Evaluate(x)
+	want := g.NewScalar()
+	want.SetUint64(255)
+	if !got.IsEqual(want) {
+		test.ReportError(t, got, want)
+	}
+}
+
+func TestLagrange(t *testing.T) {
+	g := group.P256
+	c := []group.Scalar{
+		g.NewScalar(),
+		g.NewScalar(),
+		g.NewScalar(),
+	}
+	c[0].SetUint64(1234)
+	c[1].SetUint64(166)
+	c[2].SetUint64(94)
+	p := polynomial.New(c)
+
+	x := []group.Scalar{g.NewScalar(), g.NewScalar(), g.NewScalar()}
+	x[0].SetUint64(2)
+	x[1].SetUint64(4)
+	x[2].SetUint64(5)
+
+	y := []group.Scalar{}
+	for i := range x {
+		y = append(y, p.Evaluate(x[i]))
+	}
+
+	zero := g.NewScalar()
+	l := polynomial.NewLagrangePolynomial(g, x, y)
+	test.CheckOk(l.Degree() == p.Degree(), "bad degree", t)
+
+	got := l.Evaluate(zero)
+	want := p.Evaluate(zero)
+
+	if !got.IsEqual(want) {
+		test.ReportError(t, got, want)
+	}
+
+	// Test Kronecker's delta of LagrangeBase.
+	// Thus:
+	//    L_j(x[i]) = { 1, if i == j;
+	//                { 0, otherwise.
+	one := g.NewScalar()
+	one.SetUint64(1)
+	for j := range x {
+		for i := range x {
+			got := polynomial.LagrangeBase(uint(j), x, x[i])
+
+			if i == j {
+				want = one
+			} else {
+				want = zero
+			}
+
+			if !got.IsEqual(want) {
+				test.ReportError(t, got, want)
+			}
+		}
+	}
+
+	// Test that inputs are different length
+	err := test.CheckPanic(func() { polynomial.NewLagrangePolynomial(g, x, y[:1]) })
+	test.CheckNoErr(t, err, "should panic")
+
+	// Test that nodes must be different.
+	x[0].Set(x[1])
+	err = test.CheckPanic(func() { polynomial.NewLagrangePolynomial(g, x, y) })
+	test.CheckNoErr(t, err, "should panic")
+
+	// Test LagrangeBase wrong index
+	err = test.CheckPanic(func() { polynomial.LagrangeBase(10, x, zero) })
+	test.CheckNoErr(t, err, "should panic")
+}