diff --git a/fastcrypto-tbls/src/polynomial.rs b/fastcrypto-tbls/src/polynomial.rs
index 710bbc815c..7f21fe7e7a 100644
--- a/fastcrypto-tbls/src/polynomial.rs
+++ b/fastcrypto-tbls/src/polynomial.rs
@@ -107,11 +107,32 @@ impl<C: GroupElement> Poly<C> {
             .iter()
             .fold(C::ScalarType::generator(), |acc, i| acc * i);
         let mut coeffs = Vec::new();
-        for i in &indices {
-            let denominator = indices
-                .iter()
-                .filter(|j| *j != i)
-                .fold(*i, |acc, j| acc * (*j - i));
+
+        // Keep track of differences that are already computed: deltas[i][j] = indices[j] - indices[n - 1 - i]
+        let n = indices.len();
+        let mut deltas: Vec<Vec<<C as GroupElement>::ScalarType>> = vec![Vec::new(); n];
+        for (i_idx, i) in indices.iter().enumerate() {
+            // Compute product with cached deltas
+            let mut denominator = (0..i_idx)
+                .map(|k| {
+                    deltas[k]
+                        .pop()
+                        .expect("deltas[k] is popped k-1 times and contains exactly k-1 elements")
+                })
+                .fold(*i, |acc, delta| acc * delta);
+
+            // Adjust sign since the cached values are negations of the factors we need.
+            if i_idx % 2 == 1 {
+                denominator = -denominator;
+            }
+
+            // Compute remaining deltas
+            for j in indices[i_idx + 1..].iter().rev() {
+                let delta = *j - i;
+                denominator = denominator * delta;
+                deltas[i_idx].push(delta);
+            }
+
             let coeff = full_numerator / denominator;
             coeffs.push(coeff.expect("safe since i != j"));
         }