diff --git a/ALFI/ALFI/dist.h b/ALFI/ALFI/dist.h
index d6ba31b..b2afd30 100644
--- a/ALFI/ALFI/dist.h
+++ b/ALFI/ALFI/dist.h
@@ -34,7 +34,7 @@ namespace alfi::dist {
 	};
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> uniform(SizeT n, Number a, Number b) {
+	Container<Number> uniform(SizeT n, const Number& a, const Number& b) {
 		if (n == 1)
 			return {(a+b)/2};
 		Container<Number> points(n);
@@ -63,7 +63,7 @@ namespace alfi::dist {
 		@return a container with \p n points distributed on the segment `[a, b]` according to the transform function
 	*/
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> quadratic(SizeT n, Number a, Number b) {
+	Container<Number> quadratic(SizeT n, const Number& a, const Number& b) {
 		if (n == 1)
 			return {(a+b)/2};
 		Container<Number> points(n);
@@ -91,7 +91,7 @@ namespace alfi::dist {
 		@return a container with \p n points distributed on the segment `[a, b]` according to the transform function
 	*/
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> cubic(SizeT n, Number a, Number b) {
+	Container<Number> cubic(SizeT n, const Number& a, const Number& b) {
 		if (n == 1)
 			return {(a+b)/2};
 		Container<Number> points(n);
@@ -104,7 +104,7 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev(SizeT n, Number a, Number b) {
+	Container<Number> chebyshev(SizeT n, const Number& a, const Number& b) {
 		if (n == 1)
 			return {(a+b)/2};
 		Container<Number> points(n);
@@ -116,12 +116,12 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_stretched(SizeT n, Number a, Number b) {
+	Container<Number> chebyshev_stretched(SizeT n, const Number& a, const Number& b) {
 		return points::stretched<Number,Container>(chebyshev(n, a, b), a, b);
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_augmented(SizeT n, Number a, Number b) {
+	Container<Number> chebyshev_augmented(SizeT n, const Number& a, const Number& b) {
 		if (n == 0) {
 			return {};
 		}
@@ -137,7 +137,7 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_2(SizeT n, Number a, Number b) {
+	Container<Number> chebyshev_2(SizeT n, const Number& a, const Number& b) {
 		if (n == 1)
 			return {(a+b)/2};
 		Container<Number> points(n);
@@ -149,7 +149,7 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_3(SizeT n, Number a, Number b) {
+	Container<Number> chebyshev_3(SizeT n, const Number& a, const Number& b) {
 		Container<Number> points(n);
 		for (SizeT i = 0; i < n; ++i) {
 			const Number x = 1 - std::cos(M_PI * static_cast<Number>(2*i) / static_cast<Number>(2*n - 1));
@@ -159,12 +159,12 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_3_stretched(SizeT n, Number a, Number b) {
+	Container<Number> chebyshev_3_stretched(SizeT n, const Number& a, const Number& b) {
 		return points::stretched<Number,Container>(chebyshev_3(n, a, b), a, b);
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_4(SizeT n, Number a, Number b) {
+	Container<Number> chebyshev_4(SizeT n, const Number& a, const Number& b) {
 		Container<Number> points(n);
 		for (SizeT i = 0; i < n; ++i) {
 			const Number x = 1 - std::cos(M_PI * static_cast<Number>(2*i + 1) / static_cast<Number>(2*n - 1));
@@ -174,12 +174,12 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_4_stretched(SizeT n, Number a, Number b) {
+	Container<Number> chebyshev_4_stretched(SizeT n, const Number& a, const Number& b) {
 		return points::stretched<Number,Container>(chebyshev_4(n, a, b), a, b);
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_ellipse(SizeT n, Number a, Number b, Number ratio) {
+	Container<Number> chebyshev_ellipse(SizeT n, const Number& a, const Number& b, const Number& ratio) {
 		Container<Number> points(n);
 		for (SizeT i = 0; i < n / 2; ++i) {
 			const Number theta = M_PI * (2 * static_cast<Number>(i) + 1) / (2 * static_cast<Number>(n));
@@ -193,12 +193,12 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_ellipse_stretched(SizeT n, Number a, Number b, Number ratio) {
+	Container<Number> chebyshev_ellipse_stretched(SizeT n, const Number& a, const Number& b, const Number& ratio) {
 		return points::stretched<Number,Container>(chebyshev_ellipse(n, a, b, ratio), a, b);
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_ellipse_augmented(SizeT n, Number a, Number b, Number ratio) {
+	Container<Number> chebyshev_ellipse_augmented(SizeT n, const Number& a, const Number& b, const Number& ratio) {
 		if (n == 0) {
 			return {};
 		}
@@ -214,7 +214,7 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_ellipse_2(SizeT n, Number a, Number b, Number ratio) {
+	Container<Number> chebyshev_ellipse_2(SizeT n, const Number& a, const Number& b, const Number& ratio) {
 		Container<Number> points(n);
 		for (SizeT i = 0; i < n / 2; ++i) {
 			const Number theta = M_PI * static_cast<Number>(i) / (static_cast<Number>(n) - 1);
@@ -228,7 +228,7 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_ellipse_3(SizeT n, Number a, Number b, Number ratio) {
+	Container<Number> chebyshev_ellipse_3(SizeT n, const Number& a, const Number& b, const Number& ratio) {
 		Container<Number> points(n);
 		for (SizeT i = 0; i < n; ++i) {
 			const Number theta = M_PI * static_cast<Number>(2*i) / static_cast<Number>(2*n - 1);
@@ -239,12 +239,12 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_ellipse_3_stretched(SizeT n, Number a, Number b, Number ratio) {
+	Container<Number> chebyshev_ellipse_3_stretched(SizeT n, const Number& a, const Number& b, const Number& ratio) {
 		return points::stretched<Number,Container>(chebyshev_ellipse_3(n, a, b, ratio), a, b);
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_ellipse_4(SizeT n, Number a, Number b, Number ratio) {
+	Container<Number> chebyshev_ellipse_4(SizeT n, const Number& a, const Number& b, const Number& ratio) {
 		Container<Number> points(n);
 		for (SizeT i = 0; i < n; ++i) {
 			const Number theta = M_PI * static_cast<Number>(2*i + 1) / static_cast<Number>(2*n - 1);
@@ -255,7 +255,7 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> chebyshev_ellipse_4_stretched(SizeT n, Number a, Number b, Number ratio) {
+	Container<Number> chebyshev_ellipse_4_stretched(SizeT n, const Number& a, const Number& b, const Number& ratio) {
 		return points::stretched<Number,Container>(chebyshev_ellipse_4(n, a, b, ratio), a, b);
 	}
 
@@ -283,7 +283,7 @@ namespace alfi::dist {
 		@return a container with \p n points distributed on the interval `(a, b)` according to the transform function
 	*/
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> logistic(SizeT n, Number a, Number b, Number steepness) {
+	Container<Number> logistic(SizeT n, const Number& a, const Number& b, const Number& steepness) {
 		if (n == 1)
 			return {(a+b)/2};
 		Container<Number> points(n);
@@ -296,7 +296,7 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> logistic_stretched(SizeT n, Number a, Number b, Number steepness) {
+	Container<Number> logistic_stretched(SizeT n, const Number& a, const Number& b, const Number& steepness) {
 		if (n == 0)
 			return {};
 		if (n == 1)
@@ -338,7 +338,7 @@ namespace alfi::dist {
 		@return a container with \p n points distributed on the interval `(a, b)` according to the transform function
 	*/
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> erf(SizeT n, Number a, Number b, Number steepness) {
+	Container<Number> erf(SizeT n, const Number& a, const Number& b, const Number& steepness) {
 		if (n == 0)
 			return {};
 		if (n == 1)
@@ -353,12 +353,12 @@ namespace alfi::dist {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> erf_stretched(SizeT n, Number a, Number b, Number steepness) {
+	Container<Number> erf_stretched(SizeT n, const Number& a, const Number& b, const Number& steepness) {
 		return points::stretched<Number,Container>(erf(n, a, b, steepness), a, b);
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> of_type(Type type, SizeT n, Number a, Number b, Number parameter = NAN) {
+	Container<Number> of_type(Type type, SizeT n, const Number& a, const Number& b, const Number& parameter = NAN) {
 		switch (type) {
 		case Type::QUADRATIC:
 			return quadratic(n, a, b);
diff --git a/ALFI/ALFI/misc.h b/ALFI/ALFI/misc.h
index 831b040..37a3ded 100644
--- a/ALFI/ALFI/misc.h
+++ b/ALFI/ALFI/misc.h
@@ -13,7 +13,7 @@ namespace alfi::misc {
 			const Container<Number>& Y,
 			const Container<Number>& xx,
 			dist::Type dist_type = dist::Type::GENERAL,
-			Number epsilon = std::numeric_limits<Number>::epsilon()
+			const Number& epsilon = std::numeric_limits<Number>::epsilon()
 	) {
 		if (X.size() != Y.size()) {
 			std::cerr << "Error in function " << __FUNCTION__
diff --git a/ALFI/ALFI/poly.h b/ALFI/ALFI/poly.h
index 793fe92..b920402 100644
--- a/ALFI/ALFI/poly.h
+++ b/ALFI/ALFI/poly.h
@@ -9,7 +9,7 @@
 namespace alfi::poly {
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
 	std::enable_if_t<!traits::has_size<Number>::value, Number>
-	val(const Container<Number>& coeffs, Number x) {
+	val(const Container<Number>& coeffs, const Number& x) {
 		Number result = 0;
 		for (const Number& c : coeffs) {
 			result = result * x + c;
@@ -187,7 +187,7 @@ namespace alfi::poly {
 			const Container<Number>& X,
 			const Container<Number>& Y,
 			const Container<Number>& xx,
-			Number epsilon = std::numeric_limits<Number>::epsilon()
+			const Number& epsilon = std::numeric_limits<Number>::epsilon()
 	) {
 		const auto nn = xx.size();
 
diff --git a/ALFI/ALFI/ratf.h b/ALFI/ALFI/ratf.h
index c13102c..3c857e3 100644
--- a/ALFI/ALFI/ratf.h
+++ b/ALFI/ALFI/ratf.h
@@ -35,7 +35,7 @@ namespace alfi::ratf {
 	*/
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
 	std::enable_if_t<!traits::has_size<Number>::value, Number>
-	val_mul(const RationalFunction<Number, Container>& rf, Number x) {
+	val_mul(const RationalFunction<Number, Container>& rf, const Number& x) {
 		Number n = 0;
 		for (const auto& c : rf.first) {
 			n = n * x + c;
@@ -73,7 +73,7 @@ namespace alfi::ratf {
 	*/
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
 	std::enable_if_t<!traits::has_size<Number>::value, Number>
-	val_div(const RationalFunction<Number, Container>& rf, Number x) {
+	val_div(const RationalFunction<Number, Container>& rf, const Number& x) {
 		const auto& numerator = rf.first;
 		const auto& denominator = rf.second;
 		Number n = 0;
@@ -116,7 +116,7 @@ namespace alfi::ratf {
 	*/
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
 	std::enable_if_t<!traits::has_size<Number>::value, Number>
-	val(const RationalFunction<Number, Container>& rf, Number x) {
+	val(const RationalFunction<Number, Container>& rf, const Number& x) {
 		if (std::abs(x) <= 1) {
 			return val_mul(rf, x);
 		} else {
@@ -167,7 +167,7 @@ namespace alfi::ratf {
 		@return a pair `{numerator, denominator}` representing the Pade approximant; if an approximant does not exist, an empty pair is returned
 	 */
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	RationalFunction<Number,Container> pade(Container<Number> P, SizeT n, SizeT m, Number epsilon = std::numeric_limits<Number>::epsilon()) {
+	RationalFunction<Number,Container> pade(Container<Number> P, SizeT n, SizeT m, const Number& epsilon = std::numeric_limits<Number>::epsilon()) {
 		if constexpr (std::is_signed_v<SizeT>) {
 			if (n < 0 || m < 0) {
 				return {{}, {}};
diff --git a/ALFI/ALFI/spline/cubic.h b/ALFI/ALFI/spline/cubic.h
index b743f28..cdf2f7d 100644
--- a/ALFI/ALFI/spline/cubic.h
+++ b/ALFI/ALFI/spline/cubic.h
@@ -527,25 +527,25 @@ namespace alfi::spline {
 			_coeffs = std::move(coeffs);
 		}
 
-		Number eval(Number x) const {
+		Number eval(const Number& x) const {
 			return eval(x, std::distance(_X.begin(), util::misc::first_leq_or_begin(_X.begin(), _X.end(), x)));
 		}
-		Number eval(Number x, SizeT segment) const {
+		Number eval(const Number& x, SizeT segment) const {
 			if (_coeffs.empty()) {
 				return NAN;
 			} else if (_coeffs.size() == 1) {
 				return _coeffs[0];
 			}
 			segment = std::clamp(segment, static_cast<SizeT>(0), static_cast<SizeT>(_X.size() - 2));
-			x = x - _X[segment];
-			return ((_coeffs[4*segment] * x + _coeffs[4*segment+1]) * x + _coeffs[4*segment+2]) * x + _coeffs[4*segment+3];
+			const Number x_seg = x - _X[segment];
+			return ((_coeffs[4*segment] * x_seg + _coeffs[4*segment+1]) * x_seg + _coeffs[4*segment+2]) * x_seg + _coeffs[4*segment+3];
 		}
 
 		Container<Number> eval(const Container<Number>& xx, bool sorted = true) const {
 			Container<Number> result(xx.size());
 			if (sorted) {
 				for (SizeT i = 0, i_x = 0; i < xx.size(); ++i) {
-					const Number x = xx[i];
+					const Number& x = xx[i];
 					while (i_x + 1 < _X.size() && x >= _X[i_x+1])
 						++i_x;
 					result[i] = eval(x, i_x);
@@ -558,7 +558,7 @@ namespace alfi::spline {
 			return result;
 		}
 
-		Number operator()(Number x) const {
+		Number operator()(const Number& x) const {
 			return eval(x);
 		}
 		Container<Number> operator()(const Container<Number>& xx) const {
diff --git a/ALFI/ALFI/spline/linear.h b/ALFI/ALFI/spline/linear.h
index 7de70ac..4fb5f94 100644
--- a/ALFI/ALFI/spline/linear.h
+++ b/ALFI/ALFI/spline/linear.h
@@ -76,25 +76,25 @@ namespace alfi::spline {
 			_coeffs = std::move(coeffs);
 		}
 
-		Number eval(Number x) const {
+		Number eval(const Number& x) const {
 			return eval(x, std::distance(_X.begin(), util::misc::first_leq_or_begin(_X.begin(), _X.end(), x)));
 		}
-		Number eval(Number x, SizeT segment) const {
+		Number eval(const Number& x, SizeT segment) const {
 			if (_coeffs.empty()) {
 				return NAN;
 			} else if (_coeffs.size() == 1) {
 				return _coeffs[0];
 			}
 			segment = std::clamp(segment, static_cast<SizeT>(0), static_cast<SizeT>(_X.size() - 2));
-			x = x - _X[segment];
-			return _coeffs[2*segment] * x + _coeffs[2*segment+1];
+			const Number x_seg = x - _X[segment];
+			return _coeffs[2*segment] * x_seg + _coeffs[2*segment+1];
 		}
 
 		Container<Number> eval(const Container<Number>& xx, bool sorted = true) const {
 			Container<Number> result(xx.size());
 			if (sorted) {
 				for (SizeT i = 0, i_x = 0; i < xx.size(); ++i) {
-					const Number x = xx[i];
+					const Number& x = xx[i];
 					while (i_x + 1 < _X.size() && x >= _X[i_x+1])
 						++i_x;
 					result[i] = eval(x, i_x);
@@ -107,7 +107,7 @@ namespace alfi::spline {
 			return result;
 		}
 
-		Number operator()(Number x) const {
+		Number operator()(const Number& x) const {
 			return eval(x);
 		}
 		Container<Number> operator()(const Container<Number>& xx) const {
diff --git a/ALFI/ALFI/spline/polyeqv.h b/ALFI/ALFI/spline/polyeqv.h
index 44913ad..389a7c8 100644
--- a/ALFI/ALFI/spline/polyeqv.h
+++ b/ALFI/ALFI/spline/polyeqv.h
@@ -161,10 +161,10 @@ namespace alfi::spline {
 			_coeffs = std::move(coeffs);
 		}
 
-		Number eval(Number x) const {
+		Number eval(const Number& x) const {
 			return eval(x, std::distance(_X.begin(), util::misc::first_leq_or_begin(_X.begin(), _X.end(), x)));
 		}
-		Number eval(Number x, SizeT segment) const {
+		Number eval(const Number& x, SizeT segment) const {
 			if (_coeffs.empty()) {
 				return NAN;
 			} else if (_coeffs.size() == 1) {
@@ -173,7 +173,7 @@ namespace alfi::spline {
 
 			segment = std::clamp(segment, static_cast<SizeT>(0), static_cast<SizeT>(_X.size() - 2));
 
-			x = x - _X[segment];
+			const Number x_seg = x - _X[segment];
 
 			const SizeT n = _X.size();
 
@@ -191,7 +191,7 @@ namespace alfi::spline {
 			}
 
 			for (SizeT i = 1; i < n; ++i) {
-				result *= x;
+				result *= x_seg;
 				Number current = _coeffs[segment*n+i];
 				if (std::isnan(current)) {
 					switch (_evaluation_type) {
@@ -215,7 +215,7 @@ namespace alfi::spline {
 			Container<Number> result(xx.size());
 			if (sorted) {
 				for (SizeT i = 0, i_x = 0; i < xx.size(); ++i) {
-					const Number x = xx[i];
+					const Number& x = xx[i];
 					while (i_x + 1 < _X.size() && x >= _X[i_x+1])
 						++i_x;
 					result[i] = eval(x, i_x);
@@ -228,7 +228,7 @@ namespace alfi::spline {
 			return result;
 		}
 
-		Number operator()(Number x) const {
+		Number operator()(const Number& x) const {
 			return eval(x);
 		}
 		Container<Number> operator()(const Container<Number>& xx) const {
diff --git a/ALFI/ALFI/spline/quadratic.h b/ALFI/ALFI/spline/quadratic.h
index 2b1e6c0..0ea7a96 100644
--- a/ALFI/ALFI/spline/quadratic.h
+++ b/ALFI/ALFI/spline/quadratic.h
@@ -298,25 +298,25 @@ namespace alfi::spline {
 			_coeffs = std::move(coeffs);
 		}
 
-		Number eval(Number x) const {
+		Number eval(const Number& x) const {
 			return eval(x, std::distance(_X.begin(), util::misc::first_leq_or_begin(_X.begin(), _X.end(), x)));
 		}
-		Number eval(Number x, SizeT segment) const {
+		Number eval(const Number& x, SizeT segment) const {
 			if (_coeffs.empty()) {
 				return NAN;
 			} else if (_coeffs.size() == 1) {
 				return _coeffs[0];
 			}
 			segment = std::clamp(segment, static_cast<SizeT>(0), static_cast<SizeT>(_X.size() - 2));
-			x = x - _X[segment];
-			return (_coeffs[3*segment] * x + _coeffs[3*segment+1]) * x + _coeffs[3*segment+2];
+			const Number x_seg = x - _X[segment];
+			return (_coeffs[3*segment] * x_seg + _coeffs[3*segment+1]) * x_seg + _coeffs[3*segment+2];
 		}
 
 		Container<Number> eval(const Container<Number>& xx, bool sorted = true) const {
 			Container<Number> result(xx.size());
 			if (sorted) {
 				for (SizeT i = 0, i_x = 0; i < xx.size(); ++i) {
-					const Number x = xx[i];
+					const Number& x = xx[i];
 					while (i_x + 1 < _X.size() && x >= _X[i_x+1])
 						++i_x;
 					result[i] = eval(x, i_x);
@@ -329,7 +329,7 @@ namespace alfi::spline {
 			return result;
 		}
 
-		Number operator()(Number x) const {
+		Number operator()(const Number& x) const {
 			return eval(x);
 		}
 		Container<Number> operator()(const Container<Number>& xx) const {
diff --git a/ALFI/ALFI/spline/step.h b/ALFI/ALFI/spline/step.h
index 3a1f456..a8963ea 100644
--- a/ALFI/ALFI/spline/step.h
+++ b/ALFI/ALFI/spline/step.h
@@ -60,10 +60,10 @@ namespace alfi::spline {
 			_Y = std::forward<ContainerYType>(Y);
 		}
 
-		Number eval(Number x) const {
+		Number eval(const Number& x) const {
 			return eval(x, std::distance(_X.begin(), util::misc::first_leq_or_begin(_X.begin(), _X.end(), x)));
 		}
-		Number eval(Number x, SizeT segment) const {
+		Number eval(const Number& x, SizeT segment) const {
 			if (_Y.empty()) {
 				return NAN;
 			} else if (_Y.size() == 1) {
@@ -89,7 +89,7 @@ namespace alfi::spline {
 			Container<Number> result(xx.size());
 			if (sorted) {
 				for (SizeT i = 0, i_x = 0; i < xx.size(); ++i) {
-					const Number x = xx[i];
+					const Number& x = xx[i];
 					while (i_x + 1 < _X.size() && x >= _X[i_x+1])
 						++i_x;
 					result[i] = eval(x, i_x);
@@ -102,7 +102,7 @@ namespace alfi::spline {
 			return result;
 		}
 
-		Number operator()(Number x) const {
+		Number operator()(const Number& x) const {
 			return eval(x);
 		}
 		Container<Number> operator()(const Container<Number>& xx) const {
diff --git a/ALFI/ALFI/util/linalg.h b/ALFI/ALFI/util/linalg.h
index 522c0e0..071de9d 100644
--- a/ALFI/ALFI/util/linalg.h
+++ b/ALFI/ALFI/util/linalg.h
@@ -31,7 +31,7 @@ namespace alfi::util::linalg {
 		@return the solution vector or an empty container if the matrix is degenerate
 	 */
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> lup_solve(Container<Container<Number>>&& A, Container<Number>&& B, Number epsilon = std::numeric_limits<Number>::epsilon()) {
+	Container<Number> lup_solve(Container<Container<Number>>&& A, Container<Number>&& B, const Number& epsilon = std::numeric_limits<Number>::epsilon()) {
 		const auto n = B.size();
 		assert(n == A.size());
 		for (SizeT i = 0; i < n; ++i) {
diff --git a/ALFI/ALFI/util/numeric.h b/ALFI/ALFI/util/numeric.h
index 033bae3..9248b62 100644
--- a/ALFI/ALFI/util/numeric.h
+++ b/ALFI/ALFI/util/numeric.h
@@ -6,7 +6,7 @@
 
 namespace alfi::util::numeric {
 	template <typename Number = DefaultNumber>
-	bool are_equal(Number a, Number b, Number epsilon = std::numeric_limits<Number>::epsilon()) {
+	bool are_equal(const Number& a, const Number& b, const Number& epsilon = std::numeric_limits<Number>::epsilon()) {
 		return std::abs(a - b) <= epsilon || std::abs(a - b) <= std::max(std::abs(a), std::abs(b)) * epsilon;
 	}
 
diff --git a/ALFI/ALFI/util/points.h b/ALFI/ALFI/util/points.h
index 7cb5916..e86855f 100644
--- a/ALFI/ALFI/util/points.h
+++ b/ALFI/ALFI/util/points.h
@@ -4,7 +4,7 @@
 
 namespace alfi::points {
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	void lin_map(Container<Number>& points, Number a, Number b, Number c, Number d) {
+	void lin_map(Container<Number>& points, const Number& a, const Number& b, const Number& c, const Number& d) {
 		const auto mid1 = (a + b) / 2;
 		const auto mid2 = (c + d) / 2;
 		const auto scale = (d - c) / (b - a);
@@ -14,14 +14,14 @@ namespace alfi::points {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> lin_mapped(const Container<Number>& points, Number a, Number b, Number c, Number d) {
+	Container<Number> lin_mapped(const Container<Number>& points, const Number& a, const Number& b, const Number& c, const Number& d) {
 		auto mapped_points = points;
 		lin_map<Number,Container>(mapped_points, a, b, c, d);
 		return mapped_points;
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	void stretch(Container<Number>& points, Number a, Number b) {
+	void stretch(Container<Number>& points, const Number& a, const Number& b) {
 		if (points.empty()) {
 			return;
 		}
@@ -35,7 +35,7 @@ namespace alfi::points {
 	}
 
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	Container<Number> stretched(const Container<Number>& points, Number a, Number b) {
+	Container<Number> stretched(const Container<Number>& points, const Number& a, const Number& b) {
 		auto stretched_points = points;
 		stretch<Number,Container>(stretched_points, a, b);
 		return stretched_points;
diff --git a/ALFI/ALFI/util/poly.h b/ALFI/ALFI/util/poly.h
index 0aed7bb..7b8cea4 100644
--- a/ALFI/ALFI/util/poly.h
+++ b/ALFI/ALFI/util/poly.h
@@ -24,7 +24,7 @@ namespace alfi::util::poly {
 		@param epsilon the tolerance used to determine whether a coefficient is considered zero (default is machine epsilon)
 	 */
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	void normalize(Container<Number>& p, Number epsilon = std::numeric_limits<Number>::epsilon()) {
+	void normalize(Container<Number>& p, const Number& epsilon = std::numeric_limits<Number>::epsilon()) {
 		if (p.empty()) {
 			return p.push_back(0);
 		}
@@ -87,7 +87,7 @@ namespace alfi::util::poly {
 		@return a pair `{quotient, remainder}`
 	 */
 	template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
-	std::pair<Container<Number>,Container<Number>> div(const Container<Number>& dividend, const Container<Number>& divisor, Number epsilon = std::numeric_limits<Number>::epsilon()) {
+	std::pair<Container<Number>,Container<Number>> div(const Container<Number>& dividend, const Container<Number>& divisor, const Number& epsilon = std::numeric_limits<Number>::epsilon()) {
 		const auto divisor_start = std::find_if(divisor.begin(), divisor.end(), [&epsilon](Number v) { return std::abs(v) > epsilon; });
 
 		if (divisor_start == divisor.end() || dividend.size() < divisor.size()) {
diff --git a/tests/spline/test_cubic.cpp b/tests/spline/test_cubic.cpp
index 1c1e6c8..3740834 100644
--- a/tests/spline/test_cubic.cpp
+++ b/tests/spline/test_cubic.cpp
@@ -78,7 +78,7 @@ TEST(CubicSplineTest, BasicConditions) {
 	EXPECT_NEAR(spline.coeffs()[4], spline.coeffs()[8], 1e-17);
 	spline.construct(X, Y, alfi::spline::CubicSpline<>::Types::NotAKnotEnd{});
 	EXPECT_NEAR(spline.coeffs()[4*6], spline.coeffs()[4*7], 1e-14);
-	EXPECT_NEAR(spline.coeffs()[4*7], spline.coeffs()[4*8], 1e-15);
+	EXPECT_NEAR(spline.coeffs()[4*7], spline.coeffs()[4*8], 1e-14);
 }
 
 TEST(CubicSplineTest, CustomConditions) {
diff --git a/tests/spline/test_quadratic.cpp b/tests/spline/test_quadratic.cpp
index 82d0661..69b8684 100644
--- a/tests/spline/test_quadratic.cpp
+++ b/tests/spline/test_quadratic.cpp
@@ -96,7 +96,7 @@ TEST(QuadraticSplineTest, Conditions) {
 	}
 	spline1.construct(X, Y, alfi::spline::QuadraticSpline<>::Types::Clamped{5, 10});
 	EXPECT_EQ(10, spline1.coeffs()[3*5+1]);
-	EXPECT_NEAR(10, 2*spline1.coeffs()[3*4]*dX[8] + spline1.coeffs()[3*4+1], 1e-15);
+	EXPECT_NEAR(10, 2*spline1.coeffs()[3*4]*dX[8] + spline1.coeffs()[3*4+1], 1e-14);
 	// fixed-second
 	spline1.construct(X, Y, alfi::spline::QuadraticSpline<>::Types::FixedSecondStart{10});
 	EXPECT_EQ(10, 2*spline1.coeffs()[0]);