revert: back to lagrange

internal/generator/polynomial/template/polynomial.go.tmpl

@@ -302,58 +302,63 @@ func computeLagrangeBasis(domainSize uint8) []Polynomial {
 	return res
 }
 
-
 // Interpolate fits a polynomial of degree len(X) - 1 = len(Y) - 1 to the points (X[i], Y[i])
-// Note that the runtime is O(len(X)³)
-func Interpolate(X, Y []{{.ElementType}}) (Polynomial, error) {
+// Returns the evaluation at r
+//
+// Note that this runs in O(len(X)²).
+// When possible, it is better to exploit any structure in the domain
+// (i.e. the values in X) to compute the polynomial in linear time.
+func Interpolate(X, Y []{{.ElementType}}, r {{.ElementType}}) {{.ElementType}} {
 	if len(X) != len(Y) {
-		return nil, errors.New("X and Y must have the same length")
+		panic("X and Y must have the same length")
 	}
-
-	// solve the system of equations by Gaussian elimination
-	augmentedRows := make([][]{{.ElementType}}, len(X)) // the last column is the Y values
-	for i := range augmentedRows {
-		augmentedRows[i] = make([]{{.ElementType}}, len(X)+1)
-		augmentedRows[i][0].SetOne()
-		augmentedRows[i][1].Set(&X[i])
-		for j := 2; j < len(augmentedRows[i])-1; j++ {
-			augmentedRows[i][j].Mul(&augmentedRows[i][j-1], &X[i])
-		}
-		augmentedRows[i][len(augmentedRows[i])-1].Set(&Y[i])
+	switch len(X) {
+	case 0:
+		return {{.ElementType}}{}
+	case 1:
+		return Y[0]
 	}
 
-	// make the upper triangle
-	for i := range len(augmentedRows) - 1 {
-		// use row i to eliminate the ith element in all rows below
-		var negInv {{.ElementType}}
-		if augmentedRows[i][i].IsZero() {
-			return nil, errors.New("singular matrix")
-		}
-		negInv.Inverse(&augmentedRows[i][i])
-		negInv.Neg(&negInv)
-		for j := i + 1; j < len(augmentedRows); j++ {
-			var c {{.ElementType}}
-			c.Mul(&augmentedRows[j][i], &negInv)
-			// augmentedRows[j][i].SetZero() omitted
-			for k := i + 1; k < len(augmentedRows[i]); k++ {
-				var t {{.ElementType}}
-				t.Mul(&augmentedRows[i][k], &c)
-				augmentedRows[j][k].Add(&augmentedRows[j][k], &t)
+	// weights[i] = 1 / [(X[i] - X[0]) * (X[i] - X[1]) * ... * (X[i] - X[i-1]) * (X[i] - X[i+1]) * ... * (X[i] - X[n-1])]
+	weights := make([]{{.ElementType}}, len(X))
+	for i := range weights {
+		weights[i].SetOne()
+		for j := range weights {
+			if i == j {
+				continue
 			}
-		}
-	}
-
-	// back substitution
-	res := make(Polynomial, len(X))
-	for i := len(augmentedRows) - 1; i >= 0; i-- {
-		res[i] = augmentedRows[i][len(augmentedRows[i])-1]
-		for j := i + 1; j < len(augmentedRows[i])-1; j++ {
 			var t {{.ElementType}}
-			t.Mul(&res[j], &augmentedRows[i][j])
-			res[i].Sub(&res[i], &t)
+			t.Sub(&X[i], &X[j])
+			weights[i].Mul(&weights[i], &t)
 		}
-		res[i].Div(&res[i], &augmentedRows[i][i])
 	}
+	weights = {{.FieldPackageName}}.BatchInvert(weights)
+
+	// prods[i] = (r-X[i+1]) * (r-X[i+2]) * ... (r-X[n-1])
+	prods := make([]{{.ElementType}}, len(X))
+	prods[len(X)-1].SetOne()
+	prods[len(X)-2].Sub(&r, &X[len(X)-1])
+	for i := len(X) - 3; i >= 0; i-- {
+		var t {{.ElementType}}
+		t.Sub(&r, &X[i+1])
+		prods[i].Mul(&prods[i+1], &t)
+	}
+
+	// at iteration i, prod = (r-X[0]) * (r-X[1]) * ... * (r-X[i-1])
+	prod := {{.FieldPackageName}}.One()
+	var res {{.ElementType}}
+
+	for i := range X {
+		var t {{.ElementType}}
+		t.Mul(&prod, &prods[i]) // t = (r-X[0]) * (r-X[1]) * ... * (r-X[i-1]) * (r-X[i+1]) * ... * (r-X[n-1])
+		t.Mul(&t, &weights[i])  // t = (r-X[0]) * (r-X[1]) * ... * (r-X[i-1]) * (r-X[i+1]) * ... * (r-X[n-1]) / [(X[i] - X[0]) * (X[i] - X[1]) * ... * (X[i] - X[i-1]) * (X[i] - X[i+1]) * ... * (X[i] - X[n-1])]
+		// t is now the evaluation of the i'th Lagrange basis at r
 
-	return res, nil
+		t.Mul(&t, &Y[i])
+		res.Add(&res, &t)
+
+		t.Sub(&r, &X[i])
+		prod.Mul(&prod, &t)
+	}
+	return res
 }