diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp
index 799e2ee..5d95deb 100644
--- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp
+++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp
@@ -65,7 +65,7 @@ class AdaptiveQuadratureBase
 	 * 
 	 * @return Pair (integral, estimated error).
 	 */
-	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegral(const Function& f, const Scalar& xmin, const Scalar& xmax) { return derived().estimateIntegralImpl(f, xmin, xmax); }
+	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegral(Function&& f, const Scalar& xmin, const Scalar& xmax) { return derived().estimateIntegralImpl(std::forward<Function>(f), xmin, xmax); }
 
 	/**
 	 * @brief Perform adaptive quadrature on [xmin, xmax].
@@ -75,9 +75,9 @@ class AdaptiveQuadratureBase
 	 * @param xmax Upper bound of interval.
 	 * @return Approximation of the integral.
 	 */
-	template<class Function> LongScalar integrate(const Function& f, const Scalar& xmin, const Scalar& xmax);
+	template<class Function> LongScalar integrate(Function&& f, const Scalar& xmin, const Scalar& xmax);
 
-	template<class Function> LongScalar integrateWithHints(const Function& f, const std::span<const Scalar> mu, const Scalar& sigma);
+	template<class Function> LongScalar integrateWithHints(Function&& f, const std::span<const Scalar> mu, const Scalar& sigma);
 
 	/**
 	 * @brief Perform adaptive quadrature on (-inf, xmax].
@@ -94,7 +94,7 @@ class AdaptiveQuadratureBase
 	 * @param xmax Upper bound of interval.
 	 * @return Approximation of the integral.
 	 */
-	template<class Function> LongScalar integrateLeftInfinite(const Function& f, const Scalar& xmax);
+	template<class Function> LongScalar integrateLeftInfinite(Function&& f, const Scalar& xmax);
 	
 	/**
 	 * @brief Perform adaptive quadrature on [xmin, inf).
@@ -111,7 +111,7 @@ class AdaptiveQuadratureBase
 	 * @param xmax Upper bound of interval.
 	 * @return Approximation of the integral.
 	 */
-	template<class Function> LongScalar integrateRightInfinite(const Function& f, const Scalar& xmin);
+	template<class Function> LongScalar integrateRightInfinite(Function&& f, const Scalar& xmin);
 	
 	/**
      * @brief Perform adaptive quadrature on (-inf, inf).
@@ -126,7 +126,7 @@ class AdaptiveQuadratureBase
 	 * \f]
 	 * Finally addapt the quadrature over [xmin, xmax].
      */
-	template<class Function> LongScalar integrate(const Function& f);
+	template<class Function> LongScalar integrate(Function&& f);
 	
 	/**
 	 * @brief Perform adaptive quadrature on (-inf, inf) using a coordinate-remapping technique.
@@ -141,9 +141,9 @@ class AdaptiveQuadratureBase
 	 * \int_{-1}^{1} f(x(t))\frac{1 + t^2}{(1 - t^2)^2} dt.
 	 * \f]
 	 */
-	template<class Function> LongScalar remapAndIntegrate(const Function& f);
+	template<class Function> LongScalar remapAndIntegrate(Function&& f);
 
-	template<class Function> std::invoke_result_t<Function, Scalar> integrateWithoutAdaptation(const Function& f) const;
+	template<class Function> std::invoke_result_t<Function, Scalar> integrateWithoutAdaptation(Function&& f) const;
 	
 	constexpr Size   getMaxIt()       const { return m_maxIt; }       ///<  @brief Maximum iterations allowed.
 	constexpr Size   getNits()        const { return m_it; }          ///<  @brief Number of iterations performed.
@@ -170,7 +170,7 @@ class AdaptiveQuadratureBase
 	
 	constexpr std::span<const Interval> getSubIntervals() const { return m_intervals; }
 private:
-	template<class Function> LongScalar adaptQuadrature(const Function& func);
+	template<class Function> LongScalar adaptQuadrature(Function&& func);
 
 	std::vector<Interval>   m_intervals;
 	std::vector<LongScalar> m_subIntergrals;
diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp
index 6a60cb9..f7caedf 100644
--- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp
+++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp
@@ -9,6 +9,7 @@
 #include <algorithm>
 #include <ranges>
 #include <cmath>
+#include <utility>
 
 #include <fmt/core.h>
 
@@ -28,7 +29,7 @@ AdaptiveQuadratureBase<Derived>::AdaptiveQuadratureBase(const Size& maxIt, const
 }
 
 template<class Derived> template<class Function> 
-auto AdaptiveQuadratureBase<Derived>::adaptQuadrature(const Function& f) -> LongScalar
+auto AdaptiveQuadratureBase<Derived>::adaptQuadrature(Function&& f) -> LongScalar
 {
 	using std::abs;
 	using std::isfinite;
@@ -72,7 +73,7 @@ auto AdaptiveQuadratureBase<Derived>::adaptQuadrature(const Function& f) -> Long
 }
 
 template<class Derived> template<class Function> 
-auto AdaptiveQuadratureBase<Derived>::integrate(const Function& f, const Scalar& xmin, const Scalar& xmax) -> LongScalar
+auto AdaptiveQuadratureBase<Derived>::integrate(Function&& f, const Scalar& xmin, const Scalar& xmax) -> LongScalar
 {	
 	using std::ceil;
 	
@@ -102,11 +103,11 @@ auto AdaptiveQuadratureBase<Derived>::integrate(const Function& f, const Scalar&
 	}
 
 	if (m_out) { fmt::print(m_out, "#NumericalIntegrator adapting quadrature over [{}, {}]\n", xmin, xmax); }
-	return adaptQuadrature(f);
+	return adaptQuadrature(std::forward<Function>(f));
 }
 
 template<class Derived> template<class Function> 
-auto AdaptiveQuadratureBase<Derived>::integrateWithHints(const Function& f, const std::span<const Scalar> mu, const Scalar& sigma) -> LongScalar
+auto AdaptiveQuadratureBase<Derived>::integrateWithHints(Function&& f, const std::span<const Scalar> mu, const Scalar& sigma) -> LongScalar
 {
 	using std::abs;
 	using std::isfinite;
@@ -182,11 +183,11 @@ auto AdaptiveQuadratureBase<Derived>::integrateWithHints(const Function& f, cons
 	}
 
 	if (m_out) { fmt::print(m_out, "#NumericalIntegrator adapting quadrature over {}\n", m_intervals); }
-	return adaptQuadrature(f); 
+	return adaptQuadrature(std::forward<Function>(f)); 
 }
 
 template<class Derived> template<class Function>
-auto AdaptiveQuadratureBase<Derived>::integrateLeftInfinite(const Function& f, const Scalar& xmax) -> LongScalar
+auto AdaptiveQuadratureBase<Derived>::integrateLeftInfinite(Function&& f, const Scalar& xmax) -> LongScalar
 {
 	using std::isfinite;
 	using std::abs;
@@ -200,11 +201,15 @@ auto AdaptiveQuadratureBase<Derived>::integrateLeftInfinite(const Function& f, c
 		xmin *= 2;
 		leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, xmin);
 	}
-	return isfinite(leftIntegral) ? integrate(f, xmin, xmax) : NumTraits<LongScalar>::NaN;
+	const LongScalar ret = isfinite(leftIntegral) 
+		? integrate(std::forward<Function>(f), xmin, xmax) 
+		: NumTraits<LongScalar>::NaN;
+	
+	return ret;
 }
 
 template<class Derived> template<class Function>
-auto AdaptiveQuadratureBase<Derived>::integrateRightInfinite(const Function& f, const Scalar& xmin) -> LongScalar
+auto AdaptiveQuadratureBase<Derived>::integrateRightInfinite(Function&& f, const Scalar& xmin) -> LongScalar
 {
 	using std::isfinite;
 	using std::abs;
@@ -218,11 +223,15 @@ auto AdaptiveQuadratureBase<Derived>::integrateRightInfinite(const Function& f,
 		xmax *= 2;
 		rightIntegral = gLaguerreQuad.integrateRightInfinite(f, xmax);
 	}
-	return isfinite(rightIntegral) ? integrate(f, xmin, xmax) : NumTraits<LongScalar>::NaN;
+	const LongScalar ret = isfinite(rightIntegral) 
+		? integrate(std::forward<Function>(f), xmin, xmax) 
+		: NumTraits<LongScalar>::NaN;
+	
+	return ret;
 }
 
 template<class Derived> template<class Function>
-auto AdaptiveQuadratureBase<Derived>::integrate(const Function& f) -> LongScalar
+auto AdaptiveQuadratureBase<Derived>::integrate(Function&& f) -> LongScalar
 {
 	using std::isfinite;
 	using std::abs;
@@ -237,31 +246,32 @@ auto AdaptiveQuadratureBase<Derived>::integrate(const Function& f) -> LongScalar
 		leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, xmin);
 	}
 	
-	return isfinite(leftIntegral) ? integrateRightInfinite(f, xmin) : NumTraits<LongScalar>::NaN;
+	const LongScalar ret = isfinite(leftIntegral) 
+		? integrateRightInfinite(std::forward<Function>(f), xmin) 
+		: NumTraits<LongScalar>::NaN;
+
+	return ret;
 }
 
 template<class Derived> template<class Function>
-auto AdaptiveQuadratureBase<Derived>::remapAndIntegrate(const Function& f) -> LongScalar
+auto AdaptiveQuadratureBase<Derived>::remapAndIntegrate(Function&& f) -> LongScalar
 {		
 	using std::isnan;
 	
-	constexpr Scalar eps = {};
-	
-	const auto fref = [&f](const Scalar t) -> LongScalar
+	return integrate([&f](const Scalar t) -> LongScalar
 	{			
-		const LongScalar fx = f(t / (1 - t*t));
-		const Scalar dxdt = (1 + t*t) / ((1 - t*t)*(1 - t*t));
-			
-		return isnan(fx*dxdt)
+		const LongScalar fx   = f(t / (1 - t*t));
+		const Scalar     dxdt = (1 + t*t) / ((1 - t*t)*(1 - t*t));
+		const LongScalar ret  = isnan(fx*dxdt)
 			? LongScalar{}
 			: fx*dxdt;	
-	};
-	
-	return integrate(fref, -1, 1);
+		
+		return ret;
+	}, -1, 1);
 }
 
 template<class Derived>  template<class Function> 
-auto AdaptiveQuadratureBase<Derived>::integrateWithoutAdaptation(const Function& f) const -> std::invoke_result_t<Function, Scalar>
+auto AdaptiveQuadratureBase<Derived>::integrateWithoutAdaptation(Function&& f) const -> std::invoke_result_t<Function, Scalar>
 {
 	const auto localIntegrals = m_intervals | std::views::transform([&self = derived(), &f](const Interval& interval) -> std::invoke_result_t<Function, Scalar>
 	{
diff --git a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature.hpp b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature.hpp
index ae2a1ae..e2832f3 100644
--- a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature.hpp
+++ b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature.hpp
@@ -45,9 +45,9 @@ class ClenshawCurtisAdaptiveQuadrature : public AdaptiveQuadratureBase< Clenshaw
 	 * @param xmax Upper bound.
 	 * @return Pair (integral, estimated error).
 	 */
-	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax);
+	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax);
 	
-	template<class Function> constexpr std::invoke_result_t<Function, Scalar> integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const;
+	template<class Function> constexpr std::invoke_result_t<Function, Scalar> integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const;
 	
 	constexpr Scalar getMaxDeltaXImpl(const Scalar& xmin, const Scalar& xmax) const { return (xmax - xmin)*misc::maxDiff(std::span{s_xi}); } 
 private:
diff --git a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature_impl.hpp b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature_impl.hpp
index 756c74b..4bb6e06 100644
--- a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature_impl.hpp
+++ b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature_impl.hpp
@@ -21,7 +21,7 @@ extern template class ClenshawCurtisAdaptiveQuadrature<long double, long double>
 //// method implementations ////
 
 template<typename T, typename TT> template<class Function>
-constexpr auto ClenshawCurtisAdaptiveQuadrature<T,TT>::estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) -> std::pair<LongScalar, LongScalar>
+constexpr auto ClenshawCurtisAdaptiveQuadrature<T,TT>::estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) -> std::pair<LongScalar, LongScalar>
 {	
 	using std::abs;
 	
@@ -52,7 +52,7 @@ constexpr auto ClenshawCurtisAdaptiveQuadrature<T,TT>::estimateIntegralImpl(cons
 }
 
 template<typename T, typename TT> template<class Function>
-constexpr auto ClenshawCurtisAdaptiveQuadrature<T,TT>::integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t<Function, Scalar>
+constexpr auto ClenshawCurtisAdaptiveQuadrature<T,TT>::integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t<Function, Scalar>
 {
 	const auto fx = s_xi | std::views::transform([&f, &xmin, &xmax](const Scalar& xi) -> std::invoke_result_t<Function, Scalar>
 	{
diff --git a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature.hpp b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature.hpp
index 75a97ce..8ef1005 100644
--- a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature.hpp
+++ b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature.hpp
@@ -48,9 +48,9 @@ class ClenshawCurtisHybridAdaptiveQuadrature : public AdaptiveQuadratureBase< Cl
 	 * @param xmax Upper bound.
 	 * @return Pair (integral, estimated error).
 	 */
-	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax);
+	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax);
 	
-	template<class Function> constexpr std::invoke_result_t<Function, Scalar> integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const;
+	template<class Function> constexpr std::invoke_result_t<Function, Scalar> integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const;
 	
 	constexpr Scalar getMaxDeltaXImpl(const Scalar& xmin, const Scalar& xmax) const { return (xmax - xmin)*misc::maxDiff(std::span{s_xi}); } 
 private:
diff --git a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature_impl.hpp b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature_impl.hpp
index 8247edf..34bec7a 100644
--- a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature_impl.hpp
+++ b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature_impl.hpp
@@ -20,7 +20,7 @@ extern template class ClenshawCurtisHybridAdaptiveQuadrature<long double, long d
 //// method implementations ////
 
 template<typename T, typename TT> template<class Function>
-constexpr auto ClenshawCurtisHybridAdaptiveQuadrature<T,TT>::estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) -> std::pair<LongScalar, LongScalar>
+constexpr auto ClenshawCurtisHybridAdaptiveQuadrature<T,TT>::estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) -> std::pair<LongScalar, LongScalar>
 {	
 	using std::abs;
 	
@@ -38,7 +38,7 @@ constexpr auto ClenshawCurtisHybridAdaptiveQuadrature<T,TT>::estimateIntegralImp
 }
 
 template<typename T, typename TT> template<class Function>
-constexpr auto ClenshawCurtisHybridAdaptiveQuadrature<T,TT>::integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t<Function, Scalar>
+constexpr auto ClenshawCurtisHybridAdaptiveQuadrature<T,TT>::integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t<Function, Scalar>
 {
 	const auto fx = s_xi | std::views::transform([&f, &xmin, &xmax](const Scalar& xi) -> std::invoke_result_t<Function, Scalar>
 	{
diff --git a/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature.hpp b/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature.hpp
index 40d7ba5..1c58758 100644
--- a/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature.hpp
+++ b/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature.hpp
@@ -49,9 +49,9 @@ class GLCCAdaptiveQuadrature : public AdaptiveQuadratureBase< GLCCAdaptiveQuadra
 	 * @param xmax Upper bound.
 	 * @return Pair (integral, estimated error).
 	 */
-	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const;
+	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const;
 	
-	template<class Function> constexpr std::invoke_result_t<Function, Scalar> integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const;
+	template<class Function> constexpr std::invoke_result_t<Function, Scalar> integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const;
 	
 	constexpr Scalar getMaxDeltaXImpl(const Scalar xmin, const Scalar xmax) const { return (xmax - xmin)*misc::maxDiff(std::span{s_xi_cc}); } 
 private:
diff --git a/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature_impl.hpp b/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature_impl.hpp
index 906dadc..7306575 100644
--- a/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature_impl.hpp
+++ b/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature_impl.hpp
@@ -20,7 +20,7 @@ extern template class GLCCAdaptiveQuadrature<long double, long double>;
 //// method implementations ////
 
 template<typename T, typename TT> template<class Function>
-constexpr auto GLCCAdaptiveQuadrature<T,TT>::estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::pair<LongScalar, LongScalar>
+constexpr auto GLCCAdaptiveQuadrature<T,TT>::estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::pair<LongScalar, LongScalar>
 {
 	using std::abs;
 	
@@ -43,7 +43,7 @@ constexpr auto GLCCAdaptiveQuadrature<T,TT>::estimateIntegralImpl(const Function
 }
 
 template<typename T, typename TT> template<class Function>
-constexpr auto GLCCAdaptiveQuadrature<T,TT>::integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t<Function, Scalar>
+constexpr auto GLCCAdaptiveQuadrature<T,TT>::integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t<Function, Scalar>
 {
 	const auto fx = s_xi_gl | std::views::transform([&f, &xmin, &xmax](const Scalar& xi) -> std::invoke_result_t<Function, Scalar>
 	{
diff --git a/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature.hpp b/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature.hpp
index c88c134..4b15474 100644
--- a/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature.hpp
+++ b/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature.hpp
@@ -36,9 +36,9 @@ class GaussLegendreAdaptiveQuadrature : public AdaptiveQuadratureBase< GaussLege
 	 * @param xmax Upper bound.
 	 * @return Pair (integral, estimated error).
 	 */
-	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax);
+	template<class Function> constexpr std::pair<LongScalar, LongScalar> estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax);
 	
-	template<class Function> constexpr std::invoke_result_t<Function, Scalar> integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const;
+	template<class Function> constexpr std::invoke_result_t<Function, Scalar> integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const;
 	
 	inline constexpr Scalar getMaxDeltaXImpl(const Scalar xmin, const Scalar xmax) const { return (xmax - xmin)*misc::maxDiff(std::span{s_xi}); } 
 private:
diff --git a/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature_impl.hpp b/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature_impl.hpp
index bf3fdff..a348c84 100644
--- a/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature_impl.hpp
+++ b/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature_impl.hpp
@@ -20,7 +20,7 @@ extern template class GaussLegendreAdaptiveQuadrature<long double, long double>;
 //// method implementations ////
 
 template<typename T, typename TT> template<class Function>
-constexpr auto GaussLegendreAdaptiveQuadrature<T,TT>::estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) -> std::pair<LongScalar, LongScalar>
+constexpr auto GaussLegendreAdaptiveQuadrature<T,TT>::estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) -> std::pair<LongScalar, LongScalar>
 {
 	using std::abs;
 	
@@ -43,7 +43,7 @@ constexpr auto GaussLegendreAdaptiveQuadrature<T,TT>::estimateIntegralImpl(const
 }
 
 template<typename T, typename TT> template<class Function>
-constexpr auto GaussLegendreAdaptiveQuadrature<T,TT>::integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t<Function, Scalar>
+constexpr auto GaussLegendreAdaptiveQuadrature<T,TT>::integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t<Function, Scalar>
 {
 	const auto fx = s_xi | std::views::transform([&f, &xmin, &xmax](const Scalar& xi) -> std::invoke_result_t<Function, Scalar>
 	{
diff --git a/include/LNIT/GaussHermiteQuadrature.hpp b/include/LNIT/GaussHermiteQuadrature.hpp
index 9f9fafd..c0c0f6b 100644
--- a/include/LNIT/GaussHermiteQuadrature.hpp
+++ b/include/LNIT/GaussHermiteQuadrature.hpp
@@ -24,7 +24,7 @@ class GaussHermiteQuadrature
 
 	constexpr GaussHermiteQuadrature() {} ///<  @brief Default constructor.
 
-	template<class Function> constexpr LongScalar integrate(const Function& f) const;
+	template<class Function> constexpr LongScalar integrate(Function&& f) const;
 private:
 	static constexpr std::array<Scalar, 66> s_wi = {
 		Scalar(0.72500244352094499799), Scalar(0.55908539234713480356), Scalar(0.48645731979545141956), Scalar(0.44287166385270444111), Scalar(0.41289257970879246181), 
diff --git a/include/LNIT/GaussHermiteQuadrature_impl.hpp b/include/LNIT/GaussHermiteQuadrature_impl.hpp
index bd88aab..904c196 100644
--- a/include/LNIT/GaussHermiteQuadrature_impl.hpp
+++ b/include/LNIT/GaussHermiteQuadrature_impl.hpp
@@ -22,14 +22,14 @@ extern template class GaussHermiteQuadrature<long double, long double>;
 //// method implementations ////
 
 template<typename Scalar, typename LongScalar> template<class Function> 
-constexpr LongScalar GaussHermiteQuadrature<Scalar, LongScalar>::integrate(const Function& f) const
+constexpr LongScalar GaussHermiteQuadrature<Scalar, LongScalar>::integrate(Function&& f) const
 {
 	const auto fx = s_xi | std::views::transform([&f](const Scalar& x) -> LongScalar
 	{
 		return f(x); 
 	});
 
-	return std::inner_product(s_wi.begin(), s_wi.end(), fx.end(), LongScalar{});	
+	return std::inner_product(s_wi.begin(), s_wi.end(), fx.begin(), LongScalar{});	
 }
 
 } // namespace LNIT
diff --git a/include/LNIT/GaussLaguerreQuadrature.hpp b/include/LNIT/GaussLaguerreQuadrature.hpp
index 1d0522e..ed3df28 100644
--- a/include/LNIT/GaussLaguerreQuadrature.hpp
+++ b/include/LNIT/GaussLaguerreQuadrature.hpp
@@ -38,7 +38,7 @@ class GaussLaguerreQuadrature
 	 * @param a Upper bound of integration (default = 0).
 	 * @return Approximated integral value.
 	 */
-	template<class Function> constexpr LongScalar integrateLeftInfinite  (const Function& f, const Scalar& a = Scalar{}) const;
+	template<class Function> constexpr LongScalar integrateLeftInfinite  (Function&& f, const Scalar& a = Scalar{}) const;
 	/**
 	 * @brief Approximate integral over the right semi-infinite interval.
 	 *
@@ -53,7 +53,7 @@ class GaussLaguerreQuadrature
 	 * @param a Lower bound of integration (default = 0).
 	 * @return Approximated integral value.
 	 */
-	template<class Function> constexpr LongScalar integrateRightInfinite (const Function& f, const Scalar& a = Scalar{}) const;
+	template<class Function> constexpr LongScalar integrateRightInfinite (Function&& f, const Scalar& a = Scalar{}) const;
 private:
 	static constexpr std::array<Scalar, 33> s_wi = {
 		Scalar(0.11077730587320757274), Scalar(0.25810528128189475158), Scalar(0.40622176868437369247), Scalar(0.55526230959922306292), Scalar(0.70555738765958285661), 
diff --git a/include/LNIT/GaussLaguerreQuadrature_impl.hpp b/include/LNIT/GaussLaguerreQuadrature_impl.hpp
index 9fa7486..526ce53 100644
--- a/include/LNIT/GaussLaguerreQuadrature_impl.hpp
+++ b/include/LNIT/GaussLaguerreQuadrature_impl.hpp
@@ -18,25 +18,25 @@ extern template class GaussLaguerreQuadrature<long double, long double>;
 //// method implementations ////
 
 template<typename Scalar, typename LongScalar> template<class Function> 
-constexpr LongScalar GaussLaguerreQuadrature<Scalar, LongScalar>::integrateLeftInfinite(const Function& f, const Scalar& a) const
+constexpr LongScalar GaussLaguerreQuadrature<Scalar, LongScalar>::integrateLeftInfinite(Function&& f, const Scalar& a) const
 {
 	const auto fx = s_xi | std::views::transform([&f, &a](const Scalar& x) -> LongScalar
 	{
 		return f(a - x); 
 	});
 
-	return std::inner_product(s_wi.begin(), s_wi.end(), fx.end(), LongScalar{});	
+	return std::inner_product(s_wi.begin(), s_wi.end(), fx.begin(), LongScalar{});	
 }
 
 template<typename Scalar, typename LongScalar> template<class Function> 
-constexpr LongScalar GaussLaguerreQuadrature<Scalar, LongScalar>::integrateRightInfinite(const Function& f, const Scalar& a) const
+constexpr LongScalar GaussLaguerreQuadrature<Scalar, LongScalar>::integrateRightInfinite(Function&& f, const Scalar& a) const
 {
 	const auto fx = s_xi | std::views::transform([&f, &a](const Scalar& x) -> LongScalar
 	{
 		return f(x + a); 
 	});
 
-	return std::inner_product(s_wi.begin(), s_wi.end(), fx.end(), LongScalar{});	
+	return std::inner_product(s_wi.begin(), s_wi.end(), fx.begin(), LongScalar{});	
 }
 
 } // namespace LNIT