diff --git a/tests/core/pyspec/eth2spec/test/eip4844/unittests/polynomial_commitments/test_polynomial_commitments.py b/tests/core/pyspec/eth2spec/test/eip4844/unittests/polynomial_commitments/test_polynomial_commitments.py
index 75a4e822a..b503bba57 100644
--- a/tests/core/pyspec/eth2spec/test/eip4844/unittests/polynomial_commitments/test_polynomial_commitments.py
+++ b/tests/core/pyspec/eth2spec/test/eip4844/unittests/polynomial_commitments/test_polynomial_commitments.py
@@ -1,3 +1,5 @@
+import random
+
 from eth2spec.test.context import (
     spec_state_test,
     with_eip4844_and_later,
@@ -22,6 +24,37 @@ def test_verify_kzg_proof(spec, state):
     assert spec.verify_kzg_proof_impl(commitment, x, y, proof)
 
 
+@with_eip4844_and_later
+@spec_state_test
+def test_barycentric_outside_domain(spec, state):
+    """
+    Test barycentric formula correctness by using it to evaluate a polynomial at a bunch of points outside its domain
+    (the roots of unity).
+
+    Then make sure that we would get the same result if we evaluated it from coefficient form without using the
+    barycentric formula
+    """
+    rng = random.Random(5566)
+    poly_coeff, poly_eval = get_poly_in_both_forms(spec)
+    roots_of_unity_brp = spec.bit_reversal_permutation(spec.ROOTS_OF_UNITY)
+
+    assert len(poly_coeff) == len(poly_eval) == len(roots_of_unity_brp)
+    n_samples = 12
+
+    for i in range(n_samples):
+        # Get a random evaluation point
+        z = rng.randint(0, spec.BLS_MODULUS - 1)
+
+        # Get p(z) by evaluating poly in coefficient form
+        p_z_coeff = eval_poly_in_coeff_form(spec, poly_coeff, z)
+
+        # Get p(z) by evaluating poly in evaluation form
+        p_z_eval = spec.evaluate_polynomial_in_evaluation_form(poly_eval, z)
+
+        # Both evaluations should agree
+        assert p_z_coeff == p_z_eval
+
+
 @with_eip4844_and_later
 @spec_state_test
 def test_barycentric_within_domain(spec, state):