From dd7e6d0fb9b9ca9619bd7cb273b121d770367860 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Mon, 12 Oct 2026 18:21:14 -0700 Subject: [PATCH] Add the generalized extreme value distribution, and port test_negative_binomial to math_unit_test Location, scale and shape as on Wikipedia, with log1p/expm1 forms continuous through the Gumbel case and moments summed from the zeta series of lgamma where their Gamma combinations cancel. The negative binomial quantile checks now scale with the quantile's condition number, which reached 548 near zero and failed windows_gcc. The nanobind distributions get their parameters as properties and a repr that rebuilds them. SYCL tests build with -fp-model=precise, dropping three fast-math workarounds. Co-Authored-By: Claude Opus 5.5 --- doc/distributions/dist_reference.qbk | 1 + .../generalized_extreme_value.qbk | 116 +++ include/boost/math/distributions.hpp | 1 + include/boost/math/distributions/fwd.hpp | 4 + .../generalized_extreme_value.hpp | 647 ++++++++++++ nanobind/src/distributions.hpp | 32 + nanobind/src/distributions_continuous.cpp | 30 + nanobind/src/distributions_discrete.cpp | 6 + nanobind/src/distributions_heavy.cpp | 7 + nanobind/src/distributions_noncentral.cpp | 4 + nanobind/tests/test_smoke.py | 18 + .../generalized_extreme_value_performance.cpp | 109 +++ test/CMakeLists.txt | 6 +- test/Jamfile.v2 | 10 +- .../distribution_concept_check.cpp | 1 + test/compile_test/instantiate.hpp | 3 + ...st_dist_generalized_extreme_value_incl.cpp | 25 + test/math_unit_test.hpp | 11 +- test/test_dist_deduction_guides.cpp | 5 + test/test_generalized_extreme_value.cpp | 443 +++++++++ test/test_negative_binomial.cpp | 918 +++++------------- test/test_out_of_range.hpp | 7 +- test/test_pow1p.cpp | 5 - 23 files changed, 1713 insertions(+), 696 deletions(-) create mode 100644 doc/distributions/generalized_extreme_value.qbk create mode 100644 include/boost/math/distributions/generalized_extreme_value.hpp create mode 100644 reporting/performance/generalized_extreme_value_performance.cpp create mode 100644 test/compile_test/test_dist_generalized_extreme_value_incl.cpp create mode 100644 test/test_generalized_extreme_value.cpp diff --git a/doc/distributions/dist_reference.qbk b/doc/distributions/dist_reference.qbk index d8e3028ea6..de34e51d20 100644 --- a/doc/distributions/dist_reference.qbk +++ b/doc/distributions/dist_reference.qbk @@ -13,6 +13,7 @@ [include empirical_cdf.qbk] [include exponential.qbk] [include extreme_value.qbk] +[include generalized_extreme_value.qbk] [include fisher.qbk] [include gamma.qbk] [include geometric.qbk] diff --git a/doc/distributions/generalized_extreme_value.qbk b/doc/distributions/generalized_extreme_value.qbk new file mode 100644 index 0000000000..4b6980ae97 --- /dev/null +++ b/doc/distributions/generalized_extreme_value.qbk @@ -0,0 +1,116 @@ +[section:gev_dist Generalized Extreme Value Distribution] + +``#include `` + + template + class generalized_extreme_value_distribution; + + typedef generalized_extreme_value_distribution<> generalized_extreme_value; + + template + class generalized_extreme_value_distribution + { + public: + typedef RealType value_type; + + generalized_extreme_value_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0); + + RealType location()const; + RealType scale()const; + RealType shape()const; + }; + +The [@https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution generalized extreme value distribution] +is the family of limit distributions of the normalized maxima of independent, identically distributed +random variables. With location [mu], scale [sigma] and shape [xi], its CDF is + +[expression F(x) = e[super -t(x)], t(x) = (1 + [xi](x - [mu])/[sigma])[super -1/[xi]]] + +on the support where 1 + [xi](x - [mu])/[sigma] > 0, and its PDF is + +[expression f(x) = t(x)[super [xi]+1] e[super -t(x)] / [sigma]] + +The shape selects the type: + +* [xi] > 0, the Fr[eacute]chet type: heavy tailed to the right, and bounded below at [mu] - [sigma]/[xi]. +* [xi] = 0, the Gumbel type, where t(x) = e[super -(x-[mu])/[sigma]]: this is the __extreme_distrib. +* [xi] < 0, the reversed Weibull type: bounded above at [mu] - [sigma]/[xi]. + +Note that the sign convention for the shape varies: SciPy's `genextreme`, for example, takes c = -[xi]. + +[h4 Member Functions] + + generalized_extreme_value_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0); + +Constructs a generalized extreme value distribution with the specified location, scale and shape. + +Requires that location and shape are finite and that scale is finite and > 0, otherwise calls __domain_error. + + RealType location()const; + RealType scale()const; + RealType shape()const; + +Return the parameters of the distribution. + +[h4 Non-member Accessors] + +All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] +that are generic to all distributions are supported: __usual_accessors. + +The domain of the random variable is \[-[infin], +[infin]\]: outside the support the PDF is zero +and the CDF is 0 or 1. `support` returns the bounded end for [xi] [ne] 0. + +The mean is infinite for [xi] [ge] 1 and the variance for [xi] [ge] 1/2, and both return +[infin] there. +The skewness is undefined for [xi] [ge] 1/3 and the kurtosis for [xi] [ge] 1/4, where they call __domain_error. The quantile at the unbounded ends calls __overflow_error, +as does the PDF at the upper end point when [xi] < -1, where the density is infinite. + +`logpdf`, `logcdf` and the log of the complement of the CDF are specialized for accuracy deep in the tails. + +[h4 Accuracy] + +Against references computed with mpmath at the double values of the arguments, the PDF, CDF, their +complements and logs, and the quantiles are within a few ulps, except where the function is itself +ill-conditioned: deep in the tails of the CDF, where e[super -t] magnifies the rounding of (x - [mu])/[sigma] +(up to about 30 ulps), and where the quantile is close to zero and [mu] + [sigma] q cancels. +The mean, variance, mode, median and entropy are within a few ulps. Skewness and kurtosis are accurate to +an absolute error of a few tens of epsilon for |[xi]| < 1/4, and of a few hundred for larger |[xi]|, +where they are third and fourth differences of log[Gamma](1 - k[xi]); both change sign +(the skewness near [xi] = -0.28, the excess kurtosis near [xi] = -0.15), so their relative error is unbounded there. + +[h4 Implementation] + +In the following table [mu] is the location, [sigma] the scale, [xi] the shape, z = (x - [mu])/[sigma], +/p/ the probability and /q = 1 - p/; g[sub k] = [Gamma](1 - k[xi]). + +[table +[[Function][Implementation Notes]] +[[pdf][exp(-(1 + [xi])L - t) / [sigma], where L = -log t = log1p([xi]z)/[xi], and t = e[super -L].]] +[[logpdf][-log [sigma] - (1 + [xi])L - t]] +[[cdf][exp(-t)]] +[[logcdf][-t]] +[[cdf complement][-expm1(-t)]] +[[log of the cdf complement][log(-expm1(-t)) for t < log 2, otherwise log1p(-e[super -t])]] +[[quantile][[mu] + [sigma] expm1(-[xi] log y)/[xi] with y = -log p, or [mu] + [sigma](y[super -[xi]] - 1)/[xi] when |[xi] log y| > 1.]] +[[quantile from the complement][As the quantile, with y = -log1p(-q).]] +[[mean][[mu] + [sigma] tgamma1pm1(-[xi])/[xi]]] +[[variance][([sigma]g[sub 1]/[xi])[super 2] expm1(A[sub 2]), with A[sub k] = log g[sub k] - k log g[sub 1]]] +[[skewness][sign([xi]) (g[sub 3] - 3g[sub 1]g[sub 2] + 2g[sub 1][super 3]) / (g[sub 2] - g[sub 1][super 2])[super 3/2], written in terms of the A[sub k] + so that the cancellation at small [xi] is exact: the combinations of the A[sub k] are summed from the Taylor series + log[Gamma](1 - z) = [gamma]z + [sum] [zeta](j) z[super j]/j, in which their low orders vanish identically, while the series converges quickly.]] +[[kurtosis excess][(g[sub 4] - 4g[sub 1]g[sub 3] + 6g[sub 1][super 2]g[sub 2] - 3g[sub 1][super 4]) / (g[sub 2] - g[sub 1][super 2])[super 2] - 3, as for the skewness.]] +[[mode][[mu] + [sigma]((1 + [xi])[super -[xi]] - 1)/[xi], or the upper end point [mu] - [sigma]/[xi] for [xi] [le] -1.]] +[[median][[mu] + [sigma]((log 2)[super -[xi]] - 1)/[xi]]] +[[entropy][log [sigma] + [gamma](1 + [xi]) + 1]] +] + +At [xi] = 0, and for |[xi]| below machine epsilon in the moments, the Gumbel formulas are used. + +[endsect] [/section:gev_dist Generalized Extreme Value] + +[/ generalized_extreme_value.qbk + Copyright 2026 Nick Thompson. + Distributed under 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). +] diff --git a/include/boost/math/distributions.hpp b/include/boost/math/distributions.hpp index ca09602466..4ab64adb67 100644 --- a/include/boost/math/distributions.hpp +++ b/include/boost/math/distributions.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include diff --git a/include/boost/math/distributions/fwd.hpp b/include/boost/math/distributions/fwd.hpp index f166e57917..0be403a763 100644 --- a/include/boost/math/distributions/fwd.hpp +++ b/include/boost/math/distributions/fwd.hpp @@ -42,6 +42,9 @@ class extreme_value_distribution; BOOST_MATH_EXPORT template class fisher_f_distribution; +BOOST_MATH_EXPORT template +class generalized_extreme_value_distribution; + BOOST_MATH_EXPORT template class gamma_distribution; @@ -145,6 +148,7 @@ BOOST_MATH_NAMESPACE_END typedef BOOST_MATH_NAMESPACE::extreme_value_distribution extreme_value;\ typedef BOOST_MATH_NAMESPACE::fisher_f_distribution fisher_f;\ typedef BOOST_MATH_NAMESPACE::gamma_distribution gamma;\ + typedef BOOST_MATH_NAMESPACE::generalized_extreme_value_distribution generalized_extreme_value;\ typedef BOOST_MATH_NAMESPACE::geometric_distribution geometric;\ typedef BOOST_MATH_NAMESPACE::hypergeometric_distribution hypergeometric;\ typedef BOOST_MATH_NAMESPACE::kolmogorov_smirnov_distribution kolmogorov_smirnov;\ diff --git a/include/boost/math/distributions/generalized_extreme_value.hpp b/include/boost/math/distributions/generalized_extreme_value.hpp new file mode 100644 index 0000000000..722ac01718 --- /dev/null +++ b/include/boost/math/distributions/generalized_extreme_value.hpp @@ -0,0 +1,647 @@ +// 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_DISTRIBUTIONS_GENERALIZED_EXTREME_VALUE_HPP +#define BOOST_MATH_DISTRIBUTIONS_GENERALIZED_EXTREME_VALUE_HPP + +// The generalized extreme value distribution, with location mu, scale sigma and shape xi as in +// https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution: +// +// cdf(x) = exp(-t(x)), t(x) = (1 + xi (x - mu)/sigma)^(-1/xi), or exp(-(x - mu)/sigma) for xi = 0. +// +// xi > 0 is the Frechet type, bounded below at mu - sigma/xi; xi < 0 the reversed Weibull type, bounded +// above there; xi = 0 the Gumbel distribution, extreme_value_distribution. (SciPy's genextreme takes c = -xi.) +// log t = -log1p(xi z)/xi and the quantile's expm1(-xi log y)/xi tend to their Gumbel limits as xi -> 0 +// without cancellation. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#endif + +BOOST_MATH_NAMESPACE_BEGIN + +BOOST_MATH_EXPORT template > +class generalized_extreme_value_distribution +{ +public: + using value_type = RealType; + using policy_type = Policy; + + explicit generalized_extreme_value_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0) + : m_location(location), m_scale(scale), m_shape(shape) + { + constexpr auto function = "boost::math::generalized_extreme_value_distribution<%1%>::generalized_extreme_value_distribution"; + RealType result; + detail::check_location(function, location, &result, Policy()); + detail::check_scale(function, scale, &result, Policy()); + detail::check_finite(function, shape, &result, Policy()); + } + + RealType location() const { return m_location; } + RealType scale() const { return m_scale; } + RealType shape() const { return m_shape; } + +private: + RealType m_location; + RealType m_scale; + RealType m_shape; +}; + +BOOST_MATH_EXPORT using generalized_extreme_value = generalized_extreme_value_distribution; + +#ifdef __cpp_deduction_guides +BOOST_MATH_EXPORT template +generalized_extreme_value_distribution(RealType)->generalized_extreme_value_distribution::type>; +BOOST_MATH_EXPORT template +generalized_extreme_value_distribution(RealType, RealType)->generalized_extreme_value_distribution::type>; +BOOST_MATH_EXPORT template +generalized_extreme_value_distribution(RealType, RealType, RealType)->generalized_extreme_value_distribution::type>; +#endif + +namespace detail { + +template +inline bool check_gev(const char* function, const generalized_extreme_value_distribution& dist, RealType* result) +{ + return check_location(function, dist.location(), result, Policy()) + && check_scale(function, dist.scale(), result, Policy()) + && check_finite(function, dist.shape(), result, Policy()); +} + +// Where x lies relative to the support, and inside it L = -log t(x) = log1p(xi z)/xi, so that cdf(x) = exp(-exp(-L)). +enum class gev_region { below, inside, endpoint, above }; + +// The support is decided against the end point mu - sigma/xi that support() reports, before x is standardized, +// whose rounding could otherwise carry a point outside onto it. Inside, 1 + xi z can still round to 0 or below +// next to the end point, where the end point's limit is the answer to working precision. +template +gev_region gev_log_t(RealType x, RealType mu, RealType sigma, RealType xi, RealType* L) +{ + BOOST_MATH_STD_USING + RealType z = (x - mu) / sigma; + if (xi == 0) + { + *L = z; + return gev_region::inside; + } + RealType endpoint = mu - sigma / xi; + if (x == endpoint) + { + return gev_region::endpoint; + } + if ((xi > 0) == (x < endpoint)) + { + // Beyond the lower end point for xi > 0, the upper for xi < 0: + return xi > 0 ? gev_region::below : gev_region::above; + } + RealType s = xi * z; + if (s <= -1) + { + return gev_region::endpoint; + } + // For tiny xi z, xi z may have lost bits to underflow; log1p(s)/xi = z (1 - s/2 + ...): + *L = abs(s) < tools::epsilon() ? z * (1 - s / 2) : RealType(BOOST_MATH_NAMESPACE::log1p(s) / xi); + return gev_region::inside; +} + +// (y^(-xi) - 1)/xi, the standardized quantile at t = y; -log y for xi = 0. expm1(-xi log y) carries +// log's absolute error times xi, so where |xi log y| > 1 pow, which then has nothing to cancel, is better. +template +RealType gev_standard_quantile(RealType y, RealType log_y, RealType xi, const Policy& pol) +{ + BOOST_MATH_STD_USING + RealType a = -xi * log_y; + if (xi == 0 || abs(a) < tools::epsilon()) + { + return -log_y * (1 + a / 2); + } + if (abs(a) > 1) + { + return (pow(y, -xi) - 1) / xi; + } + return BOOST_MATH_NAMESPACE::expm1(a, pol) / xi; +} + +// expm1(x) - x, without the cancellation for small x. +template +RealType gev_expm1mx(RealType x, const Policy& pol) +{ + BOOST_MATH_STD_USING + if (abs(x) > RealType(0.5)) + { + return BOOST_MATH_NAMESPACE::expm1(x, pol) - x; + } + RealType term = x * x / 2; + RealType sum = term; + for (unsigned n = 3; abs(term) > tools::epsilon() * abs(sum); ++n) + { + term *= x / n; + sum += term; + } + return sum; +} + +// sum_k c[k] A_k for k = 2, 3, 4, where A_k = lgamma(1 - k xi) - k lgamma(1 - xi) is the log of +// E[W^k]/E[W]^k for W = t^(-xi), the moments of which are Gamma(1 - k xi). The combinations the +// moments need cancel at leading order(s) in xi, so for small xi this sums the Taylor series +// lgamma(1 - z) = euler z + sum_{j >= 2} zeta(j) z^j / j term by term, in which they cancel exactly: +// sum_k c[k] A_k = sum_{j >= 2} zeta(j)/j sum_k c[k] (k^j - k) xi^j. +// The series converges for |k xi| < 1 with ratio |k xi|; it is used to |k xi| = 3/4 for the largest k, +// beyond which the lgamma differences lose at most a digit or so. +template +RealType gev_log_moment_combination(const int (&c)[3], RealType xi, const Policy& pol) +{ + BOOST_MATH_STD_USING + int k_max = c[2] != 0 ? 4 : (c[1] != 0 ? 3 : 2); + if (k_max * abs(xi) > RealType(0.75)) + { + RealType l1 = BOOST_MATH_NAMESPACE::lgamma(1 - xi, pol); + RealType sum = 0; + for (int k = 2; k <= 4; ++k) + { + if (c[k - 2] != 0) + { + sum += c[k - 2] * (BOOST_MATH_NAMESPACE::lgamma(1 - k * xi, pol) - k * l1); + } + } + return sum; + } + RealType power[3] = { 2, 3, 4 }; // k^j + RealType xi_j = xi; // xi^j + RealType sum = 0; + const std::uintmax_t max_terms = policies::get_max_series_iterations(); + for (std::uintmax_t j = 2; ; ++j) + { + RealType weight = 0; + for (int k = 2; k <= 4; ++k) + { + power[k - 2] *= k; + weight += c[k - 2] * (power[k - 2] - k); + } + xi_j *= xi; + RealType term = BOOST_MATH_NAMESPACE::zeta(RealType(j), pol) / j * weight * xi_j; + sum += term; + // The weights of the leading orders are zero, so only stop on a nonzero term: + if (weight != 0 && abs(term) <= tools::epsilon() * abs(sum)) + { + return sum; + } + if (j > max_terms) + { + return policies::raise_evaluation_error("boost::math::generalized_extreme_value moments", "Series did not converge for shape %1%.", xi, pol); + } + } +} + +// The moments of xi this small are their Gumbel limits to working precision. +template +bool gev_is_gumbel(RealType xi) +{ + BOOST_MATH_STD_USING + return abs(xi) < tools::epsilon(); +} + +// The mean for xi >= 1 and the variance for xi >= 1/2, which diverge to +infinity (max_value for a type without one). +template +RealType gev_infinity() +{ + return std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : tools::max_value(); +} + +template +RealType gev_undefined_moment(const char* function, const char* message, RealType xi) +{ + return policies::raise_domain_error(function, message, xi, Policy()); +} + +} // namespace detail + +BOOST_MATH_EXPORT template +inline std::pair range(const generalized_extreme_value_distribution& /*dist*/) +{ // Range of permissible values for random variable x; the pdf is zero outside the support. + using BOOST_MATH_NAMESPACE::tools::max_value; + return std::pair( + std::numeric_limits::has_infinity ? -std::numeric_limits::infinity() : -max_value(), + std::numeric_limits::has_infinity ? std::numeric_limits::infinity() : max_value()); +} + +BOOST_MATH_EXPORT template +inline std::pair support(const generalized_extreme_value_distribution& dist) +{ // Range where the cdf rises from 0 to 1: bounded at mu - sigma/xi below for xi > 0, above for xi < 0. + using BOOST_MATH_NAMESPACE::tools::max_value; + RealType xi = dist.shape(); + RealType endpoint = xi == 0 ? RealType(0) : RealType(dist.location() - dist.scale() / xi); + if (xi > 0) + { + return std::pair(endpoint, max_value()); + } + if (xi < 0) + { + return std::pair(-max_value(), endpoint); + } + return std::pair(-max_value(), max_value()); +} + +BOOST_MATH_EXPORT template +inline RealType logpdf(const generalized_extreme_value_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::logpdf(const generalized_extreme_value_distribution<%1%>&, %1%)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result) || !detail::check_x_not_NaN(function, x, &result, Policy())) + { + return result; + } + RealType xi = dist.shape(); + if ((BOOST_MATH_NAMESPACE::isinf)(x)) + { + return -std::numeric_limits::infinity(); + } + RealType L; + switch (detail::gev_log_t(x, dist.location(), dist.scale(), xi, &L)) + { + case detail::gev_region::inside: + break; + case detail::gev_region::endpoint: + // The bounded end, where t = 0 for xi < 0 and the density is t^(1 + xi)/sigma; t = infinity for xi > 0. + if (xi == -1) + { + return -log(dist.scale()); + } + if (xi < -1) + { + return policies::raise_overflow_error(function, nullptr, Policy()); + } + return -std::numeric_limits::infinity(); + default: + return -std::numeric_limits::infinity(); + } + // log pdf = -log sigma + (1 + xi) log t - t, with log t = -L: + RealType t = exp(-L); + if ((BOOST_MATH_NAMESPACE::isinf)(t)) + { + return -std::numeric_limits::infinity(); + } + return -log(dist.scale()) - (1 + xi) * L - t; +} + +BOOST_MATH_EXPORT template +inline RealType pdf(const generalized_extreme_value_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::pdf(const generalized_extreme_value_distribution<%1%>&, %1%)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result) || !detail::check_x_not_NaN(function, x, &result, Policy())) + { + return result; + } + RealType xi = dist.shape(); + if ((BOOST_MATH_NAMESPACE::isinf)(x)) + { + return 0; + } + RealType L; + switch (detail::gev_log_t(x, dist.location(), dist.scale(), xi, &L)) + { + case detail::gev_region::inside: + break; + case detail::gev_region::endpoint: + if (xi == -1) + { + return 1 / dist.scale(); + } + if (xi < -1) + { + return policies::raise_overflow_error(function, nullptr, Policy()); + } + return 0; + default: + return 0; + } + RealType t = exp(-L); + if ((BOOST_MATH_NAMESPACE::isinf)(t)) + { + return 0; + } + // t^(1 + xi) exp(-t) / sigma: + return exp(-(1 + xi) * L - t) / dist.scale(); +} + +namespace detail { + +// cdf, its complement, and their logs; region and t decide them all. +template +RealType gev_cdf_imp(const generalized_extreme_value_distribution& dist, const RealType& x, bool complement, bool take_log, const char* function) +{ + BOOST_MATH_STD_USING + RealType result = 0; + if (!check_gev(function, dist, &result) || !check_x_not_NaN(function, x, &result, Policy())) + { + return result; + } + const RealType zero = take_log ? RealType(-std::numeric_limits::infinity()) : RealType(0); + const RealType one = take_log ? RealType(0) : RealType(1); + if ((BOOST_MATH_NAMESPACE::isinf)(x)) + { + return (x < 0) != complement ? zero : one; + } + RealType xi = dist.shape(); + RealType L; + switch (gev_log_t(x, dist.location(), dist.scale(), xi, &L)) + { + case gev_region::inside: + break; + case gev_region::endpoint: + return (xi > 0) != complement ? zero : one; + case gev_region::below: + return complement ? one : zero; + default: + return complement ? zero : one; + } + RealType t = exp(-L); + if (!complement) + { + return take_log ? RealType(-t) : RealType(exp(-t)); + } + RealType q = -BOOST_MATH_NAMESPACE::expm1(-t, Policy()); + if (!take_log) + { + return q; + } + // log(1 - exp(-t)): log t - t/2 = -L - t/2 once t is below epsilon (t may have underflowed), through q + // while it is small, else through log1p: + if (t < tools::epsilon()) + { + return -L - t / 2; + } + return t () ? RealType(log(q)) : RealType(BOOST_MATH_NAMESPACE::log1p(-exp(-t), Policy())); +} + +template +RealType gev_quantile_imp(const generalized_extreme_value_distribution& dist, const RealType& p, bool complement, const char* function) +{ + BOOST_MATH_STD_USING + RealType result = 0; + if (!check_gev(function, dist, &result) || !check_probability(function, p, &result, Policy())) + { + return result; + } + RealType xi = dist.shape(); + RealType mu = dist.location(); + RealType sigma = dist.scale(); + // P = cdf(x): + bool lower_end = complement ? p == 1 : p == 0; + bool upper_end = complement ? p == 0 : p == 1; + if (lower_end) + { + return xi > 0 ? RealType(mu - sigma / xi) : RealType(-policies::raise_overflow_error(function, nullptr, Policy())); + } + if (upper_end) + { + return xi < 0 ? RealType(mu - sigma / xi) : policies::raise_overflow_error(function, nullptr, Policy()); + } + // t = -log P: + RealType y = complement ? RealType(-BOOST_MATH_NAMESPACE::log1p(-p, Policy())) : RealType(-log(p)); + return mu + sigma * gev_standard_quantile(y, RealType(log(y)), xi, Policy()); +} + +} // namespace detail + +BOOST_MATH_EXPORT template +inline RealType cdf(const generalized_extreme_value_distribution& dist, const RealType& x) +{ + return detail::gev_cdf_imp(dist, x, false, false, "boost::math::cdf(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType cdf(const complemented2_type, RealType>& c) +{ + return detail::gev_cdf_imp(c.dist, c.param, true, false, "boost::math::cdf(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType logcdf(const generalized_extreme_value_distribution& dist, const RealType& x) +{ + return detail::gev_cdf_imp(dist, x, false, true, "boost::math::logcdf(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType logcdf(const complemented2_type, RealType>& c) +{ + return detail::gev_cdf_imp(c.dist, c.param, true, true, "boost::math::logcdf(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType quantile(const generalized_extreme_value_distribution& dist, const RealType& p) +{ + return detail::gev_quantile_imp(dist, p, false, "boost::math::quantile(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType quantile(const complemented2_type, RealType>& c) +{ + return detail::gev_quantile_imp(c.dist, c.param, true, "boost::math::quantile(const generalized_extreme_value_distribution<%1%>&, %1%)"); +} + +BOOST_MATH_EXPORT template +inline RealType mean(const generalized_extreme_value_distribution& dist) +{ + constexpr auto function = "boost::math::mean(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (xi >= 1) + { + // The right tail decays like x^(-1/xi), so the mean is +infinity: + return detail::gev_infinity(); + } + if (detail::gev_is_gumbel(xi)) + { + return mean(extreme_value_distribution(dist.location(), dist.scale())); + } + // mu + sigma (Gamma(1 - xi) - 1)/xi: + return dist.location() + dist.scale() * BOOST_MATH_NAMESPACE::tgamma1pm1(-xi, Policy()) / xi; +} + +BOOST_MATH_EXPORT template +inline RealType variance(const generalized_extreme_value_distribution& dist) +{ + constexpr auto function = "boost::math::variance(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (2 * xi >= 1) + { + // The second moment diverges, so the variance is +infinity: + return detail::gev_infinity(); + } + RealType sigma = dist.scale(); + if (detail::gev_is_gumbel(xi)) + { + return sigma * sigma * constants::pi_sqr() / 6; + } + // sigma^2 (Gamma(1 - 2 xi) - Gamma(1 - xi)^2)/xi^2 = (sigma Gamma(1 - xi)/xi)^2 expm1(A_2): + const int a2[3] = { 1, 0, 0 }; + RealType g1 = BOOST_MATH_NAMESPACE::tgamma(1 - xi, Policy()); + RealType s = sigma * g1 / xi; + return s * s * BOOST_MATH_NAMESPACE::expm1(detail::gev_log_moment_combination(a2, xi, Policy()), Policy()); +} + +BOOST_MATH_EXPORT template +inline RealType standard_deviation(const generalized_extreme_value_distribution& dist) +{ + BOOST_MATH_STD_USING + return sqrt(variance(dist)); +} + +BOOST_MATH_EXPORT template +inline RealType skewness(const generalized_extreme_value_distribution& dist) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::skewness(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (3 * xi >= 1) + { + return detail::gev_undefined_moment(function, "The skewness is undefined for shape >= 1/3, but got %1%.", xi); + } + if (detail::gev_is_gumbel(xi)) + { + return skewness(extreme_value_distribution()); + } + // sign(xi) (g3 - 3 g1 g2 + 2 g1^3)/(g2 - g1^2)^(3/2) with g_k = Gamma(1 - k xi); over g1^3 the numerator + // is expm1(A_3) - 3 expm1(A_2) = (A_3 - 3 A_2) + expm1mx(A_3) - 3 expm1mx(A_2), which is O(xi^3): + const int a2[3] = { 1, 0, 0 }; + const int a3[3] = { 0, 1, 0 }; + const int b[3] = { -3, 1, 0 }; + RealType A2 = detail::gev_log_moment_combination(a2, xi, Policy()); + RealType A3 = detail::gev_log_moment_combination(a3, xi, Policy()); + RealType numerator = detail::gev_log_moment_combination(b, xi, Policy()) + detail::gev_expm1mx(A3, Policy()) - 3 * detail::gev_expm1mx(A2, Policy()); + RealType v = BOOST_MATH_NAMESPACE::expm1(A2, Policy()); + RealType skew = numerator / (v * sqrt(v)); + return xi > 0 ? skew : RealType(-skew); +} + +BOOST_MATH_EXPORT template +inline RealType kurtosis_excess(const generalized_extreme_value_distribution& dist) +{ + constexpr auto function = "boost::math::kurtosis_excess(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (4 * xi >= 1) + { + return detail::gev_undefined_moment(function, "The kurtosis is undefined for shape >= 1/4, but got %1%.", xi); + } + if (detail::gev_is_gumbel(xi)) + { + return kurtosis_excess(extreme_value_distribution()); + } + // (g4 - 4 g1 g3 + 6 g1^2 g2 - 3 g1^4)/(g2 - g1^2)^2 - 3; over g1^4 the numerator is + // expm1(A_4) - 4 expm1(A_3) + 6 expm1(A_2), whose part linear in the A_k is O(xi^4): + const int a2[3] = { 1, 0, 0 }; + const int a3[3] = { 0, 1, 0 }; + const int a4[3] = { 0, 0, 1 }; + const int c[3] = { 6, -4, 1 }; + RealType A2 = detail::gev_log_moment_combination(a2, xi, Policy()); + RealType A3 = detail::gev_log_moment_combination(a3, xi, Policy()); + RealType A4 = detail::gev_log_moment_combination(a4, xi, Policy()); + RealType numerator = detail::gev_log_moment_combination(c, xi, Policy()) + + detail::gev_expm1mx(A4, Policy()) - 4 * detail::gev_expm1mx(A3, Policy()) + 6 * detail::gev_expm1mx(A2, Policy()); + RealType v = BOOST_MATH_NAMESPACE::expm1(A2, Policy()); + return numerator / (v * v) - 3; +} + +BOOST_MATH_EXPORT template +inline RealType kurtosis(const generalized_extreme_value_distribution& dist) +{ + return kurtosis_excess(dist) + 3; +} + +BOOST_MATH_EXPORT template +inline RealType mode(const generalized_extreme_value_distribution& dist) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::mode(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + RealType xi = dist.shape(); + if (xi <= -1) + { + // The density increases to the upper endpoint: + return dist.location() - dist.scale() / xi; + } + // mu + sigma ((1 + xi)^(-xi) - 1)/xi, where t = 1 + xi: + return dist.location() + dist.scale() * detail::gev_standard_quantile(RealType(1 + xi), RealType(BOOST_MATH_NAMESPACE::log1p(xi, Policy())), xi, Policy()); +} + +BOOST_MATH_EXPORT template +inline RealType median(const generalized_extreme_value_distribution& dist) +{ + constexpr auto function = "boost::math::median(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + // t = log 2: + return dist.location() + dist.scale() * detail::gev_standard_quantile(constants::ln_two(), constants::ln_ln_two(), dist.shape(), Policy()); +} + +BOOST_MATH_EXPORT template +inline RealType entropy(const generalized_extreme_value_distribution& dist) +{ + BOOST_MATH_STD_USING + constexpr auto function = "boost::math::entropy(const generalized_extreme_value_distribution<%1%>&)"; + RealType result = 0; + if (!detail::check_gev(function, dist, &result)) + { + return result; + } + return log(dist.scale()) + constants::euler() * (1 + dist.shape()) + 1; +} + +BOOST_MATH_NAMESPACE_END + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include + +#endif // BOOST_MATH_DISTRIBUTIONS_GENERALIZED_EXTREME_VALUE_HPP diff --git a/nanobind/src/distributions.hpp b/nanobind/src/distributions.hpp index 6f6e5f65fe..44852970c1 100644 --- a/nanobind/src/distributions.hpp +++ b/nanobind/src/distributions.hpp @@ -14,7 +14,10 @@ #include +#include +#include #include +#include #include "vectorize.hpp" @@ -70,6 +73,35 @@ nb::class_& add_dist_methods(nb::class_& cls) { return cls; } +// A parameter: the constructor's argument name, and its value read back from the distribution. +template +using dist_param = std::pair; + +#define BMPY_PARAM(name, accessor) \ + bmpy::dist_param{name, +[](const D& d) { return nb::cast(d.accessor()); }} + +// Exposes the parameters as read-only properties, under the constructor's argument names, and a repr that +// constructs the same distribution: normal_distribution(mean=0.0, sd=1.0). A parameter named like a method +// (the mean of a normal, the mode of a triangular) stays the method, which returns the same value. +template +nb::class_& add_dist_params(nb::class_& cls, std::initializer_list> params) { + std::vector> ps(params); + for (auto const& [name, get] : ps) { + if (!nb::hasattr(cls, name)) { + cls.def_prop_ro(name, get); + } + } + std::string type = nb::cast(cls.attr("__name__")); + cls.def("__repr__", [type, ps](const D& d) { + std::string s = type + "("; + for (std::size_t i = 0; i < ps.size(); ++i) { + s += (i ? ", " : "") + std::string(ps[i].first) + "=" + nb::cast(nb::repr(ps[i].second(d))); + } + return s + ")"; + }); + return cls; +} + } // namespace bmpy #endif diff --git a/nanobind/src/distributions_continuous.cpp b/nanobind/src/distributions_continuous.cpp index 87757c0927..4c9a6c51fb 100644 --- a/nanobind/src/distributions_continuous.cpp +++ b/nanobind/src/distributions_continuous.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -40,48 +41,63 @@ void bmpy::bind_distributions_continuous(nb::module_& m) { nb::class_ c(m, "arcsine_distribution"); c.def(nb::init(), arg("x_min") = 0.0, arg("x_max") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("x_min", x_min), BMPY_PARAM("x_max", x_max)}); } { using D = bm::beta_distribution; nb::class_ c(m, "beta_distribution"); c.def(nb::init(), arg("alpha") = 1.0, arg("beta") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("alpha", alpha), BMPY_PARAM("beta", beta)}); } { using D = bm::cauchy_distribution; nb::class_ c(m, "cauchy_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::chi_squared_distribution; nb::class_ c(m, "chi_squared_distribution"); c.def(nb::init(), arg("df") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom)}); } { using D = bm::exponential_distribution; nb::class_ c(m, "exponential_distribution"); c.def(nb::init(), arg("lambda_") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("lambda_", lambda)}); } { using D = bm::extreme_value_distribution; nb::class_ c(m, "extreme_value_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); + } + { + using D = bm::generalized_extreme_value_distribution; + nb::class_ c(m, "generalized_extreme_value_distribution"); + c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0, arg("shape") = 0.0); + add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale), BMPY_PARAM("shape", shape)}); } { using D = bm::fisher_f_distribution; nb::class_ c(m, "fisher_f_distribution"); c.def(nb::init(), arg("df1"), arg("df2")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df1", degrees_of_freedom1), BMPY_PARAM("df2", degrees_of_freedom2)}); } { using D = bm::gamma_distribution; nb::class_ c(m, "gamma_distribution"); c.def(nb::init(), arg("shape"), arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("shape", shape), BMPY_PARAM("scale", scale)}); } { using D = bm::inverse_chi_squared_distribution; @@ -89,83 +105,97 @@ void bmpy::bind_distributions_continuous(nb::module_& m) { c.def(nb::init(), arg("df") = 1.0); c.def(nb::init(), arg("df"), arg("scale")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom), BMPY_PARAM("scale", scale)}); } { using D = bm::inverse_gamma_distribution; nb::class_ c(m, "inverse_gamma_distribution"); c.def(nb::init(), arg("shape") = 1.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("shape", shape), BMPY_PARAM("scale", scale)}); } { using D = bm::inverse_gaussian_distribution; nb::class_ c(m, "inverse_gaussian_distribution"); c.def(nb::init(), arg("mean") = 1.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("mean", mean), BMPY_PARAM("scale", scale)}); } { using D = bm::laplace_distribution; nb::class_ c(m, "laplace_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::logistic_distribution; nb::class_ c(m, "logistic_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::lognormal_distribution; nb::class_ c(m, "lognormal_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::normal_distribution; nb::class_ c(m, "normal_distribution"); c.def(nb::init(), arg("mean") = 0.0, arg("sd") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("mean", mean), BMPY_PARAM("sd", standard_deviation)}); } { using D = bm::pareto_distribution; nb::class_ c(m, "pareto_distribution"); c.def(nb::init(), arg("scale") = 1.0, arg("shape") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("scale", scale), BMPY_PARAM("shape", shape)}); } { using D = bm::rayleigh_distribution; nb::class_ c(m, "rayleigh_distribution"); c.def(nb::init(), arg("sigma") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("sigma", sigma)}); } { using D = bm::skew_normal_distribution; nb::class_ c(m, "skew_normal_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0, arg("shape") = 0.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale), BMPY_PARAM("shape", shape)}); } { using D = bm::students_t_distribution; nb::class_ c(m, "students_t_distribution"); c.def(nb::init(), arg("df") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom)}); } { using D = bm::triangular_distribution; nb::class_ c(m, "triangular_distribution"); c.def(nb::init(), arg("lower") = -1.0, arg("mode") = 0.0, arg("upper") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("lower", lower), BMPY_PARAM("mode", mode), BMPY_PARAM("upper", upper)}); } { using D = bm::uniform_distribution; nb::class_ c(m, "uniform_distribution"); c.def(nb::init(), arg("lower") = 0.0, arg("upper") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("lower", lower), BMPY_PARAM("upper", upper)}); } { using D = bm::weibull_distribution; nb::class_ c(m, "weibull_distribution"); c.def(nb::init(), arg("shape"), arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("shape", shape), BMPY_PARAM("scale", scale)}); } } diff --git a/nanobind/src/distributions_discrete.cpp b/nanobind/src/distributions_discrete.cpp index 6aa7ad61ec..bb8cf5354a 100644 --- a/nanobind/src/distributions_discrete.cpp +++ b/nanobind/src/distributions_discrete.cpp @@ -28,35 +28,41 @@ void bmpy::bind_distributions_discrete(nb::module_& m) { nb::class_ c(m, "bernoulli_distribution"); c.def(nb::init(), arg("p") = 0.5); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("p", success_fraction)}); } { using D = bm::binomial_distribution; nb::class_ c(m, "binomial_distribution"); c.def(nb::init(), arg("n") = 1.0, arg("p") = 0.5); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("n", trials), BMPY_PARAM("p", success_fraction)}); } { using D = bm::geometric_distribution; nb::class_ c(m, "geometric_distribution"); c.def(nb::init(), arg("p")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("p", success_fraction)}); } { using D = bm::hypergeometric_distribution; nb::class_ c(m, "hypergeometric_distribution"); c.def(nb::init(), arg("r"), arg("n"), arg("N")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("r", defective), BMPY_PARAM("n", sample_count), BMPY_PARAM("N", total)}); } { using D = bm::negative_binomial_distribution; nb::class_ c(m, "negative_binomial_distribution"); c.def(nb::init(), arg("r"), arg("p")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("r", successes), BMPY_PARAM("p", success_fraction)}); } { using D = bm::poisson_distribution; nb::class_ c(m, "poisson_distribution"); c.def(nb::init(), arg("mean") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("mean", mean)}); } } diff --git a/nanobind/src/distributions_heavy.cpp b/nanobind/src/distributions_heavy.cpp index 45764e7dd1..3687b8c569 100644 --- a/nanobind/src/distributions_heavy.cpp +++ b/nanobind/src/distributions_heavy.cpp @@ -31,36 +31,42 @@ void bmpy::bind_distributions_heavy(nb::module_& m) { nb::class_ c(m, "holtsmark_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::landau_distribution; nb::class_ c(m, "landau_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::mapairy_distribution; nb::class_ c(m, "mapairy_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::saspoint5_distribution; nb::class_ c(m, "saspoint5_distribution"); c.def(nb::init(), arg("location") = 0.0, arg("scale") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("location", location), BMPY_PARAM("scale", scale)}); } { using D = bm::von_mises_distribution; nb::class_ c(m, "von_mises_distribution"); c.def(nb::init(), arg("mean") = 0.0, arg("concentration") = 1.0); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("mean", mean), BMPY_PARAM("concentration", concentration)}); } { using D = bm::kolmogorov_smirnov_distribution; nb::class_ c(m, "kolmogorov_smirnov_distribution"); c.def(nb::init(), arg("n")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("n", number_of_observations)}); } { using D = bm::hyperexponential_distribution; @@ -71,5 +77,6 @@ void bmpy::bind_distributions_heavy(nb::module_& m) { [](D* self, const std::vector& probs, const std::vector& rates) { new (self) D(probs, rates); }, arg("probabilities"), arg("rates")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("probabilities", probabilities), BMPY_PARAM("rates", rates)}); } } diff --git a/nanobind/src/distributions_noncentral.cpp b/nanobind/src/distributions_noncentral.cpp index 0e86d5c363..1aa2d335d4 100644 --- a/nanobind/src/distributions_noncentral.cpp +++ b/nanobind/src/distributions_noncentral.cpp @@ -22,23 +22,27 @@ void bmpy::bind_distributions_noncentral(nb::module_& m) { nb::class_ c(m, "non_central_beta_distribution"); c.def(nb::init(), arg("alpha"), arg("beta"), arg("lambda_")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("alpha", alpha), BMPY_PARAM("beta", beta), BMPY_PARAM("lambda_", non_centrality)}); } { using D = bm::non_central_chi_squared_distribution; nb::class_ c(m, "non_central_chi_squared_distribution"); c.def(nb::init(), arg("df"), arg("lambda_")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom), BMPY_PARAM("lambda_", non_centrality)}); } { using D = bm::non_central_f_distribution; nb::class_ c(m, "non_central_f_distribution"); c.def(nb::init(), arg("df1"), arg("df2"), arg("lambda_")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df1", degrees_of_freedom1), BMPY_PARAM("df2", degrees_of_freedom2), BMPY_PARAM("lambda_", non_centrality)}); } { using D = bm::non_central_t_distribution; nb::class_ c(m, "non_central_t_distribution"); c.def(nb::init(), arg("df"), arg("delta")); add_dist_methods(c); + bmpy::add_dist_params(c, {BMPY_PARAM("df", degrees_of_freedom), BMPY_PARAM("delta", non_centrality)}); } } diff --git a/nanobind/tests/test_smoke.py b/nanobind/tests/test_smoke.py index 5dbb5eb9a9..a1b6d4bf6e 100644 --- a/nanobind/tests/test_smoke.py +++ b/nanobind/tests/test_smoke.py @@ -421,6 +421,7 @@ def test_array_errors_warn_and_give_nan(): ("chi_squared_distribution", (3.5,)), ("exponential_distribution", (2.0,)), ("extreme_value_distribution", ()), + ("generalized_extreme_value_distribution", (0.5, 2.0, 0.1)), ("fisher_f_distribution", (3.0, 7.0)), ("gamma_distribution", (2.5, 1.5)), ("inverse_chi_squared_distribution", (4.0,)), @@ -503,6 +504,23 @@ def test_discrete_distribution(case): assert math.isfinite(dist.mean()) and math.isfinite(dist.variance()) +@pytest.mark.parametrize("case", CONTINUOUS + DISCRETE, ids=lambda c: c[0]) +def test_distribution_repr(case): + # The repr names the constructor's arguments, each a property, and rebuilds the distribution. + name, params = case + dist = getattr(bm, name)(*params) + r = repr(dist) + assert r.startswith(name + "(") + again = eval(r, vars(bm)) + assert repr(again) == r + assert np.allclose(again.cdf(dist.quantile(0.3)), 0.3, rtol=1e-12) or name in {c[0] for c in DISCRETE} + for i, value in enumerate(params): + arg = r[len(name) + 1:-1].split(", ")[i].split("=")[0] if name != "hyperexponential_distribution" else None + if arg is not None and not isinstance(value, (list, tuple)): + got = getattr(dist, arg) + assert (got() if callable(got) else got) == value + + def test_normal_against_mpmath(): dist = bm.normal_distribution(1.0, 2.0) for x in (-3.0, 0.5, 4.0): diff --git a/reporting/performance/generalized_extreme_value_performance.cpp b/reporting/performance/generalized_extreme_value_performance.cpp new file mode 100644 index 0000000000..a6b48a9b7c --- /dev/null +++ b/reporting/performance/generalized_extreme_value_performance.cpp @@ -0,0 +1,109 @@ +// (C) 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) + +// +// The generalized extreme value distribution's pdf, cdf and quantile at shape xi = 0.2, against the +// Gumbel distribution (extreme_value, which is xi = 0) as the baseline, and its skewness, whose +// moment series is the expensive part of the moments. +// +// Promotion is disabled so that each type is timed in its own arithmetic: by default, +// float is evaluated in double, and double in long double where that is wider. +// + +#include +#include +#include +#include +#include + +using no_promotion = boost::math::policies::policy, boost::math::policies::promote_double>; + +template +std::vector abscissas() +{ + std::mt19937_64 mt {42}; + std::uniform_real_distribution dis(-2, 20); + std::vector v(1024); + for (auto& x : v) + { + x = static_cast(dis(mt)); + } + return v; +} + +template +std::vector probabilities() +{ + std::mt19937_64 mt {42}; + std::uniform_real_distribution dis(0.001, 0.999); + std::vector v(1024); + for (auto& p : v) + { + p = static_cast(dis(mt)); + } + return v; +} + +template +auto gev() +{ + return boost::math::generalized_extreme_value_distribution(Real(0.5), Real(2), Real(0.2)); +} + +template +auto gumbel() +{ + return boost::math::extreme_value_distribution(Real(0.5), Real(2)); +} + +#define BENCH(name, make, fn, points) \ +template \ +void name(benchmark::State& state) \ +{ \ + auto d = make(); \ + auto v = points(); \ + std::size_t i = 0; \ + for (auto _ : state) \ + { \ + benchmark::DoNotOptimize(fn(d, v[i++ & 1023])); \ + } \ +} + +BENCH(gev_pdf, gev, boost::math::pdf, abscissas) +BENCH(gumbel_pdf, gumbel, boost::math::pdf, abscissas) +BENCH(gev_cdf, gev, boost::math::cdf, abscissas) +BENCH(gumbel_cdf, gumbel, boost::math::cdf, abscissas) +BENCH(gev_quantile, gev, boost::math::quantile, probabilities) +BENCH(gumbel_quantile, gumbel, boost::math::quantile, probabilities) + +template +void gev_skewness(benchmark::State& state) +{ + std::mt19937_64 mt {42}; + std::uniform_real_distribution dis(-0.5, 0.3); + std::vector> v; + for (int i = 0; i < 1024; ++i) + { + v.emplace_back(Real(0), Real(1), static_cast(dis(mt))); + } + std::size_t i = 0; + for (auto _ : state) + { + benchmark::DoNotOptimize(boost::math::skewness(v[i++ & 1023])); + } +} + +BENCHMARK_TEMPLATE(gev_pdf, float); +BENCHMARK_TEMPLATE(gev_pdf, double); +BENCHMARK_TEMPLATE(gumbel_pdf, double); +BENCHMARK_TEMPLATE(gev_cdf, float); +BENCHMARK_TEMPLATE(gev_cdf, double); +BENCHMARK_TEMPLATE(gumbel_cdf, double); +BENCHMARK_TEMPLATE(gev_quantile, float); +BENCHMARK_TEMPLATE(gev_quantile, double); +BENCHMARK_TEMPLATE(gumbel_quantile, double); +BENCHMARK_TEMPLATE(gev_skewness, double); + +BENCHMARK_MAIN(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 29bcf179aa..f03362f7a9 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -38,12 +38,16 @@ if(HAVE_BOOST_TEST) message(STATUS "Building boost.math with SYCL") + # icpx defaults to -fp-model=fast, which flushes denormals to zero and reassociates arithmetic, + # defeating the error-free transforms and compensated sums the library relies on; test as + # the other compilers do (contraction stays on). + set(CMAKE_CXX_COMPILER "icpx") set(CMAKE_C_COMPILER "icx") enable_testing() - boost_test_jamfile(FILE sycl_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework sycl COMPILE_DEFINITIONS BOOST_MATH_ENABLE_SYCL=1 COMPILE_OPTIONS -fsycl ) + boost_test_jamfile(FILE sycl_jamfile LINK_LIBRARIES Boost::math Boost::assert Boost::concept_check Boost::config Boost::core Boost::integer Boost::lexical_cast Boost::multiprecision Boost::predef Boost::random Boost::throw_exception Boost::unit_test_framework sycl COMPILE_DEFINITIONS BOOST_MATH_ENABLE_SYCL=1 COMPILE_OPTIONS -fsycl -fp-model=precise ) elseif (BOOST_MATH_BUILD_MODULE) message(STATUS "Building boost.math as a C++20 module") diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index c558f5c321..359080af69 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -693,6 +693,7 @@ test-suite distribution_tests : [ run test_dist_overloads.cpp /boost/test//boost_unit_test_framework ] [ run test_exponential_dist.cpp /boost/test//boost_unit_test_framework ] [ run test_extreme_value.cpp /boost/test//boost_unit_test_framework ] + [ run test_generalized_extreme_value.cpp ] [ run test_find_location.cpp pch /boost/test//boost_unit_test_framework ] [ run test_find_scale.cpp pch /boost/test//boost_unit_test_framework ] [ run test_fisher_f.cpp /boost/test//boost_unit_test_framework ] @@ -750,28 +751,28 @@ test-suite distribution_tests : [ run test_logistic_dist.cpp /boost/test//boost_unit_test_framework ] [ run test_logcdf.cpp ] [ run test_lognormal.cpp /boost/test//boost_unit_test_framework ] - [ run test_negative_binomial.cpp /boost/test//boost_unit_test_framework + [ run test_negative_binomial.cpp : # command line : # input files : # requirements TEST_FLOAT intel:off : test_negative_binomial_float ] - [ run test_negative_binomial.cpp /boost/test//boost_unit_test_framework + [ run test_negative_binomial.cpp : # command line : # input files : # requirements TEST_DOUBLE intel:off : test_negative_binomial_double ] - [ run test_negative_binomial.cpp /boost/test//boost_unit_test_framework + [ run test_negative_binomial.cpp : # command line : # input files : # requirements TEST_LDOUBLE intel:off : test_negative_binomial_long_double ] - [ run test_negative_binomial.cpp /boost/test//boost_unit_test_framework + [ run test_negative_binomial.cpp : # command line : # input files : # requirements @@ -1650,6 +1651,7 @@ test-suite concepts : [ run compile_test/test_dist_complement_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_exponential_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_extreme_val_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] + [ run compile_test/test_dist_generalized_extreme_value_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_find_location_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_find_scale_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/test_dist_fisher_f_incl.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] diff --git a/test/compile_test/distribution_concept_check.cpp b/test/compile_test/distribution_concept_check.cpp index 3f5253cd71..0102eaa3fc 100644 --- a/test/compile_test/distribution_concept_check.cpp +++ b/test/compile_test/distribution_concept_check.cpp @@ -27,6 +27,7 @@ void instantiate(RealType) function_requires > >(); function_requires > >(); function_requires > >(); + function_requires > >(); function_requires > >(); function_requires > >(); function_requires > >(); diff --git a/test/compile_test/instantiate.hpp b/test/compile_test/instantiate.hpp index 5508268903..89999e7947 100644 --- a/test/compile_test/instantiate.hpp +++ b/test/compile_test/instantiate.hpp @@ -120,6 +120,7 @@ void instantiate(RealType) function_requires > >(); function_requires > >(); function_requires > >(); + function_requires > >(); function_requires > >(); function_requires > >(); function_requires > >(); @@ -162,6 +163,7 @@ void instantiate(RealType) function_requires > >(); function_requires > >(); function_requires > >(); + function_requires > >(); function_requires > >(); function_requires > >(); function_requires > >(); @@ -199,6 +201,7 @@ void instantiate(RealType) function_requires >(); function_requires >(); function_requires >(); + function_requires >(); function_requires >(); function_requires >(); function_requires >(); diff --git a/test/compile_test/test_dist_generalized_extreme_value_incl.cpp b/test/compile_test/test_dist_generalized_extreme_value_incl.cpp new file mode 100644 index 0000000000..313f88611e --- /dev/null +++ b/test/compile_test/test_dist_generalized_extreme_value_incl.cpp @@ -0,0 +1,25 @@ +// 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) +// +// Basic sanity check that header +// #includes all the files that it needs to. +// +#include +// +// Note this header includes no other headers, this is +// important if this test is to be meaningful: +// +#include "test_compile_result.hpp" + +void compile_and_link_test() +{ + TEST_DIST_FUNC(generalized_extreme_value) +} + +template class boost::math::generalized_extreme_value_distribution >; +template class boost::math::generalized_extreme_value_distribution >; +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS +template class boost::math::generalized_extreme_value_distribution >; +#endif diff --git a/test/math_unit_test.hpp b/test/math_unit_test.hpp index db3776077b..b165a53414 100644 --- a/test/math_unit_test.hpp +++ b/test/math_unit_test.hpp @@ -307,11 +307,11 @@ bool check_absolute_error(PreciseReal expected1, Real computed, Real acceptable_ if (!std::is_integral::value) { BOOST_MATH_ASSERT_MSG(sizeof(PreciseReal) >= sizeof(Real), "The expected number must be computed in higher (or equal) precision than the number being tested."); - BOOST_MATH_ASSERT_MSG(!isnan(expected1), "Expected value cannot be a nan (use CHECK_NAN if this is your intention)."); + BOOST_MATH_ASSERT_MSG(!boost::math::isnan(expected1), "Expected value cannot be a nan (use CHECK_NAN if this is your intention)."); } BOOST_MATH_ASSERT_MSG(acceptable_error > 0, "Error must be > 0."); - if (isnan(computed)) + if (boost::math::isnan(computed)) { std::ios_base::fmtflags f(std::cerr.flags()); std::cerr << std::setprecision(3); @@ -499,8 +499,15 @@ inline int report_errors() #define CHECK_TRUE(X) boost::math::test::check_true((X), __FILE__, __func__, __LINE__) +// Without exceptions (GPU and SYCL builds, -fno-exceptions) errors don't throw, so there is nothing to check; +// as BOOST_MATH_CHECK_THROW does there, the expressions aren't evaluated. +#ifdef BOOST_MATH_NO_EXCEPTIONS +#define CHECK_THROW(x, what) +#define CHECK_NO_THROW(x) x; +#else #define CHECK_THROW(x, what) try{ x; boost::math::test::report_non_throw(__FILE__, __LINE__); }catch(const what&){} catch(...){ boost::math::test::report_incorrect_throw(__FILE__, __LINE__); } #define CHECK_NO_THROW(x) try{ x; }catch(...){ boost::math::test::report_unexpected_throw(__FILE__, __LINE__); } +#endif #endif diff --git a/test/test_dist_deduction_guides.cpp b/test/test_dist_deduction_guides.cpp index b814ee1503..a65d6a1e0e 100644 --- a/test/test_dist_deduction_guides.cpp +++ b/test/test_dist_deduction_guides.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include //#include //#include @@ -84,6 +85,10 @@ int main() test_deduction_guide(1); test_deduction_guide(1, 1); + test_deduction_guide(0); + test_deduction_guide(0, 1); + test_deduction_guide(0, 1, 0); + test_deduction_guide(1); test_deduction_guide(1); diff --git a/test/test_generalized_extreme_value.cpp b/test/test_generalized_extreme_value.cpp new file mode 100644 index 0000000000..3101ed0fb3 --- /dev/null +++ b/test/test_generalized_extreme_value.cpp @@ -0,0 +1,443 @@ +// Copyright Nick Thompson 2026. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// https://www.boost.org/LICENSE_1_0.txt) +// +// References: mpmath at 150 digits, printed to 40, at the double values of the parameters and arguments +// (the generator is in the PR that added this file). Where the tests need more than a few ulps, it is the +// condition number of the function, not the method: +// * the cdf deep in its tails, exp(-t) for large t, magnifies the rounding of (x - mu)/sigma, so the +// tolerances scale with each function's condition number; +// * mu + sigma q cancels where the quantile is near 0 (-0.016 at xi = 1.5, p = 0.25), so the quantile's tolerance +// scales with (|mu| + |x - mu|)/|x|; +// * skewness and kurtosis excess pass through 0 (at xi = -0.28 and -0.15), so they are checked to an +// absolute error, and for |xi| >= 1/4 the kurtosis is a fourth difference of lgamma(1 - k xi). + +#include "math_unit_test.hpp" +#include +#include +#include "test_out_of_range.hpp" +#include +#include +#if __has_include() +#include +#define BOOST_MATH_TEST_GEV_MULTIPRECISION +#endif + +using boost::math::generalized_extreme_value_distribution; + +// mu, sigma, xi, x, pdf, logpdf, cdf, ccdf, logcdf, logccdf +const long double cdf_data[][10] = { + { 0.5, 2, 0, -3, 0.009115765822384945444346020789780975286465L, -4.697749856565675746283731826300868938998L, 0.003168165149053243848719669162880521440672L, 0.9968318348509467561512803308371194785593L, -5.754602676005730436866499704842692370923L, -0.003173194409419889751857546828551054431638L }, + { 0.5, 2, 0, -0.5, 0.1585209605389710880473291448008405244162L, -1.841868451260073456265882909272340139729L, 0.1922956455479649280742866860378975816894L, 0.8077043544520350719257133139621024183106L, -1.648721270700128146848650787814163571654L, -0.2135591853734341440449422035125653857933L }, + { 0.5, 2, 0, 0.25, 0.1824469864020708676088842652042828343706L, -1.701295633626771626246239349269970440641L, 0.3220178007714426731329932585931980462668L, 0.6779821992285573268670067414068019537332L, -1.133148453066826316829007227811793872566L, -0.3886342462115421229163077265709873030708L }, + { 0.5, 2, 0, 1, 0.1787176730860912736876052717082097864081L, -1.721947963631350177662402388436497215372L, 0.4589560693076638184581480769101370726848L, 0.5410439306923361815418519230898629273152L, -0.7788007830714048682451702669783206472968L, -0.614254800673943822670724491688423891635L }, + { 0.5, 2, 0, 4, 0.07302735923472654790400422830811197399479L, -2.61692112401039043609794938012454758409L, 0.8404868737475783667364934360395455344152L, 0.1595131262524216332635065639604544655848L, -0.1737739434504451266807172586663710160147L, -1.835629063889961509423681941292199838984L }, + { 0.5, 2, 0, 30, 1.963931001335069415616541086995321911745e-7L, -15.44314757334629985752113467810391568766L, 0.9999996072137225924461571189333918174662L, 3.927862774075538428810666081825337687368e-7L, -3.927863545481039025566457391195863032792e-7L, -14.75000019639317084567193797879849698995L }, + { 0.5, 2, 1e-20, -3, 0.009115765822384945443841874696649368440024L, -4.697749856565675746339036679777206683281L, 0.003168165149053243848440498834569982501757L, 0.9968318348509467561515595011654300174982L, -5.754602676005730436954617058319030113406L, -0.003173194409419889751577489231521880065215L }, + { 0.5, 2, 1e-20, -0.5, 0.1585209605389710880479932047048435720762L, -1.841868451260073456261693810860715300143L, 0.192295645547964928073890383636550153991L, 0.807704354452035071926109616363449846009L, -1.648721270700128146850711689402538731724L, -0.2135591853734341440444515507147154989699L }, + { 0.5, 2, 1e-20, 0.25, 0.1824469864020708676091104260830661932527L, -1.701295633626771626244999751492866286515L, 0.3220178007714426731329647512515727226953L, 0.6779821992285573268670352487484272773047L, -1.133148453066826316829095755034689718367L, -0.3886342462115421229162656792404559695394L }, + { 0.5, 2, 1e-20, 1, 0.1787176730860912736871708313409126657638L, -1.72194796363135017766483326368120702926L, 0.4589560693076638184580363783644582656449L, 0.5410439306923361815419636216355417343551L, -0.7788007830714048682454136422230304613048L, -0.6142548006739438226705180416400223243136L }, + { 0.5, 2, 1e-20, 4, 0.07302735923472654790365016147306552690034L, -2.616921124010390436102797793633632524827L, 0.8404868737475783667342569731629820340084L, 0.1595131262524216332657430268370179659916L, -0.1737739434504451266833781721754559568711L, -1.83562906388996150940966138430666117859L }, + { 0.5, 2, 1e-20, 30, 1.963931001335069417463249117543877889429e-7L, -15.44314757334629985658082260538182204619L, 0.9999996072137225924461566916556533395052L, 3.927862774075538433083443466604947852275e-7L, -3.927863545481039029839236454264455565347e-7L, -14.75000019639317084458412569243743621741L }, + { 0.5, 2, -1e-9, -3, 0.009115765872799554771946080968262620026985L, -4.697749851035190412359091169626025457627L, 0.003168165176970276748828468916096131826271L, 0.9968318348230297232511715310839038681737L, -5.754602667193995106259567173008525775103L, -0.003173194437425649524441523294602566123376L }, + { 0.5, 2, -1e-9, -0.5, 0.1585209604725650976744218719492888455378L, -1.841868451678983297343277604902220403243L, 0.1922956455875951681974261917104534818515L, 0.8077043544124048318025738082895465181485L, -1.648721270494037988092712111127424194838L, -0.2135591854224994238167688548412075966868L }, + { 0.5, 2, -1e-9, 0.25, 0.182446986379454779733353359398335601163L, -1.701295633750731403948963542528035519475L, 0.3220178007742934072952895375955548586645L, 0.6779821992257065927047104624044451413355L, -1.133148453057973594540194954464039676811L, -0.3886342462157468559691017126026166415997L }, + { 0.5, 2, -1e-9, 1, 0.178717673129535310427908761610173438766L, -1.721947963388262653161703324233783748909L, 0.4589560693188336730278518189664221151632L, 0.5410439306811663269721481810335778848368L, -0.7788007830470673437705128830706037529678L, -0.6142548006945888275144638144559958666686L }, + { 0.5, 2, -1e-9, 4, 0.07302735927013323140638808443280455603325L, -2.616921123525549085337949864699606090224L, 0.8404868739712246545122901543310216909817L, 0.1595131260287753454877098456689783090183L, -0.17377394318435377566550942297392605831L, -1.835629065292017209709248439735100776865L }, + { 0.5, 2, -1e-9, 30, 1.963930816664273156023001573027157430863e-7L, -15.4431476673775080906334661135206673916L, 0.9999996072137653202181010821144661558863L, 3.927862346797818989178855338441136576943e-7L, -3.927863118203151757053016200034791262688e-7L, -14.7500003051744005514612834025113567628L }, + { 0.5, 2, 1e-9, -3, 0.009115765771970336145624700217464420508732L, -4.697749862096161107626641378518741212887L, 0.003168165121136211086720674613104760872243L, 0.9968318348788637889132793253868952391278L, -5.75460268481746580152711779888655526371L, -0.003173194381414130118606643465339249650594L }, + { 0.5, 2, 1e-9, -0.5, 0.1585209606053770784750314195231078862963L, -1.84186845084116361501830958938882670287L, 0.1922956455083346879279406427428083575554L, 0.8077043544916653120720593572571916424446L, -1.648721270906218305767744173580603150022L, -0.213559185324368864246791478672148497165L }, + { 0.5, 2, 1e-9, 0.25, 0.1824469864246869554900251411162551697432L, -1.701295633502811848528132688286390169526L, 0.3220178007685919389702248226276367206498L, 0.6779821992314080610297751773723632793502L, -1.133148453075679039119364116767366211481L, -0.388634246207337389862835005334260806701L }, + { 0.5, 2, 1e-9, 1, 0.1787176730426472369684846225640671376264L, -1.721947963874437702103666158788601709202L, 0.4589560692964939638920904083634490935538L, 0.5410439307035060361079095916365509064462L, -0.7788007830957423927124756903687619088155L, -0.6142548006532988178341503453350002029345L }, + { 0.5, 2, 1e-9, 4, 0.07302735919931986439709865733298120907751L, -2.616921124495231787154938124508301140614L, 0.8404868735239320791995900204472480338069L, 0.1595131264760679208004099795527519661931L, -0.1737739437165364774824976566509613042549L, -1.835629062487905812601516027788899559118L }, + { 0.5, 2, 1e-9, 30, 1.963931186005879265734471889093366625678e-7L, -15.44314747931509354621469424527979978518L, 0.999999607213679864670405489814773606136L, 3.927863201353295945101852263938640014138e-7L, -3.927863972758964370772486657174537095069e-7L, -14.75000008761194327924907972277214799304L }, + { 0.5, 2, 0.05, -3, 0.006654207329959662526934006477968915074622L, -5.012505943175803363223445387272214509505L, 0.001945434707953234552092815540172502183168L, 0.9980545652920467654479071844598274978168L, -6.242269826651158210658094099245287635833L, -0.00194732952394737760557329759096809338509L }, + { 0.5, 2, 0.05, -0.5, 0.1619111752996693505827179422710667011729L, -1.820707394740831011982784453049470024416L, 0.1902846477536663758749416549401691954994L, 0.8097153522463336241250583450598308045006L, -1.659234181850973087825942269518598513973L, -0.2110725100569042567305917116110946846935L }, + { 0.5, 2, 0.05, 0.25, 0.1835848470423718787177487583081779818812L, -1.695078336650691841412114949055366012756L, 0.3218746711043940370184449592529980726989L, 0.6781253288956059629815550407470019273011L, -1.1335930293762498514571924436629684101L, -0.3884231572460850972272443880540950963379L }, + { 0.5, 2, 0.05, 1, 0.1765716300692055618376724509591067746702L, -1.734028648851825495277856801880003561455L, 0.4584020927196983646983517538704278540592L, 0.5415979072803016353016482461295721459408L, -0.7800085483221799657234567565864361429962L, -0.6132314213496210496449898360512151050468L }, + { 0.5, 2, 0.05, 4, 0.07125611704849486727158939803010069722299L, -2.641474610206353167989633676116143473334L, 0.8295961473591191058350015746177976544793L, 0.1704038526408808941649984253822023455207L, -0.1868162660516636254194341504841769437842L, -1.769584055447669281973206403773798709673L }, + { 0.5, 2, 0.05, 30, 4.576590121622874632575130895371787826148e-6L, -12.29455635205656301223688087145681232746L, 0.9999840962228616264473566472665571532576L, 1.590377713837355264335273344284674240573e-5L, -1.590390360477805028881394339502371816681e-5L, -11.04895392107746622240188923278841501568L }, + { 0.5, 2, -0.05, -3, 0.0116522972401220260452493356780315964923L, -4.452251930416117493176344205665701826687L, 0.004734624088379052327935505548823451317078L, 0.9952653759116209476720644944511765486829L, -5.352852945489512195361591654164601106664L, -0.004745867925273727842992985544106665896138L }, + { 0.5, 2, -0.05, -0.5, 0.1552677890261829157254418743604834157109L, -1.86260398163328395446559383287104030657L, 0.1942486110094597440911993454966736347911L, 0.8057513889905402559088006545033263652089L, -1.638616440290397162630580347606087362705L, -0.2159800344400579524535708917838999447462L }, + { 0.5, 2, -0.05, 0.25, 0.181323151227319643223971569411501005637L, -1.707474473690776936044334692444143202217L, 0.3221597500187632053101022398452834952805L, 0.6778402499812367946898977601547165047195L, -1.132707738392917036539335266553895759585L, -0.3888436383010770192104506812337850194648L }, + { 0.5, 2, -0.05, 1, 0.1809166811361984369704735047807703523638L, -1.709718678912861502344432957419012997998L, 0.4595191619005060828100913912866177786501L, 0.5404808380994939171899086087133822213499L, -0.7775746364225739497035679822722962193043L, -0.6152960947061666546825125553326396886438L }, + { 0.5, 2, -0.05, 4, 0.07478619311646186860222396651672830589823L, -2.593121995043949955032873487064120053411L, 0.8519749640523185531180141938235932211261L, 0.1480250359476814468819858061764067788739L, -0.1601981374996854657648926641497742859958L, -1.910373857715608960778425274430996387585L }, + { 0.5, 2, -0.05, 30, 4.596495628665318784510202472933638098048e-12L, -26.10572692262107364754515967239754052211L, 0.9999999999975868397949477963434007444284L, 2.413160205052203656599255571631276795675e-12L, -2.413160205055115327686884052649638212485e-12L, -26.7500839390103805297543906676252635749L }, + { 0.5, 2, 0.125, -3, 0.003423246589701594419900482066542693821116L, -5.677165882969006728455454653911608516139L, 0.0007422982779016125799840396566571464186797L, 0.9992577017220983874200159603433428535813L, -7.2057594037927936L, -0.0007425739176813640686958440054640087028603L }, + { 0.5, 2, 0.125, -0.5, 0.1672727062307847691187462917250953461827L, -1.788129826971705694097237563124980650824L, 0.1871529023412981934236701852128733036644L, 0.8128470976587018065763298147871266963356L, -1.675829336649900929736320682822740435909L, -0.2072122588825521358688092544051501334343L }, + { 0.5, 2, 0.125, 0.25, 0.1853185158707079340417764264213095687752L, -1.685679226933769823085801112302159614852L, 0.3216577267813451408584717373817422758756L, 0.6783422732186548591415282626182577241244L, -1.134267259087077031136514593991437472465L, -0.3881032906682759843935688114900746473114L }, + { 0.5, 2, 0.125, 1, 0.1734477527072901158512380207091809130925L, -1.751878862079889173325343037941231293253L, 0.4575876840435851385758078160443173339454L, 0.5424123159564148614241921839556826660546L, -0.7817867535191606685688570481121052483518L, -0.6117288362044184134725997865669490731972L }, + { 0.5, 2, 0.125, 4, 0.06863052259648018182869424109167513063265L, -2.679017907396134254294770921115472364825L, 0.8142897347074998526170434672704840171996L, 0.1857102652925001473829565327295159828004L, -0.2054390368669100216143901588846425750722L, -1.683567533247151993873221105493973708099L }, + { 0.5, 2, 0.125, 30, 4.110037028923576548359531029915430181689e-5L, -10.09949342703221363549386646401175817482L, 0.9997662143139766060126797142254146386006L, 0.0002337856860233939873202857745853613993746L, -0.0002338130181568794225181430950421601604372L, -8.361105733968211082512536218011388359181L }, + { 0.5, 2, 0.13, -3, 0.003238608780817604308032961692976858258538L, -5.732611430028139465850306873793424785492L, 0.0006870172638822821523678523743578983473903L, 0.9993129827361177178476321476256421016526L, -7.283151136680250121578402982594014422121L, -0.0006872533683874924330563580122977629701224L }, + { 0.5, 2, 0.13, -0.5, 0.167642582979750668189820404169765257474L, -1.785921048147577612224991852149429323638L, 0.186938984860475730721935780976793185975L, 0.813061015139524269278064219023206814025L, -1.676972999538390460021589330236193440741L, -0.2069491228782811202391501773295148965376L }, + { 0.5, 2, 0.13, 0.25, 0.1854352593863591939810156732461367276423L, -1.685049463854870883476988933156621950072L, 0.3216431666550400036284121905812256791238L, 0.6783568333449599963715878094187743208762L, -1.134312526011019406511383155627553356037L, -0.3880818266218485114822986459015728619994L }, + { 0.5, 2, 0.13, 1, 0.173243430462994115239555769316327097926L, -1.753057561036600643732551237372863172629L, 0.4575340801157133252464811627294380131082L, 0.5424659198842866747535188372705619868918L, -0.7819039049847524964377866345172480098276L, -0.611630016030606195189945313598941525237L }, + { 0.5, 2, 0.13, 4, 0.06845798596980971827317031124342798230274L, -2.681535065392460035884664670868056662852L, 0.813310520058397836922121976181850860999L, 0.186689479941602163077878023818149139001L, -0.2066422988649100421747385705977984967315L, -1.678308577415702260218934747197293895355L }, + { 0.5, 2, 0.13, 30, 4.537717623727485757201078903810861975631e-5L, -10.00050130537729555640686691536980221853L, 0.9997351891111566976157126505655734167769L, 0.0002648108888433023842873494344265832231021L, -0.0002648459574378944841479534482292109445275L, -8.236494613737416399497571797131769469261L }, + { 0.5, 2, -0.3, -3, 0.02257883218224408674428728479627418398383L, -3.79074244071926629082213152250139356226L, 0.01686948545248619375262787356543226767867L, 0.9831305145475138062473721264345677323213L, -4.082248883631195877498228653947753613791L, -0.01701339598251505476486047576402773046538L }, + { 0.5, 2, -0.3, -0.5, 0.1407963398662055892525447031129239136436L, -1.960440830864504084642685787505986146457L, 0.2032324564965004859201866072448043529658L, 0.7967675435034995140798133927551956470342L, -1.59340484917992907506999272775344169187L, -0.227192307096227335984562664303698179324L }, + { 0.5, 2, -0.3, 0.25, 0.1759055476293873263072756029607300024377L, -1.737808089194809348779680153454023765511L, 0.3228523505680641335944189777181274753936L, 0.6771476494319358664055810222818725246064L, -1.130560179254535433596138358277503983987L, -0.3898659361532131300854343785533963854326L }, + { 0.5, 2, -0.3, 1, 0.1927795672983714006835847556170475107805L, -1.646207881164722508291043440915108869867L, 0.4624807071474759097728591303018553747857L, 0.5375192928525240902271408696981446252143L, -0.7711504371754495264611890107350469005594L, -0.6207906259453672883184693821818092690651L }, + { 0.5, 2, -0.3, 4, 0.08096038522226132802083470239305868591372L, -2.513795315274196404958072065994092601609L, 0.9197803753929363396011351962134600433416L, 0.08021962460706366039886480378653995665844L, -0.08362035983676079918029272072320952064178L, -2.522987098194140438626166073892907539789L }, + { 0.5, 2, 0.3, -3, 8.059704992362896995119811007517873715506e-5L, -9.426048510565349992020177118882759848275L, 6.402576600575467769883962668923564384817e-6L, 0.9999935974233994245322301160373310764356L, -11.95881005477788658284444167501990557427L, -6.402597097156518221636218314392215847589e-6L }, + { 0.5, 2, 0.3, -0.5, 0.1812503225698147173179925232617782002554L, -1.707876205553407075855624216073082303607L, 0.1792492973447090255388108136039045090739L, 0.8207507026552909744611891863960954909261L, -1.71897771948381971532460359928267808018L, -0.1975358664871359474170686040679475679825L }, + { 0.5, 2, 0.3, 0.25, 0.1894940236404238306364162765021104566183L, -1.663397792474526321406233515749467943149L, 0.3211407214070816380060702435413827912124L, 0.6788592785929183619939297564586172087876L, -1.135875867468771317232919576205998713349L, -0.3873414209244915061172579286398172204468L }, + { 0.5, 2, 0.3, 1, 0.1665725681295165996385848051549330081597L, -1.79232421989259561589974139903168544228L, 0.4557606462631520925261687698121302605051L, 0.5442393537368479074738312301878697394949L, -0.7857875058209371193944533846408445480066L, -0.6083661404058928633051700276338946733242L }, + { 0.5, 2, 0.3, 4, 0.06286582155301994052150842939996142150539L, -2.766752640494270141607384901313505234844L, 0.7827334860527895458826706624083477548845L, 0.2172665139472104541173293375916522451155L, -0.2449630163437000697794851493209835780066L, -1.526630504159534364420939159431703465001L }, + { 0.5, 2, 0.3, 30, 0.0003273580834484173214067780354743419595302L, -8.024455929786768161735886630208929083997L, 0.996441826975464168416324204254454132273L, 0.003558173024535831583675795745545867727038L, -0.003564518378555365734021672304781503586974L, -5.638508061203108378546491617990135356678L }, + { 0.5, 2, -0.7, -3, 0.03065313768572185299711007995410933359305L, -3.485020250364091045199580492854250520331L, 0.04351602431599833910256362151992301784644L, 0.9564839756840016608974363784800769821536L, -3.134626033640518737066247611835648581181L, -0.04449124331691377264478415694517319717866L }, + { 0.5, 2, -0.7, -0.5, 0.1224780659178903405081370654853399188579L, -2.099823318760895943245804892832088841208L, 0.2153927028350637145526002290527971174053L, 0.7846072971649362854473997709472028825947L, -1.535292392108238402871228881400685712766L, -0.2425719447593194681492008503702743807618L }, + { 0.5, 2, -0.7, 0.25, 0.167880318057502716578366363708646221341L, -1.784503945965868755929679208375181402125L, 0.3239046900990341450574093655119228912695L, 0.6760953099009658549425906344880771087305L, -1.127305972826224354630816284263039890692L, -0.3914212218997065966819085426931104755151L }, + { 0.5, 2, -0.7, 1, 0.2153903268677851220270401478970179674622L, -1.535303423031532923288919312610236116888L, 0.4678015340166633970763864726217391180225L, 0.5321984659833366029236135273782608819775L, -0.7597111456226778590349280458738389255287L, -0.6307388028490597649920765035511791995061L }, + { 0.5, 2, 0.9, -0.5, 0.2530658734606315277479524224191294277061L, -1.374105454735238240559972359301174861478L, 0.1432650870655338913562264826490839940841L, 0.8567349129344661086437735173509160059159L, -1.943058609103825017156942536489297314252L, -0.1546267280300151672391999106680229972773L }, + { 0.5, 2, 0.9, 0.25, 0.2053591919323515073717253479034976085078L, -1.582994677310361952163458363526016508605L, 0.3192433025902280900330646960560710881853L, 0.6807566974097719099669353039439289118147L, -1.141801762863611938061786002862994555557L, -0.3845503091956913171670906007042623815087L }, + { 0.5, 2, 0.9, 1, 0.1466505062859526950907301644183770990309L, -1.919703031222647272046695408322428608797L, 0.450172165222532504053136402084222020012L, 0.549827834777467495946863597915777979988L, -0.7981251800030224208698105134477822631855L, -0.59815007743634220782812536170051942518L }, + { 0.5, 2, 0.9, 4, 0.04785650470556758324311619753289467322999L, -3.039548230811817410079180781740187256459L, 0.7049646710996157775249081280210048415821L, 0.2950353289003842224750918719789951584179L, -0.3496075893409509957988211565573886561168L, -1.220660170828498612057196401644921029177L }, + { 0.5, 2, 0.9, 30, 0.00173336156054610425360097856463650573284L, -6.357692657265018351117128996627085044189L, 0.9491997737464478260487518121425344402256L, 0.05080022625355217395124818785746555977442L, -0.05213599278291702573343607369694041681873L, -2.979854470597536862470282274060356423935L }, + { 0.5, 2, 2, 0.25, 0.2426040447120236162512332757790964753174L, -1.416324610261525447276701173471350332115L, 0.3151518986722023896141958639310454513001L, 0.6848481013277976103858041360689545486999L, -1.154700538379251529018297561003914911295L, -0.3785582151964799813383119465458937224332L }, + { 0.5, 2, 2, 1, 0.1202910057881380141120158758405166072322L, -2.117841423649917915116679819556664070255L, 0.4419773772406739579865398211578859656963L, 0.5580226227593260420134601788421140343037L, -0.816496580927726032732428024901963797322L, -0.5833557748432169237414957712355888063703L }, + { 0.5, 2, 2, 4, 0.03269059698142522532567882310927362868032L, -3.420667796515388102411015891441721522761L, 0.6241250557782609608783975539903324977094L, 0.3758749442217390391216024460096675022906L, -0.4714045207910316829338962414032326928566L, -0.9784987860692273495936509450915494075046L }, + { 0.5, 2, 2, 30, 0.002476748424901686938892812180090879881433L, -6.000808698065031280290777657414978377417L, 0.8343757053045658053987031428930575316221L, 0.1656242946954341946012968571069424683779L, -0.1810714920850370618723100630056445419812L, -1.798033340693379056209812357545670028983L }, + { 0.5, 2, -2, -3, 0.02825439601041371508207652025113584539374L, -3.56650622250772491930639438396616110654L, 0.1198732501037620330747323068407013287011L, 0.8801267498962379669252676931592986712989L, -2.121320343559642573202533086314547117855L, -0.1276893479087839697383530966932791042295L }, + { 0.5, 2, -2, -0.5, 0.08595474576918094590790089441973243385282L, -2.453934333213013012927536906396962930683L, 0.2431167344342142108048623204999460644609L, 0.7568832655657857891951376795000539355391L, -1.41421356237309504880168872420969807857L, -0.278546244085973154673030816903783655471L }, + { 0.5, 2, -2, 0.25, 0.1462039162679206421758576956058319087134L, -1.922752944966945035504966500978731937483L, 0.3269218953517579112720873007606364511319L, 0.6730781046482420887279126992393635488681L, -1.11803398874989484820458683436563811772L, -0.3958939016161662776988365565973799879904L }, + { 0.5, 2, -2, 1, 0.3486522152763511488462237401551361482114L, -1.053680371466520179109460422833937323323L, 0.4930686913952397878459573806151774752358L, 0.5069313086047602121540426193848225247642L, -0.7071067811865475244008443621048490392848L, -0.6793797705566755265109604011331633483493L }, +}; + +// mu, sigma, xi, p, quantile(p), quantile(complement(p)) +const long double quantile_data[][6] = { + { 0.5, 2, 0, 1e-300, -12.5756298398083137183963288862606664562L, 1382.051055796427410360676689140201005817L }, + { 0.5, 2, 0, 1e-10, -5.773235076484002964477948553537579912083L, 46.55170185978091360749126779238166592248L }, + { 0.5, 2, 0, 0.001, -3.365289467832130976369869903602085211958L, 14.31451014104743295818180028793985985517L }, + { 0.5, 2, 0, 0.25, -0.1532685199565619648095859264510141972425L, 2.991798647414476396761561605799278127912L }, + { 0.5, 2, 0, 0.5, 1.233025841163328654024878316465338938909L, 1.233025841163328654024878316465338938909L }, + { 0.5, 2, 0, 0.75, 2.991798647414476396761561605799278127912L, -0.1532685199565619648095859264510141972425L }, + { 0.5, 2, 0, 0.999, 14.31451014104743122259023912057020651812L, -3.365289467832130725242925488683664408872L }, + { 0.5, 2, 0, 0.99999999, 37.34136146785521232696297750383352136488L, -5.326947973310027150864118532956906440812L }, + { 0.5, 2, 1e-9, 1e-300, -12.57562979706528988462358895414425749839L, 1382.05153296736722607948550779053034854L }, + { 0.5, 2, 1e-9, 1e-10, -5.77323506664563339355703516054520922164L, 46.55170238997072872238572308430781231717L }, + { 0.5, 2, 1e-9, 0.001, -3.365289464097015311242587306799330065802L, 14.31451018875760567730617218391453783919L }, + { 0.5, 2, 1e-9, 0.25, -0.1532685198498720250296428667953639746252L, 2.991798648966741522220370573709271182235L }, + { 0.5, 2, 1e-9, 0.5, 1.23302584129766037499459114420234923096L, 1.23302584129766037499459114420234923096L }, + { 0.5, 2, 1e-9, 0.75, 2.991798648966741522220370573709271182235L, -0.1532685198498720250296428667953639746252L }, + { 0.5, 2, 1e-9, 0.999, 14.31451018875760394171459902837123186796L, -3.365289464097015060115643377220075461646L }, + { 0.5, 2, 1e-9, 0.99999999, 37.34136180717669311176557459591262533638L, -5.326947964821696488192169976806530503024L }, + { 0.5, 2, -1e-9, 1e-300, -12.57562988255133773846638851920848108813L, 1382.050578625707339929835571189279739489L }, + { 0.5, 2, -1e-9, 1e-10, -5.773235086322372555971663663944269955237L, 46.55170132959110663131118162134251665025L }, + { 0.5, 2, -1e-9, 0.001, -3.36528947156724664630958691735551269533L, 14.31451009333726045875498283928558450044L }, + { 0.5, 2, -1e-9, 0.25, -0.1532685200632519046127613791272932954428L, 2.991798645862211272592063350702808937606L }, + { 0.5, 2, -1e-9, 0.5, 1.233025841028996933087988362977232413154L, 1.233025841028996933087988362977232413154L }, + { 0.5, 2, -1e-9, 0.75, 2.991798645862211272592063350702808937606L, -0.1532685200632519046127613791272932954428L }, + { 0.5, 2, -1e-9, 0.999, 14.31451009333725872316343366008950099219L, -3.365289471567246395182642017097924755325L }, + { 0.5, 2, -1e-9, 0.99999999, 37.34136112853373570918213062556195007658L, -5.326947981798357830023087488999341876338L }, + { 0.5, 2, 0.1, 1e-300, -9.098492146372804608864205359432952262754L, 2.00000000000000755311483837239514148787e+31L }, + { 0.5, 2, 0.1, 1e-10, -4.884676695574487760921707924416759610688L, 180.4999999990000148431395292428273467397L }, + { 0.5, 2, 0.1, 0.001, -3.014721011455450203818984787388703685023L, 20.40325025512648903682365669916080180916L }, + { 0.5, 2, 0.1, 0.25, -0.1427147455436507519235484201451701256752L, 3.153677616531757077682086383344611767531L }, + { 0.5, 2, 0.1, 0.5, 1.246624642471412121321222799267411548132L, 1.246624642471412121321222799267411548132L }, + { 0.5, 2, 0.1, 0.75, 3.153677616531757077682086383344611767531L, -0.1427147455436507519235484201451701256752L }, + { 0.5, 2, 0.1, 0.999, 20.40325025512648557403643640179963007554L, -3.014721011455449996824097776356218705496L }, + { 0.5, 2, 0.1, 0.99999999, 106.6914687695347467102837941353157503273L, -4.554879256576813897134283068908618375508L }, + { 0.5, 2, -0.1, 1e-300, -17.95596288816075050927756797340028646439L, 20.49999999999999888977697537482352120613L }, + { 0.5, 2, -0.1, 1e-10, -6.868535862554625238918763399218445721122L, 18.49999999998999924915069534853500923024L }, + { 0.5, 2, -0.1, 0.001, -3.764072223342769747986519854910516424594L, 10.47575391872718881130071752413829357759L }, + { 0.5, 2, -0.1, 0.25, -0.164054630693307353519068190091098963657L, 2.842822795884769008582543984252495089731L }, + { 0.5, 2, -0.1, 0.5, 1.219755290644205266101074759104846547934L, 1.219755290644205266101074759104846547934L }, + { 0.5, 2, -0.1, 0.75, 2.842822795884769008582543984252495089731L, -0.164054630693307353519068190091098963657L }, + { 0.5, 2, -0.1, 0.999, 10.47575391872718794140087223803001241633L, -3.764072223342769443318404029362227341775L }, + { 0.5, 2, -0.1, 0.99999999, 17.330213611900137860993975842357316173L, -6.264588046304391117894908529045757023257L }, + { 0.5, 2, 0.5, 1e-300, -3.347808106759349922301748343690987107433L, 3.999999999999999949881816329582481570545e+150L }, + { 0.5, 2, 0.5, 1e-10, -2.666410670026779358283931425752205794977L, 3.999964999899999927132688700516086050025e+5L }, + { 0.5, 2, 0.5, 0.001, -1.978081067593499220751919467145709898194L, 122.9594744018793442987725313740319406105L }, + { 0.5, 2, 0.5, 0.25, -0.1027127988479238291139886635884415211394L, 3.957677382973556235894851336295208920213L }, + { 0.5, 2, 0.5, 0.5, 1.304489635145799179431213144380886890267L, 1.304489635145799179431213144380886890267L }, + { 0.5, 2, 0.5, 0.75, 3.957677382973556235894851336295208920213L, -0.1027127988479238291139886635884415211394L }, + { 0.5, 2, 0.5, 0.999, 122.9594744018792894282733809833461204596L, -1.978081067593499125203206681530880286184L }, + { 0.5, 2, 0.5, 0.99999999, 39996.49979950481432923013743236994190557L, -2.568018796303575793907240515617076686281L }, + { 0.5, 2, -0.5, 1e-300, -100.6304353951386395706955318473740513189L, 4.5L }, + { 0.5, 2, -0.5, 1e-10, -14.69410364875232481508473000078024131598L, 4.499959999999998999999271301886968726567L }, + { 0.5, 2, -0.5, 0.001, -6.013043539513863941419591712727294057535L, 4.373477253676121397844184180029057779858L }, + { 0.5, 2, -0.5, 0.25, -0.2096400900618987640462773058387985509895L, 2.354559914789393416226967493531125213178L }, + { 0.5, 2, -0.5, 0.5, 1.169781555369208974587341420419195809478L, 1.169781555369208974587341420419195809478L }, + { 0.5, 2, -0.5, 0.75, 2.354559914789393416226967493531125213178L, -0.2096400900618987640462773058387985509895L }, + { 0.5, 2, -0.5, 0.999, 4.373477253676121342946231476168160149408L, -6.013043539513863281392466567948484653939L }, + { 0.5, 2, -0.5, 0.99999999, 4.499599999997995048133242721454662269415L, -12.66772820797328732427871819361313767496L }, + { 0.5, 2, 1.5, 1e-200, -0.8331984152353507655021985957119811372087L, 1.333333333333333369132808535167814942226e+300L }, + { 0.5, 2, 1.5, 1e-10, -0.8212658917939733770375279133898950683426L, 1.333333333233332427135188704159432109032e+15L }, + { 0.5, 2, 1.5, 0.001, -0.7598932228613076226472790168569414439187L, 42131.24470703879474650865711481124332923L }, + { 0.5, 2, 1.5, 0.25, -0.01645843373705890274069126157240972584411L, 7.807776763674368551664252384128167313212L }, + { 0.5, 2, 1.5, 0.5, 1.477137790209089876562044068057611286388L, 1.477137790209089876562044068057611286388L }, + { 0.5, 2, 1.5, 0.75, 7.807776763674368551664252384128167313212L, -0.01645843373705890274069126157240972584411L }, + { 0.5, 2, 1.5, 0.999, 42131.24470703873990344933112866136953218L, -0.759893222861307608815186111606538485309L }, + { 0.5, 2, 1.5, 0.99999999, 1.333333313282981533418812268220497771132e+12L, -0.81646857107861184640351984211780469705L }, + { 0.5, 2, -1.5, 1e-300, -24205.34400274864829699997496298744631796L, 1.833333333333333333333333333333333333333L }, + { 0.5, 2, -1.5, 1e-10, -145.4868564499990778008788306250011297411L, 1.833333333333331999999999899999927127689L }, + { 0.5, 2, -1.5, 0.001, -22.37384400274864818889387098884254541468L, 1.83329113798336568488968167869093474847L }, + { 0.5, 2, -1.5, 0.25, -0.3429824999189946874784489084545169083132L, 1.627598449988495743380535845506960621981L }, + { 0.5, 2, -1.5, 0.5, 1.06388949148514695405338186368306265767L, 1.06388949148514695405338186368306265767L }, + { 0.5, 2, -1.5, 0.75, 1.627598449988495743380535845506960621981L, -0.3429824999189946874784489084545169083132L }, + { 0.5, 2, -1.5, 0.999, 1.833291137983365684834756258697668850538L, -22.3738440027486436295880129985934156298L }, + { 0.5, 2, -1.5, 0.99999999, 1.833333333331999999979950481231931414678L, -103.5804134439192277068503656272328464734L }, +}; + +// xi, mean, variance, skewness, kurtosis_excess, mode, median, entropy, for mu = 0.5, sigma = 2 +const long double moment_data[][8] = { + { -3, -2.833333333333333333333333333333333333333L, 304.0L, -19.58485896516058864744704924603960276765L, 1002.531855955678670360110803324099722992L, 1.166666666666666666666666666666666666667L, 0.9446502320073803468540976115921796372292L, 0.5387158507568795882042079412933717059912L }, + { -1, 0.5L, 4.0L, -2.0L, 6.0L, 2.5L, 1.113705638880109381165535757083646863849L, 1.693147180559945309417232121458176568076L }, + { -0.7, 0.7610321918477416120986094092934992813022L, 3.400367327861812827006952124923533645419L, -1.160392728202874330200516009502937526774L, 1.698453497221796334383449480790908683616L, 2.127109656428759269757252685159577858255L, 1.146544828562143746224956954046140776713L, 1.866311880030405193232710602402922212832L }, + { -0.3, 1.183528691291485424098573010753824634699L, 3.913853269298499701383800112689952627786L, -0.06874209942096706080288294973113119974362L, -0.289399162107949706627675781165848166877L, 1.176510388062401562284821150299102160269L, 1.194163629688140530678864252817002609567L, 2.097198145991018318250171797995864498666L }, + { -0.13, 1.424131542591410291946229744699725656592L, 4.952896477265949926405670179598615678555L, 0.5136871089846431160319333561076186998656L, 0.2858096471132615891399975041498497692222L, 0.7760180707317563468194075597259850477742L, 1.215836798220134140775267400701285930257L, 2.195324809024278895581545154437864191538L }, + { -0.12, 1.440163573101605578076943181400254978206L, 5.044061098077296378047764932023288965371L, 0.5540021135248004886617811879635306797331L, 0.3716301414765079931716820994962564984678L, 0.7537157672471984893263687435540211023558L, 1.217139783552452963126246626368292989835L, 2.201096965673294229314315246122693198937L }, + { -0.01, 1.634829761698792572529680260212485936139L, 6.412196648010445039890101186741231605832L, 1.081073759807266273858956528275016526025L, 2.125445886586571376929642583102013399481L, 0.5200996616483335806178703717697134542996L, 1.231684163594858762816687075275816643953L, 2.264590688812462841297521942887004858016L }, + { -1e-6, 1.654429351692890021456710465409222852749L, 6.57971905513086479986738955938859978589L, 1.139541132804515740836640856317379221379L, 2.399971082464274566936674468358194583232L, 0.5000020000009999996665756627996313747545L, 1.2330257068316241121591228516200509969L, 2.27036226824581326849090972512720938759L }, + { -1e-12, 1.654431329801087609222370050051952694087L, 6.579736267375693436975510933807444342468L, 1.139547099398682045080138208052414725792L, 2.399999999971082311659352991890256472674L, 0.5000000000020000000000009999597732949251L, 1.233025841163194322303925031497233263303L, 2.27036284546090095435884267869158222909L }, + { 1e-12, 1.654431329805043833203681940193961353626L, 6.579736267410118054803904146074174922329L, 1.139547099410615269905472390300595465804L, 2.400000000028917688340952240175754611781L, 0.4999999999980000000000010000402267050748L, 1.233025841163462985745831634256318863418L, 2.270362845462055385688645744389575769146L }, + { 1e-6, 1.654433307916871337273665367728207595283L, 6.579753479748693405330144223354347365232L, 1.139553066029341147467124756032196087804L, 2.400028917840957499080165724640479681121L, 0.4999980000010000003334233370193610725493L, 1.233025975495066018764882685290321447262L, 2.270363422677143071556578697953948610645L }, + { 0.01, 1.674395928821558387244769681014319615333L, 6.756655174188248780404314727298438625328L, 1.200478519762556269331078472919269293561L, 2.705134978371804552738435844656986437332L, 0.48010032835166642105707054463640603916L, 1.234370801021428012906777726473035063745L, 2.276135002110493498749966480194153140219L }, + { 0.12, 1.921796457795325141504033129127197829462L, 9.567328041748161760308999329794170737815L, 2.129206380223521923200713361070088227233L, 10.4127143737989297783691126735759442322L, 0.2748768734662167365566969712751014132572L, 1.249384593879594267985829474900117552262L, 2.339628725249662110733173176958464799298L }, + { 0.13, 1.94721813462737892879360520664898861468L, 9.932951514939585997564689348668914277938L, 2.251782857706110263308481491192789580094L, 11.99260069837783378364765088349075208532L, 0.2574963206206130134352559080015524785188L, 1.250769653649645211837358227649886079707L, 2.345400881898677444465943268643293806698L }, + { 0.2, 2.142297137253033769930034117143080376893L, 13.37614224919152628643986765706565611003L, 3.535071604621394887878195503476186178003L, 45.09151212581577318204971476397941468765L, 0.1419250400262719916570182230117110839103L, 1.260560851390051221516189461818868534605L, 2.385805978441784748553427843037065714187L }, + { 0.24, 2.269461454069749995037297391358303819376L, 16.39491346167414058557668895449269328167L, 5.024462495135007919842276150900259265379L, 309.6084113035859890772945909440607406133L, 0.08069406600135331847379173643045407240867L, 1.266231911965916800406745667483278841138L, 2.408894605037846051442602142376350599479L }, +}; + +template +bool representable(long double v) +{ + return std::abs(v) <= (std::numeric_limits::max)() && (v == 0 || std::abs(v) >= (std::numeric_limits::min)()); +} + +// The condition numbers |x f'(x)/f(x)| of the row's functions, from the reference values: with +// s = 1 + xi (x - mu)/sigma and t = -logcdf, d(log t)/dx = -1/(sigma s). +struct cdf_conditions +{ + long double pdf, logpdf, cdf, ccdf, logcdf, logccdf; +}; + +inline cdf_conditions conditions(long double const (&r)[10]) +{ + long double mu = r[0], sigma = r[1], xi = r[2], x = r[3]; + long double s = 1 + xi * (x - mu) / sigma, t = -r[8], P = r[6], Q = r[7]; + long double g = std::abs(x / (sigma * s)); + cdf_conditions c; + c.pdf = g * std::abs(1 + xi - t); + c.logpdf = c.pdf / std::abs(r[5]); + c.cdf = g * t; + c.ccdf = c.cdf * P / Q; + c.logcdf = g; + c.logccdf = c.ccdf / std::abs(r[9]); + return c; +} + +template +void test_cdf() +{ + // 16 ulps times each function's condition number, which is large deep in the tails (38 for the cdf at + // xi = 0.3, x = -3, where t = 12); the condition numbers come from the long double references. + auto ulps = [](long double kappa) { return static_cast(16 * (std::max)(1.0L, kappa)); }; + for (auto const& r : cdf_data) + { + generalized_extreme_value_distribution d(static_cast(r[0]), static_cast(r[1]), static_cast(r[2])); + T x = static_cast(r[3]); + if (static_cast(x) != r[3] || static_cast(d.shape()) != r[2]) + { + continue; + } + cdf_conditions k = conditions(r); + if (representable(r[4])) { CHECK_ULP_CLOSE(static_cast(r[4]), pdf(d, x), ulps(k.pdf)); } + CHECK_ULP_CLOSE(static_cast(r[5]), logpdf(d, x), ulps(k.logpdf)); + if (representable(r[6])) { CHECK_ULP_CLOSE(static_cast(r[6]), cdf(d, x), ulps(k.cdf)); } + if (representable(r[7])) { CHECK_ULP_CLOSE(static_cast(r[7]), cdf(complement(d, x)), ulps(k.ccdf)); } + if (representable(r[8])) { CHECK_ULP_CLOSE(static_cast(r[8]), logcdf(d, x), ulps(k.logcdf)); } + if (representable(r[9])) { CHECK_ULP_CLOSE(static_cast(r[9]), logcdf(complement(d, x)), ulps(k.logccdf)); } + } +} + +template +void test_quantile() +{ + for (auto const& r : quantile_data) + { + generalized_extreme_value_distribution d(static_cast(r[0]), static_cast(r[1]), static_cast(r[2])); + T p = static_cast(r[3]); + if (static_cast(p) != r[3] || static_cast(d.shape()) != r[2]) + { + continue; + } + // mu + sigma q cancels where the quantile is near 0, multiplying the error by (|mu| + |x - mu|)/|x|: + for (int c = 0; c < 2; ++c) + { + long double x = r[4 + c]; + if (!representable(x)) + { + continue; + } + double cancellation = static_cast((std::abs(r[0]) + std::abs(x - r[0])) / std::abs(x)); + double scale = (std::max)(1.0, cancellation); +#if defined(__MINGW32__) || defined(__CYGWIN__) + // MinGW's long double pow loses about |log x| ulps for huge results: + if (std::numeric_limits::digits > 53) + { + scale = (std::max)(scale, static_cast(std::abs(std::log(std::abs(x))))); + } +#endif + unsigned ulps = static_cast(16 * scale); + CHECK_ULP_CLOSE(static_cast(x), c ? quantile(complement(d, p)) : quantile(d, p), ulps); + } + } +} + +template +void test_moments() +{ + const T eps = std::numeric_limits::epsilon(); + for (auto const& r : moment_data) + { + T xi = static_cast(r[0]); + if (static_cast(xi) != r[0]) + { + continue; + } + generalized_extreme_value_distribution d(T(0.5), T(2), xi); + CHECK_ULP_CLOSE(static_cast(r[1]), mean(d), 8); + // Gamma(1 - xi)^2 doubles tgamma's error: + CHECK_ULP_CLOSE(static_cast(r[2]), variance(d), 32); + CHECK_ULP_CLOSE(static_cast(r[2]), standard_deviation(d) * standard_deviation(d), 32); + CHECK_MOLLIFIED_CLOSE(static_cast(r[3]), skewness(d), 64 * eps); + CHECK_MOLLIFIED_CLOSE(static_cast(r[4]), kurtosis_excess(d), 512 * eps); + CHECK_MOLLIFIED_CLOSE(static_cast(r[4] + 3), kurtosis(d), 512 * eps); + CHECK_ULP_CLOSE(static_cast(r[5]), mode(d), 8); + CHECK_ULP_CLOSE(static_cast(r[6]), median(d), 8); + CHECK_ULP_CLOSE(static_cast(r[7]), entropy(d), 8); + } +} + +// The support, its endpoints, the moments' domains, and xi = 0 against extreme_value_distribution. +template +void test_edges() +{ + using boost::math::extreme_value_distribution; + const T inf = std::numeric_limits::infinity(); + generalized_extreme_value_distribution frechet(T(0.5), T(2), T(0.25)); + generalized_extreme_value_distribution weibull(T(0.5), T(2), T(-0.25)); + // mu - sigma/xi: + CHECK_EQUAL(support(frechet).first, T(-7.5)); + CHECK_EQUAL(support(weibull).second, T(8.5)); + CHECK_EQUAL(cdf(frechet, T(-7.5)), T(0)); + CHECK_EQUAL(cdf(frechet, T(-100)), T(0)); + CHECK_EQUAL(cdf(complement(frechet, T(-100))), T(1)); + CHECK_EQUAL(pdf(frechet, T(-7.5)), T(0)); + CHECK_EQUAL(logcdf(frechet, T(-100)), -inf); + CHECK_EQUAL(cdf(weibull, T(8.5)), T(1)); + CHECK_EQUAL(cdf(weibull, T(100)), T(1)); + CHECK_EQUAL(cdf(complement(weibull, T(100))), T(0)); + CHECK_EQUAL(pdf(weibull, T(8.5)), T(0)); + CHECK_EQUAL(logpdf(weibull, T(100)), -inf); + CHECK_EQUAL(quantile(frechet, T(0)), T(-7.5)); + CHECK_EQUAL(quantile(complement(weibull, T(0))), T(8.5)); + CHECK_EQUAL(cdf(frechet, inf), T(1)); + CHECK_EQUAL(cdf(weibull, -inf), T(0)); + CHECK_EQUAL(pdf(frechet, inf), T(0)); + // Just past an end point, where (x - mu)/sigma rounds onto it: + generalized_extreme_value_distribution unit_end(T(-1), T(1), T(-1)); + CHECK_EQUAL(pdf(unit_end, T(1e-17)), T(0)); + CHECK_EQUAL(logpdf(unit_end, T(1e-17)), -inf); + CHECK_EQUAL(cdf(unit_end, T(1e-17)), T(1)); + // log(1 - cdf) = -L once t = exp(-L) underflows: + generalized_extreme_value_distribution gumbel(T(0), T(1), T(0)); + if (std::numeric_limits::max_exponent10 < 300) + { + CHECK_ULP_CLOSE(T(-100), logcdf(complement(gumbel, T(100))), 1); + } + else + { + CHECK_ULP_CLOSE(T(-750), logcdf(complement(gumbel, T(750))), 1); + } + // At the upper endpoint the density is t^(1 + xi)/sigma with t = 0: + CHECK_EQUAL(pdf(generalized_extreme_value_distribution(T(0), T(2), T(-1)), T(2)), T(0.5)); + CHECK_THROW(pdf(generalized_extreme_value_distribution(T(0), T(2), T(-2)), T(1)), std::overflow_error); + // Unbounded ends: + CHECK_THROW(quantile(frechet, T(1)), std::overflow_error); + CHECK_THROW(quantile(weibull, T(0)), std::overflow_error); + // Moments that do not exist: + generalized_extreme_value_distribution heavy(T(0), T(1), T(1)); + // The mean and variance are infinite there, the skewness and kurtosis undefined: + CHECK_EQUAL(mean(heavy), inf); + CHECK_EQUAL(variance(generalized_extreme_value_distribution(T(0), T(1), T(0.5))), inf); + CHECK_EQUAL(standard_deviation(heavy), inf); + CHECK_THROW(skewness(generalized_extreme_value_distribution(T(0), T(1), T(1) / 3)), std::domain_error); + CHECK_THROW(kurtosis(generalized_extreme_value_distribution(T(0), T(1), T(0.25))), std::domain_error); + CHECK_THROW(generalized_extreme_value_distribution(T(0), T(0), T(0)), std::domain_error); + CHECK_THROW(generalized_extreme_value_distribution(T(0), T(1), inf), std::domain_error); + check_out_of_range>(T(0.5), T(2), T(0.25)); + check_out_of_range>(T(0.5), T(2), T(-0.25)); + + // xi = 0 is the Gumbel distribution, and tiny xi tends to it; two implementations, deep in the tail at x = 30: + generalized_extreme_value_distribution g(T(0.5), T(2), T(0)); + extreme_value_distribution e(T(0.5), T(2)); + for (T x : { T(-3), T(0.5), T(4), T(30) }) + { + CHECK_ULP_CLOSE(pdf(e, x), pdf(g, x), 16); + CHECK_ULP_CLOSE(cdf(e, x), cdf(g, x), 4); + CHECK_ULP_CLOSE(cdf(complement(e, x)), cdf(complement(g, x)), 4); + } + for (T p : { T(0.001), T(0.5), T(0.999) }) + { + CHECK_ULP_CLOSE(quantile(e, p), quantile(g, p), 4); + } + CHECK_ULP_CLOSE(skewness(e), skewness(g), 0); + CHECK_ULP_CLOSE(kurtosis(e), kurtosis(g), 0); + generalized_extreme_value_distribution tiny(T(0.5), T(2), std::numeric_limits::denorm_min()); + CHECK_ULP_CLOSE(mean(e), mean(tiny), 1); + CHECK_ULP_CLOSE(cdf(e, T(4)), cdf(tiny, T(4)), 4); + CHECK_ULP_CLOSE(quantile(e, T(0.25)), quantile(tiny, T(0.25)), 4); +} + +#ifdef BOOST_MATH_TEST_GEV_MULTIPRECISION +// The moment series in many digits, against the closed forms in more, and the quantile's round trip +// (away from p -> 0, where cdf(quantile(p)) multiplies the error by -log p). +void test_multiprecision() +{ + using T = boost::multiprecision::cpp_bin_float_50; + using U = boost::multiprecision::cpp_bin_float_100; + for (const char* s : { "-0.1", "-1e-10", "1e-10", "0.0625", "0.15", "0.2" }) + { + T xi(s); + U k(xi); + generalized_extreme_value_distribution d(T(0), T(1), xi); + U g1 = boost::math::tgamma(1 - k), g2 = boost::math::tgamma(1 - 2 * k), g3 = boost::math::tgamma(1 - 3 * k), g4 = boost::math::tgamma(1 - 4 * k); + U var = (g2 - g1 * g1) / (k * k); + U skew = (k > 0 ? 1 : -1) * (g3 - 3 * g1 * g2 + 2 * g1 * g1 * g1) / pow(g2 - g1 * g1, U(1.5)); + U kex = (g4 - 4 * g1 * g3 + 6 * g2 * g1 * g1 - 3 * g1 * g1 * g1 * g1) / ((g2 - g1 * g1) * (g2 - g1 * g1)) - 3; + // tgamma itself is ~20 ulps at 50 digits, and the variance has Gamma(1 - xi)^2: + CHECK_ULP_CLOSE(T(var), variance(d), 128); + CHECK_MOLLIFIED_CLOSE(T(skew), skewness(d), T(64 * std::numeric_limits::epsilon())); + CHECK_MOLLIFIED_CLOSE(T(kex), kurtosis_excess(d), T(512 * std::numeric_limits::epsilon())); + for (T p : { T(0.001), T(0.3), T(0.9) }) + { + CHECK_ULP_CLOSE(p, cdf(d, quantile(d, p)), 64); + } + } +} +#endif + +int main() +{ + test_cdf(); + test_cdf(); + test_cdf(); + test_quantile(); + test_quantile(); + test_quantile(); + test_moments(); + test_moments(); + test_moments(); + test_edges(); + test_edges(); + test_edges(); +#ifdef BOOST_MATH_TEST_GEV_MULTIPRECISION + test_multiprecision(); +#endif + return boost::math::test::report_errors(); +} diff --git a/test/test_negative_binomial.cpp b/test/test_negative_binomial.cpp index 271e4f0cee..827dff4da5 100644 --- a/test/test_negative_binomial.cpp +++ b/test/test_negative_binomial.cpp @@ -31,32 +31,25 @@ #endif #include -#include "../include_private/boost/math/tools/test.hpp" +#include "math_unit_test.hpp" #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS #include // for real_concept -using ::boost::math::concepts::real_concept; #endif -#include // for negative_binomial_distribution -using boost::math::negative_binomial_distribution; - -#include - using boost::math::lgamma; // log gamma - -#define BOOST_TEST_MAIN -#include // for test_main -#include // for BOOST_CHECK_CLOSE +#include +#include #include "table_type.hpp" #include "test_out_of_range.hpp" -#include -using std::cout; -using std::endl; -using std::setprecision; -using std::showpoint; +#include #include -using std::numeric_limits; +#include + +using boost::math::negative_binomial_distribution; + +// The Boost.Test checks these replaced took tolerances in percent; the tolerances below are +// fractions, the same values divided by 100. template void test_spot( // Test a single spot value against 'known good' values. @@ -67,124 +60,50 @@ void test_spot( // Test a single spot value against 'known good' values. RealType Q, // Complement of CDF. RealType tol) // Test tolerance. { - boost::math::negative_binomial_distribution bn(N, p); - BOOST_CHECK_EQUAL(N, bn.successes()); - BOOST_CHECK_EQUAL(p, bn.success_fraction()); - BOOST_CHECK_CLOSE( - cdf(bn, k), P, tol); + negative_binomial_distribution bn(N, p); + CHECK_EQUAL(N, bn.successes()); + CHECK_EQUAL(p, bn.success_fraction()); + CHECK_RELATIVE_ERROR(P, cdf(bn, k), tol); if((P < 0.99) && (Q < 0.99)) { // We can only check this if P is not too close to 1, // so that we can guarantee that Q is free of error: - // - BOOST_CHECK_CLOSE( - cdf(complement(bn, k)), Q, tol); - if(k != 0) - { - BOOST_CHECK_CLOSE( - quantile(bn, P), k, tol); - } - else - { - // Just check quantile is very small: - if((std::numeric_limits::max_exponent <= std::numeric_limits::max_exponent) - && (boost::is_floating_point::value)) - { - // Limit where this is checked: if exponent range is very large we may - // run out of iterations in our root finding algorithm. - BOOST_CHECK(quantile(bn, P) < boost::math::tools::epsilon() * 10); - } - } + CHECK_RELATIVE_ERROR(Q, cdf(complement(bn, k)), tol); + // Limit where a zero quantile is checked: if exponent range is very large we may + // run out of iterations in our root finding algorithm. + bool check_zero = (std::numeric_limits::max_exponent <= std::numeric_limits::max_exponent) + && std::is_floating_point::value; if(k != 0) { - BOOST_CHECK_CLOSE( - quantile(complement(bn, Q)), k, tol); + CHECK_RELATIVE_ERROR(k, quantile(bn, P), tol); + CHECK_RELATIVE_ERROR(k, quantile(complement(bn, Q)), tol); } - else + else if(check_zero) { // Just check quantile is very small: - if((std::numeric_limits::max_exponent <= std::numeric_limits::max_exponent) - && (boost::is_floating_point::value)) - { - // Limit where this is checked: if exponent range is very large we may - // run out of iterations in our root finding algorithm. - BOOST_CHECK(quantile(complement(bn, Q)) < boost::math::tools::epsilon() * 10); - } + CHECK_LE(quantile(bn, P), boost::math::tools::epsilon() * 10); + CHECK_LE(quantile(complement(bn, Q)), boost::math::tools::epsilon() * 10); } // estimate success ratio: - BOOST_CHECK_CLOSE( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, P), - p, tol); + CHECK_RELATIVE_ERROR(p, negative_binomial_distribution::find_lower_bound_on_p(N+k, N, P), tol); // Note we bump up the sample size here, purely for the sake of the test, // internally the function has to adjust the sample size so that we get // the right upper bound, our test undoes this, so we can verify the result. - BOOST_CHECK_CLOSE( - negative_binomial_distribution::find_upper_bound_on_p( - N+k+1, N, Q), - p, tol); + CHECK_RELATIVE_ERROR(p, negative_binomial_distribution::find_upper_bound_on_p(N+k+1, N, Q), tol); - if(Q < P) - { - // - // We check two things here, that the upper and lower bounds - // are the right way around, and that they do actually bracket - // the naive estimate of p = successes / (sample size) - // - BOOST_CHECK( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, Q) - <= - negative_binomial_distribution::find_upper_bound_on_p( - N+k, N, Q) - ); - BOOST_CHECK( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, Q) - <= - N / (N+k) - ); - BOOST_CHECK( - N / (N+k) - <= - negative_binomial_distribution::find_upper_bound_on_p( - N+k, N, Q) - ); - } - else - { - // As above but when P is small. - BOOST_CHECK( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, P) - <= - negative_binomial_distribution::find_upper_bound_on_p( - N+k, N, P) - ); - BOOST_CHECK( - negative_binomial_distribution::find_lower_bound_on_p( - N+k, N, P) - <= - N / (N+k) - ); - BOOST_CHECK( - N / (N+k) - <= - negative_binomial_distribution::find_upper_bound_on_p( - N+k, N, P) - ); - } + // Check that the upper and lower bounds are the right way around, and that they + // bracket the naive estimate of p = successes / (sample size); with Q when Q < P, else with P: + RealType alpha = (Q < P) ? Q : P; + RealType lower = negative_binomial_distribution::find_lower_bound_on_p(N+k, N, alpha); + RealType upper = negative_binomial_distribution::find_upper_bound_on_p(N+k, N, alpha); + CHECK_LE(lower, upper); + CHECK_LE(lower, RealType(N / (N+k))); + CHECK_LE(RealType(N / (N+k)), upper); // Estimate sample size: - BOOST_CHECK_CLOSE( - negative_binomial_distribution::find_minimum_number_of_trials( - k, p, P), - N+k, tol); - BOOST_CHECK_CLOSE( - negative_binomial_distribution::find_maximum_number_of_trials( - k, p, Q), - N+k, tol); + CHECK_RELATIVE_ERROR(N+k, negative_binomial_distribution::find_minimum_number_of_trials(k, p, P), tol); + CHECK_RELATIVE_ERROR(N+k, negative_binomial_distribution::find_maximum_number_of_trials(k, p, Q), tol); // Double check consistency of CDF and PDF by computing the finite sum: RealType sum = 0; @@ -192,31 +111,45 @@ void test_spot( // Test a single spot value against 'known good' values. { sum += pdf(bn, RealType(i)); } - BOOST_CHECK_CLOSE(sum, P, tol); - + CHECK_RELATIVE_ERROR(P, sum, tol); // Complement is not possible since sum is to infinity. - } // + } } // test_spot +// The relative condition number |P / (x dP/dx)| of the quantile x of P, from a central difference +// of the cdf: the quantile's relative error is the cdf's times this, which is large near x = 0. +template +RealType quantile_condition_number(Dist const& d, RealType x, RealType P, bool complemented) +{ + using std::abs; + using std::pow; + if(x == 0) + { + return 1; + } + // Not tools::cbrt_epsilon, whose function-local static a SYCL kernel can't use: + RealType h = x * pow(boost::math::tools::epsilon(), RealType(1) / 3); + RealType up = complemented ? cdf(complement(d, RealType(x + h))) : cdf(d, RealType(x + h)); + RealType down = complemented ? cdf(complement(d, RealType(x - h))) : cdf(d, RealType(x - h)); + RealType kappa = abs(P * 2 * h / (x * (up - down))); + return kappa > 1 ? kappa : RealType(1); +} + template // Any floating-point type RealType. void test_spots(RealType) { // Basic sanity checks, test data is to double precision only - // so set tolerance to 1000 eps expressed as a percent, or - // 1000 eps of type double expressed as a percent, whichever - // is the larger. - + // so set tolerance to 100000 eps, or 100000 eps of type double, + // whichever is the larger. RealType tolerance = (std::max) (boost::math::tools::epsilon(), static_cast(std::numeric_limits::epsilon())); - tolerance *= 100 * 100000.0f; + tolerance *= 100000.0f; + // BOOST_CHECK_SMALL took the percentage tolerance above as an absolute bound: + RealType small_bound = tolerance * 100; - cout << "Tolerance = " << tolerance << "%." << endl; - - RealType tol1eps = boost::math::tools::epsilon() * 2; // Very tight, suit exact values. - //RealType tol2eps = boost::math::tools::epsilon() * 2; // Tight, suit exact values. - RealType tol5eps = boost::math::tools::epsilon() * 5; // Wider 5 epsilon. - cout << "Tolerance 5 eps = " << tol5eps << "%." << endl; + RealType tol1eps = boost::math::tools::epsilon() * 2 / 100; // Very tight, suit exact values. + RealType tol5eps = boost::math::tools::epsilon() * 5 / 100; // Wider 5 epsilon. // Sources of spot test values: @@ -229,521 +162,174 @@ void test_spots(RealType) // And functions.wolfram.com - using boost::math::negative_binomial_distribution; - using ::boost::math::negative_binomial; - using ::boost::math::cdf; - using ::boost::math::pdf; + using ::boost::math::cdf; + using ::boost::math::pdf; // Test negative binomial using cdf spot values from MathCAD cdf = pnbinom(k, r, p). // These test quantiles and complements as well. + // Arguments: successes r, failures k, success fraction p, CDF P, complement Q = 1 - P. test_spot( // pnbinom(1,2,0.5) = 0.5 - static_cast(2), // successes r - static_cast(1), // Number of failures, k - static_cast(0.5), // Probability of success as fraction, p - static_cast(0.5), // Probability of result (CDF), P - static_cast(0.5), // complement CCDF Q = 1 - P - tolerance); + static_cast(2), static_cast(1), static_cast(0.5), + static_cast(0.5), static_cast(0.5), tolerance); test_spot( // pbinom(0, 2, 0.25) - static_cast(2), // successes r - static_cast(0), // Number of failures, k - static_cast(0.25), - static_cast(0.0625), // Probability of result (CDF), P - static_cast(0.9375), // Q = 1 - P - tolerance); + static_cast(2), static_cast(0), static_cast(0.25), + static_cast(0.0625), static_cast(0.9375), tolerance); test_spot( // pbinom(48,8,0.25) - static_cast(8), // successes r - static_cast(48), // Number of failures, k - static_cast(0.25), // Probability of success, p - static_cast(9.826582228110670E-1), // Probability of result (CDF), P - static_cast(1 - 9.826582228110670E-1), // Q = 1 - P - tolerance); + static_cast(8), static_cast(48), static_cast(0.25), + static_cast(9.826582228110670E-1), static_cast(1 - 9.826582228110670E-1), tolerance); test_spot( // pbinom(2,5,0.4) - static_cast(5), // successes r - static_cast(2), // Number of failures, k - static_cast(0.4), // Probability of success, p - static_cast(9.625600000000020E-2), // Probability of result (CDF), P - static_cast(1 - 9.625600000000020E-2), // Q = 1 - P - tolerance); + static_cast(5), static_cast(2), static_cast(0.4), + static_cast(9.625600000000020E-2), static_cast(1 - 9.625600000000020E-2), tolerance); test_spot( // pbinom(10,100,0.9) - static_cast(100), // successes r - static_cast(10), // Number of failures, k - static_cast(0.9), // Probability of success, p - static_cast(4.535522887695670E-1), // Probability of result (CDF), P - static_cast(1 - 4.535522887695670E-1), // Q = 1 - P - tolerance); + static_cast(100), static_cast(10), static_cast(0.9), + static_cast(4.535522887695670E-1), static_cast(1 - 4.535522887695670E-1), tolerance); test_spot( // pbinom(1,100,0.991) - static_cast(100), // successes r - static_cast(1), // Number of failures, k - static_cast(0.991), // Probability of success, p - static_cast(7.693413044217000E-1), // Probability of result (CDF), P - static_cast(1 - 7.693413044217000E-1), // Q = 1 - P - tolerance); + static_cast(100), static_cast(1), static_cast(0.991), + static_cast(7.693413044217000E-1), static_cast(1 - 7.693413044217000E-1), tolerance); test_spot( // pbinom(10,100,0.991) - static_cast(100), // successes r - static_cast(10), // Number of failures, k - static_cast(0.991), // Probability of success, p - static_cast(9.999999940939000E-1), // Probability of result (CDF), P - static_cast(1 - 9.999999940939000E-1), // Q = 1 - P - tolerance); - -if(std::numeric_limits::is_specialized) -{ // An extreme value test that takes 3 minutes using the real concept type - // for which numeric_limits::is_specialized == false, deliberately - // and for which there is no Lanczos approximation defined (also deliberately) - // giving a very slow computation, but with acceptable accuracy. - // A possible enhancement might be to use a normal approximation for - // extreme values, but this is not implemented. - test_spot( // pbinom(100000,100,0.001) - static_cast(100), // successes r - static_cast(100000), // Number of failures, k - static_cast(0.001), // Probability of success, p - static_cast(5.173047534260320E-1), // Probability of result (CDF), P - static_cast(1 - 5.173047534260320E-1), // Q = 1 - P - tolerance*1000); // *1000 is OK 0.51730475350664229 versus - - // functions.wolfram.com - // for I[0.001](100, 100000+1) gives: - // Wolfram 0.517304753506834882009032744488738352004003696396461766326713 - // JM nonLanczos 0.51730475350664229 differs at the 13th decimal digit. - // MathCAD 0.51730475342603199 differs at 10th decimal digit. - - // Error tests: - check_out_of_range >(20, 0.5); - BOOST_MATH_CHECK_THROW(negative_binomial_distribution(0, 0.5), std::domain_error); - BOOST_MATH_CHECK_THROW(negative_binomial_distribution(-2, 0.5), std::domain_error); - BOOST_MATH_CHECK_THROW(negative_binomial_distribution(20, -0.5), std::domain_error); - BOOST_MATH_CHECK_THROW(negative_binomial_distribution(20, 1.5), std::domain_error); -} - // End of single spot tests using RealType - - - // Tests on PDF: - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(2), static_cast(0.5)), - static_cast(0) ), // k = 0. - static_cast(0.25), // 0 - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(4), static_cast(0.5)), - static_cast(0)), // k = 0. - static_cast(0.0625), // exact 1/16 - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(20), static_cast(0.25)), - static_cast(0)), // k = 0 - static_cast(9.094947017729270E-13), // pbinom(0,20,0.25) = 9.094947017729270E-13 - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(20), static_cast(0.2)), - static_cast(0)), // k = 0 - static_cast(1.0485760000000003e-014), // MathCAD 1.048576000000000E-14 - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(10), static_cast(0.1)), - static_cast(0)), // k = 0. - static_cast(1e-10), // MathCAD says zero, but suffers cancellation error? - tolerance); - - BOOST_CHECK_CLOSE( - pdf(negative_binomial_distribution(static_cast(20), static_cast(0.1)), - static_cast(0)), // k = 0. - static_cast(1e-20), // MathCAD says zero, but suffers cancellation error? - tolerance); - - - BOOST_CHECK_CLOSE( // . - pdf(negative_binomial_distribution(static_cast(20), static_cast(0.9)), - static_cast(0)), // k. - static_cast(1.215766545905690E-1), // k=20 p = 0.9 - tolerance); + static_cast(100), static_cast(10), static_cast(0.991), + static_cast(9.999999940939000E-1), static_cast(1 - 9.999999940939000E-1), tolerance); + + if(std::numeric_limits::is_specialized) + { // An extreme value test that takes 3 minutes using the real concept type + // for which numeric_limits::is_specialized == false, deliberately + // and for which there is no Lanczos approximation defined (also deliberately) + // giving a very slow computation, but with acceptable accuracy. + // A possible enhancement might be to use a normal approximation for + // extreme values, but this is not implemented. + test_spot( // pbinom(100000,100,0.001) + static_cast(100), static_cast(100000), static_cast(0.001), + static_cast(5.173047534260320E-1), static_cast(1 - 5.173047534260320E-1), + tolerance*1000); // *1000 is OK 0.51730475350664229 versus + + // functions.wolfram.com + // for I[0.001](100, 100000+1) gives: + // Wolfram 0.517304753506834882009032744488738352004003696396461766326713 + // JM nonLanczos 0.51730475350664229 differs at the 13th decimal digit. + // MathCAD 0.51730475342603199 differs at 10th decimal digit. + + // Error tests: + check_out_of_range >(20, 0.5); + CHECK_THROW(negative_binomial_distribution(0, 0.5), std::domain_error); + CHECK_THROW(negative_binomial_distribution(-2, 0.5), std::domain_error); + CHECK_THROW(negative_binomial_distribution(20, -0.5), std::domain_error); + CHECK_THROW(negative_binomial_distribution(20, 1.5), std::domain_error); + } + // End of single spot tests using RealType + + auto nb = [](double r, double p) { return negative_binomial_distribution(static_cast(r), static_cast(p)); }; + + // Tests on PDF at k = 0: + CHECK_RELATIVE_ERROR(static_cast(0.25), pdf(nb(2, 0.5), static_cast(0)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(0.0625), pdf(nb(4, 0.5), static_cast(0)), tolerance); // exact 1/16 + CHECK_RELATIVE_ERROR(static_cast(9.094947017729270E-13), pdf(nb(20, 0.25), static_cast(0)), tolerance); // pbinom(0,20,0.25) + CHECK_RELATIVE_ERROR(static_cast(1.0485760000000003e-014), pdf(nb(20, 0.2), static_cast(0)), tolerance); // MathCAD 1.048576000000000E-14 + CHECK_RELATIVE_ERROR(static_cast(1e-10), pdf(nb(10, 0.1), static_cast(0)), tolerance); // MathCAD says zero, but suffers cancellation error? + CHECK_RELATIVE_ERROR(static_cast(1e-20), pdf(nb(20, 0.1), static_cast(0)), tolerance); // MathCAD says zero, but suffers cancellation error? + CHECK_RELATIVE_ERROR(static_cast(1.215766545905690E-1), pdf(nb(20, 0.9), static_cast(0)), tolerance); // Tests on cdf: // MathCAD pbinom k, r, p) == failures, successes, probability. - - BOOST_CHECK_CLOSE(cdf( - negative_binomial_distribution(static_cast(2), static_cast(0.5)), // successes = 2,prob 0.25 - static_cast(0) ), // k = 0 - static_cast(0.25), // probability 1/4 - tolerance); - - BOOST_CHECK_CLOSE(cdf(complement( - negative_binomial_distribution(static_cast(2), static_cast(0.5)), // successes = 2,prob 0.25 - static_cast(0) )), // k = 0 - static_cast(0.75), // probability 3/4 - tolerance); - BOOST_CHECK_CLOSE( // k = 1. - cdf(negative_binomial_distribution(static_cast(20), static_cast(0.25)), - static_cast(1)), // k =1. - static_cast(1.455191522836700E-11), - tolerance); - - BOOST_CHECK_SMALL( // Check within an epsilon with CHECK_SMALL - cdf(negative_binomial_distribution(static_cast(20), static_cast(0.25)), - static_cast(1)) - - static_cast(1.455191522836700E-11), - tolerance ); + CHECK_RELATIVE_ERROR(static_cast(0.25), cdf(nb(2, 0.5), static_cast(0)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(0.75), cdf(complement(nb(2, 0.5), static_cast(0))), tolerance); + CHECK_RELATIVE_ERROR(static_cast(1.455191522836700E-11), cdf(nb(20, 0.25), static_cast(1)), tolerance); + CHECK_ABSOLUTE_ERROR(static_cast(1.455191522836700E-11), cdf(nb(20, 0.25), static_cast(1)), small_bound); // Some exact (probably - judging by trailing zeros) values. - BOOST_CHECK_CLOSE( - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)), // k. - static_cast(1.525878906250000E-5), - tolerance); - - BOOST_CHECK_CLOSE( - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)), // k. - static_cast(1.525878906250000E-5), - tolerance); - - BOOST_CHECK_SMALL( - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)) - - static_cast(1.525878906250000E-5), - tolerance ); - - BOOST_CHECK_CLOSE( // k = 1. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1)), // k. - static_cast(1.068115234375010E-4), - tolerance); - - BOOST_CHECK_CLOSE( // k = 2. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(2)), // k. - static_cast(4.158020019531300E-4), - tolerance); - - BOOST_CHECK_CLOSE( // k = 3. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(3)), // k.bristow - static_cast(1.188278198242200E-3), - tolerance); - - BOOST_CHECK_CLOSE( // k = 4. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(4)), // k. - static_cast(2.781510353088410E-3), - tolerance); - - BOOST_CHECK_CLOSE( // k = 5. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(5)), // k. - static_cast(5.649328231811500E-3), - tolerance); - - BOOST_CHECK_CLOSE( // k = 6. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(6)), // k. - static_cast(1.030953228473680E-2), - tolerance); - - BOOST_CHECK_CLOSE( // k = 7. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(7)), // k. - static_cast(1.729983836412430E-2), - tolerance); - - BOOST_CHECK_CLOSE( // k = 8. - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(8)), // k = n. - static_cast(2.712995628826370E-2), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(48)), // k - static_cast(9.826582228110670E-1), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(64)), // k - static_cast(9.990295004935590E-1), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(5), static_cast(0.4)), - static_cast(26)), // k - static_cast(9.989686246611190E-1), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(5), static_cast(0.4)), - static_cast(2)), // k failures - static_cast(9.625600000000020E-2), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(50), static_cast(0.9)), - static_cast(20)), // k - static_cast(9.999970854144170E-1), - tolerance); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(500), static_cast(0.7)), - static_cast(200)), // k - static_cast(2.172846379930550E-1), - tolerance* 2); - - BOOST_CHECK_CLOSE( // - cdf(negative_binomial_distribution(static_cast(50), static_cast(0.7)), - static_cast(20)), // k - static_cast(4.550203671301790E-1), - tolerance); + CHECK_RELATIVE_ERROR(static_cast(1.525878906250000E-5), cdf(nb(8, 0.25), static_cast(0)), tolerance); + CHECK_ABSOLUTE_ERROR(static_cast(1.525878906250000E-5), cdf(nb(8, 0.25), static_cast(0)), small_bound); + const double eight_quarter[][2] = { + {1, 1.068115234375010E-4}, {2, 4.158020019531300E-4}, {3, 1.188278198242200E-3}, {4, 2.781510353088410E-3}, + {5, 5.649328231811500E-3}, {6, 1.030953228473680E-2}, {7, 1.729983836412430E-2}, {8, 2.712995628826370E-2}, + {48, 9.826582228110670E-1}, {64, 9.990295004935590E-1} }; + for(auto const& row : eight_quarter) + { + CHECK_RELATIVE_ERROR(static_cast(row[1]), cdf(nb(8, 0.25), static_cast(row[0])), tolerance); + } + CHECK_RELATIVE_ERROR(static_cast(9.989686246611190E-1), cdf(nb(5, 0.4), static_cast(26)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(9.625600000000020E-2), cdf(nb(5, 0.4), static_cast(2)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(9.999970854144170E-1), cdf(nb(50, 0.9), static_cast(20)), tolerance); + CHECK_RELATIVE_ERROR(static_cast(2.172846379930550E-1), cdf(nb(500, 0.7), static_cast(200)), tolerance * 2); + CHECK_RELATIVE_ERROR(static_cast(4.550203671301790E-1), cdf(nb(50, 0.7), static_cast(20)), tolerance); // Tests of other functions, mean and other moments ... - negative_binomial_distribution dist(static_cast(8), static_cast(0.25)); + negative_binomial_distribution dist = nb(8, 0.25); using namespace std; // ADL of std names. - // mean: - BOOST_CHECK_CLOSE( - mean(dist), static_cast(8 * (1 - 0.25) /0.25), tol5eps); - BOOST_CHECK_CLOSE( - mode(dist), static_cast(21), tol1eps); - // variance: - BOOST_CHECK_CLOSE( - variance(dist), static_cast(8 * (1 - 0.25) / (0.25 * 0.25)), tol5eps); - // std deviation: - BOOST_CHECK_CLOSE( - standard_deviation(dist), // 9.79795897113271239270 - static_cast(9.797958971132712392789136298823565567864L), // using functions.wolfram.com - // 9.79795897113271152534 == sqrt(8 * (1 - 0.25) / (0.25 * 0.25))) - tol5eps * 100); - BOOST_CHECK_CLOSE( - skewness(dist), // - static_cast(0.71443450831176036), - // using http://mathworld.wolfram.com/skewness.html - tolerance); - BOOST_CHECK_CLOSE( - kurtosis_excess(dist), // - static_cast(0.7604166666666666666666666666666666666666L), // using Wikipedia Kurtosis(excess) formula - tol5eps * 100); - BOOST_CHECK_CLOSE( - kurtosis(dist), // true - static_cast(3.76041666666666666666666666666666666666666L), // - tol5eps * 100); - // hazard: + CHECK_RELATIVE_ERROR(static_cast(8 * (1 - 0.25) /0.25), mean(dist), tol5eps); + CHECK_RELATIVE_ERROR(static_cast(21), mode(dist), tol1eps); + CHECK_RELATIVE_ERROR(static_cast(8 * (1 - 0.25) / (0.25 * 0.25)), variance(dist), tol5eps); + CHECK_RELATIVE_ERROR(static_cast(9.797958971132712392789136298823565567864L), standard_deviation(dist), tol5eps * 100); // using functions.wolfram.com + CHECK_RELATIVE_ERROR(static_cast(0.71443450831176036), skewness(dist), tolerance); // using http://mathworld.wolfram.com/skewness.html + CHECK_RELATIVE_ERROR(static_cast(0.7604166666666666666666666666666666666666L), kurtosis_excess(dist), tol5eps * 100); // using Wikipedia Kurtosis(excess) formula + CHECK_RELATIVE_ERROR(static_cast(3.76041666666666666666666666666666666666666L), kurtosis(dist), tol5eps * 100); RealType x = static_cast(0.125); - BOOST_CHECK_CLOSE( - hazard(dist, x) - , pdf(dist, x) / cdf(complement(dist, x)), tol5eps); - // cumulative hazard: - BOOST_CHECK_CLOSE( - chf(dist, x), -log(cdf(complement(dist, x))), tol5eps); - // coefficient_of_variation: - BOOST_CHECK_CLOSE( - coefficient_of_variation(dist) - , standard_deviation(dist) / mean(dist), tol5eps); + CHECK_RELATIVE_ERROR(pdf(dist, x) / cdf(complement(dist, x)), hazard(dist, x), tol5eps); + CHECK_RELATIVE_ERROR(-log(cdf(complement(dist, x))), chf(dist, x), tol5eps); + CHECK_RELATIVE_ERROR(standard_deviation(dist) / mean(dist), coefficient_of_variation(dist), tol5eps); // Special cases for PDF: - BOOST_CHECK_EQUAL( - pdf( - negative_binomial_distribution(static_cast(8), static_cast(0)), // - static_cast(0)), - static_cast(0) ); - - BOOST_CHECK_EQUAL( - pdf( - negative_binomial_distribution(static_cast(8), static_cast(0)), - static_cast(0.0001)), - static_cast(0) ); - - BOOST_CHECK_EQUAL( - pdf( - negative_binomial_distribution(static_cast(8), static_cast(1)), - static_cast(0.001)), - static_cast(0) ); - - BOOST_CHECK_EQUAL( - pdf( - negative_binomial_distribution(static_cast(8), static_cast(1)), - static_cast(8)), - static_cast(0) ); - - BOOST_CHECK_SMALL( - pdf( - negative_binomial_distribution(static_cast(2), static_cast(0.25)), - static_cast(0))- - static_cast(0.0625), - 2 * boost::math::tools::epsilon() ); // Expect exact, but not quite. - // numeric_limits::epsilon()); // Not suitable for real concept! + CHECK_EQUAL(pdf(nb(8, 0), static_cast(0)), static_cast(0)); + CHECK_EQUAL(pdf(nb(8, 0), static_cast(0.0001)), static_cast(0)); + CHECK_EQUAL(pdf(nb(8, 1), static_cast(0.001)), static_cast(0)); + CHECK_EQUAL(pdf(nb(8, 1), static_cast(8)), static_cast(0)); + // Expect exact, but not quite: + CHECK_ABSOLUTE_ERROR(static_cast(0.0625), pdf(nb(2, 0.25), static_cast(0)), RealType(2 * boost::math::tools::epsilon())); // Quantile boundary cases checks: - BOOST_CHECK_EQUAL( - quantile( // zero P < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)), - static_cast(0)); - - BOOST_CHECK_EQUAL( - quantile( // min P < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(boost::math::tools::min_value())), - static_cast(0)); - - BOOST_CHECK_CLOSE_FRACTION( - quantile( // Small P < cdf(0) so should be near zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(boost::math::tools::epsilon())), // - static_cast(0), - tol5eps); - - BOOST_CHECK_CLOSE( - quantile( // Small P < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0.0001)), - static_cast(0.95854156929288470), - tolerance); - - //BOOST_CHECK( // Fails with overflow for real_concept - //quantile( // Small P near 1 so k failures should be big. - //negative_binomial_distribution(static_cast(8), static_cast(0.25)), - //static_cast(1 - boost::math::tools::epsilon())) <= - //static_cast(189.56999032670058) // 106.462769 for float - //); + CHECK_EQUAL(quantile(dist, static_cast(0)), static_cast(0)); // zero P < cdf(0) so should be exactly zero. + CHECK_EQUAL(quantile(dist, boost::math::tools::min_value()), static_cast(0)); // min P < cdf(0) + // Small P < cdf(0) so should be near zero; a relative check against zero requires zero: + CHECK_RELATIVE_ERROR(static_cast(0), quantile(dist, boost::math::tools::epsilon()), tol5eps); + CHECK_RELATIVE_ERROR(static_cast(0.95854156929288470), quantile(dist, static_cast(0.0001)), tolerance); if(std::numeric_limits::has_infinity) - { // BOOST_CHECK tests for infinity using std::numeric_limits<>::infinity() - // Note that infinity is not implemented for real_concept, so these tests - // are only done for types, like built-in float, double.. that have infinity. - // Note that these assume that BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error. - // #define BOOST_MATH_THROW_ON_OVERFLOW_POLICY == throw_on_error would throw here. - // #define BOOST_MAT_DOMAIN_ERROR_POLICY IS defined throw_on_error, - // so the throw path of error handling is tested below with BOOST_MATH_CHECK_THROW tests. - - BOOST_CHECK( - quantile( // At P == 1 so k failures should be infinite. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1)) == - //static_cast(boost::math::tools::infinity()) - static_cast(std::numeric_limits::infinity()) ); - - BOOST_CHECK_EQUAL( - quantile( // At 1 == P so should be infinite. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1)), // - std::numeric_limits::infinity() ); - - BOOST_CHECK_EQUAL( - quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0))), - std::numeric_limits::infinity() ); - } // test for infinity using std::numeric_limits<>::infinity() + { // Infinity is not implemented for real_concept, so these tests are only done for built-in types. + // These assume that BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error. + CHECK_EQUAL(quantile(dist, static_cast(1)), std::numeric_limits::infinity()); // At P == 1 so k failures should be infinite. + CHECK_EQUAL(quantile(complement(dist, static_cast(0))), std::numeric_limits::infinity()); // Q zero so P == 1. + } else - { // real_concept case, so check it throws rather than returning infinity. - BOOST_CHECK_EQUAL( - quantile( // At P == 1 so k failures should be infinite. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1)), - boost::math::tools::max_value() ); - - BOOST_CHECK_EQUAL( - quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0))), - boost::math::tools::max_value()); + { // real_concept case, which returns max_value rather than infinity. + CHECK_EQUAL(quantile(dist, static_cast(1)), boost::math::tools::max_value()); + CHECK_EQUAL(quantile(complement(dist, static_cast(0))), boost::math::tools::max_value()); } - BOOST_CHECK( // Should work for built-in and real_concept. - quantile(complement( // Q very near to 1 so P nearly 1 < so should be large > 384. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(boost::math::tools::min_value()))) - >= static_cast(384) ); - - BOOST_CHECK_EQUAL( - quantile( // P == 0 < cdf(0) so should be zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(0)), - static_cast(0)); + // Q very near to 1 so P nearly 1 so should be large > 384: + CHECK_GE(quantile(complement(dist, boost::math::tools::min_value())), static_cast(384)); + CHECK_EQUAL(quantile(dist, static_cast(0)), static_cast(0)); // P == 0 < cdf(0) so should be zero. // Quantile Complement boundary cases: - - BOOST_CHECK_EQUAL( - quantile(complement( // Q = 1 so P = 0 < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1))), - static_cast(0) - ); - - BOOST_CHECK_EQUAL( - quantile(complement( // Q very near 1 so P == epsilon < cdf(0) so should be exactly zero. - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(1 - boost::math::tools::epsilon()))), - static_cast(0) - ); + CHECK_EQUAL(quantile(complement(dist, static_cast(1))), static_cast(0)); // Q = 1 so P = 0 < cdf(0) + CHECK_EQUAL(quantile(complement(dist, static_cast(1 - boost::math::tools::epsilon()))), static_cast(0)); // P == epsilon < cdf(0) // Check that duff arguments throw domain_error: - BOOST_MATH_CHECK_THROW( - pdf( // Negative successes! - negative_binomial_distribution(static_cast(-1), static_cast(0.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - pdf( // Negative success_fraction! - negative_binomial_distribution(static_cast(8), static_cast(-0.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - pdf( // Success_fraction > 1! - negative_binomial_distribution(static_cast(8), static_cast(1.25)), - static_cast(0)), - std::domain_error - ); - BOOST_MATH_CHECK_THROW( - pdf( // Negative k argument ! - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(-1)), - std::domain_error - ); - //BOOST_MATH_CHECK_THROW( - //pdf( // Unlike binomial there is NO limit on k (failures) - //negative_binomial_distribution(static_cast(8), static_cast(0.25)), - //static_cast(9)), std::domain_error - //); - BOOST_MATH_CHECK_THROW( - cdf( // Negative k argument ! - negative_binomial_distribution(static_cast(8), static_cast(0.25)), - static_cast(-1)), - std::domain_error - ); - BOOST_MATH_CHECK_THROW( - cdf( // Negative success_fraction! - negative_binomial_distribution(static_cast(8), static_cast(-0.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - cdf( // Success_fraction > 1! - negative_binomial_distribution(static_cast(8), static_cast(1.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - quantile( // Negative success_fraction! - negative_binomial_distribution(static_cast(8), static_cast(-0.25)), - static_cast(0)), std::domain_error - ); - BOOST_MATH_CHECK_THROW( - quantile( // Success_fraction > 1! - negative_binomial_distribution(static_cast(8), static_cast(1.25)), - static_cast(0)), std::domain_error - ); + CHECK_THROW(pdf(nb(-1, 0.25), static_cast(0)), std::domain_error); // Negative successes! + CHECK_THROW(pdf(nb(8, -0.25), static_cast(0)), std::domain_error); // Negative success_fraction! + CHECK_THROW(pdf(nb(8, 1.25), static_cast(0)), std::domain_error); // Success_fraction > 1! + CHECK_THROW(pdf(nb(8, 0.25), static_cast(-1)), std::domain_error); // Negative k argument! + // Unlike binomial there is NO limit on k (failures). + CHECK_THROW(cdf(nb(8, 0.25), static_cast(-1)), std::domain_error); // Negative k argument! + CHECK_THROW(cdf(nb(8, -0.25), static_cast(0)), std::domain_error); // Negative success_fraction! + CHECK_THROW(cdf(nb(8, 1.25), static_cast(0)), std::domain_error); // Success_fraction > 1! + CHECK_THROW(quantile(nb(8, -0.25), static_cast(0)), std::domain_error); // Negative success_fraction! + CHECK_THROW(quantile(nb(8, 1.25), static_cast(0)), std::domain_error); // Success_fraction > 1! // End of check throwing 'duff' out-of-domain values. +// The table's SC_ uses BOOST_JOIN, from Boost.Config, which a standalone build doesn't have: +#ifndef BOOST_JOIN +#define BOOST_JOIN(X, Y) BOOST_MATH_JOIN(X, Y) +#endif #define T RealType #include "negative_binomial_quantile.ipp" - for(unsigned i = 0; i < negative_binomial_quantile_data.size(); ++i) + for(auto const& row : negative_binomial_quantile_data) { using namespace boost::math::policies; typedef policy > P1; @@ -752,60 +338,53 @@ if(std::numeric_limits::is_specialized) typedef policy > P4; typedef policy > P5; typedef policy > P6; + RealType r = static_cast(row[0]); + RealType sf = static_cast(row[1]); + RealType p = static_cast(row[2]); + RealType lower = static_cast(row[3]); + RealType upper = static_cast(row[4]); RealType tol = boost::math::tools::epsilon() * 700; - if(!boost::is_floating_point::value) + if(!std::is_floating_point::value) tol *= 10; // no lanczos approximation implies less accuracy // - // Check full real value first: + // Check full real value first, to the cdf's tolerance times the quantile's condition number + // (548 for r = 18.3, p = 0.906, where the complement's quantile is 0.0061): // - negative_binomial_distribution p1(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - RealType x = quantile(p1, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_CLOSE_FRACTION(x, negative_binomial_quantile_data[i][3], tol); - x = quantile(complement(p1, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_CLOSE_FRACTION(x, negative_binomial_quantile_data[i][4], tol); + negative_binomial_distribution p1(r, sf); + CHECK_RELATIVE_ERROR(lower, quantile(p1, p), tol * quantile_condition_number(p1, lower, p, false)); + CHECK_RELATIVE_ERROR(upper, quantile(complement(p1, p)), tol * quantile_condition_number(p1, upper, p, true)); // // Now with round down to integer: // - negative_binomial_distribution p2(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p2, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][3])); - x = quantile(complement(p2, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][4])); + negative_binomial_distribution p2(r, sf); + CHECK_EQUAL(quantile(p2, p), RealType(floor(lower))); + CHECK_EQUAL(quantile(complement(p2, p)), RealType(floor(upper))); // // Now with round up to integer: // - negative_binomial_distribution p3(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p3, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, ceil(negative_binomial_quantile_data[i][3])); - x = quantile(complement(p3, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, ceil(negative_binomial_quantile_data[i][4])); + negative_binomial_distribution p3(r, sf); + CHECK_EQUAL(quantile(p3, p), RealType(ceil(lower))); + CHECK_EQUAL(quantile(complement(p3, p)), RealType(ceil(upper))); // // Now with round to integer "outside": // - negative_binomial_distribution p4(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p4, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? floor(negative_binomial_quantile_data[i][3]) : ceil(negative_binomial_quantile_data[i][3])); - x = quantile(complement(p4, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? ceil(negative_binomial_quantile_data[i][4]) : floor(negative_binomial_quantile_data[i][4])); + negative_binomial_distribution p4(r, sf); + CHECK_EQUAL(quantile(p4, p), RealType(p < 0.5f ? floor(lower) : ceil(lower))); + CHECK_EQUAL(quantile(complement(p4, p)), RealType(p < 0.5f ? ceil(upper) : floor(upper))); // // Now with round to integer "inside": // - negative_binomial_distribution p5(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p5, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? ceil(negative_binomial_quantile_data[i][3]) : floor(negative_binomial_quantile_data[i][3])); - x = quantile(complement(p5, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? floor(negative_binomial_quantile_data[i][4]) : ceil(negative_binomial_quantile_data[i][4])); + negative_binomial_distribution p5(r, sf); + CHECK_EQUAL(quantile(p5, p), RealType(p < 0.5f ? ceil(lower) : floor(lower))); + CHECK_EQUAL(quantile(complement(p5, p)), RealType(p < 0.5f ? floor(upper) : ceil(upper))); // // Now with round to nearest integer: // - negative_binomial_distribution p6(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]); - x = quantile(p6, negative_binomial_quantile_data[i][2]); - BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][3] + 0.5f)); - x = quantile(complement(p6, negative_binomial_quantile_data[i][2])); - BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][4] + 0.5f)); + negative_binomial_distribution p6(r, sf); + CHECK_EQUAL(quantile(p6, p), RealType(floor(lower + 0.5f))); + CHECK_EQUAL(quantile(complement(p6, p)), RealType(floor(upper + 0.5f))); } - - return; +#undef T } // template void test_spots(RealType) // Any floating-point type RealType. template @@ -816,80 +395,53 @@ void test_entropy(RealType) RealType tol = 16 * boost::math::tools::epsilon(); auto check = [&](RealType a, RealType b, RealType h, RealType tolerance) { - BOOST_CHECK_CLOSE_FRACTION(entropy(negative_binomial_distribution(a, b)), h, tolerance); + CHECK_RELATIVE_ERROR(h, entropy(negative_binomial_distribution(a, b)), tolerance); }; check(static_cast(5.0L), static_cast(0.25L), static_cast(3.39359238247812518456836017160705370313434987L), tol); check(static_cast(0.5L), static_cast(0.5L), static_cast(0.943370759269334533806787413559674682998205284L), tol); check(static_cast(1000.5L), static_cast(0.0625L), static_cast(7.61305197413188991442493027236853895416590232L), tol); check(static_cast(1000000.0L), static_cast(0.5L), static_cast(8.67326702746660536495437769998638212050453356L), tol); - // About 1e4 terms; the sum is compensated, which a compiler that reassociates floating-point - // arithmetic (icpx's default) undoes, leaving ~70 eps: - check(static_cast(0.125L), static_cast(0.00390625L), static_cast(3.05357381884703876385483729149700634949203959L), 128 * tol / 16); + // About 1e4 terms, summed with compensation: + check(static_cast(0.125L), static_cast(0.00390625L), static_cast(3.05357381884703876385483729149700634949203959L), tol); check(static_cast(2.5L), static_cast(0.001953125L), static_cast(7.96729439209919695998955366089632316751290691L), tol); check(static_cast(100000000.0L), static_cast(0.9990234375L), static_cast(7.16451918702377261816611719177425949753810124L), tol); - BOOST_CHECK_EQUAL(entropy(negative_binomial_distribution(3, 1)), RealType(0)); + CHECK_EQUAL(entropy(negative_binomial_distribution(3, 1)), RealType(0)); } -BOOST_AUTO_TEST_CASE( test_main ) +int main() { // Check that can generate negative_binomial distribution using the two convenience methods: using namespace boost::math; - negative_binomial mynb1(2., 0.5); // Using typedef - default type is double. - negative_binomial_distribution<> myf2(2., 0.5); // Using default RealType double. - - // Basic sanity-check spot values. - - // Test some simple double only examples. - negative_binomial_distribution my8dist(8., 0.25); - // 8 successes (r), 0.25 success fraction = 35% or 1 in 4 successes. - // Note: double values (matching the distribution definition) avoid the need for any casting. + negative_binomial mynb1(2., 0.5); // Using typedef - default type is double. + negative_binomial_distribution<> myf2(2., 0.5); // Using default RealType double. + static_cast(mynb1); + static_cast(myf2); // Check accessor functions return exact values for double at least. - BOOST_CHECK_EQUAL(my8dist.successes(), static_cast(8)); - BOOST_CHECK_EQUAL(my8dist.success_fraction(), static_cast(1./4.)); + negative_binomial_distribution my8dist(8., 0.25); + CHECK_EQUAL(my8dist.successes(), static_cast(8)); + CHECK_EQUAL(my8dist.success_fraction(), static_cast(1./4.)); // (Parameter value, arbitrarily zero, only communicates the floating point type). #ifdef TEST_FLOAT - test_spots(0.0F); // Test float. + test_spots(0.0F); test_entropy(0.0F); #endif #ifdef TEST_DOUBLE - test_spots(0.0); // Test double. + test_spots(0.0); test_entropy(0.0); #endif #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS #ifdef TEST_LDOUBLE - test_spots(0.0L); // Test long double. + test_spots(0.0L); test_entropy(0.0L); #endif #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS #ifdef TEST_REAL_CONCEPT - test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. - test_entropy(boost::math::concepts::real_concept(0.)); + test_spots(boost::math::concepts::real_concept(0.)); + test_entropy(boost::math::concepts::real_concept(0.)); #endif - #endif -#else - std::cout << "The long double tests have been disabled on this platform " - "either because the long double overloads of the usual math functions are " - "not available at all, or because they are too inaccurate for these tests " - "to pass." << std::endl; #endif - - -} // BOOST_AUTO_TEST_CASE( test_main ) - -/* - -Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_negative_binomial.exe" -Running 1 test case... -Tolerance = 0.0119209%. -Tolerance 5 eps = 5.96046e-007%. -Tolerance = 2.22045e-011%. -Tolerance 5 eps = 1.11022e-015%. -Tolerance = 2.22045e-011%. -Tolerance 5 eps = 1.11022e-015%. -Tolerance = 2.22045e-011%. -Tolerance 5 eps = 1.11022e-015%. -*** No errors detected - -*/ +#endif + return boost::math::test::report_errors(); +} diff --git a/test/test_out_of_range.hpp b/test/test_out_of_range.hpp index 8b7e723c07..2971ec9708 100644 --- a/test/test_out_of_range.hpp +++ b/test/test_out_of_range.hpp @@ -10,7 +10,12 @@ #include #include -#include +// Included after math_unit_test.hpp, the checks use its CHECK_THROW instead of Boost.Test: +#ifdef BOOST_MATH_TEST_TEST_HPP +# define BOOST_CHECK_THROW(x, y) CHECK_THROW(x, y) +#else +# include +#endif /*` check_out_of_range functions check that bad parameters passed to constructors and functions throw domain_error exceptions. diff --git a/test/test_pow1p.cpp b/test/test_pow1p.cpp index 4e07646122..330f3f4a5d 100644 --- a/test/test_pow1p.cpp +++ b/test/test_pow1p.cpp @@ -68,8 +68,6 @@ void test() CHECK_EQUAL(boost::math::isnan(boost::math::pow1p(boost::math::numeric_limits::quiet_NaN(), boost::math::numeric_limits::infinity())), true); } - // icpx defaults to -fp-model=fast, which breaks the overflow handling and error-free transforms. - #ifndef BOOST_MATH_ENABLE_SYCL // Tiny x, huge finite y: the exponent y*log1p(x) is far out of range, and the // result must overflow or underflow cleanly rather than become inf*0 = NaN. using std::ldexp; @@ -79,7 +77,6 @@ void test() CHECK_EQUAL(boost::math::pow1p(tiny, -huge), T(0)); CHECK_EQUAL(boost::math::pow1p(-tiny, huge), T(0)); CHECK_EQUAL(boost::math::pow1p(-tiny, -huge), boost::math::numeric_limits::infinity()); - #endif // pow(+/-inf, y) CHECK_EQUAL(boost::math::pow1p(boost::math::numeric_limits::infinity(), T(2)), boost::math::numeric_limits::infinity()); @@ -100,7 +97,6 @@ void test() // (1+x) < 0 CHECK_ULP_CLOSE(boost::math::pow1p(T(-3), T(2)), pow(T(-2), T(2)), 10); - #ifndef BOOST_MATH_ENABLE_SYCL // MinGW's long double exp and log lose hundreds of ulp for these huge results, so skip them there. #if defined(__MINGW32__) || defined(__CYGWIN__) const bool check_huge = boost::math::tools::digits() <= 53; @@ -137,7 +133,6 @@ void test() CHECK_ULP_CLOSE(T(1.9141657450800789253871830037475352857803e-261L), boost::math::pow1p(T(-2.1413450492090156e-22), T(2.8034969310889015e+24)), 10); CHECK_ULP_CLOSE(T(1.5108017556765851783235726673009126604833e+287L), boost::math::pow1p(T(6.3377799876470245e-17), T(1.0433536087590337e+19)), 10); } - #endif // x < 0 std::mt19937_64 rng;