Support explicit zero and norm in Gauss and Gauss-Kronrod quadrature

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<std::complex<double>, 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<double> 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

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<Real>());
 
 
+[heading Vector- and matrix-valued integrands]
+
+Additional overloads accept an explicit additive identity and a scalar-valued norm:
+
+    template <class F, class Norm>
+    static auto integrate(F f, const decltype(f(Real(0)))& zero, Norm norm,
+                          Real* pL1 = nullptr) -> decltype(f(Real(0)));
+
+    template <class F, class Norm>
+    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<std::complex<double>, 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<double, 10>::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

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<F>()(std::declval<Real>()));
    };
 
+[heading Vector- and matrix-valued integrands]
+
+An additional overload accepts an explicit additive identity and a scalar norm:
+
+    template <class F, class Norm>
+    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>(),
+                          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<double, 15>::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();

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<class F, class Real, class Norm>
+    auto trapezoidal(F f, Real a, Real b, const decltype(f(a))& zero, Norm norm,
+        Real tol = boost::math::tools::root_epsilon<Real>(),
+        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.

include/boost/math/quadrature/detail/exp_sinh_detail.hpp

@@ -10,6 +10,7 @@
 #include <boost/math/tools/config.hpp>
 
 #ifndef BOOST_MATH_HAS_NVRTC
+#include <boost/math/quadrature/detail/norm_quadrature_error.hpp>
 
 #ifndef BOOST_MATH_BUILD_MODULE
 #include <cmath>
@@ -53,6 +54,8 @@ class exp_sinh_detail
 
     template<class F>
     auto integrate(const F& f, Real* error, Real* L1, const char* function, Real tolerance, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()));
+    template<class F, class Norm>
+    auto integrate(const F& f, const decltype(std::declval<F>()(std::declval<Real>()))& zero, Norm norm, Real* error, Real* L1, const char* function, Real tolerance, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()));
 
 private:
    const std::vector<Real>& get_abscissa_row(std::size_t n)const
@@ -295,6 +298,144 @@ auto exp_sinh_detail<Real, Policy>::integrate(const F& f, Real* error, Real* L1,
     return I1;
 }
 
+template<class Real, class Policy>
+template<class F, class Norm>
+auto exp_sinh_detail<Real, Policy>::integrate(const F& f, const decltype(std::declval<F>()(std::declval<Real>()))& zero, Norm norm, Real* error, Real* L1, const char* function, Real tolerance, std::size_t* levels) const ->decltype(std::declval<F>()(std::declval<Real>()))
+{
+    const auto magnitude = [&](const decltype(std::declval<F>()(std::declval<Real>()))& value) -> Real { return static_cast<Real>(norm(value)); };
+    typedef decltype(f(static_cast<Real>(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<Real>::digits10);
+
+    // Get the party started with two estimates of the integral:
+    Real min_abscissa{ 0 }, max_abscissa{ boost::math::tools::max_value<Real>() };
+    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<Real>();
+    }
+
+    I1 *= half<Real>();
+    L1_I1 *= half<Real>();
+    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<Real>()*I0;
+        L1_I1 = half<Real>()*L1_I0;
+        Real h = static_cast<Real>(1)/static_cast<Real>(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<class Real, class Policy>
 void exp_sinh_detail<Real, Policy>::init(const std::integral_constant<int, 0>&)

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 <cstddef>
+#endif
+namespace boost { namespace math { namespace quadrature { namespace detail {
+template<class K, class Real>
+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