diff --git a/tests/test_special_functions.cc b/tests/test_special_functions.cc index 9756cc7..85e107d 100644 --- a/tests/test_special_functions.cc +++ b/tests/test_special_functions.cc @@ -207,16 +207,30 @@ TEMPLATE_TEST_CASE("special_functions", "", float, double, long double) { SECTION("GammaP") { using tuple = arg_res_tuple; - auto x_yref{ - GENERATE(tuple{T(2.0l), T(0.0l), T(0.0000000000000000000000000000000000000000e+00l)}, - tuple{T(2.0l), T(1.0l), T(2.6424111765711535680895245967707826510838e-01l)}, - tuple{T(2.0l), T(2.0l), T(5.9399415029016192431800151508254678977711e-01l)}, - tuple{T(2.0l), T(3.0l), T(8.0085172652854422808263033739975289347320e-01l)}, - tuple{T(2.0l), T(4.0l), T(9.0842180555632909853140989363379378894044e-01l)}, - tuple{T(2.0l), T(5.0l), T(9.5957231800548719742018370946110945450690e-01l)}, - tuple{T(2.0l), T(6.0l), T(9.8264873476333549103868382798428332475945e-01l)}, - tuple{T(2.0l), T(7.0l), T(9.9270494427556387033597491132472573898821e-01l)}, - tuple{T(2.0l), T(8.0l), T(9.9698083634887739345060749786797225082620e-01l)})}; + auto x_yref{GENERATE( + tuple{T(2.0l), T(0.0l), T(0.0000000000000000000000000000000000000000e+00l)}, + tuple{T(2.0l), T(1.0l), T(2.6424111765711535680895245967707826510838e-01l)}, + tuple{T(2.0l), T(2.0l), T(5.9399415029016192431800151508254678977711e-01l)}, + tuple{T(2.0l), T(3.0l), T(8.0085172652854422808263033739975289347320e-01l)}, + tuple{T(2.0l), T(4.0l), T(9.0842180555632909853140989363379378894044e-01l)}, + tuple{T(2.0l), T(5.0l), T(9.5957231800548719742018370946110945450690e-01l)}, + tuple{T(2.0l), T(6.0l), T(9.8264873476333549103868382798428332475945e-01l)}, + tuple{T(2.0l), T(7.0l), T(9.9270494427556387033597491132472573898821e-01l)}, + tuple{T(2.0l), T(8.0l), T(9.9698083634887739345060749786797225082620e-01l)}, + tuple{T(12.5l), T(9.0l), T(1.5760928441953977820468234649142238289159e-01l)}, + tuple{T(12.5l), T(11.0l), T(3.6425597236137699758564203792074202059559e-01l)}, + tuple{T(12.5l), T(13.0l), T(5.9240130858255329901102834319115488696740e-01l)}, + tuple{T(12.5l), T(15.0l), T(7.7571099516559608601853775825027997448055e-01l)}, + tuple{T(12.5l), T(17.0l), T(8.9209214616519247831579080497239130365603e-01l)}, + tuple{T(40.0l), T(32.0l), T(9.5602816630230858703771575061447476503886e-02l)}, + tuple{T(40.0l), T(34.0l), T(1.7173502161590030279066281224095312864764e-01l)}, + tuple{T(40.0l), T(36.0l), T(2.7369639320460330667999631204057987978994e-01l)}, + tuple{T(40.0l), T(38.0l), T(3.9408063211889345109193252781723965433945e-01l)}, + tuple{T(40.0l), T(40.0l), T(5.2102886106105516334095689757755319315557e-01l)}, + tuple{T(40.0l), T(42.0l), T(6.4193141067359395324846555212831384747500e-01l)}, + tuple{T(40.0l), T(44.0l), T(7.4692196604585985651799646842674920219553e-01l)}, + tuple{T(40.0l), T(46.0l), T(8.3074626260026305803119722265978017795290e-01l)}, + tuple{T(40.0l), T(48.0l), T(8.9272353393334340495239865479281817396499e-01l)})}; const T s{x_yref.x[0]}; const T x{x_yref.x[1]}; const T y{trng::math::GammaP(s, x)}; @@ -225,16 +239,30 @@ TEMPLATE_TEST_CASE("special_functions", "", float, double, long double) { SECTION("GammaQ") { using tuple = arg_res_tuple; - auto x_yref{ - GENERATE(tuple{T(2.0l), T(0.0l), T(1.0000000000000000000000000000000000000000e+00l)}, - tuple{T(2.0l), T(1.0l), T(7.3575888234288464319104754032292173489162e-01l)}, - tuple{T(2.0l), T(2.0l), T(4.0600584970983807568199848491745321022289e-01l)}, - tuple{T(2.0l), T(3.0l), T(1.9914827347145577191736966260024710652680e-01l)}, - tuple{T(2.0l), T(4.0l), T(9.1578194443670901468590106366206211059560e-02l)}, - tuple{T(2.0l), T(5.0l), T(4.0427681994512802579816290538890545493098e-02l)}, - tuple{T(2.0l), T(6.0l), T(1.7351265236664508961316172015716675240545e-02l)}, - tuple{T(2.0l), T(7.0l), T(7.2950557244361296640250886752742610117898e-03l)}, - tuple{T(2.0l), T(8.0l), T(3.0191636511226065493925021320277491737981e-03l)})}; + auto x_yref{GENERATE( + tuple{T(2.0l), T(0.0l), T(1.0000000000000000000000000000000000000000e+00l)}, + tuple{T(2.0l), T(1.0l), T(7.3575888234288464319104754032292173489162e-01l)}, + tuple{T(2.0l), T(2.0l), T(4.0600584970983807568199848491745321022289e-01l)}, + tuple{T(2.0l), T(3.0l), T(1.9914827347145577191736966260024710652680e-01l)}, + tuple{T(2.0l), T(4.0l), T(9.1578194443670901468590106366206211059560e-02l)}, + tuple{T(2.0l), T(5.0l), T(4.0427681994512802579816290538890545493098e-02l)}, + tuple{T(2.0l), T(6.0l), T(1.7351265236664508961316172015716675240545e-02l)}, + tuple{T(2.0l), T(7.0l), T(7.2950557244361296640250886752742610117898e-03l)}, + tuple{T(2.0l), T(8.0l), T(3.0191636511226065493925021320277491737981e-03l)}, + tuple{T(12.5l), T(9.0l), T(8.4239071558046022179531765350857761710841e-01l)}, + tuple{T(12.5l), T(11.0l), T(6.3574402763862300241435796207925797940441e-01l)}, + tuple{T(12.5l), T(13.0l), T(4.0759869141744670098897165680884511303260e-01l)}, + tuple{T(12.5l), T(15.0l), T(2.2428900483440391398146224174972002551945e-01l)}, + tuple{T(12.5l), T(17.0l), T(1.0790785383480752168420919502760869634397e-01l)}, + tuple{T(40.0l), T(32.0l), T(9.0439718336976914129622842493855252349611e-01l)}, + tuple{T(40.0l), T(34.0l), T(8.2826497838409969720933718775904687135236e-01l)}, + tuple{T(40.0l), T(36.0l), T(7.2630360679539669332000368795942012021006e-01l)}, + tuple{T(40.0l), T(38.0l), T(6.0591936788110654890806747218276034566055e-01l)}, + tuple{T(40.0l), T(40.0l), T(4.7897113893894483665904310242244680684443e-01l)}, + tuple{T(40.0l), T(42.0l), T(3.5806858932640604675153444787168615252500e-01l)}, + tuple{T(40.0l), T(44.0l), T(2.5307803395414014348200353157325079780447e-01l)}, + tuple{T(40.0l), T(46.0l), T(1.6925373739973694196880277734021982204710e-01l)}, + tuple{T(40.0l), T(48.0l), T(1.0727646606665659504760134520718182603501e-01l)})}; const T s{x_yref.x[0]}; const T x{x_yref.x[1]}; const T y{trng::math::GammaQ(s, x)}; diff --git a/trng/math.hpp b/trng/math.hpp index 5c2cac8..9ab9694 100644 --- a/trng/math.hpp +++ b/trng/math.hpp @@ -193,6 +193,9 @@ namespace trng { return ::std::modf(x, &dummy); } + using ::std::copysign; + using ::std::signbit; + using ::std::isfinite; using ::std::isnan; using ::std::isnormal; diff --git a/trng/special_functions.hpp b/trng/special_functions.hpp index ef198dd..9258cbe 100644 --- a/trng/special_functions.hpp +++ b/trng/special_functions.hpp @@ -47,22 +47,27 @@ namespace trng { namespace math { + // --- error function and complementary error function-------------- + + using ::std::erf; + using ::std::erfc; + // --- x - ln(1 + x) ----------------------------------------------- namespace detail { template - T mln1p(T x) { - if (std::abs(x) >= T(1) / T(32)) - return x - ::std::log1p(x); + TRNG_CUDA_ENABLE T mln1p(T x) { + if (abs(x) >= T(1) / T(32)) + return x - ln1p(x); // use Taylor expansion for small arguments T y{0}; T x_to_the_n{x * x}; T sign{1}; - for (int n{2}; n < std::numeric_limits::digits; ++n) { + for (int n{2}; n < numeric_limits::digits; ++n) { const T delta{sign * x_to_the_n / n}; y += delta; - if (std::abs(delta) < 4 * std::numeric_limits::epsilon() * y) + if (abs(delta) < 4 * numeric_limits::epsilon() * y) break; x_to_the_n *= x; sign = -sign; @@ -129,13 +134,13 @@ namespace trng { namespace detail { template - T Beta(T x, T y) { - static const T ln_max{ln(std::numeric_limits::max())}; + TRNG_CUDA_ENABLE T Beta(T x, T y) { + static const T ln_max{ln(numeric_limits::max())}; if (x <= 0 or y <= 0) { #if !(defined TRNG_CUDA) errno = EDOM; #endif - return std::numeric_limits::signaling_NaN(); + return numeric_limits::signaling_NaN(); } const T z{x + y}; if (z * ln(z) - z > ln_max) @@ -181,10 +186,12 @@ namespace trng { // --- Pochhammer function ----------------------------------------- + TRNG_CUDA_ENABLE inline float Pochhammer(float x, float a) { return exp(ln_Gamma(x + a) - ln_Gamma(x)); } + TRNG_CUDA_ENABLE inline double Pochhammer(double x, double a) { return exp(ln_Gamma(x + a) - ln_Gamma(x)); } @@ -210,7 +217,7 @@ namespace trng { // by series expansion, see "Numerical Recipes" by W. H. Press et al., 3rd edition template TRNG_CUDA_ENABLE T GammaP_ser(T a, T x) { - const int itmax{64}; + const int itmax{numeric_limits::digits}; const T eps{4 * numeric_limits::epsilon()}; if (x < eps) return T{0}; @@ -243,7 +250,7 @@ namespace trng { // by continued fraction, see "Numerical Recipes" by W. H. Press et al., 3rd edition template TRNG_CUDA_ENABLE T GammaQ_cf(T a, T x) { - const T itmax{64}; + const T itmax{numeric_limits::digits}; const T eps{4 * numeric_limits::epsilon()}; const T min{4 * numeric_limits::min()}; // set up for evaluating continued fraction by modified Lentz's method @@ -272,6 +279,192 @@ namespace trng { return exp(-x + a * ln(x)) * h; } + + template + class GammaPQ_asympt_coefficients; + + template<> + class GammaPQ_asympt_coefficients { + public: + static constexpr std::size_t n = 15; + + TRNG_CUDA_ENABLE + float operator[](std::size_t index) const { + // clang-format off + static constexpr float d[n]{ + 1.f, // 1 + -3.333333333333333333333333333333333333e-01f, // -1 / 3 + 8.333333333333333333333333333333333333e-02f, // 1 / 12 + -1.481481481481481481481481481481481481e-02f, // -2 / 135 + 1.157407407407407407407407407407407407e-03f, // 1 / 864 + 3.527336860670194003527336860670194004e-04f, // 1 / 2835 + -1.787551440329218106995884773662551440e-04f, // -139 / 777600 + 3.919263178522437781697040956300215559e-05f, // 1 / 25515 + -2.185448510679992161473642955124436606e-06f, // -571 / 261273600 + -1.854062210715159960701798836229563253e-06f, // -281 / 151559100 + 8.296711340953086005016242131664432272e-07f, // 163879 / 197522841600 + -1.766595273682607930436005424574240304e-07f, // -5221 / 29554024500 + 6.707853543401498580369397100296135722e-09f, // 5246819 / 782190452736000 + 1.026180978424030804257395732272529509e-08f, // 5459 / 531972441000 + -4.382036018453353186552974622447191234e-09f // -534703531 / 122021710626816000 + }; + // clang-format on + return d[index]; + } + }; + + template<> + class GammaPQ_asympt_coefficients { + public: + static constexpr std::size_t n = 27; + + TRNG_CUDA_ENABLE + double operator[](std::size_t index) const { + // clang-format off + static constexpr double d[n]{ + 1., // 1 + -3.333333333333333333333333333333333333e-01, // -1 / 3 + 8.333333333333333333333333333333333333e-02, // 1 / 12 + -1.481481481481481481481481481481481481e-02, // -2 / 135 + 1.157407407407407407407407407407407407e-03, // 1 / 864 + 3.527336860670194003527336860670194004e-04, // 1 / 2835 + -1.787551440329218106995884773662551440e-04, // -139 / 777600 + 3.919263178522437781697040956300215559e-05, // 1 / 25515 + -2.185448510679992161473642955124436606e-06, // -571 / 261273600 + -1.854062210715159960701798836229563253e-06, // -281 / 151559100 + 8.296711340953086005016242131664432272e-07, // 163879 / 197522841600 + -1.766595273682607930436005424574240304e-07, // -5221 / 29554024500 + 6.707853543401498580369397100296135722e-09, // 5246819 / 782190452736000 + 1.026180978424030804257395732272529509e-08, // 5459 / 531972441000 + -4.382036018453353186552974622447191234e-09, // -534703531 / 122021710626816000 + 9.147699582236790234182488176331136808e-10, // 91207079 / 99704934754425000 + -2.551419399494624976687795379938870131e-11, // -4483131259 / 175711263302615040000 + -5.830772132550425067464089450400357975e-11, // -2650986803 / 45465450248017800000 + 2.436194802066741624369406967077899429e-11, // 432261921612371 / 17743323368298066739200000 + -5.027669280114175589090549859257443655e-12, // -6171801683 / 1227567156696480600000 + 1.100439203195613477083741744972934113e-13, // 6232523202521089 / 56636688191607429031526400000 + 3.371763262400985378827698841692001848e-13, // 4283933145517 / 12705320071808574210000000 + -1.392388722418162065919366184895799799e-13, // -25834629665134204969 / 185541790515705937507280486400000 + 2.853489380704744320396690990528282989e-14, // 11963983648109 / 419275562369682948930000000 + -5.139111834242572618990645803004942055e-16, // -1579029138854919086429 / 3072572050940090325120564854784000000 + -1.975228829434944283539624015807109122e-15, // -208697624924077 / 105657441717160103130360000000 + 8.099521156704561334071156687025752553e-16 // 746590869962651602203151 / 921771615282027097536169456435200000000 + }; + // clang-format on + return d[index]; + } + }; + + template<> + class GammaPQ_asympt_coefficients { + public: + static constexpr std::size_t n = 42; + + long double operator[](std::size_t index) const { + // clang-format off + static constexpr long double d[n]{ + 1.l, // 1 + -3.333333333333333333333333333333333333e-01l, // -1 / 3 + 8.333333333333333333333333333333333333e-02l, // 1 / 12 + -1.481481481481481481481481481481481481e-02l, // -2 / 135 + 1.157407407407407407407407407407407407e-03l, // 1 / 864 + 3.527336860670194003527336860670194004e-04l, // 1 / 2835 + -1.787551440329218106995884773662551440e-04l, // -139 / 777600 + 3.919263178522437781697040956300215559e-05l, // 1 / 25515 + -2.185448510679992161473642955124436606e-06l, // -571 / 261273600 + -1.854062210715159960701798836229563253e-06l, // -281 / 151559100 + 8.296711340953086005016242131664432272e-07l, // 163879 / 197522841600 + -1.766595273682607930436005424574240304e-07l, // -5221 / 29554024500 + 6.707853543401498580369397100296135722e-09l, // 5246819 / 782190452736000 + 1.026180978424030804257395732272529509e-08l, // 5459 / 531972441000 + -4.382036018453353186552974622447191234e-09l, // -534703531 / 122021710626816000 + 9.147699582236790234182488176331136808e-10l, // 91207079 / 99704934754425000 + -2.551419399494624976687795379938870131e-11l, // -4483131259 / 175711263302615040000 + -5.830772132550425067464089450400357975e-11l, // -2650986803 / 45465450248017800000 + 2.436194802066741624369406967077899429e-11l, // 432261921612371 / 17743323368298066739200000 + -5.027669280114175589090549859257443655e-12l, // -6171801683 / 1227567156696480600000 + 1.100439203195613477083741744972934113e-13l, // 6232523202521089 / 56636688191607429031526400000 + 3.371763262400985378827698841692001848e-13l, // 4283933145517 / 12705320071808574210000000 + -1.392388722418162065919366184895799799e-13l, // -25834629665134204969 / 185541790515705937507280486400000 + 2.853489380704744320396690990528282989e-14l, // 11963983648109 / 419275562369682948930000000 + -5.139111834242572618990645803004942055e-16l, // -1579029138854919086429 / 3072572050940090325120564854784000000 + -1.975228829434944283539624015807109122e-15l, // -208697624924077 / 105657441717160103130360000000 + 8.099521156704561334071156687025752553e-16l, // 746590869962651602203151 / 921771615282027097536169456435200000000 + -1.652253121639816181915148202653511616e-16l, // -29320119130515566117 / 177455371374430493811049182600000000 + 2.530543009747888423270610900602673849e-18l, // 1511513601028097903631961 / 597308006702753559203437807770009600000000 + 1.168693973855957658882308765077934746e-17l, // 2700231121460756431181 / 231046893529508502941986035745200000000 + -4.770037049820484758221678040848165974e-18l, // -8849272268392873147705987190261 / 1855178938018082279529957487152872816640000000000 + 9.699126059056237124207096858985853544e-19l, // 10084288256532215186381 / 10397110208827882632389371608534000000000 + -1.293256553803817501044325677449629120e-20l, // -6208770108287283939483943525987 / 480088045177548944685317694066691261071360000000000 + -6.969230253185693380530546315855194520e-20l, // -6782242429223267933535073 / 97316951554628981439164518255878240000000000 + 2.835145432176936599923196187131687508e-20l, // 2355444393109967510921431436000087153 / 83080196394065199975683597593629056110920990720000000000 + -5.750982159007047500162666139520842346e-21, // -51748587106835353426330148693 / 8998217291595659510809468851493269705120000000000 + 6.792953783488914564614965664370431421e-23l, // 2346608607351903737647919577082115121863 / 34544745660652310149889239879430961530920947941376000000000000 + 4.182125426111335857807972609348226866e-22l, // 7007277101869903281324331583 / 16755301163660883227024528206228847037120000000000 + -1.697153962004760373219505695058247317e-22l, // -2603072187220373277150999431416562396331667 / 15337867073329625706550822506467346919728900885970944000000000000 + 3.436215938394319882960425642469428583e-23l, // 585302872633292617248814587726421 / 17033355386471835584177000251811214753701006400000000000 + -3.643995779628021011969396867845457286e-25l, // -73239727426811935976967471475430268695630993 / 200987410128911415258641978124748114036127517209763250176000000000000 + -2.522535663578433775881227291922899288e-24l // -110855495796575034381969281033555329 / 43946056897097335807176660649672934064548596512000000000000 + }; + // clang-format on + return d[index]; + } + }; + + + template + TRNG_CUDA_ENABLE T GammaPQ_asympt_R(T a, T eta, T eta_squard_half) { + GammaPQ_asympt_coefficients coeffs; + constexpr std::size_t n{coeffs.n - 1}; + T beta[n]; + beta[n - 1] = coeffs[n]; + beta[n - 2] = coeffs[n - 1]; + for (std::ptrdiff_t i{n - 3}; i >= 0; --i) + beta[i] = beta[i + 2] * (i + 2) / a + coeffs[i + 1]; + T eta_to_the_i{1}; + T y{0}; + for (std::size_t i{0}; i < n; ++i) { + const T y_old{y}; + y += beta[i] * eta_to_the_i; + if (y == y_old) + break; + eta_to_the_i *= eta; + } + return y / (1 + beta[1] / a) * exp(-a * eta_squard_half) / sqrt(a) * + constants::one_over_sqrt_2pi; + } + + + // calculate regularized incomplete the incomplete Gamma function by an asymptotic + // expansion, see + // SIAM Journal on Scientific Computing 34(6) (2012), A2965-A2981 + // https://doi.org/10.48550/arXiv.1306.1754 and references therein + template + TRNG_CUDA_ENABLE T GammaP_asympt(T a, T x) { + const T mu{(x - a) / a}; + const T eta_squared_half{mln1p(mu)}; + const T eta{copysign(sqrt(2 * eta_squared_half), mu)}; + const T leading{erfc(-eta * sqrt(a / 2)) / 2}; + const T correction{-GammaPQ_asympt_R(a, eta, eta_squared_half)}; + return leading + correction; + } + + + // calculate regularized complementary incomplete the incomplete Gamma function by an + // asymptotic expansion, see + // SIAM Journal on Scientific Computing 34(6) (2012), A2965-A2981 + // https://doi.org/10.48550/arXiv.1306.1754 and references therein + template + TRNG_CUDA_ENABLE T GammaQ_asympt(T a, T x) { + const T mu{(x - a) / a}; + const T eta_squared_half{mln1p(mu)}; + const T eta{copysign(sqrt(2 * eta_squared_half), mu)}; + const T leading{erfc(eta * sqrt(a / 2)) / 2}; + const T correction{GammaPQ_asympt_R(a, eta, eta_squared_half)}; + return leading + correction; + } + + // P(a, x) and gamma(a, x) template TRNG_CUDA_ENABLE T GammaP(T a, T x) { @@ -282,6 +475,8 @@ namespace trng { #else if (by_Gamma_a) { #endif + if (a > 12 and x > T(3) / T(10) * a and x < T(235) / T(100) * a) + return GammaP_asympt(a, x); if (x < a + 1) return GammaP_ser(a, x); return 1 - GammaQ_cf(a, x); @@ -302,6 +497,8 @@ namespace trng { #else if (by_Gamma_a) { #endif + if (a > 12 and x > T(3) / T(10) * a and x < T(235) / T(100) * a) + return GammaQ_asympt(a, x); if (x < a + 1) return T{1} - GammaP_ser(a, x); return GammaQ_cf(a, x); @@ -397,7 +594,7 @@ namespace trng { // initial guess if (a > T{1}) { const T pp{p < T{1} / T{2} ? p : 1 - p}; - const T t = {sqrt(-2 * ln(pp))}; + const T t{sqrt(-2 * ln(pp))}; x = (T{2.30753} + t * T{0.27061}) / (1 + t * (T{0.99229} + t * T{0.04481})) - t; x = p < T{1} / T{2} ? -x : x; x = utility::max(T{1} / T{1000}, a * pow(1 - 1 / (9 * a) - x / (3 * sqrt(a)), T{3})); @@ -532,7 +729,7 @@ namespace trng { namespace detail { template - TRNG_CUDA_ENABLE inline T inv_Beta_I(T x, T p, T q, T norm) { + TRNG_CUDA_ENABLE T inv_Beta_I(T x, T p, T q, T norm) { if (x < numeric_limits::epsilon()) return 0; if (1 - x < numeric_limits::epsilon()) @@ -601,11 +798,6 @@ namespace trng { } #endif - // --- error function and complementary error function-------------- - - using ::std::erf; - using ::std::erfc; - // --- normal distribution function ------------------------------- TRNG_CUDA_ENABLE