diff --git a/doc/quadrature/gauss_hermite.qbk b/doc/quadrature/gauss_hermite.qbk index 0d4e0fa120..11975da187 100644 --- a/doc/quadrature/gauss_hermite.qbk +++ b/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()(std::declval())); template - 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 + static auto integrate_averaged(F f, Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval()(std::declval())); + + template + 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,35 @@ the rule is symmetric about zero. template static auto integrate(F f, Real* pL1 = nullptr)->decltype(std::declval()(std::declval())); -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 - 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 + static auto integrate_averaged(F f, Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval()(std::declval())); + + template + 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 computed the first time it is needed, then cached. + +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. [heading Example] @@ -78,19 +104,20 @@ 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() * sigma * x); }; - double expectation = gauss_hermite::integrate(g) / root_pi(); + double error; + double expectation = gauss_hermite::integrate_averaged(g, &error) / root_pi(); + +and `error / root_pi()`, about 3e-16, estimates its error. [heading Choosing the number of points] `gauss_hermite` has pre-computed tables of abscissas and weights for 7, 10 and 15 points, at up to 110 decimal digits of precision. For other numbers of points, or higher precision types, the abscissas and weights are computed the first time they are needed and then cached. -Each zero is isolated by bisection with a Sturm sequence, then refined with Halley's method, all evaluated -with the three-term recurrence for orthonormal Hermite polynomials. -In double precision the computed abscissas are within a few ulps of the correctly rounded values and the weights within about one hundred, -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. +They are found in O(/Points/) operations by following the differential equation of the Hermite functions +from zero to zero in twice the working precision (Glaser, Liu and Rokhlin), so there is no limit from the overflow of the polynomials, +though the outer weights underflow; in double precision the abscissas are within an ulp and the weights within a few. +If the rule cannot be computed, `integrate` and `integrate_averaged` 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. @@ -100,6 +127,8 @@ and is easily modified for other numbers of points. [heading References] * G. H. Golub and J. H. Welsch, ['Calculation of Gauss quadrature rules], Mathematics of Computation 23 (1969), 221-230. +* A. Glaser, X. Liu and V. Rokhlin, ['A fast algorithm for the calculation of the roots of special functions], SIAM Journal on Scientific Computing 29 (2007), no. 4, 1420-1438, [@https://doi.org/10.1137/06067016X doi:10.1137/06067016X]. * 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] diff --git a/include/boost/math/quadrature/detail/quadrature_constant.hpp b/include/boost/math/quadrature/detail/quadrature_constant.hpp index ea09288a9d..a0352c876c 100644 --- a/include/boost/math/quadrature/detail/quadrature_constant.hpp +++ b/include/boost/math/quadrature/detail/quadrature_constant.hpp @@ -8,8 +8,14 @@ #define BOOST_MATH_QUADRATURE_DETAIL_QUADRATURE_CONSTANT_HPP #ifndef BOOST_MATH_BUILD_MODULE +#include #include +#include #include +#include +#ifdef BOOST_MATH_HAS_THREADS +#include +#endif #endif BOOST_MATH_NAMESPACE_BEGIN namespace quadrature{ namespace detail { @@ -54,6 +60,90 @@ struct quadrature_constant_category >; }; +template +bool overflowed(const RandomAccessContainer& abscissa, const char* function, Real* pL1, Real* error, const Policy& pol) +{ + if (!(BOOST_MATH_NAMESPACE::isfinite)(static_cast(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::quiet_NaN(), pol); + if (pL1) + *pL1 = *error; + return true; + } + return false; +} + +// A rule built by the first call of get(build), which fills its abscissa x and weights w; copies share it. +template +class lazy_rule +{ + struct rule + { +#ifdef BOOST_MATH_HAS_THREADS + std::once_flag once; +#else + bool built = false; +#endif + std::vector x, w; + }; + std::shared_ptr r = std::make_shared(); + +public: + template + const rule& get(Build build) const + { +#ifdef BOOST_MATH_HAS_THREADS + std::call_once(r->once, [&] { build(r->x, r->w); }); +#else + if (!r->built) + { + build(r->x, r->w); + r->built = true; + } +#endif + return *r; + } +}; + +// 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 +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(nullptr), &nan, pol)) + { + L1 = std::numeric_limits::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(norm(fx)); + if (symmetric && x != 0) + { + K fmx = f(-x); + n += static_cast(norm(fmx)); + fx += fmx; + } + anti_result += fx * static_cast(weights[i]); + anti_L1 += n * static_cast(weights[i]); + } + L1 = (L1 + anti_L1) / 2; + if (error) + *error = (std::max)(static_cast(norm(anti_result - result)) / 2, 2 * tools::epsilon() * L1); + anti_result += result; + return anti_result * Real(0.5); +} + } // namespace detail } // namespace quadrature BOOST_MATH_NAMESPACE_END diff --git a/include/boost/math/quadrature/gauss_hermite.hpp b/include/boost/math/quadrature/gauss_hermite.hpp index ec2d7ef2f4..20a1cdb920 100644 --- a/include/boost/math/quadrature/gauss_hermite.hpp +++ b/include/boost/math/quadrature/gauss_hermite.hpp @@ -34,189 +34,190 @@ #include #include #include +#include -BOOST_MATH_NAMESPACE_BEGIN namespace quadrature { namespace detail { - -#ifndef BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND +BOOST_MATH_NAMESPACE_BEGIN namespace detail { -template -class gauss_hermite_detail +template +class hermite_family { - static_assert(N > 0, "Gauss-Hermite quadrature needs at least one point."); +public: + hermite_family(unsigned n) : N(n), a_values(N), b_values(N) + { + using std::sqrt; + for (unsigned k = 0; k < N; ++k) + { + a_values[k] = sqrt(Real(2) / Real(k + 1)); + b_values[k] = -sqrt(Real(k) / Real(k + 1)); + } + p0_value = 1 / sqrt(sqrt(BOOST_MATH_NAMESPACE::constants::pi())); + derivative_scale = sqrt(Real(2) * Real(N)); + } - // The nonnegative zeros of the Hermite polynomial H_N, in increasing order, and their weights. - // Everything is evaluated with the orthonormal recurrence - // p_0 = pi^(-1/4), p_{k+1} = x sqrt(2/(k+1)) p_k - sqrt(k/(k+1)) p_{k-1}, - // whose values stay representable long after H_N itself overflows. Then p_N' = sqrt(2N) p_{N-1}, - // p_N'' = 2x p_N' - 2N p_N, and the weight is 2/p_N'(x)^2. Evaluating H_N from its monomial - // coefficients instead would be hopelessly ill-conditioned: the largest zero of H_100 has - // condition number ~1e13. - // - // By Sturm's theorem for orthogonal polynomials, the number of sign changes in p_0(x), ..., p_N(x) - // is the number of zeros of p_N above x. Bisecting on that count isolates each zero, and Halley's - // method then converges to it inside its bracket. Asymptotic initial guesses alone are not enough: - // for N = 215, even Halley's method from those of Numerical Recipes' gauher finds the wrong zero. - struct recurrence + const Real& a(unsigned k) const { - Real p; - Real p_prime; - unsigned zeros_above; - }; + return a_values[k]; + } - class orthonormal_hermite + const Real& b(unsigned k) const { - public: - orthonormal_hermite() : a(N), b(N) - { - using std::sqrt; - for (unsigned k = 0; k < N; ++k) - { - a[k] = sqrt(Real(2) / Real(k + 1)); - b[k] = sqrt(Real(k) / Real(k + 1)); - } - p0 = 1 / sqrt(sqrt(BOOST_MATH_NAMESPACE::constants::pi())); - derivative_scale = sqrt(Real(2 * N)); - } + return b_values[k]; + } - recurrence operator()(const Real& x) const - { - Real p = p0; - Real p_previous = 0; - bool last_negative = false; - unsigned sign_changes = 0; - for (unsigned k = 0; k < N; ++k) - { - Real p_next = x * a[k] * p - b[k] * p_previous; - p_previous = p; - p = p_next; - // A zero p_k(x) lies between values of opposite sign, so skipping it keeps the count right. - if ((p != 0) && ((p < 0) != last_negative)) - { - ++sign_changes; - last_negative = !last_negative; - } - } - return recurrence{ p, derivative_scale * p_previous, sign_changes }; - } + const Real& p0() const + { + return p0_value; + } - private: - std::vector a; - std::vector b; - Real p0; - Real derivative_scale; - }; + unsigned num_roots() const + { // Need to account for odd N having root at 0 + return N / 2 + N % 2; + } - // The weight is 2/p_N'(x)^2 at the exact zero, but d(log w)/dx is about -4x, so evaluating it at - // the rounded zero z costs far more accuracy than the rounding itself. Use p'' to move p' to the - // exact zero, which lies a Newton step delta = -p/p' away. - static Real weight(const Real& z, const recurrence& r) + Real lower_bound() const { - Real delta = -r.p / r.p_prime; - Real p_prime = r.p_prime + (2 * z * r.p_prime - 2 * Real(N) * r.p) * delta; - return 2 / (p_prime * p_prime); + return 0; } - static std::pair, std::vector > calculate_values() + Real upper_bound() const { - using std::abs; using std::sqrt; - const unsigned positive_zeros = N / 2; - const unsigned offset = N & 1; - std::vector x(positive_zeros + offset), w(positive_zeros + offset); - const orthonormal_hermite evaluate; - // All zeros lie below sqrt(2N + 1), so none lie above `upper`. Every p_k grows with x beyond - // its zeros, so if the recurrence overflows anywhere we need it, it overflows here: report - // that with NaNs, which gauss_hermite turns into an evaluation error. - Real upper = sqrt(Real(2 * N + 2)); - recurrence top = evaluate(upper); - if (!(BOOST_MATH_NAMESPACE::isfinite)(top.p) || !(BOOST_MATH_NAMESPACE::isfinite)(top.p_prime)) - { - std::fill(x.begin(), x.end(), std::numeric_limits::quiet_NaN()); - std::fill(w.begin(), w.end(), std::numeric_limits::quiet_NaN()); - return std::make_pair(x, w); - } - if (offset) - { - w[0] = weight(Real(0), evaluate(Real(0))); - } - // Find the zeros from the largest down: after each one, `upper` is a point with exactly t - 1 - // zeros above it. - for (unsigned t = 1; t <= positive_zeros; ++t) - { - Real lower = 0; - Real high = upper; - for (int i = 0; i < tools::digits(); ++i) - { - Real mid = (lower + high) / 2; - unsigned above = evaluate(mid).zeros_above; - if (above >= t) - lower = mid; - else - high = mid; - if (above == t) - break; - } - // Halley's method, falling back to bisection whenever a step would leave the bracket. - // With u = p/p' and v = p''/p' = 2z - 2N u, the step is u/(1 - uv/2); the ratios cannot - // overflow even where p'^2 would. - const bool lower_negative = evaluate(lower).p < 0; - Real z = (lower + high) / 2; - recurrence r = evaluate(z); - for (int i = 0; (r.p != 0) && (i < 4 * tools::digits()); ++i) - { - if ((r.p < 0) == lower_negative) - lower = z; - else - high = z; - Real u = r.p / r.p_prime; - Real v = 2 * z - 2 * Real(N) * u; - Real next = z - u / (1 - u * v / 2); - if (!((next > lower) && (next < high))) - next = (lower + high) / 2; - const bool converged = abs(next - z) <= 2 * tools::epsilon() * abs(next); - z = next; - r = evaluate(z); - if (converged) - break; - } - x[positive_zeros + offset - t] = z; - w[positive_zeros + offset - t] = weight(z, r); - upper = lower; - } - // Every node runs through the same rounded coefficients a and b, so the weights share a relative error - // of a few eps that grows with N (against 50 digits: -5 eps at N = 200, -10 at 500); exact coefficients - // remove it only with the recurrence in twice the precision. The exact weights sum to sqrt(pi), the - // integral of exp(-x^2), so scale them to that: weights with relative errors e_i then have e_i - e, e - // their weighted mean, plus two roundings. That removes the shared part, which is what integrals of - // smooth functions see, and moves no weight by more than |e| <= max |e_i|; measured, both the largest and - // the mean |e_i| fall for N = 200, 500 and 706. Sum from the smallest weight, with Neumaier's compensation: - BOOST_MATH_NAMESPACE::detail::compensated_sum sum; - for (std::size_t i = w.size(); i-- > 0;) - { - sum.add((offset && (i == 0)) ? w[i] : Real(2 * w[i])); - } - const Real scale = BOOST_MATH_NAMESPACE::constants::root_pi() / sum.value(); - for (Real& weight_i : w) - { - weight_i *= scale; - } - return std::make_pair(x, w); + return sqrt(Real(2) * Real(N) + Real(2)); + } + + Real area() const + { + return BOOST_MATH_NAMESPACE::constants::root_pi(); + } + + unsigned weight_multiplicity(unsigned i) const + { + return (N % 2 && i == 0) ? 1 : 2; + } + + Real next(const Real& x, const Real& p, const Real& p_previous, unsigned k) const + { + return a(k) * x * p + b(k) * p_previous; + } + + Real slope(unsigned k) const + { // The coefficient of x p in next() + return a(k); + } + + Real derivative(const Real&, const Real&, const Real& p_previous) const + { + return derivative_scale * p_previous; } - static const std::pair, std::vector >& values() + Real second_derivative(const Real& x, const Real& p, const Real& p_prime) const { - static const std::pair, std::vector > data = calculate_values(); - return data; + return Real(2) * x * p_prime - Real(2) * Real(N) * p; + } + + Real second_derivative_ratio(const Real& x, const Real& p, const Real& p_prime) const + { + return Real(2) * x - Real(2) * Real(N) * p / p_prime; + } + + Real weight(const Real& z, const Real& p, const Real& p_prime) const + { + Real delta = -p / p_prime; + Real p_prime_at_root = p_prime + second_derivative(z, p, p_prime) * delta; + return Real(2) / (p_prime_at_root * p_prime_at_root); + } + + // For glaser_liu_rokhlin: u = exp(-x^2 / 2) p_N satisfies u'' + (2N + 1 - x^2) u = 0; the coefficients of + // s(t + h) = 1 and g(t + h) as polynomials in h. Its solutions are entire. + static constexpr unsigned ends = 0; + + void normal_form(const double_word_arithmetic& dw, const double_word& t, double_word (&s)[5], double_word (&g)[3], unsigned m) const + { + s[0] = double_word{ Real(1), Real(0) }; + g[0] = dw.add(dw.negate(dw.multiply(t, t)), Real(2) * Real(m) + 1); + g[1] = double_word{ -2 * t.hi, -2 * t.lo }; + g[2] = double_word{ Real(-1), Real(0) }; + } + + Real singular_distance(const Real&) const + { + return tools::max_value(); + } + + Real weight_function(const double_word& z) const + { + using std::exp; + const double_word square = double_word_arithmetic().multiply(z, z); + return exp(-square.hi) * (1 - square.lo); + } + + Real sturm_liouville(const double_word&) const + { + return 1; + } + + // p_N = (A x + B) p_{N-1} - C p_{N-2}: A = sqrt(2 / N), B = 0, and C_2 = -2 C = -2 sqrt((N - 1) / N), in double words. + void anti_gauss_recurrence(const double_word_arithmetic& dw, double_word& A, double_word& B, double_word& C2) const + { + A = square_root(dw, dw.divide(double_word{ Real(2), Real(0) }, Real(N))); + B = double_word{ Real(0), Real(0) }; + const double_word C = square_root(dw, dw.divide(double_word{ Real(N - 1), Real(0) }, Real(N))); + C2 = double_word{ -2 * C.hi, -2 * C.lo }; + } + + // One step of Newton's method for the square root, in double words. + static double_word square_root(const double_word_arithmetic& dw, const double_word& x) + { + using std::sqrt; + const Real root = sqrt(x.hi); + Real p, e; + two_product()(root, root, p, e); + return dw.add(double_word{ root, Real(0) }, (((x.hi - p) - e) + x.lo) / (2 * root)); + } + + // For the anti-Gauss polynomial: p_{N-1}(0) / p_{N-2}'(0) = -1 / sqrt(2 (N - 1)) or p_{N-1}'(0) / p_{N-2}(0) = sqrt(2 (N - 1)), + // from p_k' = sqrt(2 k) p_{k-1} and p_{N-1}(0) = b(N - 2) p_{N-3}(0). + double_word anti_gauss_ratio(const double_word_arithmetic& dw, unsigned) const + { + const double_word r = square_root(dw, double_word{ 2 * Real(N - 1), Real(0) }); + return (N - 1) % 2 ? r : dw.divide(double_word{ Real(-1), Real(0) }, r); + } + +private: + unsigned N; + std::vector a_values; + std::vector b_values; + Real p0_value; + Real derivative_scale; +}; + +} // namespace detail +BOOST_MATH_NAMESPACE_END + +BOOST_MATH_NAMESPACE_BEGIN namespace quadrature { namespace detail { + +#ifndef BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND + +template +class gauss_hermite_detail +{ + static_assert(N > 0, "Gauss-Hermite quadrature needs at least one point."); + + static const BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial>& polynomial() + { + static const BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial> value(N); + return value; } public: static const std::vector& abscissa() { - return values().first; + return polynomial().abscissa(); } + static const std::vector& weights() { - return values().second; + return polynomial().weights(); } }; @@ -516,18 +517,17 @@ class gauss_hermite : public detail::gauss_hermite_detail::value>; - // Abscissas computed on demand are NaN when the recurrence overflows Real. - static bool overflowed() + // Laurie's anti-Gauss rule with N + 1 points, computed the first time integrate_averaged is called. + static const BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial>& anti_gauss() { - return !(BOOST_MATH_NAMESPACE::isfinite)(static_cast(base::abscissa().back())); + static const BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial> value(N + 1, BOOST_MATH_NAMESPACE::detail::hermite_family(N + 1), true); + return value; } - static Real overflow_error() + template + static K anti_gauss_average(F& f, K zero, const K& result, Real& L1, Norm norm, Real* error) { - return policies::raise_evaluation_error( - "boost::math::quadrature::gauss_hermite<%1%>::integrate", - "The Hermite recurrence overflows with this many points; use fewer points or a type with a wider exponent range.", - std::numeric_limits::quiet_NaN(), Policy()); + return detail::anti_gauss_average(f, anti_gauss().abscissa(), anti_gauss().weights(), true, zero, result, L1, norm, error, "boost::math::quadrature::gauss_hermite<%1%>::integrate_averaged", Policy()); } public: @@ -541,12 +541,10 @@ class gauss_hermite : public detail::gauss_hermite_detail::value, "The return type cannot be integral, it must be either a real or complex floating point type."); using std::abs; - if (overflowed()) + Real nan; + if (detail::overflowed(base::abscissa(), "boost::math::quadrature::gauss_hermite<%1%>::integrate", pL1, &nan, Policy())) { - Real error = overflow_error(); - if (pL1) - *pL1 = error; - return static_cast(error); + return static_cast(nan); } unsigned non_zero_start = 1; K result = Real(0); @@ -578,14 +576,10 @@ class gauss_hermite : public detail::gauss_hermite_detail::value, "The return type cannot be integral."); - if (overflowed()) + Real nan; + if (detail::overflowed(base::abscissa(), "boost::math::quadrature::gauss_hermite<%1%>::integrate", pL1, &nan, Policy())) { - // Preserve the scalar error policy, including non-throwing policies, - // and use scalar multiplication to produce a result with the right shape. - Real error = overflow_error(); - if (pL1) - *pL1 = error; - K result = zero * error; + K result = zero * nan; return result; } unsigned non_zero_start = 1; @@ -610,6 +604,31 @@ class gauss_hermite : public detail::gauss_hermite_detail + static auto integrate_averaged(F f, Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval()(std::declval())) + { + using K = decltype(f(Real(0))); + Real L1; + K result = integrate(f, &L1); + result = anti_gauss_average(f, K(Real(0)), result, L1, [](const K& z) { using std::abs; return abs(z); }, error); + if (pL1) + *pL1 = L1; + return result; + } + + template + static auto integrate_averaged(F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* error = nullptr, Real* pL1 = nullptr) + ->decltype(static_cast(norm(f(Real(0)))), f(Real(0))) + { + Real L1; + auto result = integrate(f, zero, norm, &L1); + result = anti_gauss_average(f, zero, result, L1, norm, error); + if (pL1) + *pL1 = L1; + return result; + } }; } // namespace quadrature diff --git a/include/boost/math/quadrature/gauss_jacobi.hpp b/include/boost/math/quadrature/gauss_jacobi.hpp new file mode 100644 index 0000000000..b0a9bd186d --- /dev/null +++ b/include/boost/math/quadrature/gauss_jacobi.hpp @@ -0,0 +1,468 @@ +// 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_GAUSS_JACOBI_HPP +#define BOOST_MATH_QUADRATURE_GAUSS_JACOBI_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include + +#ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR +#error "The header requires C++17 or later." +#endif + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#include +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +BOOST_MATH_NAMESPACE_BEGIN namespace detail { + +template +class jacobi_family +{ +public: + // The anti-Gauss rule can have nodes outside [-1, 1], but within 3/2: doubling the last coefficient b_n <= 1 + // of the Jacobi matrix, whose norm is at most 1, moves its eigenvalues by at most (sqrt(2) - 1) sqrt(b_n). + jacobi_family(unsigned n, Real alpha_, Real beta_, bool anti_gauss_ = false) : N(n), alpha(alpha_), beta(beta_), anti_gauss(anti_gauss_) {} + + Real p0() const + { + return 1; + } + + Real lower_bound() const + { + return anti_gauss ? Real(-1.5) : Real(-1); + } + + Real upper_bound() const + { + return anti_gauss ? Real(1.5) : Real(1); + } + + unsigned num_roots() const + { + return N; + } + + // 2^(alpha + beta + 1) B(alpha + 1, beta + 1). The sums alpha + 1, beta + 1 and their sum are rounded, and the gamma + // functions take that to first order, Gamma(s + e) = Gamma(s) (1 + psi(s) e); else the area, and with it every weight, + // moves by eps s psi(s), as by 7 eps for alpha = -0.9, beta = 10. The gamma functions are also more accurate than + // beta(), which is out by 9 eps at (11, 11), where they do not overflow. + Real area() const + { + using std::pow; + const Real a = alpha + 1; + const Real b = beta + 1; + const Real c = a + b; + const Real ea = BOOST_MATH_NAMESPACE::detail::two_sum_error(alpha, Real(1), a); + const Real eb = BOOST_MATH_NAMESPACE::detail::two_sum_error(beta, Real(1), b); + const Real ec = BOOST_MATH_NAMESPACE::detail::two_sum_error(a, b, c) + ea + eb; + const Real psi_c = BOOST_MATH_NAMESPACE::digamma(c); + const Real B = c < 100 + ? BOOST_MATH_NAMESPACE::tgamma(a) * BOOST_MATH_NAMESPACE::tgamma(b) / BOOST_MATH_NAMESPACE::tgamma(c) * (1 + BOOST_MATH_NAMESPACE::digamma(a) * ea + BOOST_MATH_NAMESPACE::digamma(b) * eb - psi_c * ec) + : BOOST_MATH_NAMESPACE::beta(a, b) * (1 + (BOOST_MATH_NAMESPACE::digamma(a) - psi_c) * ea + (BOOST_MATH_NAMESPACE::digamma(b) - psi_c) * eb); + return 2 * pow(Real(2), alpha) * pow(Real(2), beta) * B; + } + + unsigned weight_multiplicity(unsigned) const + { + return 1; + } + + Real next(const Real& x, const Real& p, const Real& p_previous, unsigned k) const + { // See https://dlmf.nist.gov/18.9 + if (k == 0) + { + return ((alpha - beta) + (alpha + beta + 2) * x) / 2 * p; + } + + const Real degree = static_cast(k); + const Real sum = alpha + beta; + const Real twice_degree_sum = 2 * degree + sum; + const Real first = (twice_degree_sum + 1) * + (twice_degree_sum * (twice_degree_sum + 2) * x + alpha * alpha - beta * beta) * p; + const Real second = 2 * (degree + alpha) * (degree + beta) * (twice_degree_sum + 2) * p_previous; + const Real denominator = 2 * (degree + 1) * (degree + sum + 1) * twice_degree_sum; + return (first - second) / denominator; + } + + Real slope(unsigned k) const + { // The coefficient of x p in next() + const Real degree = static_cast(k); + const Real sum = alpha + beta; + if (k == 0) + { + return (sum + 2) / 2; + } + return (2 * degree + sum + 1) * (2 * degree + sum + 2) / (2 * (degree + 1) * (degree + sum + 1)); + } + + Real derivative(const Real& x, const Real& p, const Real& p_previous) const + { // See https://dlmf.nist.gov/18.9#E17 + const Real degree = static_cast(N); + const Real sum = alpha + beta; + if (x == 1) + { + return p * degree * (degree + sum + 1) / (2 * (alpha + 1)); + } + if (x == -1) + { + return -p * degree * (degree + sum + 1) / (2 * (beta + 1)); + } + return (degree * (alpha - beta - (2 * degree + sum) * x) * p + + 2 * (degree + alpha) * (degree + beta) * p_previous) / ((2 * degree + sum) * (1 - x * x)); + } + + Real second_derivative(const Real& x, const Real& p, const Real& p_prime) const + { + const Real degree = static_cast(N); + const Real sum = alpha + beta; + return ((alpha - beta + (sum + 2) * x) * p_prime - + degree * (degree + sum + 1) * p) / (1 - x * x); + } + + Real second_derivative_ratio(const Real& x, const Real& p, const Real& p_prime) const + { + const Real degree = static_cast(N); + const Real sum = alpha + beta; + return (alpha - beta + (sum + 2) * x - degree * (degree + sum + 1) * p / p_prime) / (1 - x * x); + } + + Real weight(const Real& x, const Real& p, const Real& p_prime) const + { + Real delta = -p / p_prime; + Real p_prime_at_root = p_prime + second_derivative(x, p, p_prime) * delta; + return 1 / ((1 - x * x) * p_prime_at_root * p_prime_at_root); + } + + // For glaser_liu_rokhlin: u = (1 - x)^((alpha + 1) / 2) (1 + x)^((beta + 1) / 2) p_N satisfies s u'' + g u = 0 with + // s = (1 - x)^2 (1 + x)^2 and g = kappa (1 - x) (1 + x) + A (1 + x)^2 + B (1 - x)^2, where A = (1 - alpha^2) / 4, + // B = (1 - beta^2) / 4 and kappa = ((2N + alpha + beta + 1)^2 + 1 - alpha^2 - beta^2) / 4. The coefficients of + // s(t + h) and g(t + h) in h, from d = 1 - t and e = 1 + t, which keep their accuracy near the ends. + + void normal_form(const double_word_arithmetic& dw, const double_word& t, double_word (&s)[5], double_word (&g)[3], unsigned m) const + { + using word = double_word; + const word alpha2 = dw.multiply(word{ alpha, Real(0) }, alpha); + const word beta2 = dw.multiply(word{ beta, Real(0) }, beta); + const word A = dw.divide(dw.add(dw.negate(alpha2), Real(1)), Real(4)); + const word B = dw.divide(dw.add(dw.negate(beta2), Real(1)), Real(4)); + const word nu = dw.add(dw.two_sum(alpha, beta), Real(2) * Real(m) + 1); + const word kappa = dw.divide(dw.add(dw.add(dw.multiply(nu, nu), dw.negate(dw.add(alpha2, beta2))), Real(1)), Real(4)); + const word d = dw.add(dw.negate(t), Real(1)); + const word e = dw.add(t, Real(1)); + // (d - h)(e + h) = q0 + q1 h - h^2 + const word q0 = dw.multiply(d, e); + const word q1 = dw.add(d, dw.negate(e)); + s[0] = dw.multiply(q0, q0); + s[1] = dw.multiply(q0, word{ 2 * q1.hi, 2 * q1.lo }); + s[2] = dw.add(dw.multiply(q1, q1), word{ -2 * q0.hi, -2 * q0.lo }); + s[3] = word{ -2 * q1.hi, -2 * q1.lo }; + s[4] = word{ Real(1), Real(0) }; + const word Ae = dw.multiply(A, e); + const word Bd = dw.multiply(B, d); + g[0] = dw.add(dw.add(dw.multiply(kappa, q0), dw.multiply(Ae, e)), dw.multiply(Bd, d)); + const word difference = dw.add(Ae, dw.negate(Bd)); + g[1] = dw.add(dw.multiply(kappa, q1), word{ 2 * difference.hi, 2 * difference.lo }); + g[2] = dw.add(dw.add(A, B), dw.negate(kappa)); + } + + Real singular_distance(const Real& x) const + { + return (std::min)(Real(1 - x), Real(1 + x)); + } + + // (1 - z)^alpha (1 + z)^beta at z = hi + lo. + Real weight_function(const double_word& z) const + { + using std::pow; + const double_word d = double_word_arithmetic::add(double_word_arithmetic::negate(z), Real(1)); + const double_word e = double_word_arithmetic::add(z, Real(1)); + return pow(d.hi, alpha) * pow(e.hi, beta) * (1 + alpha * d.lo / d.hi + beta * e.lo / e.hi); + } + + // The series of P_N about either end, in powers of w with c_0 = 1, x = 1 - 2 lambda w or -1 + 2 lambda w: P_N(x) is a + // multiple of 2F1(-N, N + alpha + beta + 1; alpha + 1; (1 - x) / 2), and of the same with alpha and beta swapped in + // (1 + x) / 2. The first end is the one where the exponent of the weight is smaller. lambda is a power of 2 near + // 1 / (2N + alpha + beta + 1)^2, which puts the nearest zero at w of order 1. + static constexpr unsigned ends = 2; + + Real series_origin(unsigned end) const + { + return at_one(end) ? Real(1) : Real(-1); + } + + Real series_step(unsigned end) const + { + using std::frexp; + using std::ldexp; + int exponent; + frexp(2 * Real(N) + alpha + beta + 1, &exponent); + return ldexp(at_one(end) ? Real(-2) : Real(2), -2 * exponent); + } + + double_word series_ratio(const double_word_arithmetic& dw, unsigned k, unsigned end, unsigned m) const + { // c_{k+1} / c_k = lambda (k - m) (k + m + alpha + beta + 1) / ((k + 1) (k + 1 + alpha)) for P_m, alpha and beta swapped at -1 + using std::abs; + const double_word numerator = dw.multiply(dw.add(dw.two_sum(alpha, beta), Real(k + m + 1)), (Real(k) - Real(m)) * abs(series_step(end)) / 2); + const double_word denominator = dw.multiply(dw.two_sum(Real(k + 1), at_one(end) ? alpha : beta), Real(k + 1)); + return dw.divide(numerator, denominator); + } + + // P_N = (A x + B) P_{N-1} - C P_{N-2} (DLMF 18.9.2), with k = N - 1, s = alpha + beta and d = 2 (k + 1) (k + s + 1) (2k + s): + // A = (2k + s + 1)(2k + s)(2k + s + 2) / d, B = (2k + s + 1)(alpha^2 - beta^2) / d, C = 2 (k + alpha)(k + beta)(2k + s + 2) / d. + void anti_gauss_recurrence(const double_word_arithmetic& dw, double_word& A, double_word& B, double_word& C2) const + { + using word = double_word; + const Real k = Real(N - 1); + const word sum = dw.two_sum(alpha, beta); + const word m0 = dw.add(sum, 2 * k), m1 = dw.add(sum, 2 * k + 1), m2 = dw.add(sum, 2 * k + 2); + const word d = dw.multiply(dw.multiply(dw.add(sum, k + 1), m0), 2 * (k + 1)); + A = dw.divide(dw.multiply(dw.multiply(m1, m0), m2), d); + B = dw.divide(dw.multiply(m1, dw.multiply(dw.two_sum(alpha, -beta), sum)), d); + const word C = dw.divide(dw.multiply(dw.multiply(dw.two_sum(k, alpha), dw.two_sum(k, beta)), m2), d); + C2 = word{ -4 * C.hi, -4 * C.lo }; + } + + // P_{N-1}(e) / P_{N-2}(e) = (N - 1 + alpha) / (N - 1) at e = 1, -(N - 1 + beta) / (N - 1) at -1 + double_word anti_gauss_ratio(const double_word_arithmetic& dw, unsigned end) const + { + const double_word r = dw.divide(dw.two_sum(Real(N - 1), at_one(end) ? alpha : beta), Real(N - 1)); + return at_one(end) ? r : dw.negate(r); + } + + // sigma' / sigma = -(alpha + 1) / (2 (1 - x)) + (beta + 1) / (2 (1 + x)) + double_word scale_log_derivative(const double_word_arithmetic& dw, const double_word& x) const + { + const double_word d = dw.add(dw.negate(x), Real(1)), e = dw.add(x, Real(1)); + return dw.add(dw.divide(dw.two_sum(beta, Real(1)), double_word{ 2 * e.hi, 2 * e.lo }), + dw.negate(dw.divide(dw.two_sum(alpha, Real(1)), double_word{ 2 * d.hi, 2 * d.lo }))); + } + + // The leading coefficient 1 - x^2 of the Sturm-Liouville form ((1 - x)^(alpha + 1) (1 + x)^(beta + 1) p')' = -lambda w p. + Real sturm_liouville(const double_word& z) const + { + const double_word_arithmetic dw; + return dw.multiply(dw.add(dw.negate(z), Real(1)), dw.add(z, Real(1))).hi; + } + +private: + bool at_one(unsigned end) const + { + return (alpha <= beta) == (end == 0); + } + + unsigned N; + Real alpha; + Real beta; + bool anti_gauss; +}; + +} // namespace detail +BOOST_MATH_NAMESPACE_END + +BOOST_MATH_NAMESPACE_BEGIN namespace quadrature { namespace detail { + +template +class gauss_jacobi_detail +{ + unsigned N; + Real alpha; + Real beta; + std::vector abscissa_vals; + std::vector weight_vals; + bool anti_gauss_internal; + lazy_rule anti; + using polynomial = BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial >; +public: + gauss_jacobi_detail(unsigned N_, Real alpha_, Real beta_) : N(N_), alpha(alpha_), beta(beta_) + { + // gauss_jacobi raises the domain error. + if (!in_domain()) + { + abscissa_vals.assign(N, std::numeric_limits::quiet_NaN()); + weight_vals.assign(N, std::numeric_limits::quiet_NaN()); + anti_gauss_internal = false; + return; + } + BOOST_MATH_NAMESPACE::detail::jacobi_family family(N, alpha, beta); + polynomial value(N, family); + abscissa_vals = value.abscissa(); + weight_vals = value.weights(); + // Laurie 1996, Theorem 3: the anti-Gauss rule has a node outside [-1, 1] when alpha or beta is below -1/2 + // and N is large, and for some N when it is near -1/2. At -1/2 it can have nodes on -1 and 1. + anti_gauss_internal = polynomial::anti_gauss_is_internal(N + 1, family, Real(1)) && polynomial::anti_gauss_is_internal(N + 1, family, Real(-1)); + } + const std::vector& abscissa() const + { + return abscissa_vals; + } + + const std::vector& weights() const + { + return weight_vals; + } + + bool anti_gauss_is_internal() const + { + return anti_gauss_internal; + } + + // Laurie's anti-Gauss rule with N + 1 points, for integrate_averaged; empty unless it is internal. + const std::vector& anti_gauss_abscissa() const + { + return anti_gauss().x; + } + + const std::vector& anti_gauss_weights() const + { + return anti_gauss().w; + } + +protected: + const auto& anti_gauss() const + { + return anti.get([this](std::vector& x, std::vector& w) + { + if (!anti_gauss_internal) + return; + const polynomial rule(N + 1, BOOST_MATH_NAMESPACE::detail::jacobi_family(N + 1, alpha, beta, true), true); + x = rule.abscissa(); + w = rule.weights(); + // Nodes on the endpoints come out a few eps either side. + for (Real& z : x) + z = (std::min)((std::max)(z, Real(-1)), Real(1)); + }); + } + + bool in_domain() const + { + return N > 0 && alpha > -1 && beta > -1 && (BOOST_MATH_NAMESPACE::isfinite)(alpha) && (BOOST_MATH_NAMESPACE::isfinite)(beta); + } +}; + +} // namespace detail + +template > +class gauss_jacobi : public detail::gauss_jacobi_detail::value> +{ + template + K anti_gauss_average(F& f, K zero, const K& result, Real& L1, Norm norm, Real* error) const + { + if (this->in_domain() && !this->anti_gauss_is_internal()) + { + L1 = policies::raise_evaluation_error("boost::math::quadrature::gauss_jacobi<%1%>::integrate_averaged", + "The anti-Gauss rule has a node outside [-1, 1] for these alpha and beta, so there is no averaged rule.", + std::numeric_limits::quiet_NaN(), Policy()); + if (error) + *error = L1; + return zero * L1; + } + return detail::anti_gauss_average(f, this->anti_gauss_abscissa(), this->anti_gauss_weights(), false, zero, result, L1, norm, error, "boost::math::quadrature::gauss_jacobi<%1%>::integrate_averaged", Policy()); + } + +public: + gauss_jacobi(unsigned N_, Real alpha_, Real beta_) : detail::gauss_jacobi_detail::value>(N_, alpha_, beta_) + { + if (N_ == 0) + { + policies::raise_domain_error("boost::math::quadrature::gauss_jacobi<%1%>::gauss_jacobi", + "The number of points must be > 0, but got %1%.", Real(0), Policy()); + } + if (!(alpha_ > -1) || !(BOOST_MATH_NAMESPACE::isfinite)(alpha_)) + { + policies::raise_domain_error("boost::math::quadrature::gauss_jacobi<%1%>::gauss_jacobi", + "The parameter alpha must be finite and > -1, but got %1%.", alpha_, Policy()); + } + if (!(beta_ > -1) || !(BOOST_MATH_NAMESPACE::isfinite)(beta_)) + { + policies::raise_domain_error("boost::math::quadrature::gauss_jacobi<%1%>::gauss_jacobi", + "The parameter beta must be finite and > -1, but got %1%.", beta_, Policy()); + } + } + + // Integrates f(x) (1 - x)^alpha * (1 + x)^beta over [-1, 1]. + template + auto integrate(F f, Real* pL1 = nullptr) const->decltype(std::declval()(std::declval())) + { + using K = decltype(f(Real(0))); + return integrate(f, K(Real(0)), [](const K& z) { using std::abs; return abs(z); }, pL1); + } + + // Explicit zero and norm for vector- and matrix-valued integrands. + template + auto integrate(F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* pL1 = nullptr) const + ->decltype(static_cast(norm(f(Real(0)))), f(Real(0))) + { + using K = decltype(f(Real(0))); + static_assert(!std::is_integral::value, + "The return type cannot be integral, it must be either a real or complex floating point type."); + Real nan = std::numeric_limits::quiet_NaN(); + if (!this->in_domain() || detail::overflowed(this->abscissa(), "boost::math::quadrature::gauss_jacobi<%1%>::integrate", pL1, &nan, Policy())) + { + if (pL1) + *pL1 = nan; + K result = zero * nan; + return result; + } + K result = zero; + Real L1 = static_cast(norm(result)); + for (unsigned i = 0; i < this->abscissa().size(); ++i) + { + K fp = f(static_cast(this->abscissa()[i])); + result += fp * static_cast(this->weights()[i]); + L1 += static_cast(norm(fp)) * static_cast(this->weights()[i]); + } + if (pL1) + *pL1 = L1; + return result; + } + // Laurie's averaged rule (G + A) / 2: N + 1 more evaluations of f than integrate, exact for polynomials of degree + // 2N + 1. The error is half the difference of G and A, a bound on the error when they bracket the integral. + template + auto integrate_averaged(F f, Real* error = nullptr, Real* pL1 = nullptr) const->decltype(std::declval()(std::declval())) + { + using K = decltype(f(Real(0))); + return integrate_averaged(f, K(Real(0)), [](const K& z) { using std::abs; return abs(z); }, error, pL1); + } + + template + auto integrate_averaged(F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* error = nullptr, Real* pL1 = nullptr) const + ->decltype(static_cast(norm(f(Real(0)))), f(Real(0))) + { + Real L1; + auto result = integrate(f, zero, norm, &L1); + result = anti_gauss_average(f, zero, result, L1, norm, error); + if (pL1) + *pL1 = L1; + return result; + } +}; +} // namespace quadrature + +BOOST_MATH_NAMESPACE_END + +#endif // BOOST_MATH_QUADRATURE_GAUSS_JACOBI_HPP diff --git a/include/boost/math/quadrature/gauss_laguerre.hpp b/include/boost/math/quadrature/gauss_laguerre.hpp new file mode 100644 index 0000000000..3079f259d9 --- /dev/null +++ b/include/boost/math/quadrature/gauss_laguerre.hpp @@ -0,0 +1,807 @@ +// 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_GAUSS_LAGUERRE_HPP +#define BOOST_MATH_QUADRATURE_GAUSS_LAGUERRE_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include + +#ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR +#error "The header requires C++17 or later." +#endif + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#include +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +BOOST_MATH_NAMESPACE_BEGIN namespace detail { + +// The generalized Laguerre polynomials L_n^(alpha), orthogonal for x^alpha exp(-x) on [0, infinity). +template +class laguerre_family +{ +public: + laguerre_family(unsigned n, Real alpha_ = 0) : N(n), alpha(alpha_) {} + + Real p0() const + { + return 1; + } + + Real lower_bound() const + { + return 0; + } + + Real upper_bound() const + { // Gershgorin's theorem for the Jacobi matrix, diagonal 2k + alpha + 1 and off-diagonal sqrt(k (k + alpha)), + // bounds the zeros by 4N + 2 alpha; (3 + sqrt 2) N + 3 alpha + 3 also covers the anti-Gauss rule. + using std::sqrt; + return (3 + sqrt(Real(2))) * static_cast(N) + 3 * alpha + 3; + } + + unsigned num_roots() const + { + return N; + } + + // Gamma(alpha + 1), with the rounding error e of alpha + 1 = s + e to first order: Gamma(s) (1 + psi(s) e). + Real area() const + { + const Real s = alpha + 1; + return BOOST_MATH_NAMESPACE::tgamma(s) * (1 + BOOST_MATH_NAMESPACE::digamma(s) * BOOST_MATH_NAMESPACE::detail::two_sum_error(alpha, Real(1), s)); + } + + unsigned weight_multiplicity(unsigned) const + { + return 1; + } + + Real next(const Real& x, const Real& p, const Real& p_previous, unsigned k) const + { + const Real degree = static_cast(k); + // x enters as x p, not through 2k + 1 + alpha - x, which would round away the low bits of the small zeros. + return ((Real(2) * degree + 1 + alpha) * p - (degree + alpha) * p_previous - x * p) / (degree + 1); + } + + // p_n(x), or Laurie's p~_n(x) with the coefficient of p_{n-2} doubled in the last step, to about twice the + // working precision: each step carries the rounding errors of its products, sums and quotient forward (a + // compensated recurrence). Near the small zeros the terms of the recurrence are about n times p_n, so in + // working precision those zeros, and the weights from p_{n-1} there, are only good to about n eps / x. + Real compensated_value(const Real& x, unsigned n, bool anti, Real& previous) const + { + const BOOST_MATH_NAMESPACE::detail::two_product product; + Real p = p0(), e = 0, p_previous = 0, e_previous = 0; + for (unsigned k = 0; k < n; ++k) + { + const Real c = (anti && (k + 1 == n)) ? Real(2) : Real(1); + const Real d = static_cast(k + 1); + const Real a = static_cast(2 * k + 1) + alpha; + const Real ea = BOOST_MATH_NAMESPACE::detail::two_sum_error(static_cast(2 * k + 1), alpha, a); + Real b = static_cast(k) + alpha; + const Real eb = c * BOOST_MATH_NAMESPACE::detail::two_sum_error(static_cast(k), alpha, b); + b *= c; + Real ap, eap, bp, ebp, xp, exp; + product(a, p, ap, eap); + product(b, p_previous, bp, ebp); + product(x, p, xp, exp); + const Real s1 = ap - bp; + const Real s = s1 - xp; + const Real err = BOOST_MATH_NAMESPACE::detail::two_sum_error(ap, -bp, s1) + BOOST_MATH_NAMESPACE::detail::two_sum_error(s1, -xp, s) + + eap - ebp - exp + (a - x) * e + ea * p - b * e_previous - eb * p_previous; + const Real q = s / d; + Real qd, eqd; + product(q, d, qd, eqd); + p_previous = p; + e_previous = e; + p = q; + e = (((s - qd) - eqd) + err) / d; + } + previous = p_previous + e_previous; + return p + e; + } + + Real slope(unsigned k) const + { // The coefficient of x p in next() + return -1 / static_cast(k + 1); + } + + Real derivative(const Real& x, const Real& p, const Real& p_previous) const + { + return (static_cast(N) * p - (static_cast(N) + alpha) * p_previous) / x; + } + + Real second_derivative(const Real& x, const Real& p, const Real& p_prime) const + { + return ((x - alpha - 1) * p_prime - static_cast(N) * p) / x; + } + + Real second_derivative_ratio(const Real& x, const Real& p, const Real& p_prime) const + { + return (x - alpha - 1) / x - static_cast(N) / x * p / p_prime; + } + + Real weight(const Real& x, const Real& p, const Real& p_prime) const + { // Up to the factor Gamma(N + alpha + 1) / N!, which the scaling to area() supplies. + Real delta = -p / p_prime; + Real p_prime_at_root = p_prime + second_derivative(x, p, p_prime) * delta; + return Real(1) / (p_prime_at_root * p_prime_at_root) / x; + } + + // For glaser_liu_rokhlin: u = x^((alpha + 1) / 2) exp(-x / 2) p_N satisfies x^2 u'' + g(x) u = 0 with + // g = (1 - alpha^2) / 4 + (N + (alpha + 1) / 2) x - x^2 / 4; the coefficients of s(t + h) and g(t + h) in h. + static constexpr unsigned ends = 1; + + void normal_form(const double_word_arithmetic& dw, const double_word& t, double_word (&s)[5], double_word (&g)[3], unsigned m) const + { + using word = double_word; + const word c = dw.add(word{ Real(m), Real(0) }, dw.divide(dw.two_sum(alpha, Real(1)), Real(2))); + const word a = dw.divide(dw.add(dw.negate(dw.multiply(word{ alpha, Real(0) }, alpha)), Real(1)), Real(4)); + const word square = dw.multiply(t, t); + s[0] = square; + s[1] = word{ 2 * t.hi, 2 * t.lo }; + s[2] = word{ Real(1), Real(0) }; + g[0] = dw.add(dw.add(word{ -square.hi / 4, -square.lo / 4 }, dw.multiply(c, t)), a); + g[1] = dw.add(c, word{ -t.hi / 2, -t.lo / 2 }); + g[2] = word{ Real(-0.25), Real(0) }; + } + + Real singular_distance(const Real& x) const + { + return x; + } + + // x^alpha exp(-x) at z = hi + lo; past where exp(-x) alone underflows, the square of its square root. + Real weight_function(const double_word& z) const + { + using std::exp; + using std::pow; + const Real correction = 1 + (alpha / z.hi - 1) * z.lo; + if (z.hi < -tools::log_min_value() / 2) + return pow(z.hi, alpha) * exp(-z.hi) * correction; + const Real root = pow(z.hi, alpha / 2) * exp(-z.hi / 2); + return root * root * correction; + } + + // L_N^alpha(lambda w) = L_N^alpha(0) sum_k c_k w^k with c_0 = 1, lambda a power of 2 near 1 / (4N), so that the + // smallest zero is at w of order 1 (Gatteschi: x_1 ~ j_{alpha,1}^2 / (4N + 2 alpha + 2)). + Real series_origin(unsigned) const + { + return 0; + } + + Real series_step(unsigned) const + { + using std::frexp; + using std::ldexp; + int exponent; + frexp(4 * Real(N) + 2 * alpha + 2, &exponent); + return ldexp(Real(1), -exponent); + } + + double_word series_ratio(const double_word_arithmetic& dw, unsigned k, unsigned end, unsigned m) const + { // c_{k+1} / c_k = lambda (k - m) / ((k + 1) (k + 1 + alpha)) for L_m^alpha + const double_word denominator = dw.multiply(dw.two_sum(Real(k + 1), alpha), Real(k + 1)); + return dw.divide(double_word{ (Real(k) - Real(m)) * series_step(end), Real(0) }, denominator); + } + + // L_N = (A x + B) L_{N-1} - C L_{N-2} with A = -1 / N, B = (2N - 1 + alpha) / N, C = (N - 1 + alpha) / N; C_2 = -2 C. + void anti_gauss_recurrence(const double_word_arithmetic& dw, double_word& A, double_word& B, double_word& C2) const + { + A = dw.divide(double_word{ Real(-1), Real(0) }, Real(N)); + B = dw.divide(dw.two_sum(2 * Real(N) - 1, alpha), Real(N)); + C2 = dw.divide(dw.two_sum(-2 * Real(N - 1), -2 * alpha), Real(N)); + } + + // L_{N-1}^alpha(0) / L_{N-2}^alpha(0) = (N - 1 + alpha) / (N - 1) + double_word anti_gauss_ratio(const double_word_arithmetic& dw, unsigned) const + { + return dw.divide(dw.two_sum(Real(N - 1), alpha), Real(N - 1)); + } + + // sigma' / sigma = (alpha + 1) / (2 x) - 1 / 2 + double_word scale_log_derivative(const double_word_arithmetic& dw, const double_word& x) const + { + return dw.add(dw.divide(dw.two_sum(alpha, Real(1)), double_word{ 2 * x.hi, 2 * x.lo }), Real(-0.5)); + } + + // The leading coefficient x of the Sturm-Liouville form (x^(alpha + 1) exp(-x) p')' = -N x^alpha exp(-x) p. + Real sturm_liouville(const double_word& z) const + { + return z.hi + z.lo; + } + +private: + unsigned N; + Real alpha; +}; + +} // namespace detail +BOOST_MATH_NAMESPACE_END + +BOOST_MATH_NAMESPACE_BEGIN namespace quadrature { namespace detail { + +#ifndef BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND + +template +class gauss_laguerre_detail +{ + static_assert(N > 1, "gauss_laguerre requires N > 1"); + static const BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial>& polynomial() + { + static const BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial> value(N); + return value; + } + +public: + static const std::vector& abscissa() + { + return polynomial().abscissa(); + } + + static const std::vector& weights() + { + return polynomial().weights(); + } +}; + +#else + +template +class gauss_laguerre_detail; + +#endif + +#ifndef BOOST_HAS_FLOAT128 +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(1.93043676560362413838247885003823725e-01L), + static_cast(1.02666489533919195034519944317483354e+00L), + static_cast(2.56787674495074620690778622665959020e+00L), + static_cast(4.90035308452648456810171437810242722e+00L), + static_cast(8.18215344456286079108182755123396163e+00L), + static_cast(1.27341802917978137580126424581951091e+01L), + static_cast(1.93957278622625403117125820576302546e+01L), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(4.09318951701273902130432880017778795e-01L), + static_cast(4.21831277861719779929281005417498934e-01L), + static_cast(1.47126348657505278395374184636545717e-01L), + static_cast(2.06335144687169398657056149642012957e-02L), + static_cast(1.07401014328074552213195962843031379e-03L), + static_cast(1.58654643485642012687326223234055954e-05L), + static_cast(3.17031547899558056227132215385324202e-08L), + }; + return data; + } + }; + +#else +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(1.93043676560362413838247885003823725e-01Q), + static_cast(1.02666489533919195034519944317483354e+00Q), + static_cast(2.56787674495074620690778622665959020e+00Q), + static_cast(4.90035308452648456810171437810242722e+00Q), + static_cast(8.18215344456286079108182755123396163e+00Q), + static_cast(1.27341802917978137580126424581951091e+01Q), + static_cast(1.93957278622625403117125820576302546e+01Q), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(4.09318951701273902130432880017778795e-01Q), + static_cast(4.21831277861719779929281005417498934e-01Q), + static_cast(1.47126348657505278395374184636545717e-01Q), + static_cast(2.06335144687169398657056149642012957e-02Q), + static_cast(1.07401014328074552213195962843031379e-03Q), + static_cast(1.58654643485642012687326223234055954e-05Q), + static_cast(3.17031547899558056227132215385324202e-08Q), + }; + return data; + } + }; + +#endif +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.9304367656036241383824788500382372459120984891356750620225136201006415345782559871118003823365115612373429871679738e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0266648953391919503451994431748335351021966391136589333607107349300045456686682315951892100561068930429806129371487e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.5678767449507462069077862266595901968369581616444700992866063303240503321974487821383346200011905724664254743082485e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.9003530845264845681017143781024272182326822536691323824442733012926988370938012768358608940248312092759732409230278e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 8.1821534445628607910818275512339616296407668517524296643048342450041449025202649459386747106022315685557229645415476e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.2734180291797813758012642458195109090853336141446261917460200935582863612586646759328152143354140682014199320963541e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.9395727862262540311712582057630254604742850103460479496941123090856173616475344405452608383727847918520964087609689e+01), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.0931895170127390213043288001777879476694318320743384936496071108439327612531393650526756975329973583687149037045411e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.2183127786171977992928100541749893425900270593200299919010224347881378577765054361914173904057093683899914455920747e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.4712634865750527839537418463654571735447659494559067156700082178865140336466631595376166029179058453234935126953389e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.0633514468716939865705614964201295701345516595854480283021401106613334093662834206930986586836486163163369652374680e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0740101432807455221319596284303137904257521979552018605203061608685037526387650167316790636375375082365847554380883e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.5865464348564201268732622323405595386094807050091717842696789736074648974951001982867672006063759990800667569423195e-05), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.1703154789955805622713221538532420152314112706016551819590923622237092653696183497591858655360389258725422336556501e-08), + }; + return data; + } + }; + +#ifndef BOOST_HAS_FLOAT128 +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(1.37793470540492430830772505652711188e-01L), + static_cast(7.29454549503170498160373121676078781e-01L), + static_cast(1.80834290174031604823292007575060883e+00L), + static_cast(3.40143369785489951448253222140839068e+00L), + static_cast(5.55249614006380363241755848686876286e+00L), + static_cast(8.33015274676449670023876719727452218e+00L), + static_cast(1.18437858379000655649185389191416140e+01L), + static_cast(1.62792578313781020995326539358336223e+01L), + static_cast(2.19965858119807619512770901955944940e+01L), + static_cast(2.99206970122738915599087933407991952e+01L), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(3.08441115765020141547470834677860696e-01L), + static_cast(4.01119929155273551515780309912819515e-01L), + static_cast(2.18068287611809421588648523474646727e-01L), + static_cast(6.20874560986777473929021293135179537e-02L), + static_cast(9.50151697518110055383907219417199123e-03L), + static_cast(7.53008388587538775455964353675663902e-04L), + static_cast(2.82592334959956556742256382685002128e-05L), + static_cast(4.24931398496268637258657665974712355e-07L), + static_cast(1.83956482397963078092153522435593825e-09L), + static_cast(9.91182721960900855837754728324473606e-13L), + }; + return data; + } + }; + +#else +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(1.37793470540492430830772505652711188e-01Q), + static_cast(7.29454549503170498160373121676078781e-01Q), + static_cast(1.80834290174031604823292007575060883e+00Q), + static_cast(3.40143369785489951448253222140839068e+00Q), + static_cast(5.55249614006380363241755848686876286e+00Q), + static_cast(8.33015274676449670023876719727452218e+00Q), + static_cast(1.18437858379000655649185389191416140e+01Q), + static_cast(1.62792578313781020995326539358336223e+01Q), + static_cast(2.19965858119807619512770901955944940e+01Q), + static_cast(2.99206970122738915599087933407991952e+01Q), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(3.08441115765020141547470834677860696e-01Q), + static_cast(4.01119929155273551515780309912819515e-01Q), + static_cast(2.18068287611809421588648523474646727e-01Q), + static_cast(6.20874560986777473929021293135179537e-02Q), + static_cast(9.50151697518110055383907219417199123e-03Q), + static_cast(7.53008388587538775455964353675663902e-04Q), + static_cast(2.82592334959956556742256382685002128e-05Q), + static_cast(4.24931398496268637258657665974712355e-07Q), + static_cast(1.83956482397963078092153522435593825e-09Q), + static_cast(9.91182721960900855837754728324473606e-13Q), + }; + return data; + } + }; + +#endif +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3779347054049243083077250565271118810799168074578329416336651137344596447646208656254375241734116529355248965683739e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 7.2945454950317049816037312167607878107607273331224924600781104880626782307305157204652499178629631939583813469724208e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.8083429017403160482329200757506088332830602823714481516932097683732967664043428239153383606124012239270279810042947e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4014336978548995144825322214083906792731566142034732294267661021823269228666946410395095870847780590663426380062014e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 5.5524961400638036324175584868687628579740642873178138284904788336252636084290435460571874259935064862774598817665864e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 8.3301527467644967002387671972745221827094389720309893447446294283573954320532476951665824491858289576075525949818506e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.1843785837900065564918538919141613985802816909465335873339591495276220150697108232717212527603241705173109616546015e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.6279257831378102099532653935833622335255995603033059077788392165002827350932597940225841507596519165747859331720868e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.1996585811980761951277090195594493976806732340001887788241503511693106375798998823067686533346980361167369020416848e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.9920697012273891559908793340799195179710670577517960166104251135309849605268452639201572864373106556343888311203255e+01), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.0844111576502014154747083467786069562872888653833744211555718020018556202932557003040685460579742269270912496733327e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.0111992915527355151578030991281951479548361696211301757496775516766091500023676907078271471115370136487330141071994e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.1806828761180942158864852347464672674277853841218894056650398132078664598285192349000126323603209982063165104569642e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 6.2087456098677747392902129313517953695909065683802092368676760303923866462990276943718400176674956216762717600777500e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 9.5015169751811005538390721941719912258624504015797533646318970848680322451340750784240877901570876971156524042858250e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 7.5300838858753877545596435367566390179203914014362886506164014367761532732553587715637849496678837150424361526704126e-04), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.8259233495995655674225638268500212828033164744374677117086795259322508664275187881955981278223738401217175272290547e-05), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.2493139849626863725865766597471235464810801986441570081254604829021791344429695503912683040809284918604347723905036e-07), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.8395648239796307809215352243559382479826127765906584095132767065586955845890346516805867002792829341205430551248512e-09), + BOOST_MATH_HUGE_CONSTANT(T, 0, 9.9118272196090085583775472832447360645810946112447692431304366686254543563865419652261172586997161662742217301111110e-13), + }; + return data; + } + }; + +#ifndef BOOST_HAS_FLOAT128 +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(9.33078120172818047629030383672077003e-02L), + static_cast(4.92691740301883908960101791412401731e-01L), + static_cast(1.21559541207094946372992716487900297e+00L), + static_cast(2.26994952620374320247421741374817219e+00L), + static_cast(3.66762272175143727724905959436028773e+00L), + static_cast(5.42533662741355316534358132596053452e+00L), + static_cast(7.56591622661306786049739555811855455e+00L), + static_cast(1.01202285680191127347927394568171697e+01L), + static_cast(1.31302824821757235640991204175965313e+01L), + static_cast(1.66544077083299578225202408429835997e+01L), + static_cast(2.07764788994487667729157175675934362e+01L), + static_cast(2.56238942267287801445868285976698414e+01L), + static_cast(3.14075191697539385152432196202214042e+01L), + static_cast(3.85306833064860094162515167595161938e+01L), + static_cast(4.80260855726857943465734308507556623e+01L), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(2.18234885940086889856413236448110800e-01L), + static_cast(3.42210177922883329638948956806686032e-01L), + static_cast(2.63027577941680097414812275021951666e-01L), + static_cast(1.26425818105930535843030549378389129e-01L), + static_cast(4.02068649210009148415854789871163157e-02L), + static_cast(8.56387780361183836391575987649400017e-03L), + static_cast(1.21243614721425207621920522466788622e-03L), + static_cast(1.11674392344251941992578595518123419e-04L), + static_cast(6.45992676202290092465319025312414635e-06L), + static_cast(2.22631690709627263033182809179377262e-07L), + static_cast(4.22743038497936500735127949331007487e-09L), + static_cast(3.92189726704108929038460981948602319e-11L), + static_cast(1.45651526407312640633273963455309983e-13L), + static_cast(1.48302705111330133546164737187127769e-16L), + static_cast(1.60059490621113323104997812369887952e-20L), + }; + return data; + } + }; + +#else +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + static_cast(9.33078120172818047629030383672077003e-02Q), + static_cast(4.92691740301883908960101791412401731e-01Q), + static_cast(1.21559541207094946372992716487900297e+00Q), + static_cast(2.26994952620374320247421741374817219e+00Q), + static_cast(3.66762272175143727724905959436028773e+00Q), + static_cast(5.42533662741355316534358132596053452e+00Q), + static_cast(7.56591622661306786049739555811855455e+00Q), + static_cast(1.01202285680191127347927394568171697e+01Q), + static_cast(1.31302824821757235640991204175965313e+01Q), + static_cast(1.66544077083299578225202408429835997e+01Q), + static_cast(2.07764788994487667729157175675934362e+01Q), + static_cast(2.56238942267287801445868285976698414e+01Q), + static_cast(3.14075191697539385152432196202214042e+01Q), + static_cast(3.85306833064860094162515167595161938e+01Q), + static_cast(4.80260855726857943465734308507556623e+01Q), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + static_cast(2.18234885940086889856413236448110800e-01Q), + static_cast(3.42210177922883329638948956806686032e-01Q), + static_cast(2.63027577941680097414812275021951666e-01Q), + static_cast(1.26425818105930535843030549378389129e-01Q), + static_cast(4.02068649210009148415854789871163157e-02Q), + static_cast(8.56387780361183836391575987649400017e-03Q), + static_cast(1.21243614721425207621920522466788622e-03Q), + static_cast(1.11674392344251941992578595518123419e-04Q), + static_cast(6.45992676202290092465319025312414635e-06Q), + static_cast(2.22631690709627263033182809179377262e-07Q), + static_cast(4.22743038497936500735127949331007487e-09Q), + static_cast(3.92189726704108929038460981948602319e-11Q), + static_cast(1.45651526407312640633273963455309983e-13Q), + static_cast(1.48302705111330133546164737187127769e-16Q), + static_cast(1.60059490621113323104997812369887952e-20Q), + }; + return data; + } + }; + +#endif +template +class gauss_laguerre_detail +{ + using storage_type = typename quadrature_constant_category::storage_type; + public: + static std::array const & abscissa() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 9.3307812017281804762903038367207700278197076330453493009360137058347064015776121761668303610899058588907097107699871e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.9269174030188390896010179141240173099215246375371941601988508193506134664433011461933813698167861160364523544308066e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.2155954120709494637299271648790029749872681589377057427014003173216015106603771651065605856143208373528444729551797e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.2699495262037432024742174137481721885680048416478717117778472877455733945904832989178549502018189830165872698672666e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.6676227217514372772490595943602877290357430909410717845566428732064149379102861748582946304744040432223255333590429e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 5.4253366274135531653435813259605345208492623922927783354709435107483423740280199191227091479154617337151884199685196e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 7.5659162266130678604973955581185545482393108008999314491628694391139408365724029944489180550403629640393078390436542e+00), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.0120228568019112734792739456817169743410911967058744917457454601072487825047671876216088650721887243391350486173212e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.3130282482175723564099120417596531310187270022548453696047974546863746229290362972280038996276387001030969282397023e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.6654407708329957822520240842983599653951575341030608785764448525487749797539465693536916048016362945551853788626249e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.0776478899448766772915717567593436247705141256327195360002221322945959767932746815884396953509961946757428840574146e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.5623894226728780144586828597669841397123940322411476950317672315339824238931352351127917412081280502785586672176735e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.1407519169753938515243219620221404196345622809889669925953438870665440073910962255366662433543120684101302171932893e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.8530683306486009416251516759516193798932119744594784284050364282888975230612867789684690042672596606317674205933510e+01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.8026085572685794346573430850755662259393479711335534147707476887606535372312894457067945653339456838525028684441788e+01), + }; + return data; + } + static std::array const & weights() + { + static std::array data = { + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.1823488594008688985641323644811080021343324068445151647893418259657190410748792541872504691025036395044680815196202e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.4221017792288332963894895680668603182827472147567140620457276907638708199822947378530621505254574834996032558277546e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.6302757794168009741481227502195166616147477737942289338398907943112508686382961577835440203774163021418786314618248e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.2642581810593053584303054937838912877184803937971685684567321902219954509817127653766637180913751240893191067737467e-01), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.0206864921000914841585478987116315744354255151588087542482770026083846642284268861374498795123424591421789047455108e-02), + BOOST_MATH_HUGE_CONSTANT(T, 0, 8.5638778036118383639157598764940001725843892261220703970182442186170667853270532099641265249737915428749887172527299e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.2124361472142520762192052246678862218825774912216089252896089351714176877531295775371012761187312001078014473388155e-03), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.1167439234425194199257859551812341890362625693267480181637850779341432852535113805714378050258665226088107162969510e-04), + BOOST_MATH_HUGE_CONSTANT(T, 0, 6.4599267620229009246531902531241463548444527335439555949894778198068019741285302492718083846291056795408895968042879e-06), + BOOST_MATH_HUGE_CONSTANT(T, 0, 2.2263169070962726303318280917937726215071887584934068919788512991159987070128917425821426819036682039019642914030592e-07), + BOOST_MATH_HUGE_CONSTANT(T, 0, 4.2274303849793650073512794933100748675498159876301045574722416288222464623638254108884085708408168629608639640555277e-09), + BOOST_MATH_HUGE_CONSTANT(T, 0, 3.9218972670410892903846098194860231885836676205596680875788417907215759905742825209015082901147543134852484805575028e-11), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.4565152640731264063327396345530998285311241361150773694115515272365278527396127670028396473753138295863865512588769e-13), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.4830270511133013354616473718712776892067670532770945176556036443926719007518190902345311172254087343813449801330193e-16), + BOOST_MATH_HUGE_CONSTANT(T, 0, 1.6005949062111332310499781236988795197921487028677950643996761551559291234920430457494041052454552370257733201426544e-20), + }; + return data; + } + }; + +} // namespace detail + +template > +class gauss_laguerre : public detail::gauss_laguerre_detail::value> +{ + using base = detail::gauss_laguerre_detail::value>; + using family = BOOST_MATH_NAMESPACE::detail::laguerre_family; + using polynomial = BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial; + + Real alpha; + std::vector abscissa_vals; + std::vector weight_vals; + detail::lazy_rule anti; + + // Laurie's anti-Gauss rule with N + 1 points, for integrate_averaged. Its nodes are all positive (Laurie 1996, Theorem 2). + const auto& anti_gauss() const + { + return anti.get([this](std::vector& x, std::vector& w) + { + if (!(alpha > -1)) + { + x.assign(N + 1, std::numeric_limits::quiet_NaN()); + w.assign(N + 1, std::numeric_limits::quiet_NaN()); + return; + } + if (alpha == 0) + { + static const polynomial standard(N + 1, family(N + 1), true); + x = standard.abscissa(); + w = standard.weights(); + } + else + { + const polynomial rule(N + 1, family(N + 1, alpha), true); + x = rule.abscissa(); + w = rule.weights(); + } + BOOST_MATH_ASSERT(!(x[0] <= 0)); + }); + } + + template + K anti_gauss_average(F& f, K zero, const K& result, Real& L1, Norm norm, Real* error) const + { + return detail::anti_gauss_average(f, anti_gauss_abscissa(), anti_gauss_weights(), false, zero, result, L1, norm, error, "boost::math::quadrature::gauss_laguerre<%1%>::integrate_averaged", Policy()); + } + +public: + // The rule for x^alpha exp(-x); alpha = 0 uses the tables, or computes the rule once and caches it. + explicit gauss_laguerre(Real alpha_ = 0) : alpha(alpha_) + { + if (!(alpha > -1) || !(BOOST_MATH_NAMESPACE::isfinite)(alpha)) + { + alpha = policies::raise_domain_error("boost::math::quadrature::gauss_laguerre<%1%>::gauss_laguerre", + "The parameter alpha must be finite and > -1, but got %1%.", alpha_, Policy()); + abscissa_vals.assign(N, std::numeric_limits::quiet_NaN()); + weight_vals.assign(N, std::numeric_limits::quiet_NaN()); + } + else if (alpha == 0) + { + for (unsigned i = 0; i < base::abscissa().size(); ++i) + { + abscissa_vals.push_back(static_cast(base::abscissa()[i])); + weight_vals.push_back(static_cast(base::weights()[i])); + } + } + else + { + const polynomial rule(N, family(N, alpha)); + abscissa_vals = rule.abscissa(); + weight_vals = rule.weights(); + } + } + + const std::vector& abscissa() const + { + return abscissa_vals; + } + + const std::vector& weights() const + { + return weight_vals; + } + + // Laurie's anti-Gauss rule with N + 1 points. + const std::vector& anti_gauss_abscissa() const + { + return anti_gauss().x; + } + + const std::vector& anti_gauss_weights() const + { + return anti_gauss().w; + } + + // Integrates f(x) x^alpha exp(-x) over [0, infinity). + template + auto integrate(F f, Real* pL1 = nullptr) const->decltype(std::declval()(std::declval())) + { + using K = decltype(f(Real(0))); + return integrate(f, K(Real(0)), [](const K& z) { using std::abs; return abs(z); }, pL1); + } + + // Explicit zero and norm for vector- and matrix-valued integrands. + template + auto integrate(F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* pL1 = nullptr) const + ->decltype(static_cast(norm(f(Real(0)))), f(Real(0))) + { + using K = decltype(f(Real(0))); + static_assert(!std::is_integral::value, + "The return type cannot be integral, it must be either a real or complex floating point type."); + // The constructor has already raised the domain error for an invalid alpha. + Real nan = std::numeric_limits::quiet_NaN(); + if (!(alpha > -1) || detail::overflowed(abscissa_vals, "boost::math::quadrature::gauss_laguerre<%1%>::integrate", pL1, &nan, Policy())) + { + if (pL1) + *pL1 = nan; + K result = zero * nan; + return result; + } + K result = zero; + Real L1 = static_cast(norm(result)); + for (unsigned i = 0; i < abscissa_vals.size(); ++i) + { + K fp = f(abscissa_vals[i]); + result += fp * weight_vals[i]; + L1 += static_cast(norm(fp)) * weight_vals[i]; + } + if (pL1) + *pL1 = L1; + return result; + } + // Laurie's averaged rule (G + A) / 2: N + 1 more evaluations of f than integrate, exact for polynomials of degree + // 2N + 1. The error is half the difference of G and A, a bound on the error when they bracket the integral. + template + auto integrate_averaged(F f, Real* error = nullptr, Real* pL1 = nullptr) const->decltype(std::declval()(std::declval())) + { + using K = decltype(f(Real(0))); + return integrate_averaged(f, K(Real(0)), [](const K& z) { using std::abs; return abs(z); }, error, pL1); + } + + template + auto integrate_averaged(F f, const decltype(std::declval()(std::declval()))& zero, Norm norm, Real* error = nullptr, Real* pL1 = nullptr) const + ->decltype(static_cast(norm(f(Real(0)))), f(Real(0))) + { + Real L1; + auto result = integrate(f, zero, norm, &L1); + result = anti_gauss_average(f, zero, result, L1, norm, error); + if (pL1) + *pL1 = L1; + return result; + } +}; +} // namespace quadrature + +BOOST_MATH_NAMESPACE_END + +#endif // BOOST_MATH_QUADRATURE_GAUSS_LAGUERRE_HPP diff --git a/include/boost/math/special_functions/detail/orthogonal_polynomial.hpp b/include/boost/math/special_functions/detail/orthogonal_polynomial.hpp new file mode 100644 index 0000000000..12388e3307 --- /dev/null +++ b/include/boost/math/special_functions/detail/orthogonal_polynomial.hpp @@ -0,0 +1,898 @@ +// Copyright Jacob Hass 2026. +// Copyright Nick Thompson 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_DETAIL_ORTHOGONAL_POLYNOMIAL_HPP +#define BOOST_MATH_DETAIL_ORTHOGONAL_POLYNOMIAL_HPP + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#endif +#include +#include +#include +#include +#include +#include + +BOOST_MATH_NAMESPACE_BEGIN namespace detail { + +// The zeros of p_N and their Gauss weights in O(N) operations, after Glaser, Liu & Rokhlin, SIAM J. Sci. Comput. 29 +// (2007), or those of Laurie's anti-Gauss polynomial p~_N = (A x + B) p_{N-1} - 2 C p_{N-2} (where p_N = (A x + B) p_{N-1} +// - C p_{N-2}) and its weights. The family supplies the normal form s u'' + g u = 0 of the equation of u = sigma p_m, with +// sigma^2 the weight function W times the leading coefficient A_SL of the Sturm-Liouville form: u is bounded where p_m +// oscillates, unlike p_m. The zeros are those of F = (B + A x) u_0 + C_2 u_1: F = u_0 = sigma p_N for the Gauss rule, and +// u_0 = sigma p_{N-1}, u_1 = sigma p_{N-2}, C_2 = -2 C for the anti-Gauss rule. +// +// From a zero, a Pruefer transform of the equation of u_0 predicts the next zero of p_N, or brackets that of p~_N between +// the next two of p_{N-1}, which interlace with them (Laurie 1996). Newton's method on the Taylor series of F, whose +// coefficients follow from the equations, finds it, and the series carry u and u' there, where the Gauss weight is +// W / F'^2 and the anti-Gauss weight W A_SL / (u_0 F'), up to a constant. The series are summed in double words, as the +// rounding errors of the steps add up, and through W move the weights by thousands of eps. They converge out to the ends +// of a finite interval, and near those the zeros come instead from the series of the polynomials about the end. +// +// solve() returns false where it cannot vouch for the result, and the caller falls back on bisection. +template +class glaser_liu_rokhlin +{ + using arithmetic = double_word_arithmetic; + using word = double_word; + // u = sigma p_m at the anchor t, and its scaled Taylor coefficients a_k = H^k u^(k)(t) / k!, b_k = k (k - 1) a_k there. + struct solution + { + unsigned degree; + word u, u_prime, g[3]; + std::vector a, b; + std::size_t terms; + }; + +public: + glaser_liu_rokhlin(const Family& family_, unsigned n, bool anti_gauss) : family(family_), N(n), anti(anti_gauss) + { + using std::ldexp; + tolerance = ldexp(tools::epsilon() * tools::epsilon(), -2); + max_terms = 8 * static_cast(tools::digits()) + 64; + u[0].degree = anti ? N - 1 : N; + u[1].degree = N - 2; + u[1].u = u[1].u_prime = word{ Real(0), Real(0) }; + if (anti) + family.anti_gauss_recurrence(dw, A, B, C2); + } + + // The zeros in increasing order, the nonnegative ones only for a symmetric family, and their weights up to a constant factor. + bool solve(std::vector& x, std::vector& w) + { + BOOST_MATH_STD_USING + const unsigned num_roots = family.num_roots(); + x.clear(); + w.clear(); + word t; + int direction = 1; + bool from_zero = true; + // The zeros found from the far end, from the end inwards, and their weights up to another constant factor. + std::vector far_zeros; + std::vector far_weights; + if (!begin(std::integral_constant(), t, direction, from_zero, x, w, far_zeros, far_weights)) + return false; + // The innermost of the far zeros is found twice, and its two weights give the ratio of the constants. + unsigned count = num_roots - static_cast(x.size() + far_zeros.size()) + (far_zeros.empty() ? 0 : 1); + for (; count > 0; --count) + { + if (!step(t, direction, from_zero)) + return false; + x.push_back(t.hi); + w.push_back(weight(t)); + } + if (!far_zeros.empty()) + { + if (!(abs(x.back() - far_zeros.back().hi) <= 4 * tools::epsilon() * (std::max)(Real(abs(x.back())), Real(1)))) + return false; + const Real scale = w.back() / far_weights.back(); + for (std::size_t i = far_zeros.size() - 1; i-- > 0;) + { + x.push_back(far_zeros[i].hi); + w.push_back(far_weights[i] * scale); + } + } + if (direction < 0) + { + std::reverse(x.begin(), x.end()); + std::reverse(w.begin(), w.end()); + } + for (std::size_t i = 0; i < x.size(); ++i) + { + if (!(BOOST_MATH_NAMESPACE::isfinite)(x[i]) || !(BOOST_MATH_NAMESPACE::isfinite)(w[i]) || !(w[i] >= 0) + || ((i > 0) && !(x[i] > x[i - 1])) || (x[i] < family.lower_bound()) || (x[i] > family.upper_bound())) + return false; + } + return true; + } + +private: + const Family& family; + unsigned N; + bool anti; + word A{ Real(0), Real(0) }, B{ Real(1), Real(0) }, C2{ Real(0), Real(0) }; + arithmetic dw; + Real tolerance; + unsigned max_terms; + solution u[2]; + std::vector reciprocals; + // s = the product of (x - e)^2 over the ends e of the interval + static constexpr unsigned s_degree = 2 * Family::ends; + word s[5]; + Real last_step = 0; + + unsigned solutions() const + { + return anti ? 2 : 1; + } + + // B + A x + word factor(const word& x) const + { + return arithmetic::add(dw.multiply(x, A), B); + } + + // The weight at a zero z of F, where u and u' are known. + Real weight(const word& z) const + { + BOOST_MATH_STD_USING + const Real F_prime = A.hi * u[0].u.hi + factor(z).hi * u[0].u_prime.hi + C2.hi * u[1].u_prime.hi; + // The anti-Gauss weights share the sign of a constant factor. + return anti ? family.weight_function(z) * family.sturm_liouville(z) / abs(u[0].u.hi * F_prime) + : family.weight_function(z) / (F_prime * F_prime); + } + + // A symmetric family: p_m is even or odd, and so is u = sigma p_m. Start at 0, where F = 0 when it is odd. + bool begin(std::true_type, word& t, int& direction, bool& from_zero, std::vector& x, std::vector& w, + std::vector&, std::vector&) + { + t = word{ Real(0), Real(0) }; + direction = 1; + // p_{N-1} and p_{N-2} with the ratio of their first nonzero derivatives at 0. + const word rho = anti ? family.anti_gauss_ratio(dw, 0) : word{ Real(1), Real(0) }; + for (unsigned i = 0; i < solutions(); ++i) + { + const word one = i == 0 ? rho : word{ Real(1), Real(0) }; + const bool odd = u[i].degree % 2; + u[i].u = odd ? word{ Real(0), Real(0) } : one; + u[i].u_prime = odd ? one : word{ Real(0), Real(0) }; + } + from_zero = N % 2; + if (from_zero) + { + x.push_back(Real(0)); + w.push_back(weight(t)); + } + return true; + } + + // Otherwise the zeros nearest each end come from the series about it, and the march starts from the last of those at + // the first end, with sigma = 1 there. Then the weights of the zeros from the series, a constant over A_SL p_N'^2 or + // p_{N-1} p~_N', are W A_SL at the last zero times those over A_SL u_0'^2 or u_0 F'. + bool begin(std::false_type, word& t, int& direction, bool& from_zero, std::vector& x, std::vector& w, + std::vector& far_zeros, std::vector& far_weights) + { + std::vector zeros; + std::vector weights; + if (!series_zeros(0, N, zeros, weights, true)) + return false; + t = zeros.back(); + const Real scale = zeros.size() < N ? family.weight_function(t) * family.sturm_liouville(t) : Real(1); + for (std::size_t i = 0; i < zeros.size(); ++i) + { + x.push_back(zeros[i].hi); + w.push_back(scale * weights[i]); + } + direction = family.series_step(0) > 0 ? 1 : -1; + from_zero = true; + // At least one zero from the far end, and at most up to the last from the first; where they overlap, the march has + // nothing to add, and the last zero from the first end is the one found twice. + if ((Family::ends == 2) && (zeros.size() < N)) + series_zeros(1, N - static_cast(zeros.size()) + 1, far_zeros, far_weights, false); + return true; + } + + // The zeros of F nearest the given end, from the series of the p_m about it, in double words, and their weights up to a + // constant factor, at least one and at most max_count. Laguerre's method from the end, with the zeros already found + // divided out, climbs monotonically to the next, as F has N real zeros, all on one side; Newton's method finishes it. + // The series cancel more the further out the zero: the march takes over where they would cost more than a few eps, + // or where its steps need not be shortened, the zeros spaced by less than a third of their distance from the end (a + // quarter for the anti-Gauss rule, whose steps span about one and a half spacings). + // With anchor, u and u' at the last zero are set for the march. + bool series_zeros(unsigned end, unsigned max_count, std::vector& zeros, std::vector& weights, bool anchor) + { + BOOST_MATH_STD_USING + const Real eps = tools::epsilon(); + const Real c_step = family.series_step(end); + const Real origin = family.series_origin(end); + // F = (B + A x) rho y_0 + C_2 y_1 with y_m = p_m / p_m(end), x = origin + c_step v. + const word rho = anti ? family.anti_gauss_ratio(dw, end) : word{ Real(1), Real(0) }; + // The ratios r_k = c_{k+1} / c_k of the coefficients of each y in v: the terms t_k = c_k v^k are summed from them, + // as the coefficients themselves fall out of range long before the terms do at the outer zeros. + std::vector r[2]; + auto ratio = [&](unsigned i, std::size_t k) -> word + { + while (r[i].size() <= k) + r[i].push_back(family.series_ratio(dw, static_cast(r[i].size()), end, u[i].degree)); + return r[i][k]; + }; + // y, y' and y'' in v at v, until two terms in a row fall below the double-word precision, then F, dF/dv and d2F/dv2, + // and the sum of the magnitudes of the terms of F; in double words, or for words = false in working precision, which + // is enough for Laguerre's method except where F, as for the anti-Gauss rule, is a small difference of its terms. + word y[2][3], F[3]; + Real magnitude; + auto evaluate_series = [&](const word& v, bool words) -> bool + { + magnitude = 0; + for (unsigned i = 0; i < solutions(); ++i) + { + const word scale = i == 0 ? rho : C2; + word term{ Real(1), Real(0) }, sum[3] = {}; + Real big = 0; + for (std::size_t k = 0, small = 0; ; ++k) + { + for (unsigned d = 0; d < 3; ++d) + { + const Real c = d == 0 ? Real(1) : d == 1 ? Real(k) : Real(k) * (Real(k) - 1); + // y'' only steers Laguerre's method, which needs it in double words only for the anti-Gauss rule. + sum[d] = (words && ((d < 2) || anti)) ? arithmetic::add(sum[d], dw.multiply(term, c)) : word{ sum[d].hi + c * term.hi, Real(0) }; + } + magnitude += abs(term.hi * scale.hi); + big = (std::max)(big, abs(term.hi)); + small = ((k > 1) && (abs(term.hi) <= tolerance * big)) ? small + 1 : 0; + if ((small == 2) || (k > u[i].degree)) + break; + if (k >= max_terms) + return false; + term = words ? dw.multiply(dw.multiply(term, ratio(i, k)), v) : word{ term.hi * ratio(i, k).hi * v.hi, Real(0) }; + } + // y' and y'' in v from v y' and v^2 y'', or at v = 0 from the coefficients; y_0 times rho. + y[i][0] = sum[0]; + y[i][1] = v.hi == 0 ? ratio(i, 0) : dw.divide(sum[1], v); + y[i][2] = v.hi == 0 ? dw.multiply(dw.multiply(ratio(i, 0), ratio(i, 1)), Real(2)) : dw.divide(dw.divide(sum[2], v), v); + for (unsigned d = 0; (i == 0) && (d < 3); ++d) + y[0][d] = dw.multiply(y[0][d], rho); + } + const word b = factor(arithmetic::add(dw.multiply(v, c_step), origin)); + const word a{ A.hi * c_step, A.lo * c_step }; + F[0] = dw.multiply(y[0][0], b); + F[1] = arithmetic::add(dw.multiply(y[0][0], a), dw.multiply(y[0][1], b)); + F[2] = arithmetic::add(dw.multiply(y[0][1], word{ 2 * a.hi, 2 * a.lo }), dw.multiply(y[0][2], b)); + for (unsigned d = 0; anti && (d < 3); ++d) + F[d] = arithmetic::add(F[d], dw.multiply(y[1][d], C2)); + return true; + }; + zeros.clear(); + weights.clear(); + std::vector found; + for (unsigned count = 0; count < max_count; ++count) + { + // From the end, then from just past the last zero, clear of the next and of the zeros divided out (v is scaled to + // put the first zeros of p_N about 1 apart). + const Real gap = count == 0 ? Real(0) : found.back() - (count > 1 ? found[count - 2] : Real(0)); + word v{ count == 0 ? Real(0) : found.back() + (std::max)(gap, Real(1)) / 64, Real(0) }; + const Real n = static_cast(N - count); + // To about sqrt(eps), and one step more, which would take it to eps but for rounding errors where the terms + // cancel, as for the anti-Gauss rule; Newton's method finishes. + for (int iteration = 0, close = 0; (iteration < 100) && (close < (anti ? 1 : 2)); ++iteration) + { + // Just past a zero, where v starts, the terms of F cancel, and the steps need it to a few digits whatever the + // precision: in float they would jump past the next zero. Use double words where those digits are lost. + if (!evaluate_series(v, anti) || (!anti && (iteration == 0) && !(magnitude * eps <= abs(F[0].hi) / 65536) && !evaluate_series(v, true))) + return !zeros.empty(); + if (F[0].hi == 0) + break; + Real G = F[1].hi / F[0].hi; + Real G2 = G * G - F[2].hi / F[0].hi; + for (const Real& z : found) + { + G -= 1 / (v.hi - z); + G2 -= 1 / ((v.hi - z) * (v.hi - z)); + } + const Real d = (std::max)(Real((n - 1) * (n * G2 - G * G)), Real(0)); + const Real next = v.hi - n / (G < 0 ? G - sqrt(d) : G + sqrt(d)); + close += (close > 0) || (abs(next - v.hi) <= sqrt(eps) * abs(next)); + v = word{ next, Real(0) }; + } + // A zero on the end itself, as anti-Gauss rules have for an exponent -1/2 there. + const bool on_end = (count == 0) && (abs(v.hi) <= 16 * eps); + if (on_end) + v = word{ Real(0), Real(0) }; + else if (!(v.hi > (found.empty() ? Real(0) : found.back()))) + return !zeros.empty(); + // Newton's method in double words, until a step of a few eps, after which y, y' and F' follow to first order. + bool converged = on_end && evaluate_series(v, true); + for (int iteration = 0; (iteration < 6) && !converged; ++iteration) + { + // The terms of the series cancel; their rounding errors must not reach the working precision. + if (!evaluate_series(v, true) || !(magnitude * eps <= abs(F[1].hi * v.hi) / 64)) + return !zeros.empty(); + const Real delta = -(F[0].hi + F[0].lo) / F[1].hi; + v = arithmetic::add(v, delta); + converged = abs(delta) <= 16 * eps * v.hi; + for (unsigned i = 0; converged && (i < solutions()); ++i) + { + y[i][0] = arithmetic::add(y[i][0], dw.multiply(y[i][1], delta)); + y[i][1] = arithmetic::add(y[i][1], dw.multiply(y[i][2], delta)); + } + if (converged) + F[1] = arithmetic::add(F[1], dw.multiply(F[2], delta)); + } + if (!converged) + return !zeros.empty(); + const word z = arithmetic::add(word{ c_step * v.hi, c_step * v.lo }, origin); + // F' and y' in x. + const Real F_prime = F[1].hi / c_step; + zeros.push_back(z); + weights.push_back(1 / abs(F_prime * (anti ? y[0][0].hi : family.sturm_liouville(z) * F_prime))); + found.push_back(v.hi); + if (anchor) + { + // u = sigma y with sigma = 1 at z: u' = y' + (sigma' / sigma) y. + const word r_sigma = family.scale_log_derivative(dw, z); + for (unsigned i = 0; i < solutions(); ++i) + { + u[i].u = y[i][0]; + u[i].u_prime = arithmetic::add(dw.divide(y[i][1], c_step), dw.multiply(y[i][0], r_sigma)); + } + last_step = abs(c_step * (v.hi - (count > 0 ? found[count - 1] : Real(0)))); + } + if ((count > 0) && ((anti ? 4 : 3) * (v.hi - found[count - 1]) <= v.hi)) + break; + } + return true; + } + + // f = g / s for u_0 and f' at t + h, in working precision. + void f(const Real& h, Real& value, Real& derivative) const + { + Real sv = 0, sd = 0, gv = 0, gd = 0; + for (unsigned i = s_degree + 1; i-- > 0;) + { + sd = sd * h + sv; + sv = sv * h + s[i].hi; + } + for (unsigned i = 3; i-- > 0;) + { + gd = gd * h + gv; + gv = gv * h + u[0].g[i].hi; + } + value = gv / sv; + derivative = (gd - value * sd) / sv; + } + + // The step h from the anchor to where the Pruefer angle of u_0, tan(theta) = sqrt(f) u / u' (f = g / s), reaches + // theta1 from theta0: dx/dtheta = 1 / (sqrt(f) + f' sin(2 theta) / (4 f)) is nearly constant where f varies slowly on the + // scale of the oscillation, and one step of Runge-Kutta per half turn predicts the zero to a fraction of a percent. + // NaN past a turning point, where there is no next zero. + Real predict(const Real& theta0, const Real& theta1) const + { + BOOST_MATH_STD_USING + const int steps = (std::max)(BOOST_MATH_NAMESPACE::itrunc(ceil(abs(theta1 - theta0) / constants::pi() - Real(0.01))), 1); + const Real d = (theta1 - theta0) / steps; + auto slope = [&](const Real& h, const Real& theta) -> Real + { + Real value, f_prime; + f(h, value, f_prime); + const Real denominator = sqrt(value) + f_prime / (4 * value) * sin(2 * theta); + return denominator > 0 ? 1 / denominator : std::numeric_limits::quiet_NaN(); + }; + Real h = 0, theta = theta0; + for (int i = 0; i < steps; ++i) + { + const Real k1 = slope(h, theta); + const Real k2 = slope(h + d / 2 * k1, theta + d / 2); + const Real k3 = slope(h + d / 2 * k2, theta + d / 2); + const Real k4 = slope(h + d * k3, theta + d); + h += d / 6 * (k1 + 2 * k2 + 2 * k3 + k4); + theta += d; + } + return h; + } + + // The scaled Taylor coefficients of a solution, until the terms fall below the double-word precision for + // |tau| <= tau_max: from s u'' + g u = 0, + // b_{m+2} = -(sum_{i>=1} s_i H^i b_{m+2-i} + sum_i g_i H^{i+2} a_{m-i}) / s_0, a_{m+2} = b_{m+2} / ((m + 2) (m + 1)). + // Returns false if they do not fall. + bool taylor(solution& v, const Real& H, const Real& tau_max) + { + BOOST_MATH_STD_USING + word st[5], gt[3]; + const word inverse = dw.divide(word{ Real(1), Real(0) }, s[0]); + Real power = H; + for (unsigned i = 1; i <= s_degree; ++i, power *= H) + st[i] = dw.multiply(word{ s[i].hi * power, s[i].lo * power }, inverse); + power = H * H; + for (unsigned i = 0; i < 3; ++i, power *= H) + gt[i] = dw.multiply(word{ v.g[i].hi * power, v.g[i].lo * power }, inverse); + std::vector& a = v.a; + std::vector& b = v.b; + a.resize((std::max)(a.size(), std::size_t(64))); + b.resize(a.size()); + a[0] = v.u; + a[1] = word{ v.u_prime.hi * H, v.u_prime.lo * H }; + b[0] = b[1] = word{ Real(0), Real(0) }; + Real big = (std::max)(abs(a[0].hi), abs(a[1].hi) * tau_max); + power = tau_max; + unsigned small = 0; + for (unsigned m = 0; m + 2 < max_terms; ++m) + { + if (m + 2 >= a.size()) + { + a.resize(2 * (m + 2)); + b.resize(2 * (m + 2)); + } + if (m >= reciprocals.size()) + reciprocals.push_back(dw.divide(word{ Real(1), Real(0) }, static_cast((m + 2) * (m + 1)))); + // From the oldest coefficients to the newest, which is then the only one the sum waits for. + Real sum = 0, error = 0; + for (unsigned i = (std::min)(2u, m) + 1; i-- > 0;) + dw.multiply_add(sum, error, gt[i], a[m - i]); + for (unsigned i = (std::min)(s_degree, m + 2); i >= 1; --i) + dw.multiply_add(sum, error, st[i], b[m + 2 - i]); + // Left unnormalized, which the products above take in their stride, to keep the renormalizations off the + // dependency chain from one coefficient to the next. + b[m + 2] = word{ -sum, -error }; + Real product = 0, product_error = 0; + dw.multiply_add(product, product_error, b[m + 2], reciprocals[m]); + a[m + 2] = word{ product, product_error }; + power *= tau_max; + const Real term = abs(a[m + 2].hi) * power; + big = (std::max)(big, term); + // Coefficients can vanish, as the odd ones do for an even u about 0: two in a row must be small. + small = (term <= tolerance * big) ? small + 1 : 0; + if (small == 2) + { + v.terms = m + 3; + return true; + } + } + return false; + } + + // A series at tau and its first derivative in double words, by compensated Horner, its second in working precision. + void evaluate(const solution& v, const Real& tau, word& value, word& slope, Real& curvature) const + { + Real y = 0, d = 0, y_error = 0, d_error = 0; + curvature = 0; + for (std::size_t k = v.terms; k > 0; --k) + { + curvature = curvature * tau + 2 * d; + dw.horner_step(d, d_error, tau, word{ y, y_error }); + dw.horner_step(y, y_error, tau, v.a[k - 1]); + } + value = arithmetic::two_sum(y, y_error); + slope = arithmetic::two_sum(d, d_error); + } + + // F and F' = dF/dtau at t + H tau in working precision. + void combination(const word& t, const Real& H, const Real& tau, Real& F, Real& F_prime) const + { + F = F_prime = 0; + for (unsigned i = 0; i < solutions(); ++i) + { + Real y = 0, d = 0; + for (std::size_t k = u[i].terms; k-- > 0;) + { + d = d * tau + y; + y = y * tau + u[i].a[k].hi; + } + const Real c = i == 0 ? B.hi + A.hi * (t.hi + H * tau) : C2.hi; + F += c * y; + F_prime += c * d + (i == 0 ? A.hi * H * y : Real(0)); + } + } + + // From the anchor t, a zero of F (or for a symmetric family, possibly 0), to the next zero in the direction. + bool step(word& t, int direction, bool& from_zero) + { + BOOST_MATH_STD_USING + const Real eps = tools::epsilon(); + const Real pi = constants::pi(); + for (unsigned i = 0; i < solutions(); ++i) + family.normal_form(dw, t, s, u[i].g, u[i].degree); + // The bracket [lo, hi] in h for the next zero, and a first guess. + Real f0, f0_prime; + f(Real(0), f0, f0_prime); + // At a zero of the Gauss F = u_0, theta = 0. + const Real theta0 = (from_zero && !anti) ? Real(0) : (u[0].u_prime.hi == 0 ? pi / 2 : atan(sqrt(f0) * u[0].u.hi / u[0].u_prime.hi)); + Real lo, hi, guess; + if (!anti) + { + // The next multiple of pi, past t if t is a zero. + guess = predict(theta0, (from_zero || (theta0 * direction > 0)) ? direction * pi : Real(0)); + lo = 3 * guess / 4; + hi = 5 * guess / 4; + } + else + { + // Between the next two zeros of p_{N-1}, the first of them t itself where it is one; past the last of them, as + // far again as the last step. + const Real theta_lo = ((u[0].u.hi == 0) || (theta0 * direction <= 0)) ? Real(0) : direction * pi; + lo = u[0].u.hi == 0 ? Real(0) : predict(theta0, theta_lo); + hi = predict(theta0, theta_lo + direction * pi); + guess = predict(theta0, theta_lo + direction * pi / 2); + // Past the last zero of p_{N-1}, and near the turning point before it, the prediction fails; the last zero of + // p~_N is then within about a step. + if (!(BOOST_MATH_NAMESPACE::isfinite)(lo * hi * guess)) + { + const Real step = last_step > 0 ? last_step : pi / sqrt(f0); + if (!(BOOST_MATH_NAMESPACE::isfinite)(lo)) + lo = direction * step / 8; + hi = lo + 3 * direction * step / 2; + guess = (lo + hi) / 2; + } + // Allow for the errors of the prediction. + const Real margin = abs(hi - lo) / 64; + lo = direction * (std::max)(Real(lo * direction - margin), Real(0)); + hi += direction * margin; + } + if (!(BOOST_MATH_NAMESPACE::isfinite)(lo * hi * guess) || !(lo * direction >= 0) || !(hi * direction > lo * direction)) + return false; + // Within 5/8 of the distance to the nearest singular point the series converge at least like (5/8)^k; the zeros + // nearer the ends came from the series about them, so the steps here stay inside. + if (!(8 * abs(hi) <= 5 * family.singular_distance(t.hi))) + return false; + int exponent; + frexp(hi, &exponent); + const Real H = direction * ldexp(Real(1), exponent); + for (unsigned i = 0; i < solutions(); ++i) + if (!taylor(u[i], H, hi / H)) + return false; + // Newton's method in working precision; for the anti-Gauss rule, whose guess is rougher, falling back on bisection + // within the bracket, where F changes sign. + Real tau_lo = lo / H, tau_hi = hi / H, tau = guess / H, F, F_prime, F_lo = -1, F_hi = 1; + if (anti) + { + combination(t, H, tau_lo, F_lo, F_prime); + combination(t, H, tau_hi, F_hi, F_prime); + if (!((F_lo < 0) != (F_hi < 0))) + return false; + } + const bool lo_negative = F_lo < 0; + bool converged = false; + for (int i = 0; (i < 100) && !converged; ++i) + { + combination(t, H, tau, F, F_prime); + ((F < 0) == lo_negative ? tau_lo : tau_hi) = tau; + // Close enough for the step in double words below to square the error; F has a floor of rounding errors above + // eps where its terms cancel, as for the anti-Gauss rule. + const Real correction = F / F_prime; + converged = abs(correction) <= pow(eps, Real(0.75)) * abs(tau); + if (!converged) + { + tau -= correction; + if (anti && !((tau - tau_lo) * (tau - tau_hi) < 0)) + tau = (tau_lo + tau_hi) / 2; + } + } + // The zero must be the one predicted: the next, not this one or the one after. + if (!converged || !((tau - lo / H) * (tau - hi / H) < 0)) + return false; + // F and its slope in double words at tau give the zero to about twice the working precision, and u, u' there to + // second and first order in the correction delta, of order eps tau. + word value[2], slope[2]; + Real curvature[2]; + for (unsigned i = 0; i < solutions(); ++i) + evaluate(u[i], tau, value[i], slope[i], curvature[i]); + const word c = factor(arithmetic::add(t, word{ H * tau, Real(0) })); + word F_word = dw.multiply(value[0], c); + Real F_slope = A.hi * H * value[0].hi + c.hi * slope[0].hi; + if (anti) + { + F_word = arithmetic::add(F_word, dw.multiply(value[1], C2)); + F_slope += C2.hi * slope[1].hi; + } + const Real delta = -(F_word.hi + F_word.lo) / F_slope; + const word tz = arithmetic::fast_two_sum(tau, delta); + for (unsigned i = 0; i < solutions(); ++i) + { + u[i].u = arithmetic::add(arithmetic::add(value[i], dw.multiply(slope[i], delta)), curvature[i] * delta * delta / 2); + const word d = arithmetic::add(slope[i], curvature[i] * delta); + u[i].u_prime = word{ d.hi / H, d.lo / H }; + } + last_step = abs(tz.hi * H); + t = arithmetic::add(t, word{ tz.hi * H, tz.lo * H }); + from_zero = true; + return true; + } +}; + +template > +class orthogonal_polynomial +{ +public: + orthogonal_polynomial(unsigned n) : N(n) + { + values = calculate_values(Family(N)); + } + + // With anti_gauss, Laurie's anti-Gauss rule with n nodes, the companion of the (n - 1)-point Gauss rule. + orthogonal_polynomial(unsigned n, const Family& family, bool anti_gauss = false) : N(n), anti(anti_gauss) + { + values = calculate_values(family); + } + + // Laurie 1996, Theorem 1: the anti-Gauss rule with n nodes has no node at or beyond the end e of the support + // iff p~_n(e) has the sign of p_n(e), with p~_n defined as in strum_counter. Allow a margin, so that nodes + // on the endpoint, which Chebyshev weights put there exactly, pass despite rounding; f is defined there. + // An overflowed recurrence passes too, so that the rule reports the overflow. + static bool anti_gauss_is_internal(unsigned n, const Family& family, const Real& e) + { + Real p = family.p0(); + Real p_previous = 0; + for (unsigned k = 0; k + 1 < n; ++k) + { + Real p_next = family.next(e, p, p_previous, k); + p_previous = p; + p = p_next; + } + return !(family.next(e, p, 2 * p_previous, n - 1) / family.next(e, p, p_previous, n - 1) <= -4 * static_cast(n) * tools::epsilon()); + } + + const std::vector& abscissa() const + { + return values.first; + } + const std::vector& weights() const + { + return values.second; + } + +private: + unsigned N; + bool anti = false; + std::pair, std::vector > values; + struct recurrence + { + Real p; + Real p_prime; + unsigned sign_changes; + Real p_previous; + Real p_prime_previous; + Real p_double_prime; + }; + + class strum_counter + { + public: + strum_counter(unsigned n, const Family& family_, bool anti_) : N(n), family(family_), anti(anti_) {} + + recurrence operator()(const Real& x) const + { + Real p = family.p0(); + Real p_previous = 0; + // p~_N satisfies no differential equation, so differentiate the recurrence instead. + Real p_prime = 0, p_prime_previous = 0, p_double_prime = 0, p_double_prime_previous = 0; + bool last_negative = false; + unsigned sign_changes = 0; + for (unsigned k = 0; k < N; ++k) + { + // Laurie's anti-Gauss rule is the Gauss rule of the same recurrence with the coefficient of p_{N-2} + // in its last step doubled: p~_N = p_N - b_{N-1} p_{N-2} in monic form. The Sturm count still holds. + const Real c = (anti && (k + 1 == N)) ? Real(2) : Real(1); + Real p_next = family.next(x, p, c * p_previous, k); + if (anti) + { + const Real slope = family.slope(k); + Real p_double_prime_next = family.next(x, p_double_prime, c * p_double_prime_previous, k) + 2 * slope * p_prime; + Real p_prime_next = family.next(x, p_prime, c * p_prime_previous, k) + slope * p; + p_double_prime_previous = p_double_prime; + p_double_prime = p_double_prime_next; + p_prime_previous = p_prime; + p_prime = p_prime_next; + } + p_previous = p; + p = p_next; + // A zero p_k(x) lies between values of opposite sign, so skipping it keeps the count right. + if ((p != 0) && ((p < 0) != last_negative)) + { + ++sign_changes; + last_negative = !last_negative; + } + } + if (anti) + return recurrence{ p, p_prime, sign_changes, p_previous, p_prime_previous, p_double_prime }; + return recurrence{ p, family.derivative(x, p, p_previous), sign_changes, p_previous, 0, 0 }; + } + + Real second_derivative_ratio(const Real& x, const recurrence& r) const + { + return anti ? r.p_double_prime / r.p_prime : family.second_derivative_ratio(x, r.p, r.p_prime); + } + + Real weight(const Real& z, const recurrence& r) const + { + // The anti-Gauss rule is the Gauss rule of the functional 2I - G_{N-1}, so for monic P its weights are + // 2 ||P_{N-1}||^2 / (P_{N-1} P~'_N). The family's normalization changes them by a constant factor, + // which the scaling to area() below removes. Like the family weights, evaluate them at the exact zero + // to first order: near the ends of the interval they are sensitive to the rounding of z. + if (!anti) + return family.weight(z, r.p, r.p_prime); + const Real delta = -r.p / r.p_prime; + return 1 / ((r.p_previous + r.p_prime_previous * delta) * (r.p_prime + r.p_double_prime * delta)); + } + + // Replaces r.p, and for the Gauss rule r.p_previous and r.p_prime, with values to about twice the working + // precision where the family can provide them; returns false where it cannot. + template + auto compensated(const Real& z, recurrence& r, int) const -> decltype(std::declval().compensated_value(z, 0u, false, r.p), bool()) + { + Real p_previous; + const Real p = family.compensated_value(z, N, anti, p_previous); + if (!(BOOST_MATH_NAMESPACE::isfinite)(p) || !(BOOST_MATH_NAMESPACE::isfinite)(p_previous)) + return false; + r.p = p; + if (!anti) + { + r.p_previous = p_previous; + r.p_prime = family.derivative(z, p, p_previous); + } + return true; + } + + bool compensated(const Real&, recurrence&, long) const + { + return false; + } + + private: + unsigned N; + Family family; + bool anti; + }; + + // The zeros of p_N in increasing order, the nonnegative ones only for a family symmetric about 0, and their weights; + // where glaser_liu_rokhlin cannot vouch for them, from the family's recurrence. By Sturm's theorem the number of sign + // changes in p_0(x), ..., p_N(x) changes by one at each zero of p_N, in a direction that depends on the family: + // bisecting on that count isolates each zero, and Halley's method then converges to it inside its bracket. + std::pair, std::vector > calculate_values(const Family& family) + { + if (N == 0) + { + static const char* function = "boost::math::detail::orthogonal_polynomial<%1%>::orthogonal_polynomial"; + BOOST_MATH_NAMESPACE::policies::raise_domain_error(function, "The argument n used for Gauss quadrature must be > 0 (got n=%1%).", static_cast(N), Policy()); + std::vector x = { std::numeric_limits::quiet_NaN() }; + std::vector w = { std::numeric_limits::quiet_NaN() }; + return std::make_pair(x, w); + } + using std::abs; + using std::sqrt; + const strum_counter evaluate(N, family, anti); + + const unsigned num_roots = family.num_roots(); + unsigned interior_roots = num_roots; + std::vector x(num_roots), w(num_roots); + { + std::vector xg, wg; + if (glaser_liu_rokhlin(family, N, anti).solve(xg, wg)) + return scaled(family, xg, wg); + } + + Real upper = family.upper_bound(); + Real lower = family.lower_bound(); + // All zeros lie below upper_bound(), so none lie above `upper`. Every p_k grows with x beyond + // its zeros, so if the recurrence overflows anywhere we need it, it overflows here: report + // that with NaNs, which turns into an evaluation error. + recurrence top = evaluate(upper); + if (!(BOOST_MATH_NAMESPACE::isfinite)(top.p) || !(BOOST_MATH_NAMESPACE::isfinite)(top.p_prime)) + { + std::fill(x.begin(), x.end(), std::numeric_limits::quiet_NaN()); + std::fill(w.begin(), w.end(), std::numeric_limits::quiet_NaN()); + return std::make_pair(x, w); + } + // Check if lower bound is 0. If so, we have an odd number of roots, and the middle one is at 0. + const recurrence bottom = evaluate(lower); + if (bottom.p == 0) + { + x[0] = lower; + w[0] = evaluate.weight(lower, bottom); + interior_roots--; + } + // Need to know whether the number of sign changes increases or decreases with x, to know which direction to bisect. + const bool count_increases = top.sign_changes > bottom.sign_changes; + // Find the zeros from the largest down, following the variation-count direction of this family. + for (unsigned t = 1; t <= interior_roots; ++t) + { + Real low = lower; + Real high = upper; + const unsigned target = count_increases + ? top.sign_changes - t + 1 + : top.sign_changes + t; + for (int i = 0; i < tools::digits(); ++i) + { + Real mid = (low + high) / 2; + unsigned variations = evaluate(mid).sign_changes; + if (count_increases) + { + if (variations >= target) + high = mid; + else + low = mid; + } + else + { + if (variations >= target) + low = mid; + else + high = mid; + } + } + // Halley's method, falling back to bisection whenever a step would leave the bracket. + // With u = p/p' and v = p''/p' from the family, the step is u/(1 - uv/2); the ratios cannot + // overflow even where p'^2 would. + const bool lower_negative = evaluate(low).p < 0; + Real z = (low + high) / 2; + recurrence r = evaluate(z); + for (int i = 0; (r.p != 0) && (i < 4 * tools::digits()); ++i) + { + if ((r.p < 0) == lower_negative) + low = z; + else + high = z; + Real u = r.p / r.p_prime; + Real v = evaluate.second_derivative_ratio(z, r); + Real next = z - u / (1 - u * v / 2); + if (!((next > low) && (next < high))) + next = (low + high) / 2; + const bool converged = abs(next - z) <= 2 * tools::epsilon() * abs(next); + z = next; + r = evaluate(z); + if (converged) + break; + } + // Halley's method has the zero of p_N as rounded; where the family can evaluate p_N more accurately, + // Newton's method on that value finds the zero itself, and the weight uses those values too. + if (evaluate.compensated(z, r, 0)) + { + for (int i = 0; i < 2; ++i) + { + // Not the bracket: it came from the signs of the rounded p_N, which are wrong this close to the zero. + const Real next = z - r.p / r.p_prime; + if (!(abs(next - z) <= sqrt(tools::epsilon()) * abs(z)) || (next == z)) + break; + z = next; + r = evaluate(z); + evaluate.compensated(z, r, 0); + } + } + x[num_roots - t] = z; + w[num_roots - t] = evaluate.weight(z, r); + upper = low; + } + return scaled(family, x, w); + } + + // The weights share a relative error of a few eps from the rounded coefficients, growing with N. The exact weights + // sum to area(), so scale them to that: weights with relative errors e_i then have e_i - e, e their weighted mean, + // plus two roundings. That removes the shared part, which is what integrals of smooth functions see, and moves no + // weight by more than |e| <= max |e_i|. Sum from the smallest weight, with Neumaier's compensation: + static std::pair, std::vector > scaled(const Family& family, std::vector& x, std::vector& w) + { + BOOST_MATH_NAMESPACE::detail::compensated_sum sum; + for (std::size_t i = w.size(); i-- > 0;) + { + sum.add(Real(family.weight_multiplicity(static_cast(i))) * w[i]); + } + const Real scale = family.area() / sum.value(); + for (Real& weight_i : w) + { + weight_i *= scale; + } + return std::make_pair(x, w); + } +}; + +} // namespace detail +BOOST_MATH_NAMESPACE_END + +#endif // BOOST_MATH_DETAIL_ORTHOGONAL_POLYNOMIAL_HPP diff --git a/include/boost/math/tools/compensated_arithmetic.hpp b/include/boost/math/tools/compensated_arithmetic.hpp index 57b566cbca..fef6ea70dd 100644 --- a/include/boost/math/tools/compensated_arithmetic.hpp +++ b/include/boost/math/tools/compensated_arithmetic.hpp @@ -70,15 +70,15 @@ template class two_product { public: - two_product() : splitter(0) + two_product() : fused(is_fused()), splitter(0) { BOOST_MATH_STD_USING - if (!is_fused()) + if (!fused) splitter = ldexp(T(1), (tools::digits() + 1) / 2) + 1; } void operator()(const T& a, const T& b, T& p, T& e) const { - if (is_fused()) + if (fused) two_product_fma::evaluate(a, b, T(0), p, 0); else p = a * b; @@ -87,7 +87,7 @@ class two_product // The rounding error of ab = fl(a b): a b = ab + error(a, b, ab) exactly. T error(const T& a, const T& b, const T& ab) const { - if (is_fused()) + if (fused) { T e; two_product_fma::evaluate(a, b, T(-ab), e, 0); @@ -119,6 +119,7 @@ class two_product }(); return result; } + bool fused; T splitter; }; @@ -149,6 +150,101 @@ BOOST_MATH_GPU_ENABLED inline T two_sum_error(T a, T b, T s) return (a - (s - bb)) + (b - bb); } +#ifndef BOOST_MATH_HAS_NVRTC + +// +// Double-word numbers hi + lo, |lo| <= ulp(hi) / 2, and the "sloppy" arithmetic of Dekker (1971) on them (see Joldes, +// Muller & Popescu 2017): its error is about u^2 times the operands, u the unit roundoff, rather than times the result, +// which is as good as twice the working precision wherever the operands are not much larger than the result. +// +template +struct double_word +{ + T hi; + T lo; +}; + +template +class double_word_arithmetic +{ +public: + using value_type = double_word; + // a + b when |a| >= |b| or a = 0. + static value_type fast_two_sum(const T& a, const T& b) + { + const T s = a + b; + return value_type{ s, b - (s - a) }; + } + static value_type two_sum(const T& a, const T& b) + { + const T s = a + b; + return value_type{ s, two_sum_error(a, b, s) }; + } + static value_type negate(const value_type& a) + { + return value_type{ -a.hi, -a.lo }; + } + static value_type add(const value_type& a, const value_type& b) + { + const value_type s = two_sum(a.hi, b.hi); + return fast_two_sum(s.hi, s.lo + (a.lo + b.lo)); + } + static value_type add(const value_type& a, const T& b) + { + const value_type s = two_sum(a.hi, b); + return fast_two_sum(s.hi, s.lo + a.lo); + } + value_type multiply(const value_type& a, const value_type& b) const + { + T p, e; + product(a.hi, b.hi, p, e); + return fast_two_sum(p, e + (a.hi * b.lo + a.lo * b.hi)); + } + value_type multiply(const value_type& a, const T& b) const + { + T p, e; + product(a.hi, b, p, e); + return fast_two_sum(p, e + a.lo * b); + } + value_type divide(const value_type& a, const T& b) const + { + const T q = a.hi / b; + T p, e; + product(q, b, p, e); + return fast_two_sum(q, (((a.hi - p) - e) + a.lo) / b); + } + value_type divide(const value_type& a, const value_type& b) const + { + const T q = a.hi / b.hi; + const value_type r = add(a, negate(multiply(b, q))); + return fast_two_sum(q, r.hi / b.hi); + } + // s + e += a b, gathering the rounding errors in e unnormalized, so that only s = s + fl(a b) is a dependency between + // consecutive calls (as in Ogita, Rump & Oishi's Dot2). + void multiply_add(T& s, T& e, const value_type& a, const value_type& b) const + { + T p, pe; + product(a.hi, b.hi, p, pe); + const T sum = s + p; + e += two_sum_error(s, p, sum) + (pe + (a.hi * b.lo + a.lo * b.hi)); + s = sum; + } + // One step s + e = (s + e) x + a of compensated Horner (Graillat, Langlois & Louvet 2005): the corrections e follow + // their own recurrence, off the dependency chain of s. + void horner_step(T& s, T& e, const T& x, const value_type& a) const + { + T p, pe; + product(s, x, p, pe); + const T sum = p + a.hi; + e = e * x + (two_sum_error(p, a.hi, sum) + (pe + a.lo)); + s = sum; + } +private: + two_product product; +}; + +#endif // BOOST_MATH_HAS_NVRTC + // // Neumaier's compensated sum: sum + comp carries the total to about one rounding, however many terms are added. // diff --git a/nanobind/src/quadrature.cpp b/nanobind/src/quadrature.cpp index 93f139de0b..958c35dfce 100644 --- a/nanobind/src/quadrature.cpp +++ b/nanobind/src/quadrature.cpp @@ -12,7 +12,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -20,14 +22,18 @@ #include #include +#include #include #include +#include #include #include #include #include +#include #include +#include namespace bmpy { @@ -66,19 +72,94 @@ std::tuple gauss_kronrod_n(scalar_function const& f, dou return {Q, error, L1}; } -template -std::pair gauss_hermite_n(scalar_function const& f) +// (integral, L1 norm) of the Gauss rule, or (integral, error estimate, L1 norm) of the averaged rule. +template +std::pair gauss_rule_integrate(Rule const& rule, scalar_function const& f) { double L1 = 0; - double Q = bq::gauss_hermite::integrate(f, &L1); + double Q = rule.integrate(f, &L1); return {Q, L1}; } +template +std::tuple gauss_rule_integrate_averaged(Rule const& rule, scalar_function const& f) +{ + double error = 0, L1 = 0; + double Q = rule.integrate_averaged(f, &error, &L1); + return {Q, error, L1}; +} + +using array1d = nb::ndarray>; + +// A copy, so that writing to it leaves the rule alone. +array1d to_array(std::vector const& v) +{ + auto* data = new std::vector(v); + nb::capsule owner(data, [](void* p) noexcept { delete static_cast*>(p); }); + return array1d(data->data(), {data->size()}, owner); +} + +// gauss_jacobi keeps its parameters to itself; remember them for repr. +struct py_gauss_jacobi : bq::gauss_jacobi +{ + unsigned n; + double alpha, beta; + py_gauss_jacobi(unsigned n_, double alpha_, double beta_) : bq::gauss_jacobi(n_, alpha_, beta_), n(n_), alpha(alpha_), beta(beta_) {} +}; + +// The n-point Gauss rule of a family, computed afresh: (nodes ascending, weights). A symmetric family keeps only +// its nonnegative nodes, so mirror them. +template +std::pair gauss_rule(unsigned n, Family const& family, bool symmetric = false) +{ + BOOST_MATH_NAMESPACE::detail::orthogonal_polynomial rule(n, family); + std::vector x = rule.abscissa(), w = rule.weights(); + if (symmetric) + { + std::size_t m = x.size(), skip = n % 2; // An odd rule has 0 once. + std::vector xs(n), ws(n); + for (std::size_t i = 0; i < m; ++i) + { + xs[n - m + i] = x[i]; + ws[n - m + i] = w[i]; + if (i >= skip) + { + xs[m - 1 - i] = -x[i]; + ws[m - 1 - i] = w[i]; + } + } + x.swap(xs); + w.swap(ws); + } + return {to_array(x), to_array(w)}; +} + +void check_exponent(char const* name, double a) +{ + if (!(a > -1) || !std::isfinite(a)) + throw nb::value_error((std::string(name) + " must be finite and > -1").c_str()); +} + [[noreturn]] void bad_points(char const* name, char const* allowed, unsigned points) { throw nb::value_error((std::string(name) + ": points must be one of " + allowed + ", got " + std::to_string(points)).c_str()); } +// fn(std::integral_constant()) for the point counts gauss_hermite and gauss_laguerre are bound for. +template +auto with_points(char const* name, unsigned points, Fn fn) +{ + switch (points) + { + case 7: return fn(std::integral_constant()); + case 10: return fn(std::integral_constant()); + case 15: return fn(std::integral_constant()); + case 20: return fn(std::integral_constant()); + case 30: return fn(std::integral_constant()); + default: bad_points(name, "7, 10, 15, 20, 30", points); + } +} + } // namespace void bind_quadrature(nb::module_& m) @@ -146,18 +227,72 @@ void bind_quadrature(nb::module_& m) "quadrature::gauss_kronrod: (integral, error estimate, L1 norm), adaptive"); m.def("gauss_hermite", [](nb::callable f, unsigned points) { - scalar_function g{f}; - switch (points) - { - case 7: return gauss_hermite_n<7>(g); - case 10: return gauss_hermite_n<10>(g); - case 15: return gauss_hermite_n<15>(g); - case 20: return gauss_hermite_n<20>(g); - case 30: return gauss_hermite_n<30>(g); - default: bad_points("gauss_hermite", "7, 10, 15, 20, 30", points); - } + return with_points("gauss_hermite", points, [g = scalar_function{f}](auto n) { + return gauss_rule_integrate(bq::gauss_hermite(), g); + }); }, arg("f"), arg("points") = 10u, "quadrature::gauss_hermite: (integral of exp(-x^2) f(x) over the real line, L1 norm)"); + m.def("gauss_hermite_averaged", [](nb::callable f, unsigned points) { + return with_points("gauss_hermite_averaged", points, [g = scalar_function{f}](auto n) { + return gauss_rule_integrate_averaged(bq::gauss_hermite(), g); + }); + }, arg("f"), arg("points") = 10u, + "quadrature::gauss_hermite::integrate_averaged: (Laurie's averaged-rule integral, error estimate, L1 norm)"); + + m.def("gauss_laguerre", [](nb::callable f, unsigned points, double alpha) { + return with_points("gauss_laguerre", points, [g = scalar_function{f}, alpha](auto n) { + return gauss_rule_integrate(bq::gauss_laguerre(alpha), g); + }); + }, arg("f"), arg("points") = 10u, arg("alpha") = 0.0, + "quadrature::gauss_laguerre(alpha): (integral of x^alpha exp(-x) f(x) over [0, inf), L1 norm)"); + + m.def("gauss_laguerre_averaged", [](nb::callable f, unsigned points, double alpha) { + return with_points("gauss_laguerre_averaged", points, [g = scalar_function{f}, alpha](auto n) { + return gauss_rule_integrate_averaged(bq::gauss_laguerre(alpha), g); + }); + }, arg("f"), arg("points") = 10u, arg("alpha") = 0.0, + "quadrature::gauss_laguerre(alpha).integrate_averaged: (Laurie's averaged-rule integral, error estimate, L1 norm)"); + + nb::class_(m, "gauss_jacobi", + "quadrature::gauss_jacobi(n, alpha, beta): the n-point rule for (1 - x)^alpha (1 + x)^beta on [-1, 1], " + "and Laurie's (n + 1)-point anti-Gauss rule") + .def(nb::init(), arg("n"), arg("alpha"), arg("beta")) + .def("integrate", [](py_gauss_jacobi const& rule, nb::callable f) { + return gauss_rule_integrate(rule, scalar_function{f}); + }, arg("f"), "(integral, L1 norm)") + .def("integrate_averaged", [](py_gauss_jacobi const& rule, nb::callable f) { + return gauss_rule_integrate_averaged(rule, scalar_function{f}); + }, arg("f"), "(Laurie's averaged-rule integral, error estimate, L1 norm); raises when anti_gauss_is_internal is False") + .def_ro("n", &py_gauss_jacobi::n) + .def_ro("alpha", &py_gauss_jacobi::alpha) + .def_ro("beta", &py_gauss_jacobi::beta) + .def_prop_ro("abscissa", [](py_gauss_jacobi const& rule) { return to_array(rule.abscissa()); }, nb::rv_policy::move) + .def_prop_ro("weights", [](py_gauss_jacobi const& rule) { return to_array(rule.weights()); }, nb::rv_policy::move) + .def_prop_ro("anti_gauss_abscissa", [](py_gauss_jacobi const& rule) { return to_array(rule.anti_gauss_abscissa()); }, nb::rv_policy::move) + .def_prop_ro("anti_gauss_weights", [](py_gauss_jacobi const& rule) { return to_array(rule.anti_gauss_weights()); }, nb::rv_policy::move) + .def_prop_ro("anti_gauss_is_internal", [](py_gauss_jacobi const& rule) { return rule.anti_gauss_is_internal(); }, + "False when the anti-Gauss rule would have a node outside [-1, 1]; its nodes and weights are then empty") + .def("__repr__", [](py_gauss_jacobi const& rule) { + return "gauss_jacobi(n=" + std::to_string(rule.n) + ", alpha=" + nb::cast(nb::repr(nb::float_(rule.alpha))) + + ", beta=" + nb::cast(nb::repr(nb::float_(rule.beta))) + ")"; + }); + + m.def("gauss_hermite_rule", [](unsigned n) { + return gauss_rule(n, BOOST_MATH_NAMESPACE::detail::hermite_family(n), true); + }, arg("n"), nb::rv_policy::move, "(nodes, weights) of the n-point rule for exp(-x^2), computed on each call; as numpy's hermgauss"); + + m.def("gauss_laguerre_rule", [](unsigned n, double alpha) { + check_exponent("alpha", alpha); + return gauss_rule(n, BOOST_MATH_NAMESPACE::detail::laguerre_family(n, alpha)); + }, arg("n"), arg("alpha") = 0.0, nb::rv_policy::move, "(nodes, weights) of the n-point rule for x^alpha exp(-x), computed on each call"); + + m.def("gauss_jacobi_rule", [](unsigned n, double alpha, double beta) { + check_exponent("alpha", alpha); + check_exponent("beta", beta); + return gauss_rule(n, BOOST_MATH_NAMESPACE::detail::jacobi_family(n, alpha, beta)); + }, arg("n"), arg("alpha"), arg("beta"), nb::rv_policy::move, + "(nodes, weights) of the n-point rule for (1 - x)^alpha (1 + x)^beta, computed on each call, without the anti-Gauss rule"); + m.def("ooura_fourier_sin", [](nb::callable f, double omega, double tolerance, std::size_t levels) { bq::ooura_fourier_sin integrator(tolerance, levels); return integrator.integrate(scalar_function{f}, omega); diff --git a/nanobind/tests/test_smoke.py b/nanobind/tests/test_smoke.py index c232eb96d1..90e3e85071 100644 --- a/nanobind/tests/test_smoke.py +++ b/nanobind/tests/test_smoke.py @@ -277,6 +277,7 @@ def test_array_matches_scalar(case): # Not elementwise, so tested on their own below. SIEVE = {"prime_sieve", "prime_range", "prime_count", "polynomial_roots", "polynomial_roots_aberth", "polynomial_roots_companion_matrix"} QUADRATURE = {"tanh_sinh", "exp_sinh", "sinh_sinh", "trapezoidal", "gauss", "gauss_kronrod", "gauss_hermite", + "gauss_laguerre", "gauss_hermite_averaged", "gauss_laguerre_averaged", "gauss_hermite_rule", "gauss_laguerre_rule", "gauss_jacobi_rule", "ooura_fourier_sin", "ooura_fourier_cos", "van_den_bos_unit_square"} OPTIMIZATION = {"differential_evolution", "jso", "random_search", "cma_es"} @@ -306,6 +307,9 @@ def bad(x): def test_quadrature(): + def second_moment(rule): + x, w = rule + return (w @ x**2,) cases = [ (bm.tanh_sinh(lambda x: 1 / math.sqrt(x), 0, 1), 2.0), (bm.exp_sinh(lambda x: math.exp(-x)), 1.0), @@ -314,6 +318,14 @@ def test_quadrature(): (bm.gauss(lambda x: x * x, 0, 1, points=7), 1 / 3), (bm.gauss_kronrod(math.sin, 0, math.pi, points=61), 2.0), (bm.gauss_hermite(lambda x: x * x), math.sqrt(math.pi) / 2), + (bm.gauss_laguerre(lambda x: x * x, alpha=0.5), math.gamma(3.5)), + (bm.gauss_jacobi(10, 0.5, -0.5).integrate(lambda x: x * x), math.pi / 2), + (bm.gauss_hermite_averaged(lambda x: x * x), math.sqrt(math.pi) / 2), + (bm.gauss_laguerre_averaged(lambda x: x * x, alpha=0.5), math.gamma(3.5)), + (bm.gauss_jacobi(10, 0.5, -0.5).integrate_averaged(lambda x: x * x), math.pi / 2), + (second_moment(bm.gauss_hermite_rule(5)), math.sqrt(math.pi) / 2), + (second_moment(bm.gauss_laguerre_rule(5, alpha=0.5)), math.gamma(3.5)), + (second_moment(bm.gauss_jacobi_rule(5, 0.5, -0.5)), math.pi / 2), (bm.ooura_fourier_sin(lambda t: 1 / t, 1.0), math.pi / 2), (bm.ooura_fourier_cos(lambda t: 1 / (1 + t * t), 1.0), math.pi / (2 * math.e)), (bm.van_den_bos_unit_square(lambda x, y: x * y), 0.25), @@ -322,6 +334,12 @@ def test_quadrature(): assert math.isclose(result[0], exact, rel_tol=1e-12) with pytest.raises(ValueError): bm.gauss(math.sin, 0, 1, points=8) + # (integral, L1) from the Gauss rule; (integral, error, L1) from the averaged rule, whose error bounds its own. + assert len(bm.gauss_hermite(math.cos)) == 2 + Q, error, L1 = bm.gauss_hermite_averaged(math.cos, points=7) + assert abs(Q - math.sqrt(math.pi) * math.exp(-0.25)) <= error < 1e-6 + with pytest.raises(RuntimeError): + bm.gauss_jacobi(10, -0.9, -0.9).integrate_averaged(math.cos) with pytest.raises(ZeroDivisionError): bm.tanh_sinh(lambda x: 1 / 0, 0, 1) @@ -469,7 +487,7 @@ def test_array_errors_warn_and_give_nan(): def test_every_distribution_has_a_case(): names = {n for n in bm.__all__ if inspect.isclass(getattr(bm, n))} - assert names - TRANSFORMS == {c[0] for c in CONTINUOUS + DISCRETE} + assert names - TRANSFORMS - {"gauss_jacobi"} == {c[0] for c in CONTINUOUS + DISCRETE} assert TRANSFORMS <= names diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 359080af69..cf83a34cb0 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1284,6 +1284,12 @@ test-suite quadrature : [ run test_gauss_hermite_quadrature.cpp : : : $(float128_type) [ requires cxx17_if_constexpr ] off msvc:/bigobj release ] [ run test_gauss_hermite_quadrature_eigen.cpp : : : [ requires cxx17_if_constexpr ] [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] + [ run test_gauss_laguerre_quadrature.cpp : : : $(float128_type) [ requires cxx17_if_constexpr ] off msvc:/bigobj release ] + [ run test_gauss_laguerre_quadrature_eigen.cpp : : : [ requires cxx17_if_constexpr ] + [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] + [ run test_gauss_jacobi_quadrature.cpp : : : $(float128_type) [ requires cxx17_if_constexpr ] off msvc:/bigobj release ] + [ run test_gauss_jacobi_quadrature_eigen.cpp : : : [ requires cxx17_if_constexpr ] + [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] [ run test_naive_monte_carlo.cpp : : : msvc:/bigobj TEST=1 [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_hdr_thread cxx11_hdr_atomic cxx11_decltype cxx11_hdr_future cxx11_hdr_chrono cxx11_hdr_random cxx11_allocator ] diff --git a/test/gauss_quadrature_test_helpers.hpp b/test/gauss_quadrature_test_helpers.hpp new file mode 100644 index 0000000000..ed5a736d2f --- /dev/null +++ b/test/gauss_quadrature_test_helpers.hpp @@ -0,0 +1,90 @@ +// 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) +// +// Checks shared by the Gauss-Hermite, Gauss-Laguerre and Gauss-Jacobi tests. + +#ifndef BOOST_MATH_TEST_GAUSS_QUADRATURE_TEST_HELPERS_HPP +#define BOOST_MATH_TEST_GAUSS_QUADRATURE_TEST_HELPERS_HPP + +#include "math_unit_test.hpp" +#include +#include +#include + +// The integrals of 1, x^2 and cos(x) against the rule's weight function, as one valarray. +template +void check_vector_valued(const Rule& rule, double m0, double m2, double cos_moment, unsigned ulps0, unsigned ulps2) +{ + using vector = std::valarray; + const vector zero(0.0, 3); + auto norm = [](const vector& v) { return std::sqrt((v * v).sum()); }; + auto f = [](double x) { return vector{1.0, x * x, std::cos(x)}; }; + double L1; + vector Q = rule.integrate(f, zero, norm, &L1); + CHECK_ULP_CLOSE(m0, Q[0], ulps0); + CHECK_ULP_CLOSE(m2, Q[1], ulps2); + CHECK_ULP_CLOSE(cos_moment, Q[2], 16); + CHECK_LE(norm(Q), L1); + // The error is the norm of half the difference from the anti-Gauss rule, at least that of any component. + double error, L1_averaged, cos_error; + const vector Q_averaged = rule.integrate_averaged(f, zero, norm, &error, &L1_averaged); + CHECK_ULP_CLOSE(m0, Q_averaged[0], ulps0); + CHECK_ULP_CLOSE(m2, Q_averaged[1], ulps2); + CHECK_ABSOLUTE_ERROR(cos_moment, Q_averaged[2], error); + CHECK_LE(norm(Q_averaged), L1_averaged); + rule.integrate_averaged([](double x) { return std::cos(x); }, &cos_error); + CHECK_LE(cos_error, error); + CHECK_LE(2 * std::numeric_limits::epsilon() * L1_averaged, error); +} + +// For smooth f, G and A bracket the integral, so the error bounds that of the average and is close to that of G +// until it reaches rounding. The error may fall short of that of the average by slack eps L1. +template +void check_error_estimate(const Rule& rule, F f, const Real& exact, const Real& slack = 0) +{ + using std::abs; + Real error, L1; + const Real Q = rule.integrate_averaged(f, &error, &L1); + const Real gauss_error = abs(rule.integrate(f) - exact); + CHECK_LE(abs(Q), L1); + CHECK_LE(abs(Q - exact), error + slack * std::numeric_limits::epsilon() * L1); + CHECK_LE(abs(Q - exact), gauss_error + 2 * std::numeric_limits::epsilon() * L1); + if (gauss_error > 1000 * std::numeric_limits::epsilon() * L1) + { + CHECK_LE(error, 2 * gauss_error); + CHECK_LE(gauss_error, 2 * error); + } + CHECK_LE(2 * std::numeric_limits::epsilon() * L1, error); +} + +// The integral of a 2x2 complex matrix polynomial, with the checks every rule shares; the callers check the moments. +// Return concrete Eigen matrices so the test does not rely on the lifetime +// of expression-template temporaries. Supply a shaped zero and a scalar norm. +template +Matrix integrate_matrix(const Rule& rule, const Matrix& zero) +{ + const auto norm = [](const Matrix& m) { return m.norm(); }; + const auto f = [&zero](double x) -> Matrix { + Matrix result = zero; + result(0, 0) = {1, x * x}; + result(0, 1) = {x, 2 * x}; + result(1, 0) = {-x * x, x}; + result(1, 1) = {x * x * x, -1}; + return result; + }; + double L1; + Matrix Q = rule.integrate(f, zero, norm, &L1); + // The integrand is a polynomial that both rules integrate exactly, so the estimate is rounding. + double error; + rule.integrate_averaged(f, zero, norm, &error); + CHECK_LE(error, 64 * std::numeric_limits::epsilon() * L1); + CHECK_EQUAL(Q.rows(), decltype(Q.rows())(2)); + CHECK_EQUAL(Q.cols(), decltype(Q.cols())(2)); + CHECK_LE(norm(Q), L1); + return Q; +} + +#endif // BOOST_MATH_TEST_GAUSS_QUADRATURE_TEST_HELPERS_HPP diff --git a/test/test_gauss_hermite_quadrature.cpp b/test/test_gauss_hermite_quadrature.cpp index be46ce9318..ac347c47c3 100644 --- a/test/test_gauss_hermite_quadrature.cpp +++ b/test/test_gauss_hermite_quadrature.cpp @@ -6,10 +6,10 @@ // or copy at http://www.boost.org/LICENSE_1_0.txt) #include "math_unit_test.hpp" +#include "gauss_quadrature_test_helpers.hpp" #include #include #include -#include #include #include #include @@ -79,17 +79,8 @@ void test_complex(std::size_t ulps, typename Complex::value_type imaginary_tol) template void test_vector_valued() { - using vector = std::valarray; - const vector zero(0.0, 3); - auto norm = [](const vector& v) { return std::sqrt((v * v).sum()); }; - auto f = [](double x) { return vector{1.0, x * x, std::cos(x)}; }; - double L1; - vector Q = gauss_hermite::integrate(f, zero, norm, &L1); const double root_pi = boost::math::constants::root_pi(); - CHECK_ULP_CLOSE(root_pi, Q[0], 4); - CHECK_ULP_CLOSE(root_pi / 2, Q[1], 4); - CHECK_ULP_CLOSE(root_pi * std::exp(-0.25), Q[2], 16); - CHECK_LE(norm(Q), L1); + check_vector_valued(gauss_hermite(), root_pi, root_pi / 2, root_pi * std::exp(-0.25), 4, 4); } // Evaluating H_N directly overflows long before these N; the orthonormal recurrence does not. @@ -105,22 +96,107 @@ void test_large_N() CHECK_ULP_CLOSE(root_pi * std::exp(-0.25), rule::integrate([](double x) { return std::cos(x); }), 16); } -// Past N = 706 in double the recurrence overflows, which must be reported through the policy. +// Past N = 706 the recurrence overflows in double, but the differential equation does not: neither rule overflows. void test_overflow() { - using ignore_policy = boost::math::policies::policy >; - double L1 = 0; - CHECK_NAN((gauss_hermite::integrate([](double x) { return std::cos(x); }, &L1))); - CHECK_NAN(L1); - using vector = std::valarray; - vector Q = gauss_hermite::integrate([](double x) { return vector{1.0, x}; }, vector(0.0, 2), - [](const vector& v) { return std::sqrt((v * v).sum()); }); - CHECK_EQUAL(Q.size(), std::size_t(2)); - CHECK_NAN(Q[0]); - CHECK_NAN(Q[1]); -#ifndef BOOST_NO_EXCEPTIONS - CHECK_THROW((gauss_hermite::integrate([](double x) { return std::cos(x); })), boost::math::evaluation_error); -#endif + const double root_pi = boost::math::constants::root_pi(); + double L1 = 0, error = 0; + const double Q = gauss_hermite::integrate([](double x) { return std::cos(x); }, &L1); + CHECK_ULP_CLOSE(root_pi * std::exp(-0.25), Q, 16); + CHECK_LE(Q, L1); + CHECK_LE(L1, root_pi); + CHECK_ULP_CLOSE(root_pi * std::exp(-0.25), (gauss_hermite::integrate_averaged([](double x) { return std::cos(x); }, &error, &L1)), 16); + CHECK_LE(error, 64 * std::numeric_limits::epsilon()); +} + +// Laurie's anti-Gauss rule with N + 1 points, which the error estimate uses, as a whole rule rather than its nonnegative half. +template +std::pair, std::vector > anti_gauss() +{ + boost::math::detail::orthogonal_polynomial > rule(N + 1, boost::math::detail::hermite_family(N + 1), true); + std::vector x, w; + for (std::size_t i = rule.abscissa().size(); i-- > 0;) + { + if (rule.abscissa()[i] != 0) + { + x.push_back(-rule.abscissa()[i]); + w.push_back(rule.weights()[i]); + } + } + x.insert(x.end(), rule.abscissa().begin(), rule.abscissa().end()); + w.insert(w.end(), rule.weights().begin(), rule.weights().end()); + return std::make_pair(x, w); +} + +// The anti-Gauss rule is the Gauss rule of 2I - G_N, so on polynomials of degree 2N + 1 its error is that of G_N +// with the opposite sign; it interlaces G_N, and its weights are positive and sum to sqrt(pi). +template +void test_anti_gauss() +{ + using std::abs; + const auto anti = anti_gauss(); + CHECK_EQUAL(anti.first.size(), static_cast(N + 1)); + Real sum = 0; + for (std::size_t i = 0; i < anti.first.size(); ++i) + { + CHECK_LE(Real(0), anti.second[i]); + sum += anti.second[i]; + } + CHECK_ULP_CLOSE(boost::math::constants::root_pi(), sum, 8); + // xi_1 < x_1 < xi_2 < ... < x_N < xi_{N+1}, with x the nonnegative Gauss abscissas mirrored. + const auto& half = gauss_hermite::abscissa(); + for (std::size_t i = 0; i < half.size(); ++i) + { + const Real x = static_cast(half[i]); + const std::size_t j = N / 2 + i; + CHECK_LE(anti.first[j], x); + CHECK_LE(x, anti.first[j + 1]); + CHECK_LE(anti.first[N - j - 1], -x); + CHECK_LE(-x, anti.first[N - j]); + } + // p(x) = sum c_k x^k with arbitrary c_k, k <= 2N + 1, whose integral is sum c_{2k} Gamma(k + 1/2). + std::vector c(2 * N + 2); + Real exact = 0; + Real moment = boost::math::constants::root_pi(); + for (unsigned k = 0; k < c.size(); ++k) + { + c[k] = Real((k * 7919) % 17) / 8 - 1; + if (k % 2 == 0) + { + exact += c[k] * moment; + moment *= Real(k + 1) / 2; + } + } + auto p = [&c](const Real& x) { Real r = 0; for (std::size_t k = c.size(); k-- > 0;) r = r * x + c[k]; return r; }; + Real L1; + const Real gauss = gauss_hermite::integrate(p, &L1); + Real anti_result = 0; + for (std::size_t i = 0; i < anti.first.size(); ++i) + { + anti_result += anti.second[i] * p(anti.first[i]); + L1 += anti.second[i] * abs(p(anti.first[i])); + } + CHECK_ABSOLUTE_ERROR(Real(0), (exact - gauss) + (exact - anti_result), (N + 8) * std::numeric_limits::epsilon() * L1); +} + +template +void test_error_estimate() +{ + using std::abs; + using std::cos; + using std::exp; + using rule = gauss_hermite; + check_error_estimate(rule(), [](const Real& x) { return cos(x); }, boost::math::constants::root_pi() * exp(-Real(1) / 4)); + // 1/(1 + x^2) has poles at +-i, so the error falls slowly. + auto g = [](const Real& x) { return 1 / (1 + x * x); }; + const Real exact_g = static_cast(1.343293421646735170437123594410589778323L); + Real error; + const Real Qg = rule::integrate_averaged(g, &error); + const Real gauss_error_g = abs(rule::integrate(g) - exact_g); + CHECK_LE(abs(Qg - exact_g), error); + CHECK_LE(abs(Qg - exact_g), gauss_error_g); + CHECK_LE(error, 2 * gauss_error_g); + CHECK_LE(gauss_error_g, 2 * error); } template @@ -135,8 +211,24 @@ void test_moments_all_N() test_moments(); } +// The averaged rule is exact on x^(2N), which G is not, so its error is that of G. +template +void test_averaged() +{ + auto f = [](double x) { return std::pow(x, 2 * N); }; + double error = 0, L1 = 0; + using rule = gauss_hermite; + const double averaged = rule::integrate_averaged(f, &error, &L1); + CHECK_ULP_CLOSE(boost::math::tgamma(N + 0.5), averaged, 64); + // The integrand is nonnegative, so the L1 norm is the integral itself. + CHECK_ULP_CLOSE(averaged, L1, 64); + CHECK_ABSOLUTE_ERROR(std::abs(boost::math::tgamma(N + 0.5) - rule::integrate(f)), error, 64 * std::numeric_limits::epsilon() * L1); +} + int main() { + test_averaged<5>(); + test_averaged<10>(); test_moments_all_N(); test_moments_all_N(); #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS @@ -173,5 +265,23 @@ int main() test_large_N<706>(); test_overflow(); + test_anti_gauss(); + test_anti_gauss(); + test_anti_gauss(); + test_anti_gauss(); + test_anti_gauss(); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_anti_gauss(); + test_anti_gauss(); +#endif + test_error_estimate(); + test_error_estimate(); + test_error_estimate(); + test_error_estimate(); + test_error_estimate(); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_error_estimate(); +#endif + return boost::math::test::report_errors(); } diff --git a/test/test_gauss_hermite_quadrature_eigen.cpp b/test/test_gauss_hermite_quadrature_eigen.cpp index 632390c86a..90dfa144dd 100644 --- a/test/test_gauss_hermite_quadrature_eigen.cpp +++ b/test/test_gauss_hermite_quadrature_eigen.cpp @@ -5,6 +5,7 @@ // or copy at http://www.boost.org/LICENSE_1_0.txt) #include "math_unit_test.hpp" +#include "gauss_quadrature_test_helpers.hpp" #include #include #include @@ -13,26 +14,12 @@ using boost::math::quadrature::gauss_hermite; -// Return concrete Eigen matrices so the test does not rely on the lifetime -// of expression-template temporaries. Supply a shaped zero and a scalar norm. template void test_matrix(const Matrix& zero) { - const auto norm = [](const Matrix& m) { return m.norm(); }; - const auto f = [&zero](double x) -> Matrix { - Matrix result = zero; - result(0, 0) = {1, x * x}; - result(0, 1) = {x, 2 * x}; - result(1, 0) = {-x * x, x}; - result(1, 1) = {x * x * x, -1}; - return result; - }; // Against exp(-x^2): 1 -> sqrt(pi), x^2 -> sqrt(pi)/2, odd powers -> 0. const double root_pi = boost::math::constants::root_pi(); - double L1; - Matrix Q = gauss_hermite::integrate(f, zero, norm, &L1); - CHECK_EQUAL(Q.rows(), Eigen::Index(2)); - CHECK_EQUAL(Q.cols(), Eigen::Index(2)); + const Matrix Q = integrate_matrix(gauss_hermite(), zero); CHECK_ULP_CLOSE(root_pi, Q(0, 0).real(), 4); CHECK_ULP_CLOSE(root_pi / 2, Q(0, 0).imag(), 4); CHECK_EQUAL(Q(0, 1), std::complex(0, 0)); @@ -40,7 +27,6 @@ void test_matrix(const Matrix& zero) CHECK_EQUAL(Q(1, 0).imag(), 0.0); CHECK_EQUAL(Q(1, 1).real(), 0.0); CHECK_ULP_CLOSE(-root_pi, Q(1, 1).imag(), 4); - CHECK_LE(norm(Q), L1); } int main() diff --git a/test/test_gauss_jacobi_quadrature.cpp b/test/test_gauss_jacobi_quadrature.cpp new file mode 100644 index 0000000000..6bc4499833 --- /dev/null +++ b/test/test_gauss_jacobi_quadrature.cpp @@ -0,0 +1,293 @@ +// 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) + +#include "math_unit_test.hpp" +#include "gauss_quadrature_test_helpers.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_HAS_FLOAT128 +#include +#endif + +using boost::math::quadrature::gauss_jacobi; + +// Gauss-Jacobi integration is exact for polynomials of degree less than 2N - 1 +// We use a polynomial of degree 10, so should be exact for N >= 6 +template +void test_moments(Real tol, const char* type_name) +{ + using std::pow; + std::cout << "Testing Moments for N = " << N << " and type " << type_name << std::endl; + auto f = [](const Real& x) { return pow(x, Real(10.0)); }; + Real L1; + gauss_jacobi integrator = gauss_jacobi(N, -Real(0.5), -Real(0.25)); + Real Q = integrator.integrate(f, &L1); + + CHECK_ABSOLUTE_ERROR(Q, Real(0.604206393406616872500520394161000920557127865245), tol); + // The integrand is nonnegative, so the L1 norm is the integral itself. + CHECK_ABSOLUTE_ERROR(Q, L1, tol); +} + +template +void test_complex(typename Complex::value_type tol) +{ + using Real = typename Complex::value_type; + using std::cos; + using std::exp; + using std::sin; + auto f = [](const Real& x) { return Complex(cos(x), sin(x)); }; + gauss_jacobi integrator = gauss_jacobi(N, 1.0, 1.0); + Complex Q = integrator.integrate(f); + Real real_expected = Real(4) * sin(Real(1)) - Real(4) * cos(Real(1)); + CHECK_ABSOLUTE_ERROR(real_expected, Q.real(), tol); + CHECK_ABSOLUTE_ERROR(Real(0.0), Q.imag(), tol); +} + +template +void test_vector_valued() +{ + check_vector_valued(gauss_jacobi(N, 1.0, 1.0), 4.0 / 3.0, 4.0 / 15.0, 4.0 * sin(1.0) - 4.0 * cos(1.0), 10, 20); +} + +// The anti-Gauss rule is the Gauss rule of 2I - G_N, so on polynomials of degree 2N + 1 its error is that of G_N +// with the opposite sign; its N + 1 abscissas interlace those of G_N, and its weights are positive and sum to +// 2^(alpha + beta + 1) B(alpha + 1, beta + 1). +template +void test_anti_gauss(unsigned N, Real alpha, Real beta) +{ + using std::abs; + using std::pow; + const gauss_jacobi rule(N, alpha, beta); + CHECK_EQUAL(true, rule.anti_gauss_is_internal()); + const std::vector& x = rule.anti_gauss_abscissa(); + const std::vector& w = rule.anti_gauss_weights(); + CHECK_EQUAL(x.size(), static_cast(N + 1)); + CHECK_LE(Real(-1), x[0]); + CHECK_LE(x[N], Real(1)); + Real sum = 0; + for (std::size_t i = 0; i < x.size(); ++i) + { + CHECK_LE(Real(0), w[i]); + sum += w[i]; + } + CHECK_ULP_CLOSE(pow(Real(2), alpha + beta + 1) * boost::math::beta(alpha + 1, beta + 1), sum, 16); + for (std::size_t i = 0; i < N; ++i) + { + CHECK_LE(x[i], rule.abscissa()[i]); + CHECK_LE(rule.abscissa()[i], x[i + 1]); + } + // p(x) = sum c_k (1 + x)^k with arbitrary c_k, k <= 2N + 1, whose integral is + // sum c_k 2^(alpha + beta + k + 1) B(alpha + 1, beta + k + 1). + std::vector c(2 * N + 2); + Real exact = 0; + for (unsigned k = 0; k < c.size(); ++k) + { + c[k] = Real((k * 7919) % 17) / 8 - 1; + exact += c[k] * pow(Real(2), alpha + beta + k + 1) * boost::math::beta(alpha + 1, beta + k + 1); + } + auto p = [&c](const Real& t) { Real r = 0; for (std::size_t k = c.size(); k-- > 0;) r = r * (1 + t) + c[k]; return r; }; + Real L1; + const Real gauss = rule.integrate(p, &L1); + Real anti_result = 0; + for (std::size_t i = 0; i < x.size(); ++i) + { + anti_result += w[i] * p(x[i]); + L1 += w[i] * abs(p(x[i])); + } + CHECK_ABSOLUTE_ERROR(Real(0), (exact - gauss) + (exact - anti_result), (N + 8) * std::numeric_limits::epsilon() * L1); + // So the averaged rule is exact on p, and integrate_averaged is (G + A) / 2. + const Real averaged = rule.integrate_averaged(p); + CHECK_ABSOLUTE_ERROR(exact, averaged, (N + 8) * std::numeric_limits::epsilon() * L1); + CHECK_ABSOLUTE_ERROR((gauss + anti_result) / 2, averaged, 4 * std::numeric_limits::epsilon() * L1); +} + +// Laurie 1996, Theorem 3: the anti-Gauss rule is internal iff both of these are nonnegative. Zero puts a node on an +// endpoint, where f is still defined. +template +void test_anti_gauss_internality() +{ + using std::abs; + using family = boost::math::detail::jacobi_family; + using polynomial = boost::math::detail::orthogonal_polynomial; + const Real parameters[] = { -0.75, -0.55, -0.5, -0.25, 0, 0.5, 3 }; + for (Real a : parameters) + { + for (Real b : parameters) + { + for (unsigned n = 1; n <= 40; ++n) + { + const Real right = (2 * a + 1) * n * n + (2 * a + 1) * (a + b + 1) * n + (a + 1) * (a + b) * (a + b + 1) / 2; + const Real left = (2 * b + 1) * n * n + (2 * b + 1) * (a + b + 1) * n + (b + 1) * (a + b) * (a + b + 1) / 2; + const family f(n + 1, a, b); + CHECK_EQUAL(polynomial::anti_gauss_is_internal(n + 1, f, Real(1)), right >= 0); + CHECK_EQUAL(polynomial::anti_gauss_is_internal(n + 1, f, Real(-1)), left >= 0); + } + } + } + // Chebyshev weights of the first kind put anti-Gauss nodes exactly on +-1. + auto internal = [](unsigned n, Real a, Real b) + { + const family f(n + 1, a, b); + return polynomial::anti_gauss_is_internal(n + 1, f, Real(1)) && polynomial::anti_gauss_is_internal(n + 1, f, Real(-1)); + }; + for (unsigned n = 1; n <= 200; n += 7) + { + CHECK_EQUAL(internal(n, Real(-0.5), Real(-0.5)), true); + CHECK_EQUAL(internal(n, Real(0), Real(0)), true); + CHECK_EQUAL(internal(n, Real(0.5), Real(0.5)), true); + CHECK_EQUAL(internal(n, Real(-0.75), Real(0)), false); + } + // Chebyshev weights: the anti-Gauss rule is Gauss-Lobatto, with nodes cos(k pi / (N + 1)), so f is evaluated at +-1. + // The integral of cos(x) / sqrt(1 - x^2) is pi J0(1). + for (unsigned N : { 4u, 6u, 10u, 40u }) + { + const gauss_jacobi chebyshev(N, Real(-0.5), Real(-0.5)); + auto f = [](const Real& x) { using std::cos; return cos(x); }; + const Real exact = boost::math::constants::pi() * boost::math::cyl_bessel_j(0, Real(1)); + Real error, L1; + const Real Q = chebyshev.integrate_averaged(f, &error, &L1); + const Real gauss_error = abs(chebyshev.integrate(f) - exact); + CHECK_LE(abs(Q - exact), error); + if (gauss_error > 1000 * std::numeric_limits::epsilon() * L1) + { + CHECK_LE(error, 2 * gauss_error); + CHECK_LE(gauss_error, 2 * error); + } + } + // With a node outside [-1, 1] there is no averaged rule, but integrate is unaffected and raises nothing. + // The integral of (1 - x)^(-3/4) is 4 2^(1/4). + const gauss_jacobi outside(10, Real(-0.75), Real(0)); + CHECK_ULP_CLOSE(4 * pow(Real(2), Real(0.25)), outside.integrate([](const Real&) { return Real(1); }), 16); +} + +template +void test_error_estimate(unsigned N) +{ + using std::cos; + using std::exp; + struct example + { + Real alpha; + Real beta; + Real (*f)(const Real&); + Real expected; + }; + const example examples[] = { + { Real(0.5), Real(0.5), [](const Real& x) { return cos(x); }, static_cast(1.382459687384168525766283327693941037437L) }, + { Real(1), Real(1), [](const Real& x) { return cos(x); }, static_cast(1.204674715759027157006262856749289583561L) }, + { Real(0.3), Real(-0.4), [](const Real& x) { return exp(x); }, static_cast(2.113791650902447106944400891924615994062L) }, + { Real(0), Real(0), [](const Real& x) { return 1 / (2 - x); }, static_cast(1.098612288668109691395245236922525704647L) }, + }; + // Once converged, rounding in the nodes and weights can exceed the 2 eps L1 floor of the error. + for (const example& e : examples) + check_error_estimate(gauss_jacobi(N, e.alpha, e.beta), e.f, e.expected, Real(4)); +} + +// No higher precision because it takes a long time because it has to +// compute the nodes/weights on the fly! +// The root finder's derivative of P_n from P_n and P_{n-1}, against jacobi_derivative. +void test_family_derivative() +{ + for (unsigned n : {1u, 2u, 7u, 30u}) { + for (double alpha : {-0.5, 0.0, 0.5, 2.0}) { + for (double beta : {-0.5, 0.0, 1.0}) { + const boost::math::detail::jacobi_family family(n, alpha, beta); + for (double x : {-0.9, -0.3, 0.0, 0.4, 0.95}) { + const double p = boost::math::jacobi(n, alpha, beta, x); + const double p_previous = boost::math::jacobi(n - 1, alpha, beta, x); + // The two terms of the identity cancel, so the tolerance scales with their size, not the result's. + const double terms = (std::abs(n * (alpha - beta - (2 * n + alpha + beta) * x) * p) + + std::abs(2 * (n + alpha) * (n + beta) * p_previous)) / ((2 * n + alpha + beta) * (1 - x * x)); + CHECK_ABSOLUTE_ERROR(boost::math::jacobi_derivative(n, alpha, beta, x, 1), family.derivative(x, p, p_previous), + 64 * std::numeric_limits::epsilon() * (1 + terms)); + } + } + } + } +} + +void test_parameter_domain() +{ + using ignore_policy = boost::math::policies::policy >; + for (double bad : { -1.0, -2.0, std::numeric_limits::quiet_NaN(), std::numeric_limits::infinity() }) + { + double error = 0; + CHECK_NAN((gauss_jacobi(5, bad, 0).integrate([](double x) { return x; }))); + CHECK_NAN((gauss_jacobi(5, bad, 0).integrate_averaged([](double x) { return x; }, &error))); + CHECK_NAN(error); + CHECK_NAN((gauss_jacobi(5, 0, bad).integrate([](double x) { return x; }))); +#ifndef BOOST_NO_EXCEPTIONS + CHECK_THROW((gauss_jacobi(5, bad, 0)), std::domain_error); + CHECK_THROW((gauss_jacobi(5, 0, bad)), std::domain_error); +#endif + } + CHECK_NAN((gauss_jacobi(0, 0.5, 0.5).integrate_averaged([](double x) { return x; }))); +#ifndef BOOST_NO_EXCEPTIONS + CHECK_THROW((gauss_jacobi(0, 0.5, 0.5)), std::domain_error); +#endif + // No averaged rule without an internal anti-Gauss rule. + using ignore_evaluation = boost::math::policies::policy >; + const gauss_jacobi external(10, -0.9, -0.9); + CHECK_EQUAL(false, external.anti_gauss_is_internal()); + CHECK_NAN(external.integrate_averaged([](double x) { return x; })); +#ifndef BOOST_NO_EXCEPTIONS + CHECK_THROW((gauss_jacobi(10, -0.9, -0.9).integrate_averaged([](double x) { return x; })), boost::math::evaluation_error); +#endif +} + +int main() +{ + test_moments(8 * std::numeric_limits::epsilon(), "float"); + test_moments(5 * std::numeric_limits::epsilon(), "double"); + test_moments(5 * std::numeric_limits::epsilon(), "double"); + +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_moments(4 * std::numeric_limits::epsilon(), "long double"); +#endif +#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS + test_moments(4e-16, "real concept"); +#endif + + + test_complex, 30>(4 * std::numeric_limits::epsilon()); + + test_vector_valued<30>(); + test_vector_valued<40>(); + + for (unsigned n : { 1u, 2u, 5u, 10u, 15u }) + { + test_anti_gauss(n, 0, 0); + test_anti_gauss(n, 0.5, 0.5); + test_anti_gauss(n, 0.3, -0.4); + test_anti_gauss(n, 2.5, 1); + } +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_anti_gauss(9, 0.25, -0.25); +#endif + test_anti_gauss_internality(); + test_parameter_domain(); + test_family_derivative(); + for (unsigned n : { 2u, 3u, 5u, 8u, 12u, 20u }) + { + test_error_estimate(n); + } +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_error_estimate(6); +#endif + + return boost::math::test::report_errors(); +} diff --git a/test/test_gauss_jacobi_quadrature_eigen.cpp b/test/test_gauss_jacobi_quadrature_eigen.cpp new file mode 100644 index 0000000000..6d0eeffe2c --- /dev/null +++ b/test/test_gauss_jacobi_quadrature_eigen.cpp @@ -0,0 +1,38 @@ +// 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) + +#include "math_unit_test.hpp" +#include "gauss_quadrature_test_helpers.hpp" +#include +#include +#include +#include + +using boost::math::quadrature::gauss_jacobi; + +template +void test_matrix(const Matrix& zero) +{ + // Against 1 -> 16/15, x -> 0, x^2 -> 16/105, odd powers are 0. + const Matrix Q = integrate_matrix(gauss_jacobi(N, 2.0, 2.0), zero); + CHECK_ULP_CLOSE(16.0 / 15.0, Q(0, 0).real(), 4); + CHECK_ULP_CLOSE(16.0 / 105.0, Q(0, 0).imag(), 4); + CHECK_ABSOLUTE_ERROR(Q(0, 1).imag(), 0.0, std::numeric_limits::epsilon()); + CHECK_ABSOLUTE_ERROR(Q(0, 1).real(), 0.0, std::numeric_limits::epsilon()); + CHECK_ULP_CLOSE(-16.0 / 105.0, Q(1, 0).real(), 4); + CHECK_ABSOLUTE_ERROR(Q(1, 0).imag(), 0.0, std::numeric_limits::epsilon()); + CHECK_ABSOLUTE_ERROR(Q(1, 1).real(), 0.0, std::numeric_limits::epsilon()); + CHECK_ULP_CLOSE(-16.0 / 15.0, Q(1, 1).imag(), 4); +} + +int main() +{ + using fixed = Eigen::Matrix, 2, 2>; + test_matrix(fixed::Zero()); + test_matrix(fixed::Zero()); + test_matrix(Eigen::MatrixXcd::Zero(2, 2)); + return boost::math::test::report_errors(); +} diff --git a/test/test_gauss_laguerre_quadrature.cpp b/test/test_gauss_laguerre_quadrature.cpp new file mode 100644 index 0000000000..4c80b39068 --- /dev/null +++ b/test/test_gauss_laguerre_quadrature.cpp @@ -0,0 +1,329 @@ +// 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) + +#include "math_unit_test.hpp" +#include "gauss_quadrature_test_helpers.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef BOOST_HAS_FLOAT128 +#include +#endif + +using boost::math::quadrature::gauss_laguerre; +using boost::math::quadrature::detail::gauss_laguerre_detail; + +// An N-point rule integrates x^(n) exp(-x) exactly for n <= 2N - 1, giving Gamma(n + 1). +// The tolerance grows with n because x^(n) magnifies the rounding of the largest abscissa +// by a factor of n. +template +void test_moments() +{ + // Gamma(k + 1) by its recurrence + Real expected; + for (unsigned k = 0; k < 2 * N - 1; ++k) + { + std::cout << "Testing N = " << N << ", k = " << k << std::endl; + expected = (k == 0) ? Real(1) : boost::math::factorial(k); + auto even = [k](const Real& x) { Real r = 1; for (unsigned j = 0; j < k; ++j) r *= x; return r; }; + Real L1; + Real Q = gauss_laguerre().integrate(even, &L1); + std::cout << std::setprecision(17) << expected << " " << Q << " " << L1 << std::endl; + CHECK_ULP_CLOSE(expected, Q, 16 * (k + 1)); + // The integrand is nonnegative, so the L1 norm is the integral itself. + CHECK_ULP_CLOSE(Q, L1, 16 * (k + 1)); + } +} + +// The tabulated abscissas and weights must agree with those computed on demand. +template +void test_tables_match_on_demand() +{ + using table = gauss_laguerre_detail; + using on_demand = gauss_laguerre_detail; + CHECK_EQUAL(table::abscissa().size(), on_demand::abscissa().size()); + for (unsigned i = 0; i < table::abscissa().size(); ++i) + { + CHECK_ULP_CLOSE(static_cast(table::abscissa()[i]), on_demand::abscissa()[i], 2); + // cpp_bin_float's exp, which the weights computed on demand take at the nodes, is out by 4 ulps at x = 19. + CHECK_ULP_CLOSE(static_cast(table::weights()[i]), on_demand::weights()[i], 16); + } +} + +// Integral of exp(ix) exp(-x) over [0, infinity) is (0.5, 0.5). +template +void test_complex(typename Complex::value_type tol) +{ + using Real = typename Complex::value_type; + using std::cos; + using std::exp; + using std::sin; + auto f = [](const Real& x) { return Complex(cos(x), sin(x)); }; + Complex Q = gauss_laguerre().integrate(f); + CHECK_ABSOLUTE_ERROR(Real(0.5), Q.real(), tol); + CHECK_ABSOLUTE_ERROR(Real(0.5), Q.imag(), tol); +} + +template +void test_vector_valued() +{ + check_vector_valued(gauss_laguerre(), 1.0, 2.0, 0.5, 10, 20); +} + +// Evaluating L_N directly overflows long before these N; the orthonormal recurrence does not. +template +void test_large_N() +{ + using rule = gauss_laguerre; + CHECK_EQUAL(rule().abscissa().size(), static_cast(N)); + for (unsigned i = 1; i < rule().abscissa().size(); ++i) + CHECK_LE(rule().abscissa()[i - 1], rule().abscissa()[i]); + CHECK_ULP_CLOSE(1.0, rule().integrate([](double) { return 1.0; }), 52); + CHECK_ULP_CLOSE(0.5, rule().integrate([](double x) { return std::cos(x); }), 52); +} + +// Past N = 326 the recurrence overflows in double. The Gauss rule comes from the differential equation, which does not, +// Past N = 326 the recurrence overflows in double, but the differential equation does not: neither rule overflows. +void test_overflow() +{ + double L1 = 0, error = 0; + CHECK_ULP_CLOSE(0.5, (gauss_laguerre().integrate([](double x) { return std::cos(x); }, &L1)), 64); + CHECK_ULP_CLOSE(0.5, (gauss_laguerre().integrate_averaged([](double x) { return std::cos(x); }, &error, &L1)), 64); + CHECK_LE(error, 64 * std::numeric_limits::epsilon()); +} + +// The anti-Gauss rule is the Gauss rule of 2I - G_N, so on polynomials of degree 2N + 1 its error is that of G_N +// with the opposite sign. Its N + 1 abscissas interlace those of G_N and are positive (Laurie 1996, Theorem 2), +// and its weights are positive and sum to Gamma(alpha + 1). +template +void test_anti_gauss(unsigned N, Real alpha = 0) +{ + using std::abs; + using family = boost::math::detail::laguerre_family; + const boost::math::detail::orthogonal_polynomial gauss(N, family(N, alpha)); + const boost::math::detail::orthogonal_polynomial anti(N + 1, family(N + 1, alpha), true); + const std::vector& x = anti.abscissa(); + const std::vector& w = anti.weights(); + CHECK_EQUAL(x.size(), static_cast(N + 1)); + CHECK_LE(std::numeric_limits::min(), x[0]); + Real sum = 0; + for (std::size_t i = 0; i < x.size(); ++i) + { + CHECK_LE(Real(0), w[i]); + sum += w[i]; + } + CHECK_ULP_CLOSE(boost::math::tgamma(alpha + 1), sum, 8); + for (std::size_t i = 0; i < N; ++i) + { + CHECK_LE(x[i], gauss.abscissa()[i]); + CHECK_LE(gauss.abscissa()[i], x[i + 1]); + } + // p = sum c_k L_k^alpha with arbitrary c_k, k <= 2N + 1, whose integral is c_0 Gamma(alpha + 1). The error of G_N + // on L_2N is 1 / binomial(2N, N), so this only constrains the rule for small N; rounding grows with N. + std::vector c(2 * N + 2); + for (unsigned k = 0; k < c.size(); ++k) + { + c[k] = Real((k * 7919) % 17) / 8 - 1; + } + const Real exact = c[0] * boost::math::tgamma(alpha + 1); + auto p = [&c, alpha](const Real& t) + { + Real r = c[0], L = 1, L_previous = 0; + for (unsigned k = 1; k < c.size(); ++k) + { + Real L_next = ((2 * k - 1 + alpha - t) * L - (k - 1 + alpha) * L_previous) / k; + L_previous = L; + L = L_next; + r += c[k] * L; + } + return r; + }; + Real gauss_result = 0, anti_result = 0, L1 = 0; + for (std::size_t i = 0; i < N; ++i) + { + gauss_result += gauss.weights()[i] * p(gauss.abscissa()[i]); + L1 += gauss.weights()[i] * abs(p(gauss.abscissa()[i])); + } + for (std::size_t i = 0; i < x.size(); ++i) + { + anti_result += w[i] * p(x[i]); + L1 += w[i] * abs(p(x[i])); + } + CHECK_ABSOLUTE_ERROR(Real(0), (exact - gauss_result) + (exact - anti_result), (N + 8) * std::numeric_limits::epsilon() * L1); +} + +template +void test_error_estimate() +{ + using std::cos; + using std::sqrt; + check_error_estimate(gauss_laguerre(), [](const Real& x) { return cos(x); }, Real(0.5)); + // The integral of exp(-x) / (1 + x) is the Gompertz constant; the pole at -1 makes the error fall slowly. + check_error_estimate(gauss_laguerre(), [](const Real& x) { return 1 / (1 + x); }, static_cast(0.5963473623231940743410784993692793760742L)); + check_error_estimate(gauss_laguerre(), [](const Real& x) { return sqrt(1 + x); }, static_cast(1.378936078070656053021675619957108958174L)); +} + +// The rule for x^alpha exp(-x) integrates x^k exactly for k <= 2N - 1, giving Gamma(k + alpha + 1), and its error +// estimate tracks the true error for the integral of x^alpha exp(-x) cos(x), Gamma(alpha + 1) cos((alpha + 1) pi / 4) / 2^((alpha + 1) / 2). +template +void test_alpha(Real alpha, Real moment_tol) +{ + using std::abs; + using std::cos; + using std::pow; + const gauss_laguerre rule(alpha); + for (unsigned k = 0; k < 2 * N; ++k) + { + Real L1; + const Real Q = rule.integrate([k](const Real& x) { return pow(x, Real(k)); }, &L1); + CHECK_ULP_CLOSE(boost::math::tgamma(k + alpha + 1), Q, static_cast(moment_tol * (k + 1))); + } + Real error, L1; + const Real Q = rule.integrate_averaged([](const Real& x) { return cos(x); }, &error, &L1); + const Real expected = boost::math::tgamma(alpha + 1) * cos((alpha + 1) * boost::math::constants::pi() / 4) / pow(Real(2), (alpha + 1) / 2); + const Real gauss_error = abs(rule.integrate([](const Real& x) { return cos(x); }) - expected); + // For fewer points the estimate is preasymptotic: 8 times the error for N = 2 and alpha = 2.5, and G and A + // need not bracket the integral. + if ((N >= 12) && (gauss_error > 1000 * std::numeric_limits::epsilon() * L1)) + { + CHECK_LE(abs(Q - expected), error); + CHECK_LE(error, 2 * gauss_error); + CHECK_LE(gauss_error, 2 * error); + } + CHECK_LE(2 * std::numeric_limits::epsilon() * L1, error); +} + +template +void test_anti_gauss_accessors(double alpha) +{ + using family = boost::math::detail::laguerre_family; + const gauss_laguerre rule(alpha); + const boost::math::detail::orthogonal_polynomial anti(N + 1, family(N + 1, alpha), true); + CHECK_EQUAL(rule.anti_gauss_abscissa().size(), static_cast(N + 1)); + for (unsigned i = 0; i <= N; ++i) + { + CHECK_EQUAL(anti.abscissa()[i], rule.anti_gauss_abscissa()[i]); + CHECK_EQUAL(anti.weights()[i], rule.anti_gauss_weights()[i]); + } + // The averaged rule is exact on x^(2N + 1), whose integral is Gamma(2N + 2 + alpha). + auto f = [](double x) { return std::pow(x, 2 * N + 1); }; + CHECK_ULP_CLOSE(boost::math::tgamma(2 * N + 2 + alpha), rule.integrate_averaged(f), 64); +} + +void test_alpha_domain() +{ + using ignore_policy = boost::math::policies::policy >; + for (double alpha : { -1.0, -2.0, std::numeric_limits::quiet_NaN(), std::numeric_limits::infinity() }) + { + double error = 0; + CHECK_NAN((gauss_laguerre(alpha).integrate([](double x) { return x; }))); + CHECK_NAN((gauss_laguerre(alpha).integrate_averaged([](double x) { return x; }, &error))); + CHECK_NAN(error); +#ifndef BOOST_NO_EXCEPTIONS + CHECK_THROW((gauss_laguerre(alpha)), std::domain_error); +#endif + } +} + +template +void test_moments_all_N() +{ + test_moments(); + test_moments(); + test_moments(); + test_moments(); + test_moments(); + test_moments(); +} + +int main() +{ + // Factorial of floats overflow for larger N > 10 + test_moments(); + test_moments(); + test_moments(); + test_moments(); + + test_moments_all_N(); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_moments_all_N(); +#endif +#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS + test_moments(); + test_moments(); +#endif + test_moments_all_N(); + test_moments(); + test_moments(); +#ifdef BOOST_HAS_FLOAT128 + test_moments(); + test_moments(); +#endif + + test_tables_match_on_demand(); + test_tables_match_on_demand(); + test_tables_match_on_demand(); + using cpp_bin_float_100 = boost::multiprecision::cpp_bin_float_100; + test_tables_match_on_demand(); + test_tables_match_on_demand(); + test_tables_match_on_demand(); + + test_complex, 30>(50 * std::numeric_limits::epsilon()); + test_complex(50 * std::numeric_limits::epsilon()); + + test_vector_valued<30>(); + test_vector_valued<40>(); + + // Overflow occurs at about N = 200 + test_large_N<100>(); + test_large_N<190>(); + + test_overflow(); + + for (unsigned n = 1; n <= 50; ++n) + { + test_anti_gauss(n); + test_anti_gauss(n, -0.9); + test_anti_gauss(n, -0.99); + test_anti_gauss(n, 2.5); + } +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_anti_gauss(7); + test_anti_gauss(20); +#endif + test_error_estimate(); + test_error_estimate(); + test_error_estimate(); + test_error_estimate(); + test_error_estimate(); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_error_estimate(); +#endif + for (double alpha : { -0.9, -0.5, 0.5, 2.5 }) + { + const double tol = 4; + test_alpha(alpha, tol); + test_alpha(alpha, tol); + test_alpha(alpha, tol); + test_alpha(alpha, tol); + test_alpha(alpha, tol); + } +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_alpha(-0.5L, 4); +#endif + test_alpha_domain(); + test_anti_gauss_accessors<7>(0); + test_anti_gauss_accessors<10>(0.5); + + return boost::math::test::report_errors(); +} diff --git a/test/test_gauss_laguerre_quadrature_eigen.cpp b/test/test_gauss_laguerre_quadrature_eigen.cpp new file mode 100644 index 0000000000..33b2eb7001 --- /dev/null +++ b/test/test_gauss_laguerre_quadrature_eigen.cpp @@ -0,0 +1,37 @@ +// 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) + +#include "math_unit_test.hpp" +#include "gauss_quadrature_test_helpers.hpp" +#include +#include +#include + +using boost::math::quadrature::gauss_laguerre; + +template +void test_matrix(const Matrix& zero) +{ + // Against exp(-x): 1 -> 1, e^(-x) * x -> 1, e^(-x) x^2 -> 2, e^(-x) x^3 -> 6. + const Matrix Q = integrate_matrix(gauss_laguerre(), zero); + CHECK_ULP_CLOSE(1.0, Q(0, 0).real(), 4); + CHECK_ULP_CLOSE(2.0, Q(0, 0).imag(), 32); + CHECK_ULP_CLOSE(1.0, Q(0, 1).real(), 4); + CHECK_ULP_CLOSE(2.0, Q(0, 1).imag(), 4); + CHECK_ULP_CLOSE(-2.0, Q(1, 0).real(), 32); + CHECK_ULP_CLOSE(Q(1, 0).imag(), 1.0, 4); + CHECK_ULP_CLOSE(6.0, Q(1, 1).real(), 16); + CHECK_ULP_CLOSE(-1.0, Q(1, 1).imag(), 4); +} + +int main() +{ + using fixed = Eigen::Matrix, 2, 2>; + test_matrix(fixed::Zero()); + test_matrix(fixed::Zero()); + test_matrix(Eigen::MatrixXcd::Zero(2, 2)); + return boost::math::test::report_errors(); +} diff --git a/tools/gauss_hermite_constants.cpp b/tools/gauss_hermite_constants.cpp index 2425c3a197..8f06949d02 100644 --- a/tools/gauss_hermite_constants.cpp +++ b/tools/gauss_hermite_constants.cpp @@ -3,29 +3,33 @@ // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) // -// Prints the tables of abscissas and weights in boost/math/quadrature/gauss_hermite.hpp. +// Prints the tables of abscissas and weights in boost/math/quadrature/gauss_hermite.hpp, or with the argument +// laguerre, in boost/math/quadrature/gauss_laguerre.hpp. // Category 999 forces the values to be computed on demand, at the precision of mp_type. #include #include #include +#include #include +#include #include template -void print_gauss_hermite_constants(const char* suffix, int prec, int tag) +void print_gauss_hermite_constants(const char* suffix, int prec, int tag, const char* family) { const std::string suffix_append = std::string(suffix) + ")"; const char* type = prec > 40 ? "T" : "storage_type"; const auto& ab = T::abscissa(); const auto& w = T::weights(); std::cout << std::setprecision(prec) << std::scientific; - std::size_t order = (ab[0] == 0) ? (ab.size() * 2) - 1 : ab.size() * 2; + // The Hermite tables hold the nonnegative abscissas only. + std::size_t order = std::strcmp(family, "hermite") ? ab.size() : (ab[0] == 0) ? (ab.size() * 2) - 1 : ab.size() * 2; std::cout << "template \n" - "class gauss_hermite_detail\n" + "class gauss_" << family << "_detail\n" "{\n" - " using storage_type = typename gauss_hermite_constant_category::storage_type;\n" + " using storage_type = typename quadrature_constant_category::storage_type;\n" " public:\n" " static std::array<" << type << ", " << ab.size() << "> const & abscissa()\n" " {\n" @@ -48,23 +52,31 @@ void print_gauss_hermite_constants(const char* suffix, int prec, int tag) " };\n\n"; } -template -void print_all(void) +using mp_type = boost::multiprecision::number >; + +template +void print_all(const char* family) { - using mp_type = boost::multiprecision::number >; - using rule = boost::math::quadrature::detail::gauss_hermite_detail; std::cout << "#ifndef BOOST_HAS_FLOAT128\n"; - print_gauss_hermite_constants("L", 35, 0); + print_gauss_hermite_constants("L", 35, 0, family); std::cout << "#else\n"; - print_gauss_hermite_constants("Q", 35, 0); + print_gauss_hermite_constants("Q", 35, 0, family); std::cout << "#endif\n"; - print_gauss_hermite_constants("", 115, 4); + print_gauss_hermite_constants("", 115, 4, family); } -int main() +int main(int argc, char* argv[]) { - print_all<7>(); - print_all<10>(); - print_all<15>(); + using namespace boost::math::quadrature::detail; + if (argc > 1 && !std::strcmp(argv[1], "laguerre")) + { + print_all >("laguerre"); + print_all >("laguerre"); + print_all >("laguerre"); + return 0; + } + print_all >("hermite"); + print_all >("hermite"); + print_all >("hermite"); return 0; }