diff --git a/ALFI/ALFI/spline/polyeqv.h b/ALFI/ALFI/spline/polyeqv.h
index 389a7c8..c2fce32 100644
--- a/ALFI/ALFI/spline/polyeqv.h
+++ b/ALFI/ALFI/spline/polyeqv.h
@@ -90,33 +90,17 @@ namespace alfi::spline {
 						case PolynomialType::IMPROVED_LAGRANGE: P = poly::imp_lagrange(X, Y); break;
 						case PolynomialType::NEWTON: P = poly::newton(X, Y); break;
 					}
-
-					const static auto binomials = [](SizeT m) {
-						Container<Container<Number>> C(m + 1);
-						for (SizeT i = 0; i <= m; ++i) {
-							C[i].resize(i + 1);
-							C[i][0] = C[i][i] = 1;
-							for (SizeT j = 1; j <= i / 2; ++j) {
-								C[i][j] = C[i][i-j] = C[i-1][j-1] + C[i-1][j];
-							}
-						}
-						return C;
-					};
-
-					const auto C = binomials(n - 1);
-
 #if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
 #pragma omp parallel for
 #endif
 					for (SizeT i = 0; i < n - 1; ++i) {
-						coeffs[i*n] = P[0];
-						for (SizeT k = 1; k < n - 1; ++k) {
-							Number r = 0;
-							for (SizeT j = 0; j < k; ++j) {
-								r += ((((k - j) & 1) == 0) ? 1 : -1) * coeffs[i*n+j] * C[n-1-j][k-j];
-								r *= X[i];
+						for (SizeT j = 0; j < n; ++j) {
+							coeffs[i*n+j] = P[j];
+						}
+						for (SizeT k = 0; k < n - 1; ++k) {
+							for (SizeT j = 1; j < n - k; ++j) {
+								coeffs[i*n+j] += coeffs[i*n+j-1] * X[i];
 							}
-							coeffs[i*n+k] = P[k] - r;
 						}
 						coeffs[i*n+n-1] = Y[i];
 					}
@@ -128,12 +112,17 @@ namespace alfi::spline {
 
 		PolyEqvSpline() = default;
 
+		template <typename ContainerXType>
+		PolyEqvSpline(ContainerXType&& X, const Container<Number>& Y, OptimizationType optimization_type)
+			: PolyEqvSpline(std::forward<ContainerXType>(X), Y, PolynomialType::Default, optimization_type) {}
+
 		template <typename ContainerXType>
 		PolyEqvSpline(ContainerXType&& X,
 					  const Container<Number>& Y,
 					  PolynomialType polynomial_type = PolynomialType::Default,
-					  OptimizationType optimization_type = OptimizationType::Default) {
-			construct(std::forward<ContainerXType>(X), Y, polynomial_type, optimization_type);
+					  OptimizationType optimization_type = OptimizationType::Default,
+					  EvaluationType evaluation_type = EvaluationType::Default) {
+			construct(std::forward<ContainerXType>(X), Y, polynomial_type, optimization_type, evaluation_type);
 		}
 
 		PolyEqvSpline(const PolyEqvSpline& other) = default;
diff --git a/tests/spline/test_polyeqv.cpp b/tests/spline/test_polyeqv.cpp
index 62d48f9..928a0c7 100644
--- a/tests/spline/test_polyeqv.cpp
+++ b/tests/spline/test_polyeqv.cpp
@@ -6,7 +6,7 @@ const auto test_data_path = TEST_DATA_DIR "/poly/poly.toml";
 
 const auto test_data = toml::parse_file(test_data_path);
 
-void test_polyeqv_spline(double epsilon) {
+void test_polyeqv_spline(double epsilon_accuracy, double epsilon_speed) {
 	const auto& test_cases = test_data["test_cases"].ref<toml::array>();
 
 	test_cases.for_each([&](const toml::table& test_case) {
@@ -15,12 +15,13 @@ void test_polyeqv_spline(double epsilon) {
 		const auto& xx = to_vector<double>(test_case["xx"].ref<toml::array>());
 		const auto& yy = to_vector<double>(test_case["yy"].ref<toml::array>());
 
-		expect_eq(alfi::spline::PolyEqvSpline<>(X, Y).eval(xx), yy, epsilon);
+		expect_eq(alfi::spline::PolyEqvSpline<>(X, Y, alfi::spline::PolyEqvSpline<>::OptimizationType::ACCURACY).eval(xx), yy, epsilon_accuracy);
+		expect_eq(alfi::spline::PolyEqvSpline<>(X, Y, alfi::spline::PolyEqvSpline<>::OptimizationType::SPEED).eval(xx), yy, epsilon_speed);
 	});
 }
 
 TEST(PolyEqvSplineTest, PolynomialData) {
-	test_polyeqv_spline(1e-10);
+	test_polyeqv_spline(1e-10, 1e-2);
 }
 
 TEST(PolyEqvSplineTest, General) {