diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 8308ec4..8f4aa22 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -17,9 +17,7 @@ endmacro()
 
 add_test_executable(test_dist dist/test_dist.cpp)
 
-add_test_executable(test_lagrange poly/test_lagrange.cpp poly/test_case.h)
-add_test_executable(test_imp_lagrange poly/test_imp_lagrange.cpp poly/test_case.h)
-add_test_executable(test_newton poly/test_newton.cpp poly/test_case.h)
+add_test_executable(test_poly poly/test_poly.cpp)
 
 add_test_executable(test_barycentric misc/test_barycentric.cpp)
 
diff --git a/tests/poly/test_case.h b/tests/poly/test_case.h
deleted file mode 100644
index 6e94848..0000000
--- a/tests/poly/test_case.h
+++ /dev/null
@@ -1,47 +0,0 @@
-#pragma once
-
-#include <vector>
-
-#include <ALFI/poly.h>
-
-#include "../test_utils.h"
-
-inline void test_case(const std::vector<double>& X, const std::vector<double>& Y,
-					  const std::vector<double>& expected_coeffs,
-					  const std::vector<double>& xx, const std::vector<double>& yy,
-					  double epsilon = 1e-9) {
-	const auto P = POLYNOMIAL(X, Y);
-	expect_eq(P, expected_coeffs);
-	const auto values = alfi::poly::val(P, xx);
-	expect_eq(values, yy, epsilon);
-	const auto vals = POLYNOMIAL_VALS(X, Y, xx);
-	expect_eq(vals, yy, epsilon);
-}
-
-inline void test_case(const std::vector<double>& X, const std::vector<double>& Y,
-					  const std::vector<double>& xx, const std::vector<double>& yy,
-					  double epsilon = 1e-9) {
-	const auto P = POLYNOMIAL(X, Y);
-	const auto values = alfi::poly::val(P, xx);
-	expect_eq(values, yy, epsilon);
-	const auto vals = POLYNOMIAL_VALS(X, Y, xx);
-	expect_eq(vals, yy, epsilon);
-}
-
-inline void test_case(const std::vector<double>& X, const std::vector<double>& Y,
-					  const std::vector<double>& expected_coeffs, double epsilon = 1e-9) {
-	test_case(X, Y, expected_coeffs, X, Y, epsilon);
-}
-
-inline void test_case(const std::vector<double>& X, const std::vector<double>& Y, double epsilon = 1e-9) {
-	test_case(X, Y, X, Y, epsilon);
-}
-
-inline void test_case(const std::vector<double>& X, const std::vector<double>& xx,
-					  const std::function<double(double)>& f, double epsilon = 1e-9) {
-	std::vector<double> Y(X.size());
-	std::transform(X.begin(), X.end(), Y.begin(), f);
-	std::vector<double> yy(xx.size());
-	std::transform(xx.begin(), xx.end(), yy.begin(), f);
-	test_case(X, Y, xx, yy, epsilon);
-}
\ No newline at end of file
diff --git a/tests/poly/test_imp_lagrange.cpp b/tests/poly/test_imp_lagrange.cpp
deleted file mode 100644
index 153052b..0000000
--- a/tests/poly/test_imp_lagrange.cpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#include <gtest/gtest.h>
-
-#include <cmath>
-#include <vector>
-
-#include <ALFI/poly.h>
-#include <ALFI/dist.h>
-
-#define POLYNOMIAL alfi::poly::imp_lagrange
-#define POLYNOMIAL_VALS alfi::poly::imp_lagrange_vals
-#include "test_case.h"
-
-TEST(ImprovedLagrangeTest, Basic) {
-	test_case({0, 1, 2}, {0, 0, 0});
-	test_case({0, 1, 2}, {1, 1, 1});
-	test_case({0, 1, 2}, {1, 2, 3});
-	test_case({0, 1, 2, 3, 5}, {3, 1, 2, 5, -1});
-}
-
-TEST(ImprovedLagrangeTest, Functions) {
-	test_case(
-		alfi::dist::uniform(3, 0.0, 20.0),
-		alfi::dist::uniform(10, 0.0, 20.0),
-		[](double x) { return x * x; });
-	test_case(
-		alfi::dist::chebyshev(21, -2.0, 2.0),
-		alfi::dist::uniform(10, -2.0, 2.0),
-		[](double x) { return std::exp(x); },
-		1e-8);
-	test_case(
-		alfi::dist::chebyshev_2(31, -2*M_PI, 2*M_PI),
-		alfi::dist::uniform(10, -2*M_PI, 2*M_PI),
-		[](double x) { return std::sin(x) + std::cos(x); },
-		1e-5);
-}
\ No newline at end of file
diff --git a/tests/poly/test_lagrange.cpp b/tests/poly/test_lagrange.cpp
deleted file mode 100644
index 497aceb..0000000
--- a/tests/poly/test_lagrange.cpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#include <gtest/gtest.h>
-
-#include <cmath>
-#include <vector>
-
-#include <ALFI/poly.h>
-#include <ALFI/dist.h>
-
-#define POLYNOMIAL alfi::poly::lagrange
-#define POLYNOMIAL_VALS alfi::poly::lagrange_vals
-#include "test_case.h"
-
-TEST(LagrangeTest, Basic) {
-	test_case({0, 1, 2}, {0, 0, 0});
-	test_case({0, 1, 2}, {1, 1, 1});
-	test_case({0, 1, 2}, {1, 2, 3});
-	test_case({0, 1, 2, 3, 5}, {3, 1, 2, 5, -1});
-}
-
-TEST(LagrangeTest, Functions) {
-	test_case(
-		alfi::dist::uniform(3, 0.0, 20.0),
-		alfi::dist::uniform(10, 0.0, 20.0),
-		[](double x) { return x * x; });
-	test_case(
-		alfi::dist::chebyshev(21, -2.0, 2.0),
-		alfi::dist::uniform(10, -2.0, 2.0),
-		[](double x) { return std::exp(x); },
-		1e-6);
-	test_case(
-		alfi::dist::chebyshev_2(31, -2*M_PI, 2*M_PI),
-		alfi::dist::uniform(10, -2*M_PI, 2*M_PI),
-		[](double x) { return std::sin(x) + std::cos(x); },
-		2e-2);
-}
\ No newline at end of file
diff --git a/tests/poly/test_newton.cpp b/tests/poly/test_newton.cpp
deleted file mode 100644
index 54819cb..0000000
--- a/tests/poly/test_newton.cpp
+++ /dev/null
@@ -1,35 +0,0 @@
-#include <gtest/gtest.h>
-
-#include <cmath>
-#include <vector>
-
-#include <ALFI/poly.h>
-#include <ALFI/dist.h>
-
-#define POLYNOMIAL alfi::poly::newton
-#define POLYNOMIAL_VALS alfi::poly::newton_vals
-#include "test_case.h"
-
-TEST(NewtonTest, Basic) {
-	test_case({0, 1, 2}, {0, 0, 0});
-	test_case({0, 1, 2}, {1, 1, 1});
-	test_case({0, 1, 2}, {1, 2, 3});
-	test_case({0, 1, 2, 3, 5}, {3, 1, 2, 5, -1});
-}
-
-TEST(NewtonTest, Functions) {
-	test_case(
-		alfi::dist::uniform(3, 0.0, 20.0),
-		alfi::dist::uniform(10, 0.0, 20.0),
-		[](double x) { return x * x; });
-	test_case(
-		alfi::dist::chebyshev(21, -2.0, 2.0),
-		alfi::dist::uniform(10, -2.0, 2.0),
-		[](double x) { return std::exp(x); },
-		1e-8);
-	test_case(
-		alfi::dist::chebyshev_2(31, -2*M_PI, 2*M_PI),
-		alfi::dist::uniform(10, -2*M_PI, 2*M_PI),
-		[](double x) { return std::sin(x) + std::cos(x); },
-		1e-6);
-}
\ No newline at end of file
diff --git a/tests/poly/test_poly.cpp b/tests/poly/test_poly.cpp
new file mode 100644
index 0000000..798cc99
--- /dev/null
+++ b/tests/poly/test_poly.cpp
@@ -0,0 +1,44 @@
+#include <ALFI/poly.h>
+
+#include "../test_utils.h"
+
+const auto test_data_path = TEST_DATA_DIR "/poly/poly.toml";
+
+const auto test_data = toml::parse_file(test_data_path);
+
+void test_polynomial(const auto& function_coeffs, const auto& function_vals,
+					 double epsilon_coeffs, double epsilon_values, double epsilon_vals) {
+	const auto& test_cases = test_data["test_cases"].ref<toml::array>();
+
+	test_cases.for_each([&](const toml::table& test_case) {
+		const auto& X = to_vector<double>(test_case["X"].ref<toml::array>());
+		const auto& Y = to_vector<double>(test_case["Y"].ref<toml::array>());
+		const auto& coeffs = to_vector<double>(test_case["coeffs"].ref<toml::array>());
+		const auto& xx = to_vector<double>(test_case["xx"].ref<toml::array>());
+		const auto& yy = to_vector<double>(test_case["yy"].ref<toml::array>());
+
+		const auto P = function_coeffs(X, Y);
+		expect_eq(P, coeffs, epsilon_coeffs);
+		const auto values = alfi::poly::val(P, xx);
+		expect_eq(values, yy, epsilon_values);
+		const auto vals = function_vals(X, Y, xx);
+		expect_eq(vals, yy, epsilon_vals);
+	});
+}
+
+// need this workaround because of the fourth parameter of `alfi::poly::imp_lagrange_vals`
+static const auto& imp_lagrange_vals = [](const auto& X, const auto& Y, const auto& xx) {
+	return alfi::poly::imp_lagrange_vals<>(X, Y, xx);
+};
+
+TEST(PolynomialsTest, Lagrange) {
+	test_polynomial(alfi::poly::lagrange<>, alfi::poly::lagrange_vals<>, 1e-7, 1e-1, 1e-11);
+}
+
+TEST(PolynomialsTest, ImprovedLagrange) {
+	test_polynomial(alfi::poly::imp_lagrange<>, imp_lagrange_vals, 1e-4, 1e-4, 1e-11);
+}
+
+TEST(PolynomialsTest, Newton) {
+	test_polynomial(alfi::poly::newton<>, alfi::poly::newton_vals<>, 1e-8, 1e-4, 1e-5);
+}
\ No newline at end of file
diff --git a/tests/test_data b/tests/test_data
index 615bbf9..b743ee0 160000
--- a/tests/test_data
+++ b/tests/test_data
@@ -1 +1 @@
-Subproject commit 615bbf9b03840562e69db706b205ba8c7d8c128e
+Subproject commit b743ee09d375c7202c638cc0c17d4fdb4d5cf6c6