Van den Bos quadrature

doc/math.qbk

@@ -767,6 +767,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [endmathpart]
 
 [mathpart quadrature Quadrature and Differentiation]
+[include quadrature/van_den_bos_unit_square.qbk]
 [include quadrature/trapezoidal.qbk]
 [include quadrature/gauss.qbk]
 [include quadrature/gauss_kronrod.qbk]

doc/quadrature/van_den_bos_unit_square.qbk

@@ -0,0 +1,105 @@
+[/
+Copyright (c) 2026 Nicholas Thompson and Matt Borland
+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)
+]
+
+[section:van_den_bos_unit_square Van den Bos Unit-Square Quadrature]
+
+[heading Synopsis]
+
+    #include <boost/math/quadrature/van_den_bos_unit_square.hpp>
+    namespace boost { namespace math { namespace quadrature {
+
+    template<class F>
+    auto van_den_bos_unit_square(
+        F const& f,
+        double tolerance,
+        double* error_estimate = nullptr,
+        double* L1 = nullptr,
+        std::size_t* evaluations = nullptr);
+
+    }}} // namespaces
+
+[heading Description]
+
+`van_den_bos_unit_square` computes an adaptive positive cubature of a
+real- or complex-valued function over the reference square [0,1]^2. The
+integrand receives the two coordinates as separate `double` arguments. For a
+physical domain, map the coordinates and include the Jacobian in the
+integrand.
+
+The precomputed rules form a nested hierarchy. After the first level,
+only points newly introduced by the next level are evaluated. The routine
+requires two consecutive adjacent-level differences to satisfy `tolerance`
+before stopping, which avoids accepting an isolated cancellation. If the
+maximum level is reached first, the last rule is returned and its difference
+is used as the final `error_estimate` (subject to the roundoff floor).
+
+The optional `L1` output is the weighted sum of the absolute integrand values.
+`evaluations` reports the number of function calls, including calls reused by
+all subsequent levels. The supported result types are `float` and `double`,
+including complex values with `float` or `double` components. Higher-precision
+return types do not increase the precision of the nodes or weights and are not
+supported.
+
+The header and generated tables require C++17. Table construction is
+performed offline in `cpp_bin_float_quad`; emitted nodes and weights are
+rounded once to binary64 hexadecimal literals. The reproducible generator is
+`tools/dump_van_den_bos_unit_square.cpp`. No rules are constructed at runtime.
+If the requested tolerance is not reached by the final precomputed level, the
+routine returns that level's result and error estimate; it does not compute
+additional levels.
+
+Positive cubature rules satisfying the same moment conditions are not unique.
+The stored nodes and weights are one reproducible construction; different
+candidate sets or valid reduction paths can produce different positive nested
+rules with the same polynomial exactness.
+
+[heading Dilogarithm example]
+
+For a smooth special-function example, consider
+[$\int_0^1\!\int_0^1 (1-xy/2)^{-1}\,dx\,dy
+= \pi^2/6-\log^2 2].
+
+    using boost::math::quadrature::van_den_bos_unit_square;
+
+    double error = 0;
+    double value = van_den_bos_unit_square(
+        [](double x, double y)
+        {
+            return 1 / (1 - x * y / 2);
+        },
+        1e-10,
+        &error);
+
+[heading Square-lattice Green function]
+
+The resolvent of the nearest-neighbor square lattice gives the genuinely
+two-dimensional periodic integral
+[$G(z)=\int_0^1\!\int_0^1
+(z-\cos(2\pi x)-\cos(2\pi y))^{-1}\,dx\,dy]. For [^z > 2], its
+closed form is [$G(z)=2K(2/z)/(\pi z)], where [^K] is the complete elliptic
+integral of the first kind. For example, [^z = 3] can be computed with
+
+    double const z = 3;
+    double const two_pi = 2 * boost::math::constants::pi<double>();
+    double value = van_den_bos_unit_square(
+        [=](double x, double y)
+        {
+            return 1 / (z - std::cos(two_pi * x)
+                          - std::cos(two_pi * y));
+        },
+        1e-6);
+
+The lattice value [^z = 2], and the parameter value [^z = 1] in the
+dilogarithm family [$(1-zxy)^{-1}], give improper integrals. Singular
+integrands generally require special treatment and should not be expected to
+have reliable adjacent-level error estimates from this routine.
+
+The construction follows the positive, nested reduction approach of
+[@https://arxiv.org/abs/1809.09842 van den Bos et al.], adapted here to the
+uniform measure on the unit square.
+
+[endsect] [/section:van_den_bos_unit_square Van den Bos Unit-Square Quadrature]