From 951e4e381d04e91ca911ebfdcfe926a4a8a3a451 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 5 Sep 2026 11:05:31 -0700 Subject: [PATCH] Support explicit zero and norm across quadrature integrators --- doc/quadrature/double_exponential.qbk | 58 +++ doc/quadrature/gauss.qbk | 47 +++ doc/quadrature/gauss_kronrod.qbk | 43 ++ doc/quadrature/trapezoidal.qbk | 30 ++ .../quadrature/detail/exp_sinh_detail.hpp | 141 +++++++ .../detail/norm_quadrature_error.hpp | 20 + .../quadrature/detail/sinh_sinh_detail.hpp | 134 ++++++ .../quadrature/detail/tanh_sinh_detail.hpp | 380 ++++++++++++++++++ include/boost/math/quadrature/exp_sinh.hpp | 62 +++ include/boost/math/quadrature/gauss.hpp | 124 ++++++ .../boost/math/quadrature/gauss_kronrod.hpp | 190 +++++++++ include/boost/math/quadrature/sinh_sinh.hpp | 9 + include/boost/math/quadrature/tanh_sinh.hpp | 245 +++++++++++ include/boost/math/quadrature/trapezoidal.hpp | 98 +++++ test/Jamfile.v2 | 9 + test/double_exponential_eigen_test.cpp | 252 ++++++++++++ test/gauss_kronrod_eigen_test.cpp | 121 ++++++ test/gauss_quadrature_eigen_test.cpp | 168 ++++++++ 18 files changed, 2131 insertions(+) create mode 100644 include/boost/math/quadrature/detail/norm_quadrature_error.hpp create mode 100644 test/double_exponential_eigen_test.cpp create mode 100644 test/gauss_kronrod_eigen_test.cpp create mode 100644 test/gauss_quadrature_eigen_test.cpp diff --git a/doc/quadrature/double_exponential.qbk b/doc/quadrature/double_exponential.qbk index 2959b94cd3..cb0a779901 100644 --- a/doc/quadrature/double_exponential.qbk +++ b/doc/quadrature/double_exponential.qbk @@ -380,6 +380,64 @@ check we end up with `0 * Infinity` as the result (a NaN). [endsect] [/section:de_exp_sinh exp_sinh] +[section:de_matrix Vector- and matrix-valued integrands] + +The CPU `exp_sinh`, `sinh_sinh`, and `tanh_sinh` classes also accept an explicit +additive identity and scalar norm. Insert `zero, norm` after the bounds, or after +`f` when using the default domain: + + integrator.integrate(f, zero, norm, tolerance, &error, &L1, &levels); + integrator.integrate(f, a, b, zero, norm, tolerance, &error, &L1, &levels); + +The bounded form is available for `exp_sinh` and `tanh_sinh`; `sinh_sinh` always +integrates over the whole real line. Arguments after `norm` retain their existing +defaults. The new overloads participate only if applying `norm` to the integrand's +result can be explicitly converted to `Real`. Existing overloads remain unchanged. + +For example, using Eigen with exp-sinh: + + using Matrix = Eigen::Matrix, 2, 2>; + Matrix zero = Matrix::Zero(); + auto f = [](double x) -> Matrix { + Matrix result; + result << std::exp(-x), 0., 0., std::exp(-2*x); + return result; + }; + auto norm = [](const Matrix& m) { return m.stableNorm(); }; + boost::math::quadrature::exp_sinh integrator; + Matrix result = integrator.integrate(f, zero, norm); + +Return concrete owning values, with dimensions matching `zero`, rather than +expression templates referencing temporary values. The value type must support +addition, addition assignment, unary negation, scalar multiplication on either +side, and multiplication assignment by a scalar. The norm must be nonnegative, +vanish only at zero, and return a scalar convertible to `Real`. It should also +report non-finite values for non-finite inputs; endpoint checks use this property. +A norm implementation that avoids intermediate overflow and underflow is useful +for the large dynamic ranges encountered near endpoints. + +The selected norm is used for successive-estimate errors, the L1 integral, +tail truncation, and checks for non-finite values. Convergence compares the error +with `tolerance * L1`; it does not impose a separate relative tolerance on every +entry. Error estimates are not rigorous bounds, and exhausting refinement levels +does not guarantee that the requested tolerance was reached. + +Both forms of the two-argument tanh-sinh integrand are supported: `f(x, xc)` may +return a matrix, while `x` and the signed endpoint distance `xc` remain scalars. +For the new bounded tanh-sinh overloads, equal bounds return `zero` and set requested +error, L1, and level outputs to zero without calling `f`. Reversed finite bounds +negate the integral; the two-argument functor receives distances for the endpoints +in increasing order. Exp-sinh retains its requirement for a half-infinite domain. + +Domain errors use the scalar policy. If it returns, the result is `zero` times +its scalar error value; requested error and L1 outputs receive that value, and +levels is set to zero. Evaluation-error policies receive the scalar norm of the +problematic value; a non-throwing evaluation policy returns that concrete value +or current estimate, as appropriate. Diagnostic outputs are not specified after +an evaluation error. GPU free-function interfaces are unchanged. + +[endsect] + [section:de_tol Setting the Termination Condition for Integration] The integrate method for all three double-exponential quadratures supports ['tolerance] argument that acts as the diff --git a/doc/quadrature/gauss.qbk b/doc/quadrature/gauss.qbk index e2bdd6b478..2da438bbd2 100644 --- a/doc/quadrature/gauss.qbk +++ b/doc/quadrature/gauss.qbk @@ -85,6 +85,53 @@ so it can be effectively computed via Gaussian quadrature using the following co Complex W = integrator.integrate(lw, (Real) 0, pi()); +[heading Vector- and matrix-valued integrands] + +Additional overloads accept an explicit additive identity and a scalar-valued norm: + + template + static auto integrate(F f, const decltype(f(Real(0)))& zero, Norm norm, + Real* pL1 = nullptr) -> decltype(f(Real(0))); + + template + static auto integrate(F f, Real a, Real b, + const decltype(f(Real(0)))& zero, Norm norm, + Real* pL1 = nullptr) -> decltype(f(Real(0))); + +These overloads participate in overload resolution only when `norm(f(Real(0)))` +is convertible to `Real` by an explicit cast. Existing overloads are unchanged. + +The integrand must return a concrete value type supporting addition, addition +assignment, unary negation, and multiplication by `Real` on either side. +All returned values and `zero` must have compatible dimensions. In particular, +return an owning matrix rather than an unevaluated expression referencing temporaries. +`norm` must return a nonnegative scalar norm convertible to `Real`. +If requested, `pL1` receives the quadrature approximation to the integral of +`norm(f(x))` over the interval in increasing order. This is not an error estimate. + +For example, with Eigen available: + + using Matrix = Eigen::Matrix, 2, 2>; + Matrix zero = Matrix::Zero(); + auto f = [](double x) -> Matrix { + Matrix result; + result << x, 0., 0., x*x; + return result; + }; + auto norm = [](const Matrix& m) { return m.norm(); }; + Matrix result = boost::math::quadrature::gauss::integrate( + f, 0., 1., zero, norm); + +The explicit zero also allows dynamically sized matrices: initialize it with the required dimensions. +For equal bounds the supplied zero is returned, `*pL1` is set to zero if requested, +and the integrand is not evaluated. + +Invalid bounds invoke the existing scalar domain-error policy. If that policy +returns a scalar error value instead of throwing, the result is `zero` multiplied +by that value, and `*pL1` receives that value if requested. For floating-point Eigen +matrices and the ignore-error policy, this produces a matrix of NaNs with the +supplied dimensions. Custom value types must support this multiplication too. + [heading Choosing the number of points] Internally class `gauss` has pre-computed tables of abscissa and weights for 7, 15, 20, 25 and 30 points at up to 100-decimal diff --git a/doc/quadrature/gauss_kronrod.qbk b/doc/quadrature/gauss_kronrod.qbk index 31aeba9bba..449ab5c7d1 100644 --- a/doc/quadrature/gauss_kronrod.qbk +++ b/doc/quadrature/gauss_kronrod.qbk @@ -42,6 +42,49 @@ with no end point singularities. For difficult functions, or those with end poi Real* pL1 = nullptr)->decltype(std::declval()(std::declval())); }; +[heading Vector- and matrix-valued integrands] + +An additional overload accepts an explicit additive identity and a scalar norm: + + template + static auto integrate(F f, Real a, Real b, + const decltype(f(a))& zero, Norm norm, + unsigned max_depth = 15, + Real tol = tools::root_epsilon(), + Real* error = nullptr, Real* pL1 = nullptr) + -> decltype(f(a)); + +This overload participates only when `norm(f(a))` can be explicitly converted +to `Real`. Existing overloads and their behavior are unchanged. + +As with the [link math_toolkit.gauss Gauss overloads], `f` must return a concrete +value type supporting addition, addition assignment, unary negation, and scalar +multiplication on either side. The supplied zero and all integrand values must +have compatible dimensions. The norm must be nonnegative and convertible to `Real`. +For example, with Eigen: + + Eigen::MatrixXcd zero = Eigen::MatrixXcd::Zero(rows, cols); + auto norm = [](const Eigen::MatrixXcd& m) { return m.norm(); }; + Eigen::MatrixXcd result = boost::math::quadrature::gauss_kronrod::integrate( + f, a, b, zero, norm, 15, 1e-10); + +The supplied norm is used for the Gauss-Kronrod difference, the roundoff floor, +the relative convergence test, and the L1 integral. Tolerance therefore applies +to that norm, not separately to each matrix entry. `error` is an estimate, not a +rigorous bound, and reaching `max_depth` does not guarantee the requested tolerance. +In this overload, local error estimates are scaled by interval width, and by the +additional factor in the half-infinite substitutions, to match the units of the +returned integral. Reversing finite bounds negates the integral but leaves error +and L1 nonnegative. + +Equal bounds return the supplied zero and set requested error and L1 outputs to +zero without evaluating `f`. Invalid bounds invoke the scalar domain-error policy; +if it returns instead of throwing, the result is `zero` multiplied by the policy's +scalar error value, and requested error and L1 outputs receive that value. For +floating-point Eigen matrices the standard ignore-error policy thus produces a +matrix of NaNs with the supplied dimensions. No Eigen dependency is added to +Boost.Math headers. + [heading Description] static const RandomAccessContainer& abscissa(); diff --git a/doc/quadrature/trapezoidal.qbk b/doc/quadrature/trapezoidal.qbk index de01cb838a..dd2cb0553c 100644 --- a/doc/quadrature/trapezoidal.qbk +++ b/doc/quadrature/trapezoidal.qbk @@ -25,6 +25,36 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) }}} // namespaces +[heading Vector- and matrix-valued integrands] + +An additional overload accepts an explicit zero value and scalar norm: + + template + auto trapezoidal(F f, Real a, Real b, const decltype(f(a))& zero, Norm norm, + Real tol = boost::math::tools::root_epsilon(), + std::size_t max_refinements = 12, + Real* error_estimate = nullptr, Real* L1 = nullptr) -> decltype(f(a)); + +A policy overload accepts the same arguments without defaults, followed by +`const Policy& pol`. These overloads participate only if `norm(f(a))` can be +explicitly converted to `Real`. Existing overloads remain unchanged. + +Return a concrete owning vector or matrix, with dimensions matching `zero`. +The value type needs addition, addition assignment, unary negation, and scalar +multiplication. The supplied norm controls the successive-estimate error and +L1 integral, so tolerance is measured against `tol * L1`, not independently for +each entry. For an Eigen matrix, a typical norm is +`[](const Matrix& m) { return m.norm(); }`. + + Matrix result = boost::math::quadrature::trapezoidal( + f, a, b, zero, norm, 1e-10); + +Equal bounds return `zero` and set requested error and L1 outputs to zero without +calling `f`. Reversed bounds negate the integral and preserve nonnegative error +and L1 estimates. Invalid bounds invoke the scalar domain-error policy; if it +returns, the result is `zero` multiplied by its scalar error value, and requested +error and L1 outputs receive that value. + [heading Description] The functional `trapezoidal` calculates the integral of a function /f/ using the surprisingly simple trapezoidal rule. diff --git a/include/boost/math/quadrature/detail/exp_sinh_detail.hpp b/include/boost/math/quadrature/detail/exp_sinh_detail.hpp index cdce12c0c8..32a3718e3c 100644 --- a/include/boost/math/quadrature/detail/exp_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/exp_sinh_detail.hpp @@ -10,6 +10,7 @@ #include #ifndef BOOST_MATH_HAS_NVRTC +#include #ifndef BOOST_MATH_BUILD_MODULE #include @@ -53,6 +54,8 @@ class exp_sinh_detail template auto integrate(const F& f, Real* error, Real* L1, const char* function, Real tolerance, std::size_t* levels) const ->decltype(std::declval()(std::declval())); + template + auto integrate(const F& f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* error, Real* L1, const char* function, Real tolerance, std::size_t* levels) const ->decltype(std::declval()(std::declval())); private: const std::vector& get_abscissa_row(std::size_t n)const @@ -295,6 +298,144 @@ auto exp_sinh_detail::integrate(const F& f, Real* error, Real* L1, return I1; } +template +template +auto exp_sinh_detail::integrate(const F& f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* error, Real* L1, const char* function, Real tolerance, std::size_t* levels) const ->decltype(std::declval()(std::declval())) +{ + const auto magnitude = [&](const decltype(std::declval()(std::declval()))& value) -> Real { return static_cast(norm(value)); }; + typedef decltype(f(static_cast(0))) K; + using std::abs; + using std::floor; + using std::tanh; + using std::sinh; + using std::sqrt; + using boost::math::constants::half; + using boost::math::constants::half_pi; + + + //std::cout << std::setprecision(5*std::numeric_limits::digits10); + + // Get the party started with two estimates of the integral: + Real min_abscissa{ 0 }, max_abscissa{ boost::math::tools::max_value() }; + K I0 = zero; + Real L1_I0 = 0; + for(size_t i = 0; i < m_abscissas[0].size(); ++i) + { + K y = f(m_abscissas[0][i]); + K I0_last = I0; + I0 += y*m_weights[0][i]; + L1_I0 += magnitude(y)*m_weights[0][i]; + if ((magnitude(I0_last - I0) == 0) && (magnitude(I0) != 0)) + { + max_abscissa = m_abscissas[0][i]; + break; + } + } + + //std::cout << "First estimate : " << I0 << std::endl; + K I1 = I0; + Real L1_I1 = L1_I0; + bool have_first_j = false; + std::size_t first_j = 0; + for (size_t i = 0; (i < m_abscissas[1].size()) && (m_abscissas[1][i] < max_abscissa); ++i) + { + K y = f(m_abscissas[1][i]); + K I1_last = I1; + I1 += y*m_weights[1][i]; + L1_I1 += magnitude(y)*m_weights[1][i]; + if (!have_first_j && (magnitude(I1_last - I1) == 0)) + { + // No change to the sum, disregard these values on the LHS: + if ((i < m_abscissas[1].size() - 1) && (m_abscissas[1][i + 1] > max_abscissa)) + { + // The summit is so high, that we found nothing in this row which added to the integral!! + have_first_j = true; + } + else + { + min_abscissa = m_abscissas[1][i]; + first_j = i; + } + } + else + have_first_j = true; + } + + if (magnitude(I0) == 0) + { + // We failed to find anything, is the integral zero, or have we just not found it yet? + // We'll try one more level, if that still finds nothing then it'll terminate. + min_abscissa = 0; + max_abscissa = boost::math::tools::max_value(); + } + + I1 *= half(); + L1_I1 *= half(); + Real err = magnitude(I0 - I1); + //std::cout << "Second estimate: " << I1 << " Error estimate at level " << 1 << " = " << err << std::endl; + + size_t i = 2; + for(; i < m_abscissas.size(); ++i) + { + I0 = I1; + L1_I0 = L1_I1; + + I1 = half()*I0; + L1_I1 = half()*L1_I0; + Real h = static_cast(1)/static_cast(1 << i); + K sum = zero; + Real absum = 0; + + auto abscissas_row = get_abscissa_row(i); + auto weight_row = get_weight_row(i); + + first_j = first_j == 0 ? 0 : 2 * first_j - 1; // appoximate location to start looking for lowest meaningful abscissa value + std::size_t j = first_j; + while (abscissas_row[j] < min_abscissa) + ++j; + for(; (j < m_weights[i].size()) && (abscissas_row[j] < max_abscissa); ++j) + { + Real x = abscissas_row[j]; + K y = f(x); + sum += y*weight_row[j]; + Real abterm0 = magnitude(y)*weight_row[j]; + absum += abterm0; + } + + I1 += sum*h; + L1_I1 += absum*h; + err = magnitude(I0 - I1); + //std::cout << "Estimate: " << I1 << " Error estimate at level " << i << " = " << err << std::endl; + // Use L1_I1 here to make it work with both complex and real valued integrands: + if (!(boost::math::isfinite)(L1_I1)) + { + policies::raise_evaluation_error(function, "The exp_sinh quadrature evaluated your function at a singular point and returned %1%. Please ensure your function evaluates to a finite number over its entire domain.", magnitude(I1), Policy()); + return I1; + } + if (err <= tolerance*L1_I1) + { + break; + } + } + + if (error) + { + *error = err; + } + + if(L1) + { + *L1 = L1_I1; + } + + if (levels) + { + *levels = i; + } + + return I1; +} + template void exp_sinh_detail::init(const std::integral_constant&) diff --git a/include/boost/math/quadrature/detail/norm_quadrature_error.hpp b/include/boost/math/quadrature/detail/norm_quadrature_error.hpp new file mode 100644 index 0000000000..72050a0ce7 --- /dev/null +++ b/include/boost/math/quadrature/detail/norm_quadrature_error.hpp @@ -0,0 +1,20 @@ +// Copyright Nick Thompson, 2026 +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +#ifndef BOOST_MATH_QUADRATURE_DETAIL_NORM_QUADRATURE_ERROR_HPP +#define BOOST_MATH_QUADRATURE_DETAIL_NORM_QUADRATURE_ERROR_HPP +#ifndef BOOST_MATH_BUILD_MODULE +#include +#endif +namespace boost { namespace math { namespace quadrature { namespace detail { +template +K norm_quadrature_error(const K& zero, Real invalid, Real* error, Real* l1, std::size_t* levels = nullptr) +{ + if (error) *error = invalid; + if (l1) *l1 = invalid; + if (levels) *levels = 0; + K result = zero * invalid; + return result; +} +}}}} +#endif diff --git a/include/boost/math/quadrature/detail/sinh_sinh_detail.hpp b/include/boost/math/quadrature/detail/sinh_sinh_detail.hpp index b5c257e5af..91fdd2ecbb 100644 --- a/include/boost/math/quadrature/detail/sinh_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/sinh_sinh_detail.hpp @@ -11,6 +11,7 @@ #include #ifndef BOOST_MATH_HAS_NVRTC +#include #ifndef BOOST_MATH_BUILD_MODULE #include @@ -56,6 +57,8 @@ class sinh_sinh_detail template auto integrate(const F f, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval()(std::declval())); + template + auto integrate(const F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval()(std::declval())); private: @@ -283,6 +286,137 @@ auto sinh_sinh_detail::integrate(const F f, Real tolerance, Real* return I1; } +template +template +auto sinh_sinh_detail::integrate(const F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(std::declval()(std::declval())) +{ + const auto magnitude = [&](const decltype(std::declval()(std::declval()))& value) -> Real { return static_cast(norm(value)); }; + using std::abs; + using std::sqrt; + using boost::math::constants::half; + using boost::math::constants::half_pi; + + static const char* function = "boost::math::quadrature::sinh_sinh<%1%>::integrate"; + + typedef decltype(f(static_cast(0))) K; + static_assert(!std::is_integral::value, + "The return type cannot be integral."); + K y_max = f(boost::math::tools::max_value()); + if(magnitude(y_max) > boost::math::tools::epsilon()) + { + return norm_quadrature_error(zero, policies::raise_domain_error(function, + "The function you are trying to integrate does not go to zero at infinity, and instead evaluates to %1%", magnitude(y_max), Policy()), error, L1, levels); + } + + K y_min = f(-boost::math::tools::max_value()); + if(magnitude(y_min) > boost::math::tools::epsilon()) + { + return norm_quadrature_error(zero, policies::raise_domain_error(function, + "The function you are trying to integrate does not go to zero at -infinity, and instead evaluates to %1%", magnitude(y_min), Policy()), error, L1, levels); + } + + // Get the party started with two estimates of the integral: + K I0 = f(0)*half_pi(); + Real L1_I0 = magnitude(I0); + for(size_t i = 0; i < m_abscissas[0].size(); ++i) + { + Real x = m_abscissas[0][i]; + K yp = f(x); + K ym = f(-x); + I0 += (yp + ym)*m_weights[0][i]; + L1_I0 += (magnitude(yp)+magnitude(ym))*m_weights[0][i]; + } + + // Uncomment the estimates to work the convergence on the command line. + // std::cout << std::setprecision(std::numeric_limits::digits10); + // std::cout << "First estimate : " << I0 << std::endl; + K I1 = I0; + Real L1_I1 = L1_I0; + for (size_t i = 0; i < m_abscissas[1].size(); ++i) + { + Real x= m_abscissas[1][i]; + K yp = f(x); + K ym = f(-x); + I1 += (yp + ym)*m_weights[1][i]; + L1_I1 += (magnitude(yp) + magnitude(ym))*m_weights[1][i]; + } + + I1 *= half(); + L1_I1 *= half(); + Real err = magnitude(I0 - I1); + // std::cout << "Second estimate: " << I1 << " Error estimate at level " << 1 << " = " << err << std::endl; + + size_t i = 2; + for(; i <= m_max_refinements; ++i) + { + I0 = I1; + L1_I0 = L1_I1; + + I1 = half()*I0; + L1_I1 = half()*L1_I0; + Real h = static_cast(1) / static_cast(1 << i); + K sum = zero; + Real absum = 0; + + Real abterm1 = 1; + Real eps = boost::math::tools::epsilon()*L1_I1; + + auto abscissa_row = get_abscissa_row(i); + auto weight_row = get_weight_row(i); + + for(size_t j = 0; j < abscissa_row.size(); ++j) + { + Real x = abscissa_row[j]; + K yp = f(x); + K ym = f(-x); + sum += (yp + ym)*weight_row[j]; + Real abterm0 = (magnitude(yp) + magnitude(ym))*weight_row[j]; + absum += abterm0; + + // We require two consecutive terms to be < eps in case we hit a zero of f. + if (x > static_cast(100) && abterm0 < eps && abterm1 < eps) + { + break; + } + abterm1 = abterm0; + } + + I1 += sum*h; + L1_I1 += absum*h; + err = magnitude(I0 - I1); + // std::cout << "Estimate: " << I1 << " Error estimate at level " << i << " = " << err << std::endl; + if (!(boost::math::isfinite)(L1_I1)) + { + const char* err_msg = "The sinh_sinh quadrature evaluated your function at a singular point, leading to the value %1%.\n" + "sinh_sinh quadrature cannot handle singularities in the domain.\n" + "If you are sure your function has no singularities, please submit a bug against boost.math\n"; + policies::raise_evaluation_error(function, err_msg, magnitude(I1), Policy()); + return I1; + } + if (err <= tolerance*L1_I1) + { + break; + } + } + + if (error) + { + *error = err; + } + + if (L1) + { + *L1 = L1_I1; + } + + if (levels) + { + *levels = i; + } + + return I1; +} + template void sinh_sinh_detail::init(const std::integral_constant&) { diff --git a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp index e1e3ef625e..b51e7aa2f5 100644 --- a/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp +++ b/include/boost/math/quadrature/detail/tanh_sinh_detail.hpp @@ -16,6 +16,7 @@ #include #include #include +#include #ifdef BOOST_MATH_HAS_THREADS #ifndef BOOST_MATH_BUILD_MODULE @@ -54,6 +55,8 @@ class tanh_sinh_detail template decltype(std::declval()(std::declval(), std::declval())) integrate(const F f, Real* error, Real* L1, const char* function, Real left_min_complement, Real right_min_complement, Real tolerance, std::size_t* levels) const; + template + decltype(std::declval()(std::declval(), std::declval())) integrate(const F f, const decltype(std::declval()(std::declval(), std::declval()))& zero, Norm norm, Real* error, Real* L1, const char* function, Real left_min_complement, Real right_min_complement, Real tolerance, std::size_t* levels) const; private: const std::vector& get_abscissa_row(std::size_t n)const @@ -561,6 +564,383 @@ decltype(std::declval()(std::declval(), std::declval())) tanh_sin return I1; } +template +template +decltype(std::declval()(std::declval(), std::declval())) tanh_sinh_detail::integrate(const F f, const decltype(std::declval()(std::declval(), std::declval()))& zero, Norm norm, Real* error, Real* L1, const char* function, Real left_min_complement, Real right_min_complement, Real tolerance, std::size_t* levels) const +{ + const auto magnitude = [&](const decltype(std::declval()(std::declval(), std::declval()))& value) -> Real { return static_cast(norm(value)); }; + using std::abs; + using std::fabs; + using std::floor; + using std::tanh; + using std::sinh; + using std::sqrt; + using boost::math::constants::half; + using boost::math::constants::half_pi; + + // + // The type of the result: + typedef decltype(std::declval()(std::declval(), std::declval())) result_type; + + Real h = m_t_max / m_inital_row_length; + result_type I0 = half_pi() * f(0, 1); + Real L1_I0 = magnitude(I0); + // + // We maintain 4 integer values: + // max_left_position is the logical index of the abscissa value closest to the + // left endpoint of the range that we can call f(x_i) on without rounding error + // inside f(x_i) causing evaluation at the endpoint. + // max_left_index is the actual position in the current row that has a logical index + // no higher than max_left_position. Remember that since we only store odd numbered + // indexes in each row, this may actually be one position to the left of max_left_position + // in the case that is even. Then, if we only evaluate f(-x_i) for abscissa values + // i <= max_left_index we will never evaluate f(-x_i) at the left endpoint. + // max_right_position and max_right_index are defined similarly for the right boundary + // and are used to guard evaluation of f(x_i). + // + // max_left_position and max_right_position start off as the last element in row zero: + // + std::size_t max_left_position(m_abscissas[0].size() - 1); + std::size_t max_left_index, max_right_position(max_left_position), max_right_index; + // + // Decrement max_left_position and max_right_position until the complement + // of the abscissa value is greater than the smallest permitted (as specified + // by the function caller): + // + while ((max_left_position > 1) && fabs(m_abscissas[0][max_left_position]) < left_min_complement) + --max_left_position; + while ((max_right_position > 1) && fabs(m_abscissas[0][max_right_position]) < right_min_complement) + --max_right_position; + // + // Check for non-finite values at the end points: + // + result_type yp{ f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]) }; + result_type ym{ f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]) }; + Real tail_tolerance{ (std::max)(boost::math::tools::epsilon(), Real(tolerance * tolerance)) }; + while (max_left_position) + { + if ((boost::math::isfinite)(magnitude(yp))) + break; + --max_left_position; + yp = f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]); + } + // + // Also remove points which are insignificant or zero: + // + while (max_left_position > 1) + { + if (magnitude(yp * m_weights[0][max_left_position]) > abs(L1_I0 * tail_tolerance)) + break; + --max_left_position; + yp = f(-1 - m_abscissas[0][max_left_position], m_abscissas[0][max_left_position]); + } + // + // Over again for the right hand side: + // + while (max_right_position) + { + if ((boost::math::isfinite)(magnitude(ym))) + break; + --max_right_position; + ym = f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]); + } + while (max_right_position > 1) + { + if (magnitude(ym * m_weights[0][max_right_position]) > abs(L1_I0 * tail_tolerance)) + break; + --max_right_position; + ym = f(1 + m_abscissas[0][max_right_position], -m_abscissas[0][max_right_position]); + } + + if ((max_left_position == 0) || (max_right_position == 0)) + { + policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", magnitude(ym), Policy()); + return ym; + } + + I0 += yp * m_weights[0][max_left_position] + ym * m_weights[0][max_right_position]; + L1_I0 += magnitude(yp * m_weights[0][max_left_position]) + magnitude(ym * m_weights[0][max_right_position]); + // + // Assumption: left_min_complement/right_min_complement are sufficiently small that we only + // ever decrement through the stored values that are complements (the negative ones), and + // never ever hit the true abscissa values (positive stored values). + // + BOOST_MATH_ASSERT(m_abscissas[0][max_left_position] < 0); + BOOST_MATH_ASSERT(m_abscissas[0][max_right_position] < 0); + + for(size_t i = 1; i < m_abscissas[0].size(); ++i) + { + if ((i >= max_right_position) && (i >= max_left_position)) + break; + Real x = m_abscissas[0][i]; + Real xc = x; + Real w = m_weights[0][i]; + if ((boost::math::signbit)(x)) + { + // We have stored x - 1: + x = 1 + xc; + } + else + xc = x - 1; + yp = i < max_right_position ? f(x, -xc) : zero; + ym = i < max_left_position ? f(-x, xc) : zero; + I0 += (yp + ym)*w; + L1_I0 += (magnitude(yp) + magnitude(ym))*w; + } + // + // We have: + // k = current row. + // I0 = last integral value. + // I1 = current integral value. + // L1_I0 and L1_I1 are the absolute integral values. + // + size_t k = 1; + result_type I1 = I0; + Real L1_I1 = L1_I0; + Real err = 0; + // + // thrash_count is a heuristic - it counts how many time the error has actually increased + // rather than decreased, if this gets too high we abort... + // + unsigned thrash_count = 0; + + while (k < 4 || (k < m_weights.size() && k < m_max_refinements) ) + { + I0 = I1; + L1_I0 = L1_I1; + + I1 = half()*I0; + L1_I1 = half()*L1_I0; + h *= half(); + result_type sum = zero; + Real absum = 0; + Real endpoint_error = 0; + auto const& abscissa_row = this->get_abscissa_row(k); + auto const& weight_row = this->get_weight_row(k); + std::size_t first_complement_index = this->get_first_complement_index(k); + // + // At the start of each new row we need to update the max left/right indexes + // at which we can evaluate f(x_i). The new logical position is simply twice + // the old value. The new max index is one position to the left of the new + // logical value (remember each row contains only odd numbered positions). + // Then we have to make a single check, to see if one position to the right + // is also in bounds (this is the new abscissa value in this row which is + // known to be in between a value known to be in bounds, and one known to be + // not in bounds). + // Thus, we filter which abscissa values generate a call to f(x_i), with a single + // floating point comparison per loop. Everything else is integer logic. + // + BOOST_MATH_ASSERT(max_left_position); + max_left_index = max_left_position - 1; + max_left_position *= 2; + BOOST_MATH_ASSERT(max_right_position); + max_right_index = max_right_position - 1; + max_right_position *= 2; + if ((abscissa_row.size() > max_left_index + 1) && (fabs(abscissa_row[max_left_index + 1]) > left_min_complement)) + { + ++max_left_position; + ++max_left_index; + } + if ((abscissa_row.size() > max_right_index + 1) && (fabs(abscissa_row[max_right_index + 1]) > right_min_complement)) + { + ++max_right_position; + ++max_right_index; + } + // + // We also check that our endpoints don't hit singularities: + // + do + { + yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]); + if ((boost::math::isfinite)(magnitude(yp))) + break; + if(max_left_position <= 2) + { + policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", magnitude(ym), Policy()); + return ym; + } + max_left_position -= 2; + --max_left_index; + } while (abscissa_row[max_left_index] < 0); + bool truncate_left(false), truncate_right(false); + if (abs(L1_I1 * tail_tolerance) > magnitude(yp * weight_row[max_left_index])) + truncate_left = true; + do + { + ym = f(1 + abscissa_row[max_right_index], -abscissa_row[max_right_index]); + if ((boost::math::isfinite)(magnitude(ym))) + break; + if (max_right_position <= 2) + { + policies::raise_evaluation_error(function, "The tanh_sinh quadrature found your function to be non-finite everywhere! Please check your function for singularities.", magnitude(ym), Policy()); + return ym; + } + --max_right_index; + max_right_position -= 2; + } while (abscissa_row[max_right_index] < 0); + if (abs(L1_I1 * tail_tolerance) > magnitude(ym * weight_row[max_right_index])) + truncate_right = true; + + sum += yp * weight_row[max_left_index] + ym * weight_row[max_right_index]; + absum += magnitude(yp * weight_row[max_left_index]) + magnitude(ym * weight_row[max_right_index]); + // + // We estimate the error due to truncation as the value contributed by the two most extreme points. + // In most cases this is tiny and can be ignored, if it is significant then either the area of the + // integral is so far our in the tails that our exponent range can't reach it (example x^-p at double + // precision and p ~ 1), or our function is truncated near epsilon, and we have had to narrow our endpoints. + // In this latter case we may over-estimate the error, but this is the best we can do. + // In any event, we do not add endpoint_error to the error estimate until we terminate the main loop, + // otherwise it can make things appear to be non-converged, when in reality, they are as converged as they + // will ever be. + // + endpoint_error = absum; + + for(size_t j = 0; j < weight_row.size(); ++j) + { + // If both left and right abscissa values are out of bounds at this step + // we can just stop this loop right now: + if ((j >= max_left_index) && (j >= max_right_index)) + break; + Real x = abscissa_row[j]; + Real xc = x; + Real w = weight_row[j]; + if (j >= first_complement_index) + { + // We have stored x - 1: + BOOST_MATH_ASSERT(x < 0); + x = 1 + xc; + } + else + { + BOOST_MATH_ASSERT(x >= 0); + xc = x - 1; + } + + yp = j >= max_right_index ? zero : f(x, -xc); + ym = j >= max_left_index ? zero : f(-x, xc); + result_type term = (yp + ym)*w; + sum += term; + + // A question arises as to how accurately we actually need to estimate the L1 integral. + // For simple integrands, computing the L1 norm makes the integration 20% slower, + // but for more complicated integrands, this calculation is not noticeable. + Real abterm = (magnitude(yp) + magnitude(ym))*w; + absum += abterm; + } + + I1 += sum*h; + L1_I1 += absum*h; + + ++k; + Real last_err = err; + err = magnitude(I0 - I1); + + if (!(boost::math::isfinite)(magnitude(I1))) + { + policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Please narrow the bounds of integration or check your function for singularities.", magnitude(I1), Policy()); + return I1; + } + // + // If the error is increasing, and we're past level 4, something bad is very likely happening: + // + if ((err * 1.5 > last_err) && (k > 4)) + { + bool terminate = false; + if ((++thrash_count > 1) && (last_err < 1e-3)) + // Probably just thrashing, abort: + terminate = true; + else if(thrash_count > 2) + // OK, terrible error, but giving up anyway! + terminate = true; + else if (last_err < boost::math::tools::root_epsilon()) + // Trying to squeeze precision that probably isn't there, abort: + terminate = true; + else + { + // Take a look at the end points, if there's significant new area being + // discovered, then we're not able to get close enough to the endpoints + // to ever find the integral: + if (endpoint_error * h > err) + terminate = true; + } + + if (terminate) + { + // We could raise an evaluation_error, but since we likely have some sort of result, just return the last one + // (ie before the error started going up) + I1 = I0; + L1_I1 = L1_I0; + --k; + err = last_err + endpoint_error * h; + break; + } + // Fall through and keep going, assume we've discovered a new feature of f(x).... + } + // + // Termination condition: + // No more levels are considered once the error is less than the specified tolerance. + // Note however, that we always go down at least 4 levels, otherwise we risk missing + // features of interest in f() - imagine for example a function which flatlines, except + // for a very small "spike". An example would be the incomplete beta integral with large + // parameters. We could keep hunting until we find something, but that would handicap + // integrals which really are zero.... so a compromise then! + // + if ((err <= abs(tolerance*L1_I1)) && (k >= 4)) + { + // + // A quick sanity check: have we at some point narrowed our boundaries as a result + // of non-finite values? If so let's check that the area isn't on an increasing + // trajectory at our new end point, and increase our error estimate by the last + // good value as an estimate for what we may have discarded. + // + if (max_left_index && (max_left_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_left_index + 1]) > left_min_complement)) + { + yp = f(-1 - abscissa_row[max_left_index], abscissa_row[max_left_index]) * weight_row[max_left_index]; + ym = f(-1 - abscissa_row[max_left_index - 1], abscissa_row[max_left_index - 1]) * weight_row[max_left_index - 1]; + if (magnitude(yp) > magnitude(ym)) + { + policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Integration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", magnitude(I1), Policy()); + return I1; + } + } + if (max_right_index && (max_right_index < abscissa_row.size() - 1) && (abs(abscissa_row[max_right_index + 1]) > right_min_complement)) + { + yp = f(1 + abscissa_row[max_right_index], -abscissa_row[max_right_index]) * weight_row[max_right_index]; + ym = f(1 + abscissa_row[max_right_index - 1], -abscissa_row[max_right_index - 1]) * weight_row[max_right_index - 1]; + if (magnitude(yp) > magnitude(ym)) + { + policies::raise_evaluation_error(function, "The tanh_sinh quadrature evaluated your function at a singular point and got %1%. Integration bounds were automatically narrowed, but the integral was found to be increasing at the new endpoint. Please check your function, and consider providing a 2-argument functor.", magnitude(I1), Policy()); + return I1; + } + } + err += endpoint_error * h; + break; + } + + if ((truncate_left) && (max_left_position > 1)) + --max_left_position; + if ((truncate_right) && (max_right_position > 1)) + --max_right_position; + + } + if (error) + { + *error = err; + } + + if (L1) + { + *L1 = L1_I1; + } + + if (levels) + { + *levels = k; + } + + return I1; +} + template void tanh_sinh_detail::init(const Real& min_complement, const std::integral_constant&) { diff --git a/include/boost/math/quadrature/exp_sinh.hpp b/include/boost/math/quadrature/exp_sinh.hpp index 7eae7249c0..3e7f939394 100644 --- a/include/boost/math/quadrature/exp_sinh.hpp +++ b/include/boost/math/quadrature/exp_sinh.hpp @@ -41,6 +41,11 @@ class exp_sinh template auto integrate(const F& f, Real tol = boost::math::tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval()(std::declval())); + template + auto integrate(const F& f, Real a, Real b, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(static_cast(norm(std::declval()(std::declval()))), std::declval()(std::declval())); + template + auto integrate(const F& f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(static_cast(norm(std::declval()(std::declval()))), std::declval()(std::declval())); + private: std::shared_ptr> m_imp; }; @@ -103,6 +108,63 @@ auto exp_sinh::integrate(const F& f, Real tolerance, Real* error, } +template +template +auto exp_sinh::integrate(const F& f, Real a, Real b, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(static_cast(norm(std::declval()(std::declval()))), std::declval()(std::declval())) +{ + typedef decltype(f(a)) K; + static_assert(!std::is_integral::value, + "The return type cannot be integral, it must be either a real or complex floating point type."); + using std::abs; + using boost::math::constants::half; + using boost::math::quadrature::detail::exp_sinh_detail; + + static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate"; + + // Neither limit may be a NaN: + if((boost::math::isnan)(a) || (boost::math::isnan)(b)) + { + return detail::norm_quadrature_error(zero, policies::raise_domain_error(function, "NaN supplied as one limit of integration - sorry I don't know what to do", a, Policy()), error, L1, levels); + } + // Right limit is infinite: + if ((boost::math::isfinite)(a) && (b >= boost::math::tools::max_value())) + { + // If a = 0, don't use an additional level of indirection: + if (a == static_cast(0)) + { + return m_imp->integrate(f, zero, norm, error, L1, function, tolerance, levels); + } + const auto u = [&](Real t)->K { return f(t + a); }; + return m_imp->integrate(u, zero, norm, error, L1, function, tolerance, levels); + } + + if ((boost::math::isfinite)(b) && a <= -boost::math::tools::max_value()) + { + const auto u = [&](Real t)->K { return f(b-t);}; + return m_imp->integrate(u, zero, norm, error, L1, function, tolerance, levels); + } + + // Infinite limits: + if ((a <= -boost::math::tools::max_value()) && (b >= boost::math::tools::max_value())) + { + return detail::norm_quadrature_error(zero, policies::raise_domain_error(function, "Use sinh_sinh quadrature for integration over the whole real line; exp_sinh is for half infinite integrals.", a, Policy()), error, L1, levels); + } + // If we get to here then both ends must necessarily be finite: + return detail::norm_quadrature_error(zero, policies::raise_domain_error(function, "Use tanh_sinh quadrature for integration over finite domains; exp_sinh is for half infinite integrals.", a, Policy()), error, L1, levels); +} + +template +template +auto exp_sinh::integrate(const F& f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(static_cast(norm(std::declval()(std::declval()))), std::declval()(std::declval())) +{ + static const char* function = "boost::math::quadrature::exp_sinh<%1%>::integrate"; + using std::abs; + if (abs(tolerance) > 1) { + return detail::norm_quadrature_error(zero, policies::raise_domain_error(function, "The tolerance provided (%1%) is unusually large; did you confuse it with a domain bound?", tolerance, Policy()), error, L1, levels); + } + return m_imp->integrate(f, zero, norm, error, L1, function, tolerance, levels); +} + }}} #endif // BOOST_MATH_HAS_NVRTC diff --git a/include/boost/math/quadrature/gauss.hpp b/include/boost/math/quadrature/gauss.hpp index 414c9bef4c..a85dc2ea77 100644 --- a/include/boost/math/quadrature/gauss.hpp +++ b/include/boost/math/quadrature/gauss.hpp @@ -893,6 +893,130 @@ class gauss : public detail::gauss_detail(policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy())); } + + // Explicit zero and norm for vector- and matrix-valued integrands. + // The original overloads above retain their existing behavior. + template + static auto integrate(F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* pL1 = nullptr) + ->decltype(static_cast(norm(f(Real(0)))), f(Real(0))) + { + typedef decltype(f(Real(0))) K; + static_assert(!std::is_integral::value, + "The return type cannot be integral."); + unsigned non_zero_start = 1; + K result = zero; + if (N & 1) { + result = f(Real(0)) * static_cast(base::weights()[0]); + } + else { + result = zero; + non_zero_start = 0; + } + Real L1 = static_cast(norm(result)); + for (unsigned i = non_zero_start; i < base::abscissa().size(); ++i) + { + K fp = f(static_cast(base::abscissa()[i])); + K fm = f(static_cast(-base::abscissa()[i])); + result += (fp + fm) * static_cast(base::weights()[i]); + L1 += (static_cast(norm(fp)) + static_cast(norm(fm))) * static_cast(base::weights()[i]); + } + if (pL1) + *pL1 = L1; + return result; + } + template + static auto integrate(F f, Real a, Real b, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* pL1 = nullptr) + ->decltype(static_cast(norm(f(Real(0)))), f(Real(0))) + { + typedef decltype(f(a)) K; + static const char* function = "boost::math::quadrature::gauss<%1%>::integrate(f, %1%, %1%)"; + if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b)) + { + // Infinite limits: + Real min_inf = -tools::max_value(); + if ((a <= min_inf) && (b >= tools::max_value())) + { + auto u = [&](const Real& t)->K + { + Real t_sq = t*t; + Real inv = 1 / (1 - t_sq); + K res = f(t*inv)*(1 + t_sq)*inv*inv; + return res; + }; + return integrate(u, zero, norm, pL1); + } + + // Right limit is infinite: + if ((boost::math::isfinite)(a) && (b >= tools::max_value())) + { + auto u = [&](const Real& t)->K + { + Real z = 1 / (t + 1); + Real arg = 2 * z + a - 1; + K res = f(arg)*z*z; + return res; + }; + K Q = Real(2) * integrate(u, zero, norm, pL1); + if (pL1) + { + *pL1 *= 2; + } + return Q; + } + + if ((boost::math::isfinite)(b) && (a <= -tools::max_value())) + { + auto v = [&](const Real& t)->K + { + Real z = 1 / (t + 1); + Real arg = 2 * z - 1; + K res = f(b - arg) * z * z; + return res; + }; + K Q = Real(2) * integrate(v, zero, norm, pL1); + if (pL1) + { + *pL1 *= 2; + } + return Q; + } + + if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) + { + if (a == b) + { + if (pL1) + *pL1 = Real(0); + return zero; + } + if (b < a) + { + return -integrate(f, b, a, zero, norm, pL1); + } + Real avg = (a + b)*constants::half(); + Real scale = (b - a)*constants::half(); + + auto u = [&](Real z)->K + { + return f(avg + scale*z); + }; + K Q = scale*integrate(u, zero, norm, pL1); + + if (pL1) + { + *pL1 *= scale; + } + return Q; + } + } + // Preserve the scalar error policy, including non-throwing policies, + // and use scalar multiplication to produce a result with the right shape. + Real error = policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()); + if (pL1) + *pL1 = error; + K result = zero * error; + return result; + } }; } // namespace quadrature diff --git a/include/boost/math/quadrature/gauss_kronrod.hpp b/include/boost/math/quadrature/gauss_kronrod.hpp index b4d6f1f9f8..05a2066ac4 100644 --- a/include/boost/math/quadrature/gauss_kronrod.hpp +++ b/include/boost/math/quadrature/gauss_kronrod.hpp @@ -1306,6 +1306,196 @@ class gauss_kronrod : public detail::gauss_kronrod_detail(policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy())); } + +private: + // Explicit zero/norm path; preserve the original scalar/complex path above. + template + static auto integrate_non_adaptive_m1_1(F f, const decltype(f(Real(0)))& zero, Norm norm, Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval()(std::declval())) + { + typedef decltype(f(Real(0))) K; + unsigned gauss_start = 2; + unsigned kronrod_start = 1; + unsigned gauss_order = (N - 1) / 2; + K kronrod_result = zero; + K gauss_result = zero; + K fp = zero, fm = zero; + if (gauss_order & 1) + { + fp = f(value_type(0)); + kronrod_result = fp * static_cast(base::weights()[0]); + gauss_result += fp * static_cast(gauss::weights()[0]); + } + else + { + fp = f(value_type(0)); + kronrod_result = fp * static_cast(base::weights()[0]); + gauss_start = 1; + kronrod_start = 2; + } + Real L1 = static_cast(norm(kronrod_result)); + for (unsigned i = gauss_start; i < base::abscissa().size(); i += 2) + { + fp = f(static_cast(base::abscissa()[i])); + fm = f(static_cast(-base::abscissa()[i])); + kronrod_result += (fp + fm) * static_cast(base::weights()[i]); + L1 += (static_cast(norm(fp)) + static_cast(norm(fm))) * static_cast(base::weights()[i]); + gauss_result += (fp + fm) * static_cast(gauss::weights()[i / 2]); + } + for (unsigned i = kronrod_start; i < base::abscissa().size(); i += 2) + { + fp = f(static_cast(base::abscissa()[i])); + fm = f(static_cast(-base::abscissa()[i])); + kronrod_result += (fp + fm) * static_cast(base::weights()[i]); + L1 += (static_cast(norm(fp)) + static_cast(norm(fm))) * static_cast(base::weights()[i]); + } + if (pL1) + *pL1 = L1; + if (error) + *error = (std::max)(static_cast(norm(kronrod_result - gauss_result)), static_cast(norm(kronrod_result * tools::epsilon() * Real(2)))); + return kronrod_result; + } + + template + struct norm_recursive_info + { + F f; + Real tol; + const decltype(std::declval()(std::declval()))& zero; + Norm norm; + }; + + template + static auto recursive_adaptive_integrate(const norm_recursive_info* info, Real a, Real b, unsigned max_levels, Real abs_tol, Real* error, Real* L1)->decltype(std::declval()(std::declval())) + { + typedef decltype(info->f(Real(a))) K; + Real error_local; + Real mean = (b + a) / 2; + Real scale = (b - a) / 2; + auto ff = [&](const Real& x)->K + { + return info->f(scale * x + mean); + }; + K r1 = integrate_non_adaptive_m1_1(ff, info->zero, info->norm, &error_local, L1); + K estimate = scale * r1; + error_local *= scale; + + K tmp = estimate * info->tol; + Real abs_tol1 = static_cast(info->norm(tmp)); + if (abs_tol == 0) + abs_tol = abs_tol1; + + if (max_levels && (abs_tol1 < error_local) && (abs_tol < error_local)) + { + Real mid = (a + b) / 2; + Real L1_local; + estimate = recursive_adaptive_integrate(info, a, mid, max_levels - 1, abs_tol / 2, error, L1); + estimate += recursive_adaptive_integrate(info, mid, b, max_levels - 1, abs_tol / 2, &error_local, &L1_local); + if (error) + *error += error_local; + if (L1) + *L1 += L1_local; + return estimate; + } + if(L1) + *L1 *= scale; + if (error) + *error = error_local; + return estimate; + } + +public: + template + static auto integrate(F f, Real a, Real b, const decltype(f(a))& zero, Norm norm, unsigned max_depth = 15, Real tol = tools::root_epsilon(), Real* error = nullptr, Real* pL1 = nullptr)->decltype(static_cast(norm(f(a))), f(a)) + { + typedef decltype(f(a)) K; + static_assert(!std::is_integral::value, + "The return type cannot be integral."); + static const char* function = "boost::math::quadrature::gauss_kronrod<%1%>::integrate(f, %1%, %1%)"; + if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b)) + { + // Infinite limits: + if ((a <= -tools::max_value()) && (b >= tools::max_value())) + { + auto u = [&](const Real& t)->K + { + Real t_sq = t*t; + Real inv = 1 / (1 - t_sq); + Real w = (1 + t_sq)*inv*inv; + Real arg = t*inv; + K res = f(arg)*w; + return res; + }; + norm_recursive_info info = { u, tol, zero, norm }; + K res = recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + return res; + } + + // Right limit is infinite: + if ((boost::math::isfinite)(a) && (b >= tools::max_value())) + { + auto u = [&](const Real& t)->K + { + Real z = 1 / (t + 1); + Real arg = 2 * z + a - 1; + K res = f(arg)*z*z; + return res; + }; + norm_recursive_info info = { u, tol, zero, norm }; + K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + if (pL1) + { + *pL1 *= 2; + } + if (error) + *error *= 2; + return Q; + } + + if ((boost::math::isfinite)(b) && (a <= -tools::max_value())) + { + auto v = [&](const Real& t)->K + { + Real z = 1 / (t + 1); + Real arg = 2 * z - 1; + return f(b - arg) * z * z; + }; + norm_recursive_info info = { v, tol, zero, norm }; + K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1); + if (pL1) + { + *pL1 *= 2; + } + if (error) + *error *= 2; + return Q; + } + + if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) + { + if (a==b) + { + if (error) + *error = Real(0); + if (pL1) + *pL1 = Real(0); + return zero; + } + norm_recursive_info info = { f, tol, zero, norm }; + if (b < a) + { + return -recursive_adaptive_integrate(&info, b, a, max_depth, Real(0), error, pL1); + } + return recursive_adaptive_integrate(&info, a, b, max_depth, Real(0), error, pL1); + } + } + Real invalid = policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()); + if (error) + *error = invalid; + if (pL1) + *pL1 = invalid; + K result = zero * invalid; + return result; + } }; } // namespace quadrature diff --git a/include/boost/math/quadrature/sinh_sinh.hpp b/include/boost/math/quadrature/sinh_sinh.hpp index 12a0640222..26867f19df 100644 --- a/include/boost/math/quadrature/sinh_sinh.hpp +++ b/include/boost/math/quadrature/sinh_sinh.hpp @@ -45,6 +45,15 @@ class sinh_sinh return m_imp->integrate(f, tol, error, L1, levels); } + template + auto integrate(const F f, const decltype(f(Real(0)))& zero, Norm norm, + Real tol = tools::root_epsilon(), Real* error = nullptr, + Real* L1 = nullptr, std::size_t* levels = nullptr) const + ->decltype(static_cast(norm(f(Real(0)))), f(Real(0))) + { + return m_imp->integrate(f, zero, norm, tol, error, L1, levels); + } + private: std::shared_ptr> m_imp; }; diff --git a/include/boost/math/quadrature/tanh_sinh.hpp b/include/boost/math/quadrature/tanh_sinh.hpp index f417b36064..9ac2c156f1 100644 --- a/include/boost/math/quadrature/tanh_sinh.hpp +++ b/include/boost/math/quadrature/tanh_sinh.hpp @@ -55,6 +55,15 @@ class tanh_sinh template auto integrate(const F f, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(std::declval()(std::declval(), std::declval())); + template + auto integrate(const F f, Real a, Real b, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(static_cast(norm(std::declval()(std::declval()))), std::declval()(std::declval())); + template + auto integrate(const F f, Real a, Real b, const decltype(std::declval()(std::declval(), std::declval()))& zero, Norm norm, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(static_cast(norm(std::declval()(std::declval(), std::declval()))), std::declval()(std::declval(), std::declval())); + template + auto integrate(const F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(static_cast(norm(std::declval()(std::declval()))), std::declval()(std::declval())); + template + auto integrate(const F f, const decltype(std::declval()(std::declval(), std::declval()))& zero, Norm norm, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr, std::size_t* levels = nullptr) const ->decltype(static_cast(norm(std::declval()(std::declval(), std::declval()))), std::declval()(std::declval(), std::declval())); + private: std::shared_ptr> m_imp; }; @@ -285,6 +294,242 @@ auto tanh_sinh::integrate(const F f, Real tolerance, Real* error, return m_imp->integrate(f, error, L1, function, min_complement, min_complement, tolerance, levels); } + +template +template +auto tanh_sinh::integrate(const F f, Real a, Real b, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(static_cast(norm(std::declval()(std::declval()))), std::declval()(std::declval())) +{ + BOOST_MATH_STD_USING + using boost::math::constants::half; + using boost::math::quadrature::detail::tanh_sinh_detail; + + static const char* function = "tanh_sinh<%1%>::integrate"; + + typedef decltype(std::declval()(std::declval())) result_type; + static_assert(!std::is_integral::value, + "The return type cannot be integral, it must be either a real or complex floating point type."); + if (!(boost::math::isnan)(a) && !(boost::math::isnan)(b)) + { + + // Infinite limits: + if ((a <= -tools::max_value()) && (b >= tools::max_value())) + { + auto u = [&](const Real& t, const Real& tc)->result_type + { + Real t_sq = t*t; + Real inv; + if (t > 0.5f) + inv = 1 / ((2 - tc) * tc); + else if(t < -0.5) + inv = 1 / ((2 + tc) * -tc); + else + inv = 1 / (1 - t_sq); + return f(t*inv)*(1 + t_sq)*inv*inv; + }; + Real limit = sqrt(tools::min_value()) * 4; + return m_imp->integrate(u, zero, norm, error, L1, function, limit, limit, tolerance, levels); + } + + // Right limit is infinite: + if ((boost::math::isfinite)(a) && (b >= tools::max_value())) + { + auto u = [&](const Real& t, const Real& tc)->result_type + { + Real z, arg; + if (t > -0.5f) + z = 1 / (t + 1); + else + z = -1 / tc; + if (t < 0.5) + arg = 2 * z + a - 1; + else + arg = a + tc / (2 - tc); + return f(arg)*z*z; + }; + Real left_limit = sqrt(tools::min_value()) * 4; + result_type Q = Real(2) * m_imp->integrate(u, zero, norm, error, L1, function, left_limit, tools::min_value(), tolerance, levels); + if (L1) + { + *L1 *= 2; + } + if (error) + { + *error *= 2; + } + + return Q; + } + + if ((boost::math::isfinite)(b) && (a <= -tools::max_value())) + { + auto v = [&](const Real& t, const Real& tc)->result_type + { + Real z; + if (t > -0.5) + z = 1 / (t + 1); + else + z = -1 / tc; + Real arg; + if (t < 0.5) + arg = 2 * z - 1; + else + arg = tc / (2 - tc); + return f(b - arg) * z * z; + }; + + Real left_limit = sqrt(tools::min_value()) * 4; + result_type Q = Real(2) * m_imp->integrate(v, zero, norm, error, L1, function, left_limit, tools::min_value(), tolerance, levels); + if (L1) + { + *L1 *= 2; + } + if (error) + { + *error *= 2; + } + return Q; + } + + if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) + { + if (a == b) + { + if (error) *error = Real(0); + if (L1) *L1 = Real(0); + if (levels) *levels = 0; + return zero; + } + if (b < a) + { + return -this->integrate(f, b, a, zero, norm, tolerance, error, L1, levels); + } + Real avg = (a + b)*half(); + Real diff = (b - a)*half(); + Real avg_over_diff_m1 = a / diff; + Real avg_over_diff_p1 = b / diff; + bool have_small_left = fabs(a) < 0.5f; + bool have_small_right = fabs(b) < 0.5f; + Real left_min_complement = float_next(avg_over_diff_m1) - avg_over_diff_m1; + Real min_complement_limit = (std::max)(tools::min_value(), float_next(Real(tools::min_value() / diff))); + if (left_min_complement < min_complement_limit) + left_min_complement = min_complement_limit; + Real right_min_complement = avg_over_diff_p1 - float_prior(avg_over_diff_p1); + if (right_min_complement < min_complement_limit) + right_min_complement = min_complement_limit; + // + // These asserts will fail only if rounding errors on + // type Real have accumulated so much error that it's + // broken our internal logic. Should that prove to be + // a persistent issue, we might need to add a bit of fudge + // factor to move left_min_complement and right_min_complement + // further from the end points of the range. + // + BOOST_MATH_ASSERT((left_min_complement * diff + a) > a); + BOOST_MATH_ASSERT((b - right_min_complement * diff) < b); + auto u = [&](Real z, Real zc)->result_type + { + Real position; + if (z < -0.5) + { + if(have_small_left) + return f(diff * (avg_over_diff_m1 - zc)); + position = a - diff * zc; + } + else if (z > 0.5) + { + if(have_small_right) + return f(diff * (avg_over_diff_p1 - zc)); + position = b - diff * zc; + } + else + position = avg + diff*z; + BOOST_MATH_ASSERT(position != a); + BOOST_MATH_ASSERT(position != b); + return f(position); + }; + result_type Q = diff*m_imp->integrate(u, zero, norm, error, L1, function, left_min_complement, right_min_complement, tolerance, levels); + + if (L1) + { + *L1 *= diff; + } + if (error) + { + *error *= diff; + } + return Q; + } + } + return detail::norm_quadrature_error(zero, policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()), error, L1, levels); +} + +template +template +auto tanh_sinh::integrate(const F f, Real a, Real b, const decltype(std::declval()(std::declval(), std::declval()))& zero, Norm norm, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(static_cast(norm(std::declval()(std::declval(), std::declval()))), std::declval()(std::declval(), std::declval())) +{ + BOOST_MATH_STD_USING + using boost::math::constants::half; + using boost::math::quadrature::detail::tanh_sinh_detail; + + static const char* function = "tanh_sinh<%1%>::integrate"; + + if ((boost::math::isfinite)(a) && (boost::math::isfinite)(b)) + { + if (a == b) + { + if (error) *error = Real(0); + if (L1) *L1 = Real(0); + if (levels) *levels = 0; + return zero; + } + if (b < a) + { + return -this->integrate(f, b, a, zero, norm, tolerance, error, L1, levels); + } + auto u = [&](Real z, Real zc)->decltype(std::declval()(std::declval(), std::declval())) + { + if (z < 0) + return f((a - b) * zc / 2 + a, (b - a) * zc / 2); + else + return f((a - b) * zc / 2 + b, (b - a) * zc / 2); + }; + Real diff = (b - a)*half(); + Real left_min_complement = tools::min_value() * 4; + Real right_min_complement = tools::min_value() * 4; + decltype(std::declval()(std::declval(), std::declval())) Q = diff*m_imp->integrate(u, zero, norm, error, L1, function, left_min_complement, right_min_complement, tolerance, levels); + + if (L1) + { + *L1 *= diff; + } + if (error) + { + *error *= diff; + } + return Q; + } + return detail::norm_quadrature_error(zero, policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()), error, L1, levels); +} + +template +template +auto tanh_sinh::integrate(const F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(static_cast(norm(std::declval()(std::declval()))), std::declval()(std::declval())) +{ + using boost::math::quadrature::detail::tanh_sinh_detail; + static const char* function = "tanh_sinh<%1%>::integrate"; + Real min_complement = tools::epsilon(); + return m_imp->integrate([&](const Real& arg, const Real&) { return f(arg); }, zero, norm, error, L1, function, min_complement, min_complement, tolerance, levels); +} + +template +template +auto tanh_sinh::integrate(const F f, const decltype(std::declval()(std::declval(), std::declval()))& zero, Norm norm, Real tolerance, Real* error, Real* L1, std::size_t* levels) const ->decltype(static_cast(norm(std::declval()(std::declval(), std::declval()))), std::declval()(std::declval(), std::declval())) +{ + using boost::math::quadrature::detail::tanh_sinh_detail; + static const char* function = "tanh_sinh<%1%>::integrate"; + Real min_complement = tools::min_value() * 4; + return m_imp->integrate(f, zero, norm, error, L1, function, min_complement, min_complement, tolerance, levels); +} } } } diff --git a/include/boost/math/quadrature/trapezoidal.hpp b/include/boost/math/quadrature/trapezoidal.hpp index 7a257f85a0..e1770b6809 100644 --- a/include/boost/math/quadrature/trapezoidal.hpp +++ b/include/boost/math/quadrature/trapezoidal.hpp @@ -22,6 +22,7 @@ #include #endif #include +#include #include #include #include @@ -123,5 +124,102 @@ auto trapezoidal(F f, Real a, Real b, Real tol = boost::math::tools::root_epsilo return trapezoidal(f, a, b, tol, max_refinements, error_estimate, L1, boost::math::policies::policy<>()); } +BOOST_MATH_EXPORT template +auto trapezoidal(F f, Real a, Real b, const decltype(f(a))& zero, Norm norm, Real tol, std::size_t max_refinements, Real* error_estimate, Real* L1, const Policy& pol)->decltype(static_cast(norm(f(a))), f(a)) +{ + static const char* function = "boost::math::quadrature::trapezoidal<%1%>(F, %1%, %1%, %1%)"; + const auto magnitude = [&](const decltype(f(a))& value) -> Real { return static_cast(norm(value)); }; + using boost::math::constants::half; + // In many math texts, K represents the field of real or complex numbers. + // Too bad we can't put blackboard bold into C++ source! + typedef decltype(f(a)) K; + static_assert(!std::is_integral::value, + "The return type cannot be integral, it must be either a real or complex floating point type."); + if (!(boost::math::isfinite)(a)) + { + return detail::norm_quadrature_error(zero, boost::math::policies::raise_domain_error(function, "Left endpoint of integration must be finite for adaptive trapezoidal integration but got a = %1%.\n", a, pol), error_estimate, L1); + } + if (!(boost::math::isfinite)(b)) + { + return detail::norm_quadrature_error(zero, boost::math::policies::raise_domain_error(function, "Right endpoint of integration must be finite for adaptive trapezoidal integration but got b = %1%.\n", b, pol), error_estimate, L1); + } + + if (a == b) + { + if (error_estimate) *error_estimate = Real(0); + if (L1) *L1 = Real(0); + return zero; + } + if(a > b) + { + return -trapezoidal(f, b, a, zero, norm, tol, max_refinements, error_estimate, L1, pol); + } + + + K ya = f(a); + K yb = f(b); + Real h = (b - a)*half(); + K I0 = (ya + yb)*h; + Real IL0 = (magnitude(ya) + magnitude(yb))*h; + + K yh = f(a + h); + K I1 = zero; + I1 = I0*half() + yh*h; + Real IL1 = IL0*half() + magnitude(yh)*h; + + // The recursion is: + // I_k = 1/2 I_{k-1} + 1/2^k \sum_{j=1; j odd, j < 2^k} f(a + j(b-a)/2^k) + std::size_t k = 2; + // We want to go through at least 5 levels so we have sampled the function at least 20 times. + // Otherwise, we could terminate prematurely and miss essential features. + // This is of course possible anyway, but 20 samples seems to be a reasonable compromise. + Real error = magnitude(I0 - I1); + // I take k < 5, rather than k < 4, or some other smaller minimum number, + // because I hit a truly exceptional bug where the k = 2 and k =3 refinement were bitwise equal, + // but the quadrature had not yet converged. + while (k < 5 || (k < max_refinements && error > tol*IL1) ) + { + I0 = I1; + IL0 = IL1; + + I1 = I0*half(); + IL1 = IL0*half(); + std::size_t p = static_cast(1u) << k; + h *= half(); + K sum = zero; + Real absum = 0; + + for(std::size_t j = 1; j < p; j += 2) + { + K y = f(a + j*h); + sum += y; + absum += magnitude(y); + } + + I1 += sum*h; + IL1 += absum*h; + ++k; + error = magnitude(I0 - I1); + } + + if (error_estimate) + { + *error_estimate = error; + } + + if (L1) + { + *L1 = IL1; + } + + return static_cast(I1); +} + +BOOST_MATH_EXPORT template +auto trapezoidal(F f, Real a, Real b, const decltype(f(a))& zero, Norm norm, Real tol = boost::math::tools::root_epsilon(), std::size_t max_refinements = 12, Real* error_estimate = nullptr, Real* L1 = nullptr)->decltype(static_cast(norm(f(a))), f(a)) +{ + return trapezoidal(f, a, b, zero, norm, tol, max_refinements, error_estimate, L1, boost::math::policies::policy<>()); +} + }}} #endif diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index e795dcf17a..15062f4fa9 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1217,6 +1217,15 @@ test-suite quadrature : : : : release TEST10 $(float128_type) [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax ] [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] : exp_sinh_quadrature_test_10 ] + [ run double_exponential_eigen_test.cpp : : : + [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax ] + [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] + [ run gauss_kronrod_eigen_test.cpp : : : + [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax ] + [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] + [ run gauss_quadrature_eigen_test.cpp : : : + [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax ] + [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] [ run gauss_quadrature_test.cpp : : : TEST1 $(float128_type) [ requires cxx11_auto_declarations cxx11_lambdas cxx11_smart_ptr cxx11_unified_initialization_syntax ] off msvc:/bigobj release : gauss_quadrature_test_1 ] [ run gauss_quadrature_test.cpp : : : TEST2 $(float128_type) diff --git a/test/double_exponential_eigen_test.cpp b/test/double_exponential_eigen_test.cpp new file mode 100644 index 0000000000..8782335120 --- /dev/null +++ b/test/double_exponential_eigen_test.cpp @@ -0,0 +1,252 @@ +// Copyright Nick Thompson, 2026 +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +#define BOOST_TEST_MODULE double_exponential_eigen_test +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace q = boost::math::quadrature; +using matrix = Eigen::MatrixXcd; +const matrix zero = matrix::Zero(2, 3); +const matrix value = (matrix(2, 3) << std::complex(1, 2), 2., -3., 4., 0., std::complex(-2, 1)).finished(); +const auto norm = [](const matrix& m) { return m.stableNorm(); }; +const auto scaled_norm = [](const matrix& m) { return 8*m.stableNorm(); }; +const double inf = std::numeric_limits::infinity(); +const double pi = boost::math::constants::pi(); + +void check(const matrix& actual, const matrix& expected, double tolerance = 1e-9) +{ + BOOST_CHECK_EQUAL(actual.rows(), expected.rows()); + BOOST_CHECK_EQUAL(actual.cols(), expected.cols()); + BOOST_CHECK_SMALL((actual-expected).norm(), tolerance); +} + +BOOST_AUTO_TEST_CASE(exp_sinh_matrices) +{ + q::exp_sinh integrator; + unsigned calls = 0; + const auto f = [&](double x) -> matrix { ++calls; return value*std::exp(-std::sqrt(x)); }; + double error = 0, l1 = 0; + std::size_t levels = 0; + matrix coarse = integrator.integrate(f, zero, norm, 1e-3); + const unsigned coarse_calls = calls; + calls = 0; + matrix result = integrator.integrate(f, zero, norm, 1e-11, &error, &l1, &levels); + check(result, value*2); + BOOST_CHECK_GT(calls, coarse_calls); + BOOST_CHECK_GE(levels, 2); + BOOST_CHECK_LT(error, 1e-9); + BOOST_CHECK_CLOSE_FRACTION(l1, norm(value)*2, 1e-10); + const unsigned refined_calls = calls; + calls = 0; + double scaled_error = 0, scaled_l1 = 0; + check(integrator.integrate(f, zero, scaled_norm, 1e-11, &scaled_error, &scaled_l1), result); + BOOST_CHECK_EQUAL(calls, refined_calls); + BOOST_CHECK_CLOSE_FRACTION(scaled_error, 8*error, 1e-10); + BOOST_CHECK_CLOSE_FRACTION(scaled_l1, 8*l1, 1e-10); + check(integrator.integrate(f, 0., inf, zero, norm), value*2); + const auto shifted = [&](double x) -> matrix { return value*std::exp(-(x-2)); }; + check(integrator.integrate(shifted, 2., inf, zero, norm), value); + const auto left = [&](double x) -> matrix { return value*std::exp(x-2); }; + check(integrator.integrate(left, -inf, 2., zero, norm), value); + // Distinct entrywise decay rates have independently known integrals. + const auto mixed = [](double x) -> matrix { + matrix m(2, 3); + for (Eigen::Index i = 0; i < m.rows(); ++i) + for (Eigen::Index j = 0; j < m.cols(); ++j) + m(i, j) = value(i, j)*std::exp(-(1+i+2*j)*x); + return m; + }; + matrix expected(2, 3); + for (Eigen::Index i = 0; i < expected.rows(); ++i) + for (Eigen::Index j = 0; j < expected.cols(); ++j) + expected(i, j) = value(i, j)/double(1+i+2*j); + check(integrator.integrate(mixed, zero, norm), expected); + BOOST_CHECK_THROW(integrator.integrate(f, 0., 1., zero, norm), std::domain_error); +} + +BOOST_AUTO_TEST_CASE(sinh_sinh_matrices) +{ + q::sinh_sinh integrator; + unsigned calls = 0; + const auto f = [&](double x) -> matrix { ++calls; return value*std::exp(-x*x); }; + double error = 0, l1 = 0; + std::size_t levels = 0; + integrator.integrate(f, zero, norm, 1e-3); + const unsigned coarse_calls = calls; + calls = 0; + matrix result = integrator.integrate(f, zero, norm, 1e-11, &error, &l1, &levels); + check(result, value*std::sqrt(pi)); + BOOST_CHECK_GT(calls, coarse_calls); + BOOST_CHECK_GE(levels, 2); + BOOST_CHECK_LT(error, 1e-9); + BOOST_CHECK_CLOSE_FRACTION(l1, norm(value)*std::sqrt(pi), 1e-10); + const unsigned refined_calls = calls; + calls = 0; + double scaled_error = 0, scaled_l1 = 0; + check(integrator.integrate(f, zero, scaled_norm, 1e-11, &scaled_error, &scaled_l1), result); + BOOST_CHECK_EQUAL(calls, refined_calls); + BOOST_CHECK_CLOSE_FRACTION(scaled_error, 8*error, 1e-10); + BOOST_CHECK_CLOSE_FRACTION(scaled_l1, 8*l1, 1e-10); + const auto constant = [&](double) -> matrix { return value; }; + BOOST_CHECK_THROW(integrator.integrate(constant, zero, norm), std::domain_error); +} + +BOOST_AUTO_TEST_CASE(tanh_sinh_matrices) +{ + q::tanh_sinh integrator; + unsigned calls = 0; + // A narrow analytic peak, with a known antiderivative, needs refinement. + const auto f = [&](double x) -> matrix { ++calls; return value/(1+10000*(x-0.37)*(x-0.37)); }; + const double exact = (std::atan(63.)+std::atan(37.))/100; + integrator.integrate(f, 0., 1., zero, norm, 1e-3); + const unsigned coarse_calls = calls; + calls = 0; + double error = 0, l1 = 0; + std::size_t levels = 0; + matrix result = integrator.integrate(f, 0., 1., zero, norm, 1e-11, &error, &l1, &levels); + check(result, value*exact); + BOOST_CHECK_GT(calls, coarse_calls); + BOOST_CHECK_GT(levels, 4); + BOOST_CHECK_LT(error, 1e-9); + BOOST_CHECK_CLOSE_FRACTION(l1, norm(value)*exact, 1e-10); + const unsigned refined_calls = calls; + calls = 0; + double scaled_error = 0, scaled_l1 = 0; + check(integrator.integrate(f, 0., 1., zero, scaled_norm, 1e-11, &scaled_error, &scaled_l1), result); + BOOST_CHECK_EQUAL(calls, refined_calls); + BOOST_CHECK_CLOSE_FRACTION(scaled_error, 8*error, 1e-10); + BOOST_CHECK_CLOSE_FRACTION(scaled_l1, 8*l1, 1e-10); + check(integrator.integrate(f, 1., 0., zero, norm, 1e-11), -result); + const auto smooth = [&](double x) -> matrix { return value*std::exp(x); }; + check(integrator.integrate(smooth, zero, norm), value*(std::exp(1.)-std::exp(-1.))); + const auto gaussian = [&](double x) -> matrix { return value*std::exp(-x*x); }; + check(integrator.integrate(gaussian, -inf, inf, zero, norm), value*std::sqrt(pi)); + const auto right = [&](double x) -> matrix { return value*std::exp(-x); }; + const auto left = [&](double x) -> matrix { return value*std::exp(x); }; + check(integrator.integrate(right, 0., inf, zero, norm), value); + check(integrator.integrate(left, -inf, 0., zero, norm), value); + // Endpoint-distance functor: avoids cancellation in 1-x*x near +/-1. + const auto endpoint = [&](double, double complement) -> matrix { + const double c = std::abs(complement); + return value/std::sqrt(c*(2-c)); + }; + check(integrator.integrate(endpoint, zero, norm), value*pi); + check(integrator.integrate(endpoint, -1., 1., zero, norm), value*pi); + check(integrator.integrate(endpoint, 1., -1., zero, norm), -value*pi); + const auto shifted_endpoint = [&](double x, double complement) -> matrix { + BOOST_CHECK_GE(x, 2.); + BOOST_CHECK_LE(x, 6.); // x can round to the endpoint; complement retains distance. + const double c = std::abs(complement); + return value/std::sqrt(c*(4-c)); + }; + check(integrator.integrate(shifted_endpoint, 2., 6., zero, norm, 1e-11, &error, &l1), value*pi); + BOOST_CHECK_CLOSE_FRACTION(l1, norm(value)*pi, 1e-10); + // An endpoint singularity in the ordinary one-argument interface. + const auto logarithm = [&](double x) -> matrix { return value*std::log(x); }; + check(integrator.integrate(logarithm, 0., 1., zero, norm), -value); + calls = 0; + check(integrator.integrate(f, 2., 2., zero, norm, 1e-11, &error, &l1, &levels), zero); + BOOST_CHECK_EQUAL(calls, 0); + BOOST_CHECK_EQUAL(error, 0); + BOOST_CHECK_EQUAL(l1, 0); + BOOST_CHECK_EQUAL(levels, 0); + check(integrator.integrate(endpoint, 2., 2., zero, norm), zero); +} + +BOOST_AUTO_TEST_CASE(trapezoidal_matrices) +{ + unsigned calls = 0; + // Poisson kernel: its integral over one period is 2*pi. + const double r = 0.8; + const auto f = [&](double x) -> matrix { ++calls; return value*((1-r*r)/(1-2*r*std::cos(x)+r*r)); }; + q::trapezoidal(f, 0., 2*pi, zero, norm, 1e-3); + const unsigned coarse_calls = calls; + calls = 0; + double error = 0, l1 = 0; + matrix result = q::trapezoidal(f, 0., 2*pi, zero, norm, 1e-11, 12, &error, &l1); + check(result, value*(2*pi)); + BOOST_CHECK_GT(calls, coarse_calls); + BOOST_CHECK_LT(error, 1e-8); + BOOST_CHECK_CLOSE_FRACTION(l1, norm(value)*(2*pi), 1e-10); + const unsigned refined_calls = calls; + calls = 0; + double scaled_error = 0, scaled_l1 = 0; + check(q::trapezoidal(f, 0., 2*pi, zero, scaled_norm, 1e-11, 12, &scaled_error, &scaled_l1), result); + BOOST_CHECK_EQUAL(calls, refined_calls); + BOOST_CHECK_CLOSE_FRACTION(scaled_error, 8*error, 1e-10); + BOOST_CHECK_CLOSE_FRACTION(scaled_l1, 8*l1, 1e-10); + check(q::trapezoidal(f, 2*pi, 0., zero, norm, 1e-11), -result); + calls = 0; + check(q::trapezoidal(f, 1., 1., zero, norm, 1e-11, 12, &error, &l1), zero); + BOOST_CHECK_EQUAL(calls, 0); + BOOST_CHECK_EQUAL(error, 0); + BOOST_CHECK_EQUAL(l1, 0); +} + +BOOST_AUTO_TEST_CASE(nonthrowing_errors) +{ + using ignore = boost::math::policies::policy>; + const auto f = [&](double) -> matrix { return value; }; + const auto invalid = [](const matrix& m) { + BOOST_CHECK_EQUAL(m.rows(), 2); + BOOST_CHECK_EQUAL(m.cols(), 3); + for (Eigen::Index i = 0; i < m.size(); ++i) + BOOST_CHECK(std::isnan(m.data()[i].real())); + }; + const double nan = std::numeric_limits::quiet_NaN(); + invalid(q::exp_sinh().integrate(f, 0., 1., zero, norm)); + invalid(q::exp_sinh().integrate(f, zero, norm, 2.)); + invalid(q::sinh_sinh().integrate(f, zero, norm)); + invalid(q::tanh_sinh().integrate(f, nan, 1., zero, norm)); + invalid(q::trapezoidal(f, 0., inf, zero, norm, 1e-10, 12, static_cast(nullptr), static_cast(nullptr), ignore())); +} + +BOOST_AUTO_TEST_CASE(evaluation_errors_and_scalar_compatibility) +{ + const auto bad = [](double x) -> matrix { + if (std::abs(x) > 1e100) return zero; + matrix m = value; + m(0, 0) = std::numeric_limits::quiet_NaN(); + return m; + }; + BOOST_CHECK_THROW(q::exp_sinh().integrate(bad, zero, norm), boost::math::evaluation_error); + BOOST_CHECK_THROW(q::sinh_sinh().integrate(bad, zero, norm), boost::math::evaluation_error); + BOOST_CHECK_THROW(q::tanh_sinh().integrate(bad, zero, norm), boost::math::evaluation_error); + using ignore = boost::math::policies::policy>; + // Non-throwing evaluation policies retain the concrete estimate, as in the + // existing interface, rather than trying to construct a matrix from a scalar. + matrix result = q::exp_sinh().integrate(bad, zero, norm); + BOOST_CHECK_EQUAL(result.rows(), 2); + BOOST_CHECK_EQUAL(result.cols(), 3); + BOOST_CHECK(std::isnan(result(0, 0).real())); + result = q::sinh_sinh().integrate(bad, zero, norm); + BOOST_CHECK(std::isnan(result(0, 0).real())); + result = q::tanh_sinh().integrate(bad, zero, norm); + BOOST_CHECK(std::isnan(result(0, 0).real())); + + const auto f = [](double x) { return std::exp(-x*x); }; + const auto scalar_norm = [](double x) { return std::abs(x); }; + using function = decltype(f); + using exp_integrator = q::exp_sinh; + using sinh_integrator = q::sinh_sinh; + using tanh_integrator = q::tanh_sinh; + double (exp_integrator::*exp_original)(const function&, double, double*, double*, std::size_t*) const = &exp_integrator::integrate; + double (sinh_integrator::*sinh_original)(function, double, double*, double*, std::size_t*) const = &sinh_integrator::integrate; + double (tanh_integrator::*tanh_original)(function, double, double*, double*, std::size_t*) const = &tanh_integrator::integrate; + exp_integrator e; + sinh_integrator s; + tanh_integrator t; + BOOST_CHECK_EQUAL((e.*exp_original)(f, 1e-10, nullptr, nullptr, nullptr), e.integrate(f, 0., scalar_norm, 1e-10)); + BOOST_CHECK_EQUAL((s.*sinh_original)(f, 1e-10, nullptr, nullptr, nullptr), s.integrate(f, 0., scalar_norm, 1e-10)); + BOOST_CHECK_EQUAL((t.*tanh_original)(f, 1e-10, nullptr, nullptr, nullptr), t.integrate(f, 0., scalar_norm, 1e-10)); + double (*trap_original)(function, double, double, double, std::size_t, double*, double*) = &q::trapezoidal; + BOOST_CHECK_EQUAL(trap_original(f, -1., 1., 1e-10, 12, nullptr, nullptr), q::trapezoidal(f, -1., 1., 0., scalar_norm, 1e-10)); +} diff --git a/test/gauss_kronrod_eigen_test.cpp b/test/gauss_kronrod_eigen_test.cpp new file mode 100644 index 0000000000..f2fca33d36 --- /dev/null +++ b/test/gauss_kronrod_eigen_test.cpp @@ -0,0 +1,121 @@ +// Copyright Nick Thompson, 2026 +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt +#define BOOST_TEST_MODULE gauss_kronrod_eigen_test +#include +#include +#include +#include +#include +#include + +using boost::math::quadrature::gauss_kronrod; +using matrix = Eigen::MatrixXcd; + +// Both parities of the embedded Gauss rule, with rectangular dynamic matrices. +template +void test_adaptive() +{ + const matrix zero = matrix::Zero(2, 3); + matrix value(2, 3); + value << std::complex(1, 2), 2., -3., 4., 5., std::complex(-2, 1); + const auto norm = [](const matrix& m) { return m.norm(); }; + unsigned calls = 0; + const auto f = [&](double x) -> matrix { + ++calls; + return value / (1 + 10000*(x-0.37)*(x-0.37)); + }; + const double integral = (std::atan(100*(1-0.37)) + std::atan(100*0.37))/100; + const matrix expected = value * integral; + double error = 0, l1 = 0; + matrix coarse = gauss_kronrod::integrate(f, 0., 1., zero, norm, 0, 1e-10, &error, &l1); + BOOST_CHECK_EQUAL(calls, N); + const double coarse_error = (coarse - expected).norm(); + calls = 0; + matrix result = gauss_kronrod::integrate(f, 0., 1., zero, norm, 15, 1e-10, &error, &l1); + BOOST_CHECK_GT(calls, N); + const unsigned refined_calls = calls; + BOOST_CHECK_LT((result - expected).norm(), coarse_error); + BOOST_CHECK_SMALL((result - expected).norm(), 1e-11); + BOOST_CHECK_CLOSE_FRACTION(l1, value.norm()*integral, 1e-10); + BOOST_CHECK_GE(error, 0); + BOOST_CHECK_LT(error, 1e-10*expected.norm()); + + // Rescaling a norm rescales error/L1, but preserves refinement decisions. + const auto scaled_norm = [](const matrix& m) { return 8*m.norm(); }; + double scaled_error = 0, scaled_l1 = 0; + calls = 0; + matrix scaled = gauss_kronrod::integrate(f, 0., 1., zero, scaled_norm, 15, 1e-10, &scaled_error, &scaled_l1); + BOOST_CHECK_EQUAL(calls, refined_calls); + BOOST_CHECK_SMALL((scaled-result).norm(), 1e-14); + BOOST_CHECK_CLOSE_FRACTION(scaled_error, 8*error, 1e-12); + BOOST_CHECK_CLOSE_FRACTION(scaled_l1, 8*l1, 1e-12); + result = gauss_kronrod::integrate(f, 1., 0., zero, norm, 15, 1e-10, &error, &l1); + BOOST_CHECK_SMALL((result + expected).norm(), 1e-11); + BOOST_CHECK_CLOSE_FRACTION(l1, value.norm()*integral, 1e-10); +} + +BOOST_AUTO_TEST_CASE(adaptive_matrix_integral) +{ + test_adaptive<15>(); + test_adaptive<21>(); +} + +BOOST_AUTO_TEST_CASE(bounds_and_error_policies) +{ + using integrator = gauss_kronrod; + const matrix zero = matrix::Zero(2, 3); + const matrix value = matrix::Ones(2, 3); + const auto norm = [](const matrix& m) { return m.norm(); }; + const double inf = std::numeric_limits::infinity(); + const auto right = [&](double x) -> matrix { return value / ((1+x)*(1+x)); }; + const auto left = [&](double x) -> matrix { return value / ((1-x)*(1-x)); }; + const auto whole = [&](double x) -> matrix { return value / (boost::math::constants::pi()*(1+x*x)); }; + double error = 0, l1 = 0; + matrix result = integrator::integrate(right, 0., inf, zero, norm, 15, 1e-10, &error, &l1); + BOOST_CHECK_SMALL((result-value).norm(), 1e-12); + BOOST_CHECK_CLOSE_FRACTION(l1, value.norm(), 1e-12); + result = integrator::integrate(left, -inf, 0., zero, norm); + BOOST_CHECK_SMALL((result-value).norm(), 1e-12); + result = integrator::integrate(whole, -inf, inf, zero, norm, 15, 1e-10); + BOOST_CHECK_SMALL((result-value).norm(), 1e-10); + unsigned calls = 0; + const auto counted = [&](double) -> matrix { ++calls; return value; }; + result = integrator::integrate(counted, 1., 1., zero, norm, 15, 1e-10, &error, &l1); + BOOST_CHECK_EQUAL(result.rows(), 2); + BOOST_CHECK_EQUAL(result.cols(), 3); + BOOST_CHECK_EQUAL(result.norm(), 0); + BOOST_CHECK_EQUAL(error, 0); + BOOST_CHECK_EQUAL(l1, 0); + BOOST_CHECK_EQUAL(calls, 0); + const double nan = std::numeric_limits::quiet_NaN(); + BOOST_CHECK_THROW(integrator::integrate(counted, nan, 1., zero, norm), std::domain_error); + using ignore = boost::math::policies::policy>; + result = gauss_kronrod::integrate(counted, nan, 1., zero, norm, 15, 1e-10, &error, &l1); + BOOST_CHECK_EQUAL(result.rows(), 2); + BOOST_CHECK_EQUAL(result.cols(), 3); + BOOST_CHECK(std::isnan(error)); + BOOST_CHECK(std::isnan(l1)); + for (Eigen::Index i = 0; i < result.size(); ++i) + BOOST_CHECK(std::isnan(result.data()[i].real())); + BOOST_CHECK_EQUAL(calls, 0); +} + +BOOST_AUTO_TEST_CASE(scalar_compatibility_and_interval_scaling) +{ + using integrator = gauss_kronrod; + const auto f = [](double x) { return std::exp(x); }; + const auto norm = [](double x) { return std::abs(x); }; + using function = decltype(f); + double (*original)(function, double, double, unsigned, double, double*, double*) = &integrator::integrate; + BOOST_CHECK_EQUAL(original(f, 0., 1., 0, 1e-10, nullptr, nullptr), + integrator::integrate(f, 0., 1., 0., norm, 0, 1e-10)); + // An affine change of variables must scale the error estimate and L1. + const auto stretched = [&](double x) { return f(x/4); }; + double error1 = 0, error2 = 0, l1 = 0, l2 = 0; + double q1 = integrator::integrate(f, 0., 1., 0., norm, 0, 1e-10, &error1, &l1); + double q2 = integrator::integrate(stretched, 0., 4., 0., norm, 0, 1e-10, &error2, &l2); + BOOST_CHECK_EQUAL(q2, 4*q1); + BOOST_CHECK_EQUAL(error2, 4*error1); + BOOST_CHECK_EQUAL(l2, 4*l1); +} diff --git a/test/gauss_quadrature_eigen_test.cpp b/test/gauss_quadrature_eigen_test.cpp new file mode 100644 index 0000000000..d9d9544a84 --- /dev/null +++ b/test/gauss_quadrature_eigen_test.cpp @@ -0,0 +1,168 @@ +// Copyright Nick Thompson, 2026 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_TEST_MODULE gauss_quadrature_eigen_test + +#include +#include +#include +#include +#include +#include + +using boost::math::quadrature::gauss; +using matrix = Eigen::Matrix, 2, 2>; + +// Regression for https://github.com/boostorg/math/issues/337. +// Return concrete Eigen matrices so the test does not rely on the lifetime +// of expression-template temporaries. Supply a shaped zero and a scalar norm. +template +void test_matrix() +{ + const matrix zero = matrix::Zero(); + const auto norm = [](const matrix& m) { return m.norm(); }; + const auto f = [](double x) -> matrix { + matrix result; + result(0, 0) = {1, x * x}; + result(0, 1) = {x, 2 * x}; + result(1, 0) = {-x * x, x}; + result(1, 1) = {x * x * x, -1}; + return result; + }; + const auto primitive = [](double x) -> matrix { + matrix result; + result(0, 0) = {x, x * x * x / 3}; + result(0, 1) = {x * x / 2, x * x}; + result(1, 0) = {-x * x * x / 3, x * x / 2}; + result(1, 1) = {x * x * x * x / 4, -x}; + return result; + }; + const auto check = [](const matrix& actual, const matrix& expected) { + for (unsigned i = 0; i < 2; ++i) + for (unsigned j = 0; j < 2; ++j) + BOOST_CHECK_SMALL(std::abs(actual(i, j) - expected(i, j)), + 100 * std::numeric_limits::epsilon()); + }; + + check(gauss::integrate(f, zero, norm), primitive(1) + -primitive(-1)); + check(gauss::integrate(f, 0., 1., zero, norm), primitive(1)); + check(gauss::integrate(f, 1., 0., zero, norm), -primitive(1)); + check(gauss::integrate(f, 1., 1., zero, norm), matrix::Zero()); + + // A constant matrix has an exactly known integral and L1 norm. + const matrix value = f(1); + double l1 = 0; + check(gauss::integrate([&](double) { return value; }, 0., 3., zero, norm, &l1), + value * 3); + BOOST_CHECK_CLOSE_FRACTION(l1, 3 * value.norm(), + 100 * std::numeric_limits::epsilon()); +} + +BOOST_AUTO_TEST_CASE(matrix_valued_quadrature) +{ + test_matrix<7>(); + test_matrix<10>(); +} + +BOOST_AUTO_TEST_CASE(dynamic_matrices_and_infinite_bounds) +{ + using dynamic_matrix = Eigen::MatrixXcd; + const dynamic_matrix zero = dynamic_matrix::Zero(2, 3); + dynamic_matrix value(2, 3); + value << 1., 2., 3., 4., 5., 6.; + const auto norm = [](const dynamic_matrix& m) { return m.norm(); }; + const auto check = [&](const dynamic_matrix& actual) { + BOOST_CHECK_EQUAL(actual.rows(), 2); + BOOST_CHECK_EQUAL(actual.cols(), 3); + BOOST_CHECK_SMALL((actual - value).norm(), 1e-10); + }; + const double inf = std::numeric_limits::infinity(); + using integrator = gauss; + // These functions become constants under the half-line substitution. + const auto right = [&](double x) -> dynamic_matrix { return value / ((1+x)*(1+x)); }; + const auto left = [&](double x) -> dynamic_matrix { return value / ((1-x)*(1-x)); }; + double l1 = 0; + check(integrator::integrate(right, 0., inf, zero, norm, &l1)); + BOOST_CHECK_CLOSE_FRACTION(l1, value.norm(), 1e-12); + check(integrator::integrate(left, -inf, 0., zero, norm)); + const auto whole = [&](double x) -> dynamic_matrix { + return value / (boost::math::constants::pi() * (1+x*x)); + }; + check(integrator::integrate(whole, -inf, inf, zero, norm)); + + unsigned calls = 0; + const auto counted = [&](double) -> dynamic_matrix { ++calls; return value; }; + l1 = -1; + dynamic_matrix empty = integrator::integrate(counted, 2., 2., zero, norm, &l1); + BOOST_CHECK_EQUAL(empty.rows(), 2); + BOOST_CHECK_EQUAL(empty.cols(), 3); + BOOST_CHECK_EQUAL(empty.norm(), 0); + BOOST_CHECK_EQUAL(l1, 0); + BOOST_CHECK_EQUAL(calls, 0); + + const double nan = std::numeric_limits::quiet_NaN(); + BOOST_CHECK_THROW(integrator::integrate(counted, nan, 1., zero, norm), std::domain_error); + using ignore = boost::math::policies::policy< + boost::math::policies::domain_error>; + dynamic_matrix invalid = gauss::integrate(counted, nan, 1., zero, norm, &l1); + BOOST_CHECK_EQUAL(invalid.rows(), 2); + BOOST_CHECK_EQUAL(invalid.cols(), 3); + BOOST_CHECK(std::isnan(l1)); + for (Eigen::Index i = 0; i < invalid.size(); ++i) + BOOST_CHECK(std::isnan(invalid.data()[i].real())); + BOOST_CHECK_EQUAL(calls, 0); +} + +BOOST_AUTO_TEST_CASE(scalar_overload_compatibility) +{ + const auto f = [](double x) { return x*x; }; + const auto norm = [](double x) { return std::abs(x); }; + using integrator = gauss; + double old_l1 = 0, new_l1 = 0; + BOOST_CHECK_EQUAL(integrator::integrate(f, &old_l1), + integrator::integrate(f, 0., norm, &new_l1)); + BOOST_CHECK_EQUAL(old_l1, new_l1); + BOOST_CHECK_EQUAL(integrator::integrate(f, 0., 1., &old_l1), + integrator::integrate(f, 0., 1., 0., norm, &new_l1)); + BOOST_CHECK_EQUAL(old_l1, new_l1); + using function = decltype(f); + double (*original)(function, double*) = &integrator::integrate; + BOOST_CHECK_EQUAL(original(f, nullptr), integrator::integrate(f)); +} + +BOOST_AUTO_TEST_CASE(controllability_gramian) +{ + // The dimensions and matrix-exponential expression from issue #337. + // A*A = 0, so exp(A*t) = I + A*t and the integral is polynomial. + using gramian_matrix = Eigen::Matrix, 12, 12>; + Eigen::Matrix a = Eigen::Matrix::Zero(); + Eigen::Matrix b = Eigen::Matrix::Zero(); + a(0, 1) = 1; + b(1, 0) = 1; + const auto f = [&](double t) -> gramian_matrix { + return (a*t).exp() * b * b.transpose() * (a.transpose()*t).exp(); + }; + const gramian_matrix zero = gramian_matrix::Zero(); + const auto norm = [](const gramian_matrix& m) { return m.norm(); }; + gramian_matrix expected = zero; + expected(0, 0) = 1./3; + expected(0, 1) = expected(1, 0) = 0.5; + expected(1, 1) = 1; + gramian_matrix result = gauss::integrate(f, 0., 1., zero, norm); + BOOST_CHECK_SMALL((result - expected).norm(), 1e-14); + + // Verify that L1 follows the supplied norm, including on the canonical + // interval and when the integration bounds are reversed. + const auto scaled_norm = [](const gramian_matrix& m) { return 2*m.norm(); }; + const auto constant = [&](double) -> gramian_matrix { return expected; }; + double l1 = 0; + result = gauss::integrate(constant, zero, scaled_norm, &l1); + BOOST_CHECK_SMALL((result - 2*expected).norm(), 1e-14); + BOOST_CHECK_CLOSE_FRACTION(l1, 4*expected.norm(), 1e-14); + result = gauss::integrate(constant, 3., 0., zero, scaled_norm, &l1); + BOOST_CHECK_SMALL((result + 3*expected).norm(), 1e-14); + BOOST_CHECK_CLOSE_FRACTION(l1, 6*expected.norm(), 1e-14); +}