From 7f20e15f23c88f7da3d9bb7d4e0ab4a79725dd41 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Fri, 28 Aug 2026 10:11:53 -0700 Subject: [PATCH] Van den Bos quadrature --- include/boost/math/quadrature/van_den_bos.hpp | 513 ++++++++++++++++++ test/test_van_den_bos.cpp | 455 ++++++++++++++++ 2 files changed, 968 insertions(+) create mode 100644 include/boost/math/quadrature/van_den_bos.hpp create mode 100644 test/test_van_den_bos.cpp diff --git a/include/boost/math/quadrature/van_den_bos.hpp b/include/boost/math/quadrature/van_den_bos.hpp new file mode 100644 index 0000000000..43baaf6cf0 --- /dev/null +++ b/include/boost/math/quadrature/van_den_bos.hpp @@ -0,0 +1,513 @@ +/* + * Copyright 2026 Nick Thompson + * + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. + * + * Positive cubature reduction in the style of: + * + * L. M. M. van den Bos, B. Sanderse, W. A. A. M. Bierbooms, + * G. J. W. van Bussel, + * "Generating nested quadrature rules with positive weights based on + * arbitrary sample sets", SIAM/ASA Journal on Uncertainty Quantification. + * + * This header implements the invariant-preserving reduction step: + * + * A w = m, w >= 0 + * + * and repeatedly moves along null vectors c of A, + * + * w <- w - alpha c, + * + * choosing alpha at a positivity boundary so that a node weight becomes zero. + * Exactness is preserved because A c = 0. + * + * The implementation is deliberately header-only and has no dependency on an + * external linear-algebra package. The null vector is obtained with + * pivoted Gaussian elimination. + */ + +#ifndef BOOST_MATH_QUADRATURE_VAN_DEN_BOS_HPP +#define BOOST_MATH_QUADRATURE_VAN_DEN_BOS_HPP + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +#include + +namespace boost { namespace math { namespace quadrature { + +template +class positive_cubature_rule +{ +public: + typedef Point point_type; + typedef Real real_type; + + positive_cubature_rule() = default; + + positive_cubature_rule(std::vector points, std::vector weights) + : points_(std::move(points)), weights_(std::move(weights)) + { + if (points_.size() != weights_.size()) + throw std::invalid_argument( + "positive_cubature_rule: points and weights must have the same size."); + } + + std::size_t size() const noexcept { return points_.size(); } + + std::vector const& points() const noexcept { return points_; } + std::vector const& weights() const noexcept { return weights_; } + + Point const& point(std::size_t i) const noexcept { return points_[i]; } + Real weight(std::size_t i) const noexcept { return weights_[i]; } + + template + auto integrate(F const& f, Real* L1 = nullptr) const + -> decltype(std::declval()(std::declval())) + { + using std::abs; + typedef decltype(f(std::declval())) K; + + static_assert(!std::is_integral::value, + "The integrand return type must be a real or complex floating-point type."); + + K result = K(0); + Real l1 = Real(0); + + for (std::size_t i = 0; i < points_.size(); ++i) + { + K y = f(points_[i]); + result += weights_[i] * y; + l1 += weights_[i] * abs(y); + } + + if (L1) + *L1 = l1; + return result; + } + +private: + std::vector points_; + std::vector weights_; +}; + +namespace detail { + +template +Real max_abs(std::mdspan A) +{ + using std::abs; + Real m = 0; + for (std::size_t i = 0; i < A.extent(0); ++i) + { + for (std::size_t j = 0; j < A.extent(1); ++j) + { + m = (std::max)(m, abs(A[i, j])); + } + } + return m; +} + +/* + * Compute a null vector of A by row reduction. + * + * free_choice selects which free variable is set to one. Returning different + * free choices gives different null directions and is useful when some nodes + * are protected from removal. + */ +template +bool null_vector( + std::mdspan input, + std::size_t free_choice, + std::vector& c, + Real rank_tolerance) +{ + std::vector storage(input.extent(0) * input.extent(1)); + std::mdspan A(storage.data(), input.extent(0), input.extent(1)); + + for (std::size_t i = 0; i < input.extent(0); ++i) + { + for (std::size_t j = 0; j < input.extent(1); ++j) + { + A[i, j] = input[i, j]; + } + } + using std::abs; + + const std::size_t m = A.extent(0); + const std::size_t n = A.extent(1); + + if (n == 0) + return false; + + Real scale = max_abs(A); + if (scale == Real(0)) + { + c.assign(n, Real(0)); + c[free_choice % n] = Real(1); + return true; + } + + Real tol = rank_tolerance * scale; + + std::vector pivot_col; + std::size_t row = 0; + + for (std::size_t col = 0; col < n && row < m; ++col) + { + std::size_t p = row; + Real best = abs(A[row, col]); + for (std::size_t i = row + 1; i < m; ++i) + { + Real a = abs(A[i, col]); + if (a > best) + { + best = a; + p = i; + } + } + + if (best <= tol) + continue; + + if (p != row) + { + for (std::size_t j = 0; j < n; ++j) + { + std::swap(A[row, j], A[p, j]); + } + } + + Real pivot = A[row, col]; + for (std::size_t j = col; j < n; ++j) + A[row, j] /= pivot; + + for (std::size_t i = 0; i < m; ++i) + { + if (i == row) + continue; + Real q = A[i, col]; + if (abs(q) <= tol) + continue; + for (std::size_t j = col; j < n; ++j) + A[i, j] -= q * A[row, j]; + } + + pivot_col.push_back(col); + ++row; + } + + std::vector is_pivot(n, 0); + for (std::size_t j : pivot_col) + is_pivot[j] = 1; + + std::vector free_cols; + for (std::size_t j = 0; j < n; ++j) + if (!is_pivot[j]) + free_cols.push_back(j); + + if (free_cols.empty()) + return false; + + std::size_t f = free_cols[free_choice % free_cols.size()]; + c.assign(n, Real(0)); + c[f] = Real(1); + + // Since A is in reduced row-echelon form: + // x_pivot + sum_free A[row, free] x_free = 0. + for (std::size_t r = 0; r < pivot_col.size(); ++r) + { + std::size_t p = pivot_col[r]; + c[p] = -A[r, f]; + } + + return true; +} + +template +Real null_residual( + std::mdspan A, + std::vector const& c) +{ + using std::abs; + Real r = 0; + for (std::size_t i = 0; i < A.extent(0); ++i) + { + Real s = 0; + for (std::size_t j = 0; j < A.extent(1); ++j) + s += A[i, j] * c[j]; + r = (std::max)(r, abs(s)); + } + return r; +} + +template +struct reduction_candidate +{ + bool valid = false; + bool use_upper = true; + std::size_t removed_local = 0; + Real alpha = 0; + Real protected_margin = 0; +}; + +/* + * For w(alpha)=w-alpha*c, determine the two positivity-boundary values of + * alpha. Prefer a boundary that removes an unprotected node. Among valid + * choices, maximize the minimum protected weight. + */ +template +reduction_candidate choose_boundary( + std::vector const& w, + std::vector const& c, + std::vector const& protected_local, + Real positivity_tolerance) +{ + using std::abs; + + const Real inf = (std::numeric_limits::infinity)(); + + Real upper = inf; + Real lower = -inf; + std::size_t upper_idx = 0; + std::size_t lower_idx = 0; + bool have_upper = false; + bool have_lower = false; + + for (std::size_t i = 0; i < w.size(); ++i) + { + if (c[i] > positivity_tolerance) + { + Real a = w[i] / c[i]; + if (!have_upper || a < upper) + { + upper = a; + upper_idx = i; + have_upper = true; + } + } + else if (c[i] < -positivity_tolerance) + { + Real a = w[i] / c[i]; + if (!have_lower || a > lower) + { + lower = a; + lower_idx = i; + have_lower = true; + } + } + } + + reduction_candidate best; + + auto inspect = [&](Real alpha, std::size_t boundary_idx, bool upper_side) + { + if (protected_local[boundary_idx]) + return; + + Real margin = inf; + for (std::size_t i = 0; i < w.size(); ++i) + { + Real wi = w[i] - alpha * c[i]; + if (wi < -positivity_tolerance) + return; + if (protected_local[i]) + margin = (std::min)(margin, wi); + } + + if (!best.valid || margin > best.protected_margin) + { + best.valid = true; + best.use_upper = upper_side; + best.removed_local = boundary_idx; + best.alpha = alpha; + best.protected_margin = margin; + } + }; + + if (have_upper) + inspect(upper, upper_idx, true); + if (have_lower) + inspect(lower, lower_idx, false); + + return best; +} + +} // namespace detail + +/* + * Reduce a positive exact cubature rule while preserving the moments of a + * user-supplied basis. + * + * basis(point, k) evaluates basis function k at point. + * moments[k] is the exact desired integral of basis function k. + * + * The input weights are assumed to already satisfy the moment equations. + * + * protected_count: + * the first protected_count nodes may not be removed. This is convenient + * for constructing nested rules: place old nodes first, append new + * candidates, and protect the prefix. + * + * target_size: + * reduction stops once this many nodes remain, or when no admissible + * positivity-preserving null direction is found. + * + * The routine preserves nonnegative weights up to positivity_tolerance. + */ +template +positive_cubature_rule van_den_bos_reduce( + std::vector points, + std::vector weights, + Basis const& basis, + std::vector const& moments, + std::size_t target_size = 0, + std::size_t protected_count = 0, + Real rank_tolerance = Real(100) * boost::math::tools::epsilon(), + Real positivity_tolerance = Real(1000) * boost::math::tools::epsilon()) +{ + using std::abs; + + if (points.size() != weights.size()) + throw std::invalid_argument( + "van_den_bos_reduce: points and weights must have the same size."); + if (protected_count > points.size()) + throw std::invalid_argument( + "van_den_bos_reduce: protected_count exceeds number of points."); + if (moments.empty()) + throw std::invalid_argument( + "van_den_bos_reduce: moments must not be empty."); + + for (Real const& w : weights) + if (w < -positivity_tolerance) + throw std::domain_error( + "van_den_bos_reduce: input rule must have nonnegative weights."); + + if (target_size == 0) + target_size = moments.size(); + + target_size = (std::max)(target_size, protected_count); + + std::vector protected_flag(points.size(), 0); + for (std::size_t i = 0; i < protected_count; ++i) + protected_flag[i] = 1; + + // Remove numerically-zero unprotected input weights immediately. + for (std::size_t i = points.size(); i-- > protected_count;) + { + if (abs(weights[i]) <= positivity_tolerance) + { + points.erase(points.begin() + static_cast(i)); + weights.erase(weights.begin() + static_cast(i)); + protected_flag.erase(protected_flag.begin() + static_cast(i)); + } + } + + while (points.size() > target_size) + { + const std::size_t m = moments.size(); + const std::size_t n = points.size(); + + std::vector matrix_storage(m * n); + std::mdspan A(matrix_storage.data(), m, n); + for (std::size_t k = 0; k < m; ++k) + { + for (std::size_t j = 0; j < n; ++j) + { + A[k, j] = basis(points[j], k); + } + } + + std::vector local_protected = protected_flag; + + bool reduced = false; + + // Different free variables generate different null directions. + // Try all n choices; modulo the nullity this cycles through the + // available elementary free-variable directions. + for (std::size_t choice = 0; choice < n && !reduced; ++choice) + { + std::vector c; + if (!detail::null_vector(A, choice, c, rank_tolerance)) + break; + + Real cscale = 0; + for (Real const& x : c) + cscale = (std::max)(cscale, abs(x)); + if (cscale == Real(0)) + continue; + for (Real& x : c) + x /= cscale; + + auto candidate = detail::choose_boundary( + weights, c, local_protected, positivity_tolerance); + + if (!candidate.valid) + continue; + + for (std::size_t i = 0; i < n; ++i) + { + weights[i] -= candidate.alpha * c[i]; + if (abs(weights[i]) <= positivity_tolerance) + weights[i] = Real(0); + } + + // Remove all zero-weight unprotected nodes generated by this step. + for (std::size_t i = points.size(); i-- > 0;) + { + if (!protected_flag[i] && + abs(weights[i]) <= positivity_tolerance) + { + points.erase(points.begin() + static_cast(i)); + weights.erase(weights.begin() + static_cast(i)); + protected_flag.erase( + protected_flag.begin() + static_cast(i)); + } + } + + reduced = true; + } + + if (!reduced) + break; + } + + return positive_cubature_rule( + std::move(points), std::move(weights)); +} + +/* + * Return the maximum moment residual of a rule. + */ +template +Real van_den_bos_moment_residual( + positive_cubature_rule const& rule, + Basis const& basis, + std::vector const& moments) +{ + using std::abs; + + Real result = 0; + for (std::size_t k = 0; k < moments.size(); ++k) + { + Real q = 0; + for (std::size_t j = 0; j < rule.size(); ++j) + q += rule.weight(j) * basis(rule.point(j), k); + result = (std::max)(result, abs(q - moments[k])); + } + return result; +} + +}}} // namespaces + +#endif diff --git a/test/test_van_den_bos.cpp b/test/test_van_den_bos.cpp new file mode 100644 index 0000000000..20d72478a0 --- /dev/null +++ b/test/test_van_den_bos.cpp @@ -0,0 +1,455 @@ +/* + * Copyright 2026 Nick Thompson + * + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. + */ + +#include +#include "math_unit_test.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +template +std::pair, std::vector> gauss_legendre(std::size_t n) +{ + using std::abs; + using std::acos; + using std::cos; + + const Real pi = acos(Real(-1)); + std::vector x(n); + std::vector w(n); + const std::size_t m = (n + 1) / 2; + + for (std::size_t i = 0; i < m; ++i) + { + Real z = cos(pi * (Real(i) + Real(0.75)) / (Real(n) + Real(0.5))); + Real z1 = 0; + Real pp = 0; + do + { + Real p1 = 1; + Real p2 = 0; + for (std::size_t j = 1; j <= n; ++j) + { + Real p3 = p2; + p2 = p1; + p1 = ((Real(2 * j - 1) * z * p2) - Real(j - 1) * p3) / Real(j); + } + pp = Real(n) * (z * p1 - p2) / (z * z - Real(1)); + z1 = z; + z = z1 - p1 / pp; + } + while (abs(z - z1) > Real(8) * std::numeric_limits::epsilon()); + + x[i] = (Real(1) - z) / Real(2); + x[n - 1 - i] = (Real(1) + z) / Real(2); + const Real wi = Real(1) / ((Real(1) - z * z) * pp * pp); + w[i] = wi; + w[n - 1 - i] = wi; + } + return {x, w}; +} + +template +void append_exponents( + unsigned remaining, + std::size_t coordinate, + std::array& current, + std::vector>& result) +{ + if (coordinate + 1 == Dim) + { + current[coordinate] = remaining; + result.push_back(current); + return; + } + + for (unsigned k = 0; k <= remaining; ++k) + { + current[coordinate] = k; + append_exponents(remaining - k, coordinate + 1, current, result); + } +} + +template +std::vector> total_degree_exponents(unsigned degree) +{ + std::vector> result; + std::array current{}; + for (unsigned total = 0; total <= degree; ++total) + { + append_exponents(total, 0, current, result); + } + return result; +} + +template +Real shifted_legendre(unsigned n, Real x) +{ + const Real t = Real(2) * x - Real(1); + if (n == 0) + { + return Real(1); + } + if (n == 1) + { + return t; + } + + Real p0 = 1; + Real p1 = t; + for (unsigned k = 2; k <= n; ++k) + { + const Real p = (Real(2 * k - 1) * t * p1 - Real(k - 1) * p0) / Real(k); + p0 = p1; + p1 = p; + } + return p1; +} + +template +auto make_unit_cube_rule(unsigned degree, std::size_t tensor_order) +{ + using point = std::array; + using boost::math::quadrature::van_den_bos_moment_residual; + using boost::math::quadrature::van_den_bos_reduce; + + const auto exponents = total_degree_exponents(degree); + const auto gl = gauss_legendre(tensor_order); + + std::vector points; + std::vector weights; + std::array index{}; + bool done = false; + while (!done) + { + point p{}; + Real w = 1; + for (std::size_t d = 0; d < Dim; ++d) + { + p[d] = gl.first[index[d]]; + w *= gl.second[index[d]]; + } + points.push_back(p); + weights.push_back(w); + + for (std::size_t d = Dim; d-- > 0;) + { + ++index[d]; + if (index[d] < tensor_order) + { + break; + } + index[d] = 0; + if (d == 0) + { + done = true; + } + } + } + + auto basis = [exponents](point const& p, std::size_t k) + { + Real value = 1; + for (std::size_t d = 0; d < Dim; ++d) + { + value *= shifted_legendre(exponents[k][d], p[d]); + } + return value; + }; + + std::vector moments(exponents.size(), Real(0)); + moments[0] = Real(1); + + auto rule = van_den_bos_reduce( + std::move(points), std::move(weights), basis, moments, + exponents.size(), 0, Real(1e-12), Real(1e-13)); + + CHECK_EQUAL(rule.size(), exponents.size()); + CHECK_LE(van_den_bos_moment_residual(rule, basis, moments), Real(5e-10)); + for (Real w : rule.weights()) + { + CHECK_GE(w, Real(-1e-13)); + } + return rule; +} + +template +double genz_oscillatory_exact( + std::array const& a, + std::array const& u) +{ + using complex = std::complex; + const double pi = std::acos(-1.0); + complex value = std::exp(complex(0, 2 * pi * u[0])); + for (double ai : a) + { + if (ai == 0) + { + continue; + } + value *= (std::exp(complex(0, ai)) - complex(1, 0)) / complex(0, ai); + } + return value.real(); +} + +template +double genz_product_peak_exact( + std::array const& a, + std::array const& u) +{ + double value = 1; + for (std::size_t i = 0; i < Dim; ++i) + { + value *= a[i] * (std::atan(a[i] * (1 - u[i])) + std::atan(a[i] * u[i])); + } + return value; +} + +template +double genz_corner_peak_exact(std::array const& a) +{ + double alternating_sum = 0; + const std::size_t subsets = std::size_t(1) << Dim; + for (std::size_t mask = 0; mask < subsets; ++mask) + { + double denominator = 1; + unsigned bits = 0; + for (std::size_t i = 0; i < Dim; ++i) + { + if (mask & (std::size_t(1) << i)) + { + denominator += a[i]; + ++bits; + } + } + alternating_sum += (bits & 1 ? -1.0 : 1.0) / denominator; + } + + double scale = 1; + for (std::size_t i = 2; i <= Dim; ++i) + { + scale *= double(i); + } + for (double ai : a) + { + scale *= ai; + } + return alternating_sum / scale; +} + +template +double genz_gaussian_exact( + std::array const& a, + std::array const& u) +{ + const double sqrt_pi = std::sqrt(std::acos(-1.0)); + double value = 1; + for (std::size_t i = 0; i < Dim; ++i) + { + value *= sqrt_pi / (2 * a[i]) + * (std::erf(a[i] * u[i]) + std::erf(a[i] * (1 - u[i]))); + } + return value; +} + +template +double genz_continuous_exact( + std::array const& a, + std::array const& u) +{ + double value = 1; + for (std::size_t i = 0; i < Dim; ++i) + { + value *= (2 - std::exp(-a[i] * u[i]) - std::exp(-a[i] * (1 - u[i]))) / a[i]; + } + return value; +} + +template +double genz_discontinuous_exact( + std::array const& a, + std::array const& u) +{ + double value = 1; + for (std::size_t i = 0; i < Dim; ++i) + { + const double upper = i < 2 ? u[i] : 1.0; + value *= std::expm1(a[i] * upper) / a[i]; + } + return value; +} + +void test_basic_reduction() +{ + using point = std::array; + using boost::math::quadrature::van_den_bos_moment_residual; + using boost::math::quadrature::van_den_bos_reduce; + + const std::vector x = {-1.0, 0.0, 1.0}; + const std::vector w1 = {1.0 / 3.0, 4.0 / 3.0, 1.0 / 3.0}; + + std::vector points; + std::vector weights; + for (std::size_t i = 0; i < 3; ++i) + { + for (std::size_t j = 0; j < 3; ++j) + { + points.push_back({x[i], x[j]}); + weights.push_back(w1[i] * w1[j]); + } + } + + auto basis = [](point const& p, std::size_t k) + { + switch (k) + { + case 0: return 1.0; + case 1: return p[0]; + case 2: return p[1]; + case 3: return p[0] * p[0]; + case 4: return p[0] * p[1]; + case 5: return p[1] * p[1]; + default: return 0.0; + } + }; + + const std::vector moments = { + 4.0, 0.0, 0.0, 4.0 / 3.0, 0.0, 4.0 / 3.0 + }; + + auto rule = van_den_bos_reduce(points, weights, basis, moments, 6, 0); + CHECK_EQUAL(rule.size(), std::size_t(6)); + CHECK_LE(van_den_bos_moment_residual(rule, basis, moments), 1e-12); + + const double q = rule.integrate([](point const& p) + { + return 1.0 + p[0] * p[0] + p[1] * p[1] + p[0] * p[1]; + }); + CHECK_ABSOLUTE_ERROR(20.0 / 3.0, q, 1e-12); +} + +void test_complex_integrand() +{ + using point = std::array; + const auto rule = make_unit_cube_rule<2>(12, 10); + + std::complex exact = 0; + std::complex power = 1; + double factorial = 1; + for (unsigned k = 0; k < 100; ++k) + { + if (k > 0) + { + power *= std::complex(0, 5); + factorial *= double(k); + } + const double kp1 = double(k + 1); + const std::complex term = power / (factorial * kp1 * kp1); + exact += term; + if (std::abs(term) < 1e-18) + { + break; + } + } + + const auto computed = rule.integrate([](point const& p) + { + return std::exp(std::complex(0, 5 * p[0] * p[1])); + }); + + CHECK_ABSOLUTE_ERROR(exact.real(), computed.real(), 2e-8); + CHECK_ABSOLUTE_ERROR(exact.imag(), computed.imag(), 5e-8); +} + +void test_genz_3d() +{ + using point = std::array; + + // Mild but nontrivial parameters. Genz's a-vector controls difficulty; + // keeping it moderate makes these deterministic regression tests rather + // than convergence benchmarks. + const std::array a = {0.35, 0.45, 0.55}; + const std::array u = {0.30, 0.60, 0.75}; + const double pi = std::acos(-1.0); + + const auto rule = make_unit_cube_rule<3>(8, 6); + + const double oscillatory = rule.integrate([&](point const& x) + { + return std::cos(2 * pi * u[0] + a[0] * x[0] + a[1] * x[1] + a[2] * x[2]); + }); + CHECK_ABSOLUTE_ERROR(genz_oscillatory_exact(a, u), oscillatory, 2e-8); + + const double product_peak = rule.integrate([&](point const& x) + { + double value = 1; + for (std::size_t i = 0; i < 3; ++i) + { + const double d = x[i] - u[i]; + value /= 1 / (a[i] * a[i]) + d * d; + } + return value; + }); + CHECK_ABSOLUTE_ERROR(genz_product_peak_exact(a, u), product_peak, 2e-8); + + const double corner_peak = rule.integrate([&](point const& x) + { + const double s = 1 + a[0] * x[0] + a[1] * x[1] + a[2] * x[2]; + return 1 / (s * s * s * s); + }); + CHECK_ABSOLUTE_ERROR(genz_corner_peak_exact(a), corner_peak, 5e-7); + + const double gaussian = rule.integrate([&](point const& x) + { + double exponent = 0; + for (std::size_t i = 0; i < 3; ++i) + { + const double d = x[i] - u[i]; + exponent -= a[i] * a[i] * d * d; + } + return std::exp(exponent); + }); + CHECK_ABSOLUTE_ERROR(genz_gaussian_exact(a, u), gaussian, 2e-8); + + const double continuous = rule.integrate([&](point const& x) + { + double exponent = 0; + for (std::size_t i = 0; i < 3; ++i) + { + exponent -= a[i] * std::abs(x[i] - u[i]); + } + return std::exp(exponent); + }); + CHECK_ABSOLUTE_ERROR(genz_continuous_exact(a, u), continuous, 2e-3); + + const double discontinuous = rule.integrate([&](point const& x) + { + if (x[0] > u[0] || x[1] > u[1]) + { + return 0.0; + } + return std::exp(a[0] * x[0] + a[1] * x[1] + a[2] * x[2]); + }); + CHECK_ABSOLUTE_ERROR(genz_discontinuous_exact(a, u), discontinuous, 9e-2); +} + +} // namespace + +int main() +{ + test_basic_reduction(); + test_complex_integrand(); + test_genz_3d(); + return boost::math::test::report_errors(); +}