Multi dimensional integrator

include/boost/math/quadrature/multi_integrals.hpp

@@ -0,0 +1,84 @@
+// Copyright Jacob Hass, 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)
+
+#ifndef BOOST_MATH_QUADRATURE_MULTI_INTEGRALS_HPP
+#define BOOST_MATH_QUADRATURE_MULTI_INTEGRALS_HPP
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <functional>
+#include <cmath>
+#include <stdexcept>
+#include <boost/math/quadrature/gauss.hpp>
+#include <boost/math/quadrature/trapezoidal.hpp>
+#include <boost/math/policies/policy.hpp>
+#include <boost/math/policies/error_handling.hpp>
+
+namespace boost { namespace math { namespace quadrature { namespace detail {
+
+// Recursive engine tracking runtime depth and coordinate state
+template <typename F, typename Integrator, typename RealType>
+auto integrate_recursive(const F& f, 
+                         const std::vector<RealType>& a, 
+                         const std::vector<RealType>& b, 
+                         std::vector<RealType>& working_coords, 
+                         const size_t& depth,
+                         const Integrator& integrate) -> decltype(f(a))
+{
+    RealType low = a[depth];
+    RealType high = b[depth];
+
+    // Base Case: Innermost dimension reached
+    if (depth == a.size() - 1) {
+        auto innermost_f = [&](const RealType& x) {
+            working_coords[depth] = x; // Assign final coordinate
+            return f(working_coords);  // Evaluate the arbitrary-D function
+        };
+        return integrate(innermost_f, low, high);
+    }
+
+    // Recursive Step: Outer dimensions
+    auto next_dimension_f = [&](const RealType& x) {
+        working_coords[depth] = x; // Assign current dimension's coordinate
+        return integrate_recursive(f, a, b, working_coords, depth + 1, integrate);
+    };
+
+    return integrate(next_dimension_f, low, high);
+}
+
+} // namespace detail
+
+template <typename F, typename Integrator, typename RealType, class Policy>
+auto integrateND(const F& f, const std::vector<RealType>& a, const std::vector<RealType>& b, 
+                                   const Integrator& integrate, const Policy& pol) -> decltype(f(a))
+{
+    static const char* function = "boost::math::quadrature::integrateND(f, a, b, integrate, pol)";
+
+    // Error checking on integration bounds
+    if (a.empty() || b.empty()){
+        boost::math::policies::raise_evaluation_error(function, "Integration limits cannot be empty.\n", 0, pol);
+    }
+    if (a.size() != b.size()) {
+        std::string msg = "Integration limits must be the same size but got (" + std::to_string(a.size()) + "!=" + std::to_string(b.size()) + ")";
+        boost::math::policies::raise_evaluation_error(function, msg.c_str(), 0, pol);
+    }
+
+    // Allocate the coordinate state tracking vector dynamically based on runtime size
+    std::vector<RealType> working_coords(a.size(), 0.0);
+    return detail::integrate_recursive(f, a, b, working_coords, 0, integrate);
+}
+
+template <typename F, typename Integrator, typename RealType>
+auto integrateND(const F& f, const std::vector<RealType>& a, const std::vector<RealType>& b, 
+                 const Integrator& integrate) -> decltype(f(a))
+{
+    return integrateND(f, a, b, integrate, boost::math::policies::policy<>());
+}
+
+}}}
+
+#endif // BOOST_MATH_QUADRATURE_MULTI_INTEGRALS_HPP

test/Jamfile.v2

@@ -1345,6 +1345,7 @@ test-suite quadrature :
    [ run test_trapezoidal.cpp /boost/test//boost_unit_test_framework : : :
       release [ requires cxx11_lambdas cxx11_auto_declarations cxx11_decltype cxx11_unified_initialization_syntax cxx11_variadic_templates ]
       $(float128_type) ]
+   [ run test_multi_integrals.cpp ]
 ;
 
 test-suite autodiff :

test/test_multi_integrals.cpp

@@ -0,0 +1,146 @@
+/*
+ * Copyright Jacob Hass, 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 multi_integrals
+
+#include <algorithm>
+#include <type_traits>
+#include <vector>
+#include <boost/math/tools/config.hpp>
+#include <boost/test/included/unit_test.hpp>
+#include <boost/test/tools/floating_point_comparison.hpp>
+#include <boost/math/tools/test_value.hpp>
+#include <boost/math/concepts/real_concept.hpp>
+#include <boost/math/quadrature/multi_integrals.hpp>
+#include <boost/math/quadrature/tanh_sinh.hpp>
+#include <boost/math/quadrature/trapezoidal.hpp>
+#include <boost/math/quadrature/gauss_kronrod.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+#if __has_include(<stdfloat>)
+#  include <stdfloat>
+#endif
+
+using boost::math::quadrature::integrateND;
+using boost::math::quadrature::tanh_sinh;
+using boost::math::quadrature::gauss_kronrod;
+
+template <typename RealType>
+void test_catalan()
+{
+    auto f = [](const std::vector<RealType>& coords) {
+        return 1 / (1 + coords[0] * coords[0] * coords[1] * coords[1]);
+    };
+
+    std::vector<RealType> lower_bounds = {0.0, 0.0};
+    std::vector<RealType> upper_bounds = {1.0, 1.0}; 
+
+    typedef std::function<RealType(RealType)> func;
+
+    auto integrate = [](func f, RealType a, RealType b)
+    {
+        tanh_sinh<RealType> integrator;
+        return integrator.integrate(f, a, b);
+    }; 
+
+    RealType result = integrateND(f, lower_bounds, upper_bounds, integrate);
+    RealType tol = boost::math::tools::epsilon<RealType>();
+    BOOST_CHECK_CLOSE_FRACTION(result, boost::math::constants::catalan<RealType>(), 4 * tol);
+}
+
+template <typename RealType>
+void test_euler()
+{
+    BOOST_MATH_STD_USING
+
+    auto f = [](const std::vector<RealType>& coords) {
+        RealType x = coords[0];
+        RealType y = coords[1];
+        return exp( -x * y ) * (1 / (1 + x) - exp(-x));
+    };
+
+    std::vector<RealType> lower_bounds = {0.0, 0.0};
+    std::vector<RealType> upper_bounds = {std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()}; 
+
+    typedef std::function<RealType(RealType)> func;
+
+    auto integrate = [](func f, RealType a, RealType b)
+    {
+        return gauss_kronrod<RealType, 61>::integrate(f, a, b);
+    };
+
+    RealType result = integrateND(f, lower_bounds, upper_bounds, integrate);
+    RealType tol = boost::math::tools::epsilon<RealType>();
+    // Long double precision is not good compared to other precisions
+    RealType check_tol = std::is_same<long double, RealType>::value ? 3000 * tol : 4 * tol;
+    BOOST_CHECK_CLOSE_FRACTION(result, boost::math::constants::euler<RealType>(), check_tol);
+}
+
+template <typename RealType>
+void test_square()
+{
+    BOOST_MATH_STD_USING 
+
+    auto f = [](const std::vector<RealType>& coords)
+    {
+        return abs(coords[0] + coords[1] + coords[2]);
+    };
+
+    typedef std::function<RealType(RealType)> func;
+    auto integrate = [](func f, RealType a, RealType b)
+    {
+        return gauss_kronrod<RealType, 61>::integrate(f, a, b);
+    }; 
+
+    std::vector<RealType> lower_bounds = {-1.0, -1.0, -1.0};
+    std::vector<RealType> upper_bounds = {1.0, 1.0, 1.0};
+
+    RealType result = integrateND(f, lower_bounds, upper_bounds, integrate);
+    RealType tol = boost::math::tools::epsilon<RealType>();
+    RealType multFactor;
+
+    // Does not seem to perform well on this integral
+    if(std::is_same<RealType, float>::value){
+        multFactor = 200;
+    }
+    else if (std::is_same<RealType, double>::value){
+        multFactor = 4 * pow(10, 6);
+    }
+    else if (std::is_same<RealType, long double>::value){
+        multFactor = 2 * pow(10, 8);
+    }
+
+    BOOST_CHECK_CLOSE_FRACTION(result, RealType(13) / RealType(2), multFactor * tol);
+}
+
+BOOST_AUTO_TEST_CASE(multi_integrals)
+{
+    // Catalan constant tests
+    test_catalan<float>();
+    test_catalan<double>();
+    #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+        test_catalan<long double>();
+    #endif
+    #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
+        test_catalan<boost::math::concepts::real_concept>();
+    #endif
+
+    // Euler constant tests
+    test_euler<float>();
+    test_euler<double>();
+    #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+        test_euler<long double>();
+    #endif
+
+    // Square integral tests
+    test_square<float>();
+    test_square<double>();
+    #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+        test_square<long double>();
+    #endif
+}