diff --git a/src/Makefile b/src/Makefile
index c4b282b..20dcc24 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,6 +1,6 @@
 TESTS = bls12_381_test das_extension_test c_kzg_alloc_test fft_common_test fft_fr_test fft_g1_test \
 	fk20_proofs_test kzg_proofs_test poly_test recover_test utility_test zero_poly_test
-BENCH = fft_fr_bench fft_g1_bench recover_bench zero_poly_bench kzg_proofs_bench poly_bench
+BENCH = fft_fr_bench fft_g1_bench recover_bench zero_poly_bench kzg_proofs_bench poly_bench poly_barycentric_bench compute_proof_lagrange_bench
 TUNE = poly_mul_tune poly_div_tune
 LIB_SRC = bls12_381.c c_kzg_alloc.c das_extension.c fft_common.c fft_fr.c fft_g1.c fk20_proofs.c kzg_proofs.c poly.c recover.c utility.c zero_poly.c
 LIB_OBJ = $(LIB_SRC:.c=.o)
diff --git a/src/c_kzg.h b/src/c_kzg.h
index f699889..5c5a5e9 100644
--- a/src/c_kzg.h
+++ b/src/c_kzg.h
@@ -86,14 +86,22 @@ typedef struct {
     uint64_t length; /**< One more than the polynomial's degree */
 } poly;
 
+typedef struct {
+    fr_t *values;    /**< `values[i]` is value of the polynomial at `ω^i`. */
+    uint64_t length; /**< One more than the polynomial's degree */
+} poly_l; // Lagrange form
+
 void eval_poly(fr_t *out, const poly *p, const fr_t *x);
+C_KZG_RET eval_poly_l(fr_t *out, const poly_l *p, const fr_t *x, const FFTSettings *fs);
 C_KZG_RET poly_inverse(poly *out, poly *b);
 C_KZG_RET poly_mul(poly *out, const poly *a, const poly *b);
 C_KZG_RET poly_mul_(poly *out, const poly *a, const poly *b, FFTSettings *fs);
 C_KZG_RET new_poly_div(poly *out, const poly *dividend, const poly *divisor);
 C_KZG_RET new_poly(poly *out, uint64_t length);
+C_KZG_RET new_poly_l(poly_l *out, uint64_t length);
 C_KZG_RET new_poly_with_coeffs(poly *out, const fr_t *coeffs, uint64_t length);
 void free_poly(poly *p);
+void free_poly_l(poly_l *p);
 
 //
 // kzg_proofs.c
@@ -107,12 +115,17 @@ void free_poly(poly *p);
 typedef struct {
     const FFTSettings *fs; /**< The corresponding settings for performing FFTs */
     g1_t *secret_g1;       /**< G1 group elements from the trusted setup */
+    g1_t *secret_g1_l;     /**< secret_g1 in Lagrange form */
     g2_t *secret_g2;       /**< G2 group elements from the trusted setup */
     uint64_t length;       /**< The number of elements in secret_g1 and secret_g2 */
 } KZGSettings;
 
+C_KZG_RET new_poly_l_from_poly(poly_l *out, const poly *in, const KZGSettings *ks);
+
 C_KZG_RET commit_to_poly(g1_t *out, const poly *p, const KZGSettings *ks);
+C_KZG_RET commit_to_poly_l(g1_t *out, const poly_l *p, const KZGSettings *ks);
 C_KZG_RET compute_proof_single(g1_t *out, const poly *p, const fr_t *x0, const KZGSettings *ks);
+C_KZG_RET compute_proof_single_l(g1_t *out, const poly_l *p, const fr_t *x0, const fr_t *y, const KZGSettings *ks);
 C_KZG_RET check_proof_single(bool *out, const g1_t *commitment, const g1_t *proof, const fr_t *x, fr_t *y,
                              const KZGSettings *ks);
 C_KZG_RET compute_proof_multi(g1_t *out, const poly *p, const fr_t *x0, uint64_t n, const KZGSettings *ks);
diff --git a/src/compute_proof_lagrange_bench.c b/src/compute_proof_lagrange_bench.c
new file mode 100644
index 0000000..a77c81a
--- /dev/null
+++ b/src/compute_proof_lagrange_bench.c
@@ -0,0 +1,96 @@
+/*
+ * Copyright 2021 Benjamin Edgington
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <stdlib.h> // malloc(), free(), atoi()
+#include <stdio.h>  // printf()
+#include <assert.h> // assert()
+#include <unistd.h> // EXIT_SUCCESS/FAILURE
+#include "bench_util.h"
+#include "test_util.h"
+#include "c_kzg.h"
+
+// Run the benchmark for `max_seconds` and return the time per iteration in nanoseconds.
+long run_bench(int scale, int max_seconds) {
+    timespec_t t0, t1;
+    unsigned long total_time = 0, nits = 0;
+    FFTSettings fs;
+    KZGSettings ks;
+
+    assert(C_KZG_OK == new_fft_settings(&fs, scale));
+
+    // Allocate on the heap to avoid stack overflow for large sizes
+    g1_t *s1 = malloc(fs.max_width * sizeof(g1_t));
+    g2_t *s2 = malloc(fs.max_width * sizeof(g2_t));
+
+    generate_trusted_setup(s1, s2, &secret, fs.max_width);
+    assert(C_KZG_OK == new_kzg_settings(&ks, s1, s2, fs.max_width, &fs));
+
+    poly_l p;
+    assert(C_KZG_OK == new_poly_l(&p, fs.max_width));
+    for (int i = 0; i < fs.max_width; i++) {
+        p.values[i] = rand_fr();
+    }
+
+    fr_t x = rand_fr();
+    fr_t y;
+    assert(C_KZG_OK == eval_poly_l(&y, &p, &x, &fs));
+
+    while (total_time < max_seconds * NANO) {
+        g1_t proof;
+        clock_gettime(CLOCK_REALTIME, &t0);
+
+        assert(C_KZG_OK == compute_proof_single_l(&proof, &p, &x, &y, &ks));
+
+        clock_gettime(CLOCK_REALTIME, &t1);
+        nits++;
+        total_time += tdiff(t0, t1);
+    }
+
+    free_poly_l(&p);
+    free(s1);
+    free(s2);
+    free_kzg_settings(&ks);
+    free_fft_settings(&fs);
+
+    return total_time / nits;
+}
+
+int main(int argc, char *argv[]) {
+    int nsec = 0;
+
+    switch (argc) {
+    case 1:
+        nsec = NSEC;
+        break;
+    case 2:
+        nsec = atoi(argv[1]);
+        break;
+    default:
+        break;
+    };
+
+    if (nsec == 0) {
+        printf("Usage: %s [test time in seconds > 0]\n", argv[0]);
+        exit(EXIT_FAILURE);
+    }
+
+    printf("*** Benchmarking KZG Proof from Lagrange, %d second%s per test.\n", nsec, nsec == 1 ? "" : "s");
+    for (int scale = 1; scale <= 15; scale++) {
+        printf("compute_proof_single_l/scale_%d %lu ns/op\n", scale, run_bench(scale, nsec));
+    }
+
+    return EXIT_SUCCESS;
+}
diff --git a/src/fft_g1.c b/src/fft_g1.c
index cfe08a0..ec13255 100644
--- a/src/fft_g1.c
+++ b/src/fft_g1.c
@@ -198,4 +198,4 @@ TEST_LIST = {
     {NULL, NULL} /* zero record marks the end of the list */
 };
 
-#endif // KZGTEST
\ No newline at end of file
+#endif // KZGTEST
diff --git a/src/fk20_proofs.c b/src/fk20_proofs.c
index b692de7..c373caa 100644
--- a/src/fk20_proofs.c
+++ b/src/fk20_proofs.c
@@ -485,7 +485,7 @@ void fk_single(void) {
     FFTSettings fs;
     KZGSettings ks;
     FK20SingleSettings fk;
-    uint64_t secrets_len = n_len + 1;
+    uint64_t secrets_len = n_len;
     g1_t s1[secrets_len];
     g2_t s2[secrets_len];
     poly p;
@@ -558,7 +558,7 @@ void fk_single_strided(void) {
     FFTSettings fs;
     KZGSettings ks;
     FK20SingleSettings fk;
-    uint64_t secrets_len = n_len + 1;
+    uint64_t secrets_len = n_len;
     g1_t s1[secrets_len];
     g2_t s2[secrets_len];
     poly p;
@@ -606,7 +606,7 @@ void fk_multi_settings(void) {
     KZGSettings ks;
     FK20MultiSettings fk;
     uint64_t n = 5;
-    uint64_t secrets_len = 33;
+    uint64_t secrets_len = 32;
     g1_t s1[secrets_len];
     g2_t s2[secrets_len];
 
@@ -764,4 +764,4 @@ TEST_LIST = {
     {NULL, NULL} /* zero record marks the end of the list */
 };
 
-#endif
\ No newline at end of file
+#endif
diff --git a/src/kzg_proofs.c b/src/kzg_proofs.c
index aa5a087..0eb2285 100644
--- a/src/kzg_proofs.c
+++ b/src/kzg_proofs.c
@@ -43,6 +43,21 @@ C_KZG_RET commit_to_poly(g1_t *out, const poly *p, const KZGSettings *ks) {
     return C_KZG_OK;
 }
 
+/**
+ * Make a KZG commitment to a polynomial in Lagrange form.
+ *
+ * @param[out] out The commitment to the polynomial, in the form of a G1 group point
+ * @param[in]  p_l The polynomial to be committed to
+ * @param[in]  ks  The settings containing the secrets, previously initialised with #new_kzg_settings
+ * @retval C_CZK_OK      All is well
+ * @retval C_CZK_BADARGS Invalid parameters were supplied
+ */
+C_KZG_RET commit_to_poly_l(g1_t *out, const poly_l *p_l, const KZGSettings *ks) {
+    CHECK(p_l->length <= ks->length);
+    g1_linear_combination(out, ks->secret_g1_l, p_l->values, p_l->length);
+    return C_KZG_OK;
+}
+
 /**
  * Compute KZG proof for polynomial at position x0.
  *
@@ -85,6 +100,66 @@ C_KZG_RET check_proof_single(bool *out, const g1_t *commitment, const g1_t *proo
     return C_KZG_OK;
 }
 
+/**
+ * Compute KZG proof for polynomial in Lagrange form at position x0
+ *
+ * @param[out] out The combined proof as a single G1 element
+ * @param[in]  p   The polynomial in Lagrange form
+ * @param[in]  x  The generator x-value for the evaluation points
+ * @param[in]  y   The value of @p p at @p x
+ * @param[in]  ks  The settings containing the secrets, previously initialised with #new_kzg_settings
+ * @retval C_KZG_OK      All is well
+ * @retval C_KZG_ERROR   An internal error occurred
+ * @retval C_KZG_MALLOC  Memory allocation failed
+ */
+C_KZG_RET compute_proof_single_l(g1_t *out, const poly_l *p, const fr_t *x, const fr_t *y, const KZGSettings *ks) {
+  fr_t tmp, tmp2;
+  poly_l q;
+  uint64_t i, m = 0;
+
+  new_poly_l(&q, p->length);
+
+  fr_t *inverses_in, *inverses;
+
+  TRY(new_fr_array(&inverses_in, p->length));
+  TRY(new_fr_array(&inverses, p->length));
+
+  for (i = 0; i < q.length; i++) {
+    if (fr_equal(x, &ks->fs->expanded_roots_of_unity[i])) {
+      m = i + 1;
+      continue;
+    }
+    // (p_i - y) / (ω_i - x)
+    fr_sub(&q.values[i], &p->values[i], y);
+    fr_sub(&inverses_in[i], &ks->fs->expanded_roots_of_unity[i], x);
+  }
+
+  TRY(fr_batch_inv(inverses, inverses_in, q.length));
+
+  for (i = 0; i < q.length; i++) {
+    fr_mul(&q.values[i], &q.values[i], &inverses[i]);
+  }  
+  if (m) { // ω_m == x
+    q.values[--m] = fr_zero;
+    for (i=0; i < q.length; i++) {
+      if (i == m) continue;
+      // (p_i - y) * ω_i / (x * (x - ω_i))
+      fr_sub(&tmp, x, &ks->fs->expanded_roots_of_unity[i]);
+      fr_mul(&inverses_in[i], &tmp, x);
+    }
+    TRY(fr_batch_inv(inverses, inverses_in, q.length));
+    for (i=0; i < q.length; i++) {
+      fr_sub(&tmp2, &p->values[i], y);
+      fr_mul(&tmp, &tmp2, &inverses[i]);
+      fr_mul(&tmp, &tmp, &ks->fs->expanded_roots_of_unity[i]);
+      fr_add(&q.values[m], &q.values[m], &tmp);
+    }
+  }
+  free(inverses_in);
+  free(inverses);
+  return commit_to_poly_l(out, &q, ks);
+}
+
 /**
  * Compute KZG proof for polynomial at positions x0 * w^y where w is an n-th root of unity.
  *
@@ -96,10 +171,10 @@ C_KZG_RET check_proof_single(bool *out, const g1_t *commitment, const g1_t *proo
  * @param[in]  x0  The generator x-value for the evaluation points
  * @param[in]  n   The number of points at which to evaluate the polynomial, must be a power of two
  * @param[in]  ks  The settings containing the secrets, previously initialised with #new_kzg_settings
- * @retval C_CZK_OK      All is well
- * @retval C_CZK_BADARGS Invalid parameters were supplied
- * @retval C_CZK_ERROR   An internal error occurred
- * @retval C_CZK_MALLOC  Memory allocation failed
+ * @retval C_KZG_OK      All is well
+ * @retval C_KZG_BADARGS Invalid parameters were supplied
+ * @retval C_KZG_ERROR   An internal error occurred
+ * @retval C_KZG_MALLOC  Memory allocation failed
  */
 C_KZG_RET compute_proof_multi(g1_t *out, const poly *p, const fr_t *x0, uint64_t n, const KZGSettings *ks) {
     poly divisor, q;
@@ -216,6 +291,7 @@ C_KZG_RET new_kzg_settings(KZGSettings *ks, const g1_t *secret_g1, const g2_t *s
 
     // Allocate space for the secrets
     TRY(new_g1_array(&ks->secret_g1, ks->length));
+    TRY(new_g1_array(&ks->secret_g1_l, ks->length));
     TRY(new_g2_array(&ks->secret_g2, ks->length));
 
     // Populate the secrets
@@ -225,7 +301,8 @@ C_KZG_RET new_kzg_settings(KZGSettings *ks, const g1_t *secret_g1, const g2_t *s
     }
     ks->fs = fs;
 
-    return C_KZG_OK;
+    // Add Lagrange form (and return its success)
+    return fft_g1(ks->secret_g1_l, ks->secret_g1, true, length, fs);
 }
 
 /**
@@ -235,6 +312,7 @@ C_KZG_RET new_kzg_settings(KZGSettings *ks, const g1_t *secret_g1, const g2_t *s
  */
 void free_kzg_settings(KZGSettings *ks) {
     free(ks->secret_g1);
+    free(ks->secret_g1_l);
     free(ks->secret_g2);
     ks->length = 0;
 }
@@ -248,7 +326,7 @@ void proof_single(void) {
     // Our polynomial: degree 15, 16 coefficients
     uint64_t coeffs[] = {1, 2, 3, 4, 7, 7, 7, 7, 13, 13, 13, 13, 13, 13, 13, 13};
     int poly_len = sizeof coeffs / sizeof coeffs[0];
-    uint64_t secrets_len = poly_len + 1;
+    uint64_t secrets_len = poly_len;
 
     FFTSettings fs;
     KZGSettings ks;
@@ -291,13 +369,116 @@ void proof_single(void) {
     free_poly(&p);
 }
 
+void proof_single_l(void) {
+    // Our polynomial: degree 15, 16 coefficients
+    uint64_t coeffs[] = {1, 2, 3, 4, 7, 7, 7, 7, 13, 13, 13, 13, 13, 13, 13, 13};
+    int poly_len = sizeof coeffs / sizeof coeffs[0];
+    uint64_t secrets_len = poly_len;
+
+    FFTSettings fs;
+    KZGSettings ks;
+    g1_t s1[secrets_len];
+    g2_t s2[secrets_len];
+    poly p;
+    poly_l p_l;
+    g1_t commitment, proof;
+    fr_t x, value;
+    bool result;
+
+    // Create the polynomial
+    new_poly(&p, poly_len);
+    for (int i = 0; i < poly_len; i++) {
+        fr_from_uint64(&p.coeffs[i], coeffs[i]);
+    }
+
+    // Initialise the secrets and data structures
+    generate_trusted_setup(s1, s2, &secret, secrets_len);
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4)); // ln_2 of poly_len
+    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks, s1, s2, secrets_len, &fs));
+
+    // Create Lagrange form
+    new_poly_l_from_poly(&p_l, &p, &ks);
+
+    // Compute the proof for x = 25
+    fr_from_uint64(&x, 25);
+    TEST_CHECK(C_KZG_OK == commit_to_poly_l(&commitment, &p_l, &ks));
+    eval_poly_l(&value, &p_l, &x, &fs);
+    TEST_CHECK(C_KZG_OK == compute_proof_single_l(&proof, &p_l, &x, &value, &ks));
+
+    // Verify the proof that the (unknown) polynomial has y = value at x = 25
+    TEST_CHECK(C_KZG_OK == check_proof_single(&result, &commitment, &proof, &x, &value, &ks));
+    TEST_CHECK(true == result);
+
+    // Change the value and check that the proof fails
+    fr_add(&value, &value, &fr_one);
+    TEST_CHECK(C_KZG_OK == check_proof_single(&result, &commitment, &proof, &x, &value, &ks));
+    TEST_CHECK(false == result);
+
+    free_fft_settings(&fs);
+    free_kzg_settings(&ks);
+    free_poly(&p);
+    free_poly_l(&p_l);
+}
+
+void proof_single_l_at_root(void) {
+    // Our polynomial: degree 15, 16 coefficients
+    uint64_t coeffs[] = {3, 2, 13, 4, 7, 7, 9, 7, 9913, 13, 8813, 13, 7713, 13, 5513, 14};
+    int poly_len = sizeof coeffs / sizeof coeffs[0];
+    uint64_t secrets_len = poly_len;
+
+    FFTSettings fs;
+    KZGSettings ks;
+    g1_t s1[secrets_len];
+    g2_t s2[secrets_len];
+    poly p;
+    poly_l p_l;
+    g1_t commitment, proof;
+    fr_t value;
+    fr_t *x;
+    bool result;
+
+    // Create the polynomial
+    new_poly(&p, poly_len);
+    for (int i = 0; i < poly_len; i++) {
+        fr_from_uint64(&p.coeffs[i], coeffs[i]);
+    }
+
+    // Initialise the secrets and data structures
+    generate_trusted_setup(s1, s2, &secret, secrets_len);
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4)); // ln_2 of poly_len
+    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks, s1, s2, secrets_len, &fs));
+
+    // Create Lagrange form
+    new_poly_l_from_poly(&p_l, &p, &ks);
+
+    // Compute the proof for x = the 5th root of unity
+    x = &fs.expanded_roots_of_unity[6];
+    TEST_CHECK(C_KZG_OK == commit_to_poly_l(&commitment, &p_l, &ks));
+    eval_poly_l(&value, &p_l, x, &fs);
+    TEST_CHECK(C_KZG_OK == compute_proof_single_l(&proof, &p_l, x, &value, &ks));
+
+    // Verify the proof that the (unknown) polynomial has y = value at x = ω_5
+    TEST_CHECK(C_KZG_OK == check_proof_single(&result, &commitment, &proof, x, &value, &ks));
+    TEST_CHECK(true == result);
+
+    // Change the value and check that the proof fails
+    fr_add(&value, &value, &fr_one);
+    TEST_CHECK(C_KZG_OK == check_proof_single(&result, &commitment, &proof, x, &value, &ks));
+    TEST_CHECK(false == result);
+
+    free_fft_settings(&fs);
+    free_kzg_settings(&ks);
+    free_poly(&p);
+    free_poly_l(&p_l);
+}
+
 void proof_multi(void) {
     // Our polynomial: degree 15, 16 coefficients
     uint64_t coeffs[] = {1, 2, 3, 4, 7, 7, 7, 7, 13, 13, 13, 13, 13, 13, 13, 13};
     int poly_len = sizeof coeffs / sizeof coeffs[0];
 
-    FFTSettings fs1, fs2;
-    KZGSettings ks1, ks2;
+    FFTSettings fs;
+    KZGSettings ks;
     poly p;
     g1_t commitment, proof;
     fr_t x, tmp;
@@ -307,7 +488,7 @@ void proof_multi(void) {
     int coset_scale = 3, coset_len = (1 << coset_scale);
     fr_t y[coset_len];
 
-    uint64_t secrets_len = poly_len > coset_len ? poly_len + 1 : coset_len + 1;
+    uint64_t secrets_len = poly_len > coset_len ? poly_len : coset_len;
     g1_t s1[secrets_len];
     g2_t s2[secrets_len];
 
@@ -319,38 +500,34 @@ void proof_multi(void) {
 
     // Initialise the secrets and data structures
     generate_trusted_setup(s1, s2, &secret, secrets_len);
-    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs1, 4)); // ln_2 of poly_len
-    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks1, s1, s2, secrets_len, &fs1));
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4)); // ln_2 of poly_len
+    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks, s1, s2, secrets_len, &fs));
 
     // Commit to the polynomial
-    TEST_CHECK(C_KZG_OK == commit_to_poly(&commitment, &p, &ks1));
-
-    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs2, coset_scale));
-    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks2, s1, s2, secrets_len, &fs2));
+    TEST_CHECK(C_KZG_OK == commit_to_poly(&commitment, &p, &ks));
 
     // Compute proof at the points [x * root_i] 0 <= i < coset_len
     fr_from_uint64(&x, 5431);
-    TEST_CHECK(C_KZG_OK == compute_proof_multi(&proof, &p, &x, coset_len, &ks2));
+    TEST_CHECK(C_KZG_OK == compute_proof_multi(&proof, &p, &x, coset_len, &ks));
 
     // y_i is the value of the polynomial at each x_i
+    uint64_t stride = secrets_len / coset_len;
     for (int i = 0; i < coset_len; i++) {
-        fr_mul(&tmp, &x, &ks2.fs->expanded_roots_of_unity[i]);
+        fr_mul(&tmp, &x, &fs.expanded_roots_of_unity[i * stride]);
         eval_poly(&y[i], &p, &tmp);
     }
 
     // Verify the proof that the (unknown) polynomial has value y_i at x_i
-    TEST_CHECK(C_KZG_OK == check_proof_multi(&result, &commitment, &proof, &x, y, coset_len, &ks2));
+    TEST_CHECK(C_KZG_OK == check_proof_multi(&result, &commitment, &proof, &x, y, coset_len, &ks));
     TEST_CHECK(true == result);
 
     // Change a value and check that the proof fails
     fr_add(y + coset_len / 2, y + coset_len / 2, &fr_one);
-    TEST_CHECK(C_KZG_OK == check_proof_multi(&result, &commitment, &proof, &x, y, coset_len, &ks2));
+    TEST_CHECK(C_KZG_OK == check_proof_multi(&result, &commitment, &proof, &x, y, coset_len, &ks));
     TEST_CHECK(false == result);
 
-    free_fft_settings(&fs1);
-    free_fft_settings(&fs2);
-    free_kzg_settings(&ks1);
-    free_kzg_settings(&ks2);
+    free_fft_settings(&fs);
+    free_kzg_settings(&ks);
     free_poly(&p);
 }
 
@@ -397,12 +574,164 @@ void commit_to_too_long_poly(void) {
     free_kzg_settings(&ks);
 }
 
+void commit_to_poly_lagrange(void) {
+    // Our polynomial: degree 15, 16 coefficients
+    uint64_t coeffs[] = {12, 2, 8, 4, 7, 9, 1337, 227, 3, 13, 13, 130, 13, 13111, 13, 12223};
+    int poly_len = sizeof coeffs / sizeof coeffs[0];
+    uint64_t secrets_len = poly_len;
+
+    FFTSettings fs;
+    KZGSettings ks;
+    g1_t s1[secrets_len];
+    g2_t s2[secrets_len];
+    poly p;
+    poly_l p_l;
+    g1_t commitment, commitment_l;
+
+    // Create the polynomial
+    new_poly(&p, poly_len);
+    for (int i = 0; i < poly_len; i++) {
+        fr_from_uint64(&p.coeffs[i], coeffs[i]);
+    }
+
+    // Initialise the secrets and data structures
+    generate_trusted_setup(s1, s2, &secret, secrets_len);
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4)); // ln_2 of poly_len
+    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks, s1, s2, secrets_len, &fs));
+
+    // Create Lagrange form
+    new_poly_l_from_poly(&p_l, &p, &ks);
+
+    // Compute commitments
+    TEST_CHECK(C_KZG_OK == commit_to_poly(&commitment, &p, &ks));
+    TEST_CHECK(C_KZG_OK == commit_to_poly_l(&commitment_l, &p_l, &ks));
+
+    // Check commitments are equal
+    TEST_CHECK(g1_equal(&commitment, &commitment_l));
+
+    free_fft_settings(&fs);
+    free_kzg_settings(&ks);
+    free_poly(&p);
+    free_poly_l(&p_l);
+}
+
+void poly_eval_l_check(void) {
+    uint64_t n = 10;
+    fr_t actual, expected;
+    poly p;
+    new_poly(&p, n);
+    for (uint64_t i = 0; i < n; i++) {
+        fr_from_uint64(&p.coeffs[i], i + 1);
+    }
+    fr_t x;
+    fr_from_uint64(&x, 39);
+    // x = fr_one;
+    eval_poly(&expected, &p, &x);
+
+    poly_l p_l;
+    FFTSettings fs;
+    KZGSettings ks;
+    uint64_t secrets_len = 16;
+    g1_t s1[secrets_len];
+    g2_t s2[secrets_len];
+
+    generate_trusted_setup(s1, s2, &secret, secrets_len);
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4)); // log_2(secrets_len)
+    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks, s1, s2, secrets_len, &fs));
+
+    TEST_CHECK(C_KZG_OK == new_poly_l_from_poly(&p_l, &p, &ks));
+
+    eval_poly_l(&actual, &p_l, &x, &fs);
+
+    TEST_CHECK(fr_equal(&expected, &actual));
+
+    free_fft_settings(&fs);
+    free_kzg_settings(&ks);
+    free_poly(&p);
+    free_poly_l(&p_l);
+}
+
+void eval_poly_l_at_first_root_of_unity(void) {
+    uint64_t n = 10;
+    fr_t actual, expected;
+    poly p;
+    new_poly(&p, n);
+    for (uint64_t i = 0; i < n; i++) {
+        fr_from_uint64(&p.coeffs[i], i + 2);
+    }
+
+    poly_l p_l;
+    FFTSettings fs;
+    KZGSettings ks;
+    uint64_t secrets_len = 16;
+    g1_t s1[secrets_len];
+    g2_t s2[secrets_len];
+
+    generate_trusted_setup(s1, s2, &secret, secrets_len);
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4)); // log_2(secrets_len)
+    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks, s1, s2, secrets_len, &fs));
+
+    eval_poly(&expected, &p, &fs.expanded_roots_of_unity[0]);
+
+    TEST_CHECK(C_KZG_OK == new_poly_l_from_poly(&p_l, &p, &ks));
+
+    eval_poly_l(&actual, &p_l, &fs.expanded_roots_of_unity[0], &fs);
+
+    TEST_CHECK(fr_equal(&expected, &actual));
+
+    free_fft_settings(&fs);
+    free_kzg_settings(&ks);
+    free_poly(&p);
+    free_poly_l(&p_l);
+}
+
+void eval_poly_l_at_another_root_of_unity(void) {
+    uint64_t n = 13;
+    fr_t actual, expected;
+    poly p;
+    new_poly(&p, n);
+    for (uint64_t i = 0; i < n; i++) {
+        fr_from_uint64(&p.coeffs[i], 2 * n - i);
+    }
+
+    poly_l p_l;
+    FFTSettings fs;
+    KZGSettings ks;
+    uint64_t secrets_len = 16;
+    g1_t s1[secrets_len];
+    g2_t s2[secrets_len];
+
+    generate_trusted_setup(s1, s2, &secret, secrets_len);
+    TEST_CHECK(C_KZG_OK == new_fft_settings(&fs, 4)); // log_2(secrets_len)
+    TEST_CHECK(C_KZG_OK == new_kzg_settings(&ks, s1, s2, secrets_len, &fs));
+
+    eval_poly(&expected, &p, &fs.expanded_roots_of_unity[5]);
+
+    TEST_CHECK(C_KZG_OK == new_poly_l_from_poly(&p_l, &p, &ks));
+
+    eval_poly_l(&actual, &p_l, &fs.expanded_roots_of_unity[5], &fs);
+
+    TEST_CHECK(fr_equal(&expected, &actual));
+
+    free_fft_settings(&fs);
+    free_kzg_settings(&ks);
+    free_poly(&p);
+    free_poly_l(&p_l);
+}
+
+
 TEST_LIST = {
     {"KZG_PROOFS_TEST", title},
+    {"poly_eval_l_check", poly_eval_l_check},
+    {"eval_poly_l_at_first_root_of_unity", eval_poly_l_at_first_root_of_unity},
+    {"eval_poly_l_at_another_root_of_unity", eval_poly_l_at_another_root_of_unity},
     {"proof_single", proof_single},
+    {"proof_single_l", proof_single_l},
+    {"proof_single_l_at_root", proof_single_l_at_root},
     {"proof_multi", proof_multi},
     {"commit_to_nil_poly", commit_to_nil_poly},
     {"commit_to_too_long_poly", commit_to_too_long_poly},
+    {"commit_to_poly_lagrange", commit_to_poly_lagrange},
     {NULL, NULL} /* zero record marks the end of the list */
 };
 
diff --git a/src/poly.c b/src/poly.c
index 9d06607..9b10b95 100644
--- a/src/poly.c
+++ b/src/poly.c
@@ -102,6 +102,49 @@ void eval_poly(fr_t *out, const poly *p, const fr_t *x) {
     }
 }
 
+/**
+ * Evaluate a polynomial in Lagrange form over the finite field at a point.
+ *
+ * @param[out] out The value of the polynomial at the point @p x
+ * @param[in]  p   The polynomial in Lagrange form
+ * @param[in]  x   The x-coordinate to be evaluated
+ */
+C_KZG_RET eval_poly_l(fr_t *out, const poly_l *p, const fr_t *x, const FFTSettings *fs) {
+  fr_t tmp, *inverses_in, *inverses;
+  uint64_t i;
+  const uint64_t stride = fs->max_width / p->length;
+
+  TRY(new_fr_array(&inverses_in, p->length));
+  TRY(new_fr_array(&inverses, p->length));
+  for (i = 0; i < p->length; i++) {
+    if (fr_equal(x, &fs->expanded_roots_of_unity[i * stride])) {
+      *out = p->values[i];
+      free(inverses_in);
+      free(inverses);
+      return C_KZG_OK;
+    }
+    fr_sub(&inverses_in[i], x, &fs->expanded_roots_of_unity[i * stride]);
+  }
+  TRY(fr_batch_inv(inverses, inverses_in, p->length));
+
+  *out = fr_zero;
+  for (i = 0; i < p->length; i++) {
+    fr_mul(&tmp, &inverses[i], &fs->expanded_roots_of_unity[i * stride]);
+    fr_mul(&tmp, &tmp, &p->values[i]);
+    fr_add(out, out, &tmp);
+  }
+  fr_from_uint64(&tmp, p->length);
+  fr_div(out, out, &tmp);
+  fr_pow(&tmp, x, p->length);
+  fr_sub(&tmp, &tmp, &fr_one);
+  fr_mul(out, out, &tmp);
+
+  free(inverses_in);
+  free(inverses);
+
+  return C_KZG_OK;
+}
+
 /**
  * Polynomial division in the finite field via long division.
  *
@@ -541,6 +584,21 @@ C_KZG_RET new_poly(poly *out, uint64_t length) {
     return new_fr_array(&out->coeffs, length);
 }
 
+/**
+ * Initialise an empty polynomial in Lagrange form of the given size.
+ *
+ * @remark This allocates space for the Lagrange values that must be later reclaimed by calling #free_poly_l.
+ *
+ * @param[out] out    The initialised polynomial structure
+ * @param[in]  length The number of coefficients required, which is one more than the polynomial's degree
+ * @retval C_CZK_OK      All is well
+ * @retval C_CZK_MALLOC  Memory allocation failed
+ */
+C_KZG_RET new_poly_l(poly_l *out, uint64_t length) {
+    out->length = length;
+    return new_fr_array(&out->values, length);
+}
+
 /**
  * Initialise a polynomial of the given size with the given coefficients.
  *
@@ -562,6 +620,39 @@ C_KZG_RET new_poly_with_coeffs(poly *out, const fr_t *coeffs, uint64_t length) {
     return C_KZG_OK;
 }
 
+/**
+ * Initialise a polynomial in Lagrange form from the given polynomial in coefficient form.
+ *
+ * @remark This allocates space for the Lagrange values that must be later reclaimed by calling #free_poly_l.
+ *
+ * @param[out] out The initialised Lagrange polynomial structure
+ * @param[in]  in  The polynomial of which to compute the Lagrange form
+ * @param[in]  ks  The settings containing the roots of unity to use for DFT
+ * @retval C_CZK_OK      All is well
+ * @retval C_CZK_BADARGS Invalid settings, e.g., fft max_width too low for this polynomial
+ * @retval C_CZK_MALLOC  Memory allocation failed
+ */
+C_KZG_RET new_poly_l_from_poly(poly_l *out, const poly *in, const KZGSettings *ks) {
+  TRY(new_poly_l(out, ks->length));
+  if (out->length <= in->length) {
+    return fft_fr(out->values, in->coeffs, false, out->length, ks->fs);
+  }
+  else {
+    int i;
+    fr_t *coeffs;
+    TRY(new_fr_array(&coeffs, out->length));
+    for (i = 0; i < in->length; i++) {
+      coeffs[i] = in->coeffs[i];
+    }
+    for (; i < out->length; i++) {
+      coeffs[i] = fr_zero;
+    }
+    TRY(fft_fr(out->values, coeffs, false, out->length, ks->fs));
+    free(coeffs);
+    return C_KZG_OK;
+  }
+}
+
 /**
  * Reclaim the memory used by a polynomial.
  *
@@ -576,6 +667,20 @@ void free_poly(poly *p) {
     }
 }
 
+/**
+ * Reclaim the memory used by a polynomial in Lagrange form.
+ *
+ * @remark To avoid memory leaks, this must be called for polynomials initialised with #new_poly_l or
+ * #new_poly_l_from_poly after use.
+ *
+ * @param[in,out] p The polynomial
+ */
+void free_poly_l(poly_l *p) {
+    if (p->values != NULL) {
+        free(p->values);
+    }
+}
+
 #ifdef KZGTEST
 
 #include "../inc/acutest.h"
@@ -978,4 +1083,4 @@ TEST_LIST = {
     {NULL, NULL} /* zero record marks the end of the list */
 };
 
-#endif // KZGTEST
\ No newline at end of file
+#endif // KZGTEST
diff --git a/src/poly_barycentric_bench.c b/src/poly_barycentric_bench.c
new file mode 100644
index 0000000..33a70b4
--- /dev/null
+++ b/src/poly_barycentric_bench.c
@@ -0,0 +1,86 @@
+/*
+ * Copyright 2021 Benjamin Edgington
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <stdlib.h> // malloc(), free(), atoi()
+#include <stdio.h>  // printf()
+#include <assert.h> // assert()
+#include <unistd.h> // EXIT_SUCCESS/FAILURE
+#include "bench_util.h"
+#include "test_util.h"
+#include "c_kzg.h"
+
+// Run the benchmark for `max_seconds` and return the time per iteration in nanoseconds.
+long run_bench(int scale, int max_seconds) {
+    timespec_t t0, t1;
+    unsigned long total_time = 0, nits = 0;
+
+    uint64_t width = (uint64_t)1 << scale;
+
+    FFTSettings fs;
+
+    new_fft_settings(&fs, scale);
+
+    poly_l p;
+    fr_t x, y;
+    new_poly_l(&p, width);
+
+    for (int i = 0; i < width; i++) {
+        p.values[i] = rand_fr();
+    }
+    x = rand_fr();
+
+    while (total_time < max_seconds * NANO) {
+        clock_gettime(CLOCK_REALTIME, &t0);
+
+        eval_poly_l(&y, &p, &x, &fs);
+
+        clock_gettime(CLOCK_REALTIME, &t1);
+        nits++;
+        total_time += tdiff(t0, t1);
+    }
+
+    free_poly_l(&p);
+    free_fft_settings(&fs);
+
+    return total_time / nits;
+}
+
+int main(int argc, char *argv[]) {
+    int nsec = 0;
+
+    switch (argc) {
+    case 1:
+        nsec = NSEC;
+        break;
+    case 2:
+        nsec = atoi(argv[1]);
+        break;
+    default:
+        break;
+    };
+
+    if (nsec == 0) {
+        printf("Usage: %s [test time in seconds > 0]\n", argv[0]);
+        exit(EXIT_FAILURE);
+    }
+
+    printf("*** Benchmarking Barycentric Formula (eval_poly_l), %d second%s per test.\n", nsec, nsec == 1 ? "" : "s");
+    for (int scale = 6; scale <= 15; scale++) {
+        printf("eval_poly_l/scale_%d %lu ns/op\n", scale, run_bench(scale, nsec));
+    }
+
+    return EXIT_SUCCESS;
+}
diff --git a/src/utility.c b/src/utility.c
index 916871a..b86e9a9 100644
--- a/src/utility.c
+++ b/src/utility.c
@@ -23,6 +23,7 @@
 #include <string.h> // memcpy()
 #include "control.h"
 #include "utility.h"
+#include "c_kzg_alloc.h"
 
 /**
  * Utility function to test whether the argument is a power of two.
@@ -154,6 +155,39 @@ C_KZG_RET reverse_bit_order(void *values, size_t size, uint64_t n) {
     return C_KZG_OK;
 }
 
+/**
+ * Montgomery batch inversion in finite field
+ *
+ * @param[out] out The inverses of @p a, length @p len
+ * @param[in]  a   A vector of field elements, length @p len
+ * @param[in]  len Length
+ */
+C_KZG_RET fr_batch_inv(fr_t *out, const fr_t *a, size_t len) {
+    fr_t *prod;
+    fr_t inv;
+    size_t i;
+
+    TRY(new_fr_array(&prod, len));
+
+    prod[0] = a[0];
+
+    for(i = 1; i < len; i++) {
+        fr_mul(&prod[i], &a[i], &prod[i - 1]);
+    }
+
+    blst_fr_eucl_inverse(&inv, &prod[len - 1]);
+
+    for(i = len - 1; i > 0; i--) {
+        fr_mul(&out[i], &inv, &prod[i - 1]);
+        fr_mul(&inv, &a[i], &inv);
+    }
+    out[0] = inv;
+
+    free(prod);
+
+    return C_KZG_OK;
+}
+
 #ifdef KZGTEST
 
 #include "../inc/acutest.h"
@@ -185,6 +219,28 @@ void is_power_of_two_works(void) {
     TEST_CHECK(false == is_power_of_two(1234567));
 }
 
+void test_batch_inv(void) {
+    fr_t *inputs, *actual, *expected;
+    int i;
+
+    TEST_CHECK(C_KZG_OK == new_fr_array(&inputs, 32));
+    TEST_CHECK(C_KZG_OK == new_fr_array(&actual, 32));
+    TEST_CHECK(C_KZG_OK == new_fr_array(&expected, 32));
+
+    for (i = 0; i < 32; i++) {
+        inputs[i] = rand_fr();
+        fr_inv(&expected[i], &inputs[i]);
+    }
+    fr_batch_inv(actual, inputs, 32);
+    for (i = 0; i < 32; i++) {
+        TEST_CHECK(fr_equal(&expected[i], &actual[i]));
+    }
+
+    free(inputs);
+    free(actual);
+    free(expected);
+}
+
 void test_log2_pow2(void) {
     int actual, expected;
     for (int i = 0; i < 32; i++) {
@@ -309,6 +365,7 @@ void test_reverse_bit_order_fr_large(void) {
 
 TEST_LIST = {
     {"UTILITY_TEST", title},
+    {"test_batch_inv", test_batch_inv},
     {"is_power_of_two_works", is_power_of_two_works},
     {"test_log2_pow2", test_log2_pow2},
     {"test_next_power_of_two_powers", test_next_power_of_two_powers},
@@ -322,4 +379,4 @@ TEST_LIST = {
     {NULL, NULL} /* zero record marks the end of the list */
 };
 
-#endif // KZGTEST
\ No newline at end of file
+#endif // KZGTEST
diff --git a/src/utility.h b/src/utility.h
index e5af837..9803cb6 100644
--- a/src/utility.h
+++ b/src/utility.h
@@ -49,3 +49,4 @@ uint64_t next_power_of_two(uint64_t v);
 uint32_t reverse_bits(uint32_t a);
 uint32_t reverse_bits_limited(uint32_t n, uint32_t value);
 C_KZG_RET reverse_bit_order(void *values, size_t size, uint64_t n);
+C_KZG_RET fr_batch_inv(fr_t *out, const fr_t *a, size_t len);
diff --git a/src/zero_poly.c b/src/zero_poly.c
index 7370d23..8c405bb 100644
--- a/src/zero_poly.c
+++ b/src/zero_poly.c
@@ -516,11 +516,10 @@ void zero_poly_random(void) {
             TEST_CHECK(len_missing + 1 == zero_poly.length);
             TEST_MSG("ZeroPolyLen: expected %d, got %lu", len_missing + 1, zero_poly.length);
 
-            int ret = 0;
             for (int i = 0; i < len_missing; i++) {
                 fr_t out;
                 eval_poly(&out, &zero_poly, &fs.expanded_roots_of_unity[missing[i]]);
-                ret = TEST_CHECK(fr_is_zero(&out));
+                TEST_CHECK(fr_is_zero(&out));
                 TEST_MSG("Failed for missing[%d] = %lu", i, missing[i]);
             }
 
@@ -567,11 +566,10 @@ void zero_poly_all_but_one(void) {
     TEST_CHECK(len_missing + 1 == zero_poly.length);
     TEST_MSG("ZeroPolyLen: expected %d, got %lu", len_missing + 1, zero_poly.length);
 
-    int ret = 0;
     for (int i = 0; i < len_missing; i++) {
         fr_t out;
         eval_poly(&out, &zero_poly, &fs.expanded_roots_of_unity[missing[i]]);
-        ret = TEST_CHECK(fr_is_zero(&out));
+        TEST_CHECK(fr_is_zero(&out));
         TEST_MSG("Failed for missing[%d] = %lu", i, missing[i]);
     }
 
@@ -615,11 +613,10 @@ void zero_poly_252(void) {
     TEST_CHECK(len_missing + 1 == zero_poly.length);
     TEST_MSG("ZeroPolyLen: expected %d, got %lu", len_missing + 1, zero_poly.length);
 
-    int ret = 0;
     for (int i = 0; i < len_missing; i++) {
         fr_t out;
         eval_poly(&out, &zero_poly, &fs.expanded_roots_of_unity[missing[i]]);
-        ret = TEST_CHECK(fr_is_zero(&out));
+        TEST_CHECK(fr_is_zero(&out));
         TEST_MSG("Failed for missing[%d] = %lu", i, missing[i]);
     }