More Gaussian Quadrature Methods

doc/quadrature/gauss_hermite.qbk

@@ -24,8 +24,14 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
       static auto integrate(F f, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()));
 
       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)));
+      static auto integrate(F f, const decltype(f(Real(0)))& zero, Norm norm, Real* pL1 = nullptr) -> decltype(f(Real(0)));
+
+      template <class F>
+      static auto integrate_averaged(F f, Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()));
+
+      template <class F, class Norm>
+      static auto integrate_averaged(F f, const decltype(f(Real(0)))& zero, Norm norm,
+                                     Real* error = nullptr, Real* pL1 = nullptr) -> decltype(f(Real(0)));
    };
 
    }}} // namespaces
@@ -40,9 +46,10 @@ using /Points/ evaluations of /f/, where the abscissas /x[sub i]/ are the zeros
 The rule is exact when /f/ is a polynomial of degree at most 2/Points/ - 1.
 Note that the weight function exp(-/x/[super 2]) is not part of /f/: pass only the factor that multiplies it.
 
-Like [link math_toolkit.gauss `gauss`], this is a "one shot" non-adaptive routine.
-It obtains no estimate of the error, and is intended for smooth /f/ that grow no faster than a polynomial,
+Like [link math_toolkit.gauss `gauss`], this is a "one shot" non-adaptive routine,
+intended for smooth /f/ that grow no faster than a polynomial,
 such as expectations of smooth functions of a normally distributed variable.
+For /Points/ + 1 more evaluations of /f/, `integrate_averaged` is exact to higher degree and estimates its error.
 
    static const RandomAccessContainer& abscissa();
    static const RandomAccessContainer& weights();
@@ -53,16 +60,38 @@ the rule is symmetric about zero.
    template <class F>
    static auto integrate(F f, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()));
 
-Integrates /f/(/x/) exp(-/x/[super 2]) over the real line, and optionally sets `*pL1` to the corresponding integral of |/f/|.
+Integrates /f/(/x/) exp(-/x/[super 2]) over the real line with the /Points/-point Gauss rule /G/. /f/ may return a real or complex value.
+
+If `pL1` is not null, `*pL1` is set to the corresponding integral of |/f/|.
 If that is substantially larger than the absolute value of the result, then the sum was ill-conditioned.
-/f/ may return a real or complex value.
 
    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)));
+   static auto integrate(F f, const decltype(f(Real(0)))& zero, Norm norm, Real* pL1 = nullptr) -> decltype(f(Real(0)));
 
 As for [link math_toolkit.gauss `gauss`], this overload handles vector- and matrix-valued integrands:
 `zero` is the additive identity of the returned type, and `norm` returns a nonnegative scalar convertible to `Real`.
+The L1 norm is measured with `norm`.
+
+   template <class F>
+   static auto integrate_averaged(F f, Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()));
+
+   template <class F, class Norm>
+   static auto integrate_averaged(F f, const decltype(f(Real(0)))& zero, Norm norm,
+                                  Real* error = nullptr, Real* pL1 = nullptr) -> decltype(f(Real(0)));
+
+Returns Laurie's averaged rule (/G/(/f/) + /A/(/f/))/2, at the cost of /Points/ + 1 more evaluations of /f/.
+/A/ is the anti-Gauss rule: the rule with /Points/ + 1 nodes whose error on every polynomial of degree up to 2/Points/ + 1
+is exactly that of /G/ with the opposite sign, so the average is exact for polynomials of degree 2/Points/ + 1.
+Its nodes interlace the Gauss nodes and its weights are positive.
+It is the Gauss rule of the same recurrence with its last coefficient doubled, so it is computed (the first time it is needed, then cached)
+by the same method as the Gauss rule.
+
+If `error` is not null, `*error` is set to |/A/(/f/) - /G/(/f/)|/2, but never less than 2[epsilon] times the L1 norm.
+When /G/ and /A/ bracket the integral, as they do for smooth /f/ once the error falls steadily with /Points/, this bounds the error of the average;
+it is then also within a few percent of the error of /G/.
+Otherwise it is a heuristic, and integrands that alias onto the null space of /A/ - /G/ can fool it.
+If `pL1` is not null, `*pL1` is set to the L1 norm of the averaged rule, measured with `norm` in the second overload.
+The averaged rule uses /Points/ + 1 points, so it overflows for one fewer point than `integrate` does.
 
 [heading Example]
 
@@ -78,7 +107,10 @@ so, for example, E\[cos /X/\] = exp(-[sigma][super 2]/2) cos [mu] is computed by
    double mu = 0.5;
    double sigma = 2;
    auto g = [&](double x) { return std::cos(mu + root_two<double>() * sigma * x); };
-   double expectation = gauss_hermite<double, 30>::integrate(g) / root_pi<double>();
+   double error;
+   double expectation = gauss_hermite<double, 30>::integrate_averaged(g, &error) / root_pi<double>();
+
+and `error / root_pi<double>()`, about 3e-16, estimates its error.
 
 [heading Choosing the number of points]
 
@@ -90,7 +122,7 @@ In double precision the computed abscissas are within a few ulps of the correctl
 for up to at least 200 points.
 
 In double precision the recurrence overflows for more than 706 points (the outer weights underflow long before that).
-`integrate` then raises an __evaluation_error through the /Policy/; if the policy does not throw, the result is a NaN.
+`integrate` and `integrate_averaged` then raise an __evaluation_error through the /Policy/; if the policy does not throw, the result is a NaN.
 
 As with `gauss`, defining `BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND` turns any combination of type and number of points without
 pre-computed values into a compile-time error.
@@ -101,5 +133,6 @@ and is easily modified for other numbers of points.
 
 * G. H. Golub and J. H. Welsch, ['Calculation of Gauss quadrature rules], Mathematics of Computation 23 (1969), 221-230.
 * W. Gautschi, ['Orthogonal Polynomials: Computation and Approximation], Oxford University Press, 2004.
+* D. P. Laurie, ['Anti-Gaussian quadrature formulas], Mathematics of Computation 65 (1996), no. 214, 739-747, [@https://doi.org/10.1090/S0025-5718-96-00713-2 doi:10.1090/S0025-5718-96-00713-2].
 
 [endsect] [/section:gauss_hermite Gauss-Hermite quadrature]

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

@@ -8,6 +8,7 @@
 #define BOOST_MATH_QUADRATURE_DETAIL_QUADRATURE_CONSTANT_HPP
 
 #ifndef BOOST_MATH_BUILD_MODULE
+#include <algorithm>
 #include <limits>
 #include <type_traits>
 #endif
@@ -54,6 +55,58 @@ struct quadrature_constant_category
       >;
 };
 
+template<typename RandomAccessContainer, typename Real, typename Policy>
+bool overflowed(const RandomAccessContainer& abscissa, const char* function, Real* pL1, Real* error, const Policy& pol)
+{
+   if (!(BOOST_MATH_NAMESPACE::isfinite)(static_cast<Real>(abscissa.back())))
+   {
+      *error = policies::raise_evaluation_error(
+               function,
+               "The recurrence relation overflows with this many points; use fewer points or a type with a wider exponent range.",
+               std::numeric_limits<Real>::quiet_NaN(), pol);
+      if (pL1)
+         *pL1 = *error;
+      return true;
+   }
+   return false;
+}
+
+// The averaged rule (G + A) / 2, exact for polynomials of degree 2N + 1. Half the difference is set as its error;
+// it bounds the error when G and A bracket the integral. L1 goes in as that of G and comes out as that of the average.
+template <class F, class RandomAccessContainer, class K, class Real, class Norm, class Policy>
+K anti_gauss_average(F& f, const RandomAccessContainer& abscissa, const RandomAccessContainer& weights, bool symmetric,
+                     K anti_result, const K& result, Real& L1, Norm norm, Real* error, const char* function, const Policy& pol)
+{
+   Real nan;
+   if (!(BOOST_MATH_NAMESPACE::isfinite)(L1) || overflowed(abscissa, function, static_cast<Real*>(nullptr), &nan, pol))
+   {
+      L1 = std::numeric_limits<Real>::quiet_NaN();
+      if (error)
+         *error = L1;
+      return anti_result * L1;
+   }
+   Real anti_L1 = 0;
+   for (unsigned i = 0; i < abscissa.size(); ++i)
+   {
+      Real x = abscissa[i];
+      K fx = f(x);
+      Real n = static_cast<Real>(norm(fx));
+      if (symmetric && x != 0)
+      {
+         K fmx = f(-x);
+         n += static_cast<Real>(norm(fmx));
+         fx += fmx;
+      }
+      anti_result += fx * static_cast<Real>(weights[i]);
+      anti_L1 += n * static_cast<Real>(weights[i]);
+   }
+   L1 = (L1 + anti_L1) / 2;
+   if (error)
+      *error = (std::max)(static_cast<Real>(norm(anti_result - result)) / 2, 2 * tools::epsilon<Real>() * L1);
+   anti_result += result;
+   return anti_result * Real(0.5);
+}
+
 } // namespace detail
 } // namespace quadrature
 BOOST_MATH_NAMESPACE_END