fix: lagrange interpolation

barretenberg/cpp/src/barretenberg/polynomials/polynomial_arithmetic.cpp

@@ -1210,6 +1210,8 @@ template <typename Fr> Fr compute_sum(const Fr* src, const size_t n)
 // This function computes the polynomial (x - a)(x - b)(x - c)... given n distinct roots (a, b, c, ...).
 template <typename Fr> void compute_linear_polynomial_product(const Fr* roots, Fr* dest, const size_t n)
 {
+    for (size_t i = 0; i < n; ++i) {
+    };
     auto scratch_space_ptr = get_scratch_space<Fr>(n);
     auto scratch_space = scratch_space_ptr.get();
     memcpy((void*)scratch_space, (void*)roots, n * sizeof(Fr));
@@ -1328,15 +1330,14 @@ void compute_efficient_interpolation(const Fr* src, Fr* dest, const Fr* evaluati
     */
     Fr numerator_polynomial[n + 1];
     polynomial_arithmetic::compute_linear_polynomial_product(evaluation_points, numerator_polynomial, n);
-
+    // First half contains roots, second half contains denominators (to be inverted)
     Fr roots_and_denominators[2 * n];
     Fr temp_src[n];
     for (size_t i = 0; i < n; ++i) {
         roots_and_denominators[i] = -evaluation_points[i];
         temp_src[i] = src[i];
         dest[i] = 0;
-
-        // compute constant denominator
+        // compute constant denominators
         roots_and_denominators[n + i] = 1;
         for (size_t j = 0; j < n; ++j) {
             if (j == i) {
@@ -1345,21 +1346,65 @@ void compute_efficient_interpolation(const Fr* src, Fr* dest, const Fr* evaluati
             roots_and_denominators[n + i] *= (evaluation_points[i] - evaluation_points[j]);
         }
     }
-
+    // at this point roots_and_denominators is populated as follows
+    // (x_0,\ldots, x_{n-1}, d_0, \ldots, d_{n-1})
     Fr::batch_invert(roots_and_denominators, 2 * n);
 
     Fr z, multiplier;
     Fr temp_dest[n];
-    for (size_t i = 0; i < n; ++i) {
-        z = roots_and_denominators[i];
-        multiplier = temp_src[i] * roots_and_denominators[n + i];
-        temp_dest[0] = multiplier * numerator_polynomial[0];
-        temp_dest[0] *= z;
-        dest[0] += temp_dest[0];
-        for (size_t j = 1; j < n; ++j) {
-            temp_dest[j] = multiplier * numerator_polynomial[j] - temp_dest[j - 1];
-            temp_dest[j] *= z;
-            dest[j] += temp_dest[j];
+    size_t idx_zero = 0;
+    bool interpolation_domain_contains_zero = false;
+    if (numerator_polynomial[0] == Fr(0)) {
+        for (size_t i = 0; i < n; ++i) {
+            if (evaluation_points[i] == Fr(0)) {
+                idx_zero = i;
+                interpolation_domain_contains_zero = true;
+                break;
+            }
+        }
+    };
+    if (!interpolation_domain_contains_zero) {
+        for (size_t i = 0; i < n; ++i) {
+            // set z = - 1/x_i for x_i <> 0
+            z = roots_and_denominators[i];
+            // temp_src[i] is y_i, it gets multiplied by 1/d_i
+            multiplier = temp_src[i] * roots_and_denominators[n + i];
+            temp_dest[0] = multiplier * numerator_polynomial[0];
+            temp_dest[0] *= z;
+            dest[0] += temp_dest[0];
+            for (size_t j = 1; j < n; ++j) {
+                temp_dest[j] = multiplier * numerator_polynomial[j] - temp_dest[j - 1];
+                temp_dest[j] *= z;
+                dest[j] += temp_dest[j];
+            }
+        }
+    } else {
+        for (size_t i = 0; i < n; ++i) {
+            if (i == idx_zero) {
+                // the contribution from the term corresponding to idx_zero is computed separately
+                continue;
+            } else {
+                // get the next inverted root
+                z = roots_and_denominators[i];
+                // compute f(x_i) * d_{x_i}^{-1}
+                multiplier = temp_src[i] * roots_and_denominators[n + i];
+                // get x_i^{-1} * f(x_i) * d_{x_i}^{-1} into the "free" term
+                temp_dest[1] = multiplier * numerator_polynomial[1];
+                temp_dest[1] *= z;
+                // correct the first coefficient as it is now accumulating free terms from
+                // f(x_i) d_i^{-1} prod_(X-x_i, x_i <> 0) (X-x_i) * 1/(X-x_i)
+                dest[1] += temp_dest[1];
+                // compute the quotient N(X)/(X-x_i) f(x_i)/d_{x_i} and its contribution to the target coefficients
+                for (size_t j = 2; j < n; ++j) {
+                    temp_dest[j] = multiplier * numerator_polynomial[j] - temp_dest[j - 1];
+                    temp_dest[j] *= z;
+                    dest[j] += temp_dest[j];
+                };
+            }
+        }
+        // correct the target coefficients by the contribution from N(X)/X * d_{idx_zero}^{-1} * f(idx_zero)
+        for (size_t i = 0; i < n; ++i) {
+            dest[i] += temp_src[idx_zero] * roots_and_denominators[n + idx_zero] * numerator_polynomial[i + 1];
         }
     }
 }

barretenberg/cpp/src/barretenberg/polynomials/polynomial_arithmetic.test.cpp

@@ -999,6 +999,58 @@ TYPED_TEST(PolynomialTests, compute_efficient_interpolation)
         EXPECT_EQ(src[i], poly[i]);
     }
 }
+// Test efficient Lagrange interpolation when interpolation points contain 0
+TYPED_TEST(PolynomialTests, compute_efficient_interpolation_domain_with_zero)
+{
+    using FF = TypeParam;
+    constexpr size_t n = 15;
+    FF src[n], poly[n], x[n];
+
+    for (size_t i = 0; i < n; ++i) {
+        poly[i] = FF(i + 1);
+    }
+
+    for (size_t i = 0; i < n; ++i) {
+        x[i] = FF(i);
+        src[i] = polynomial_arithmetic::evaluate(poly, x[i], n);
+    }
+    polynomial_arithmetic::compute_efficient_interpolation(src, src, x, n);
+
+    for (size_t i = 0; i < n; ++i) {
+        EXPECT_EQ(src[i], poly[i]);
+    }
+    // Test for the domain (-n/2, ..., 0, ... , n/2)
+
+    for (size_t i = 0; i < n; ++i) {
+        poly[i] = FF(i + 54);
+    }
+
+    for (size_t i = 0; i < n; ++i) {
+        x[i] = FF(i - n / 2);
+        src[i] = polynomial_arithmetic::evaluate(poly, x[i], n);
+    }
+    polynomial_arithmetic::compute_efficient_interpolation(src, src, x, n);
+
+    for (size_t i = 0; i < n; ++i) {
+        EXPECT_EQ(src[i], poly[i]);
+    }
+
+    // Test for the domain (-n+1, ..., 0)
+
+    for (size_t i = 0; i < n; ++i) {
+        poly[i] = FF(i * i + 57);
+    }
+
+    for (size_t i = 0; i < n; ++i) {
+        x[i] = FF(i - (n - 1));
+        src[i] = polynomial_arithmetic::evaluate(poly, x[i], n);
+    }
+    polynomial_arithmetic::compute_efficient_interpolation(src, src, x, n);
+
+    for (size_t i = 0; i < n; ++i) {
+        EXPECT_EQ(src[i], poly[i]);
+    }
+}
 
 TYPED_TEST(PolynomialTests, interpolation_constructor_single)
 {