Rvalue ref

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;

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>
 	{

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:

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>
 	{

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:

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>
 	{

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:

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>
 	{

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:

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>
 	{

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), 

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

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),