diff --git a/tests/test_special_functions.cc b/tests/test_special_functions.cc index 85e107d..feaf41e 100644 --- a/tests/test_special_functions.cc +++ b/tests/test_special_functions.cc @@ -36,6 +36,8 @@ #include #include #include +#include +#include #include #include @@ -110,6 +112,13 @@ void check_function(const T_ref &x_yref, const T y) { const std::tuple y_min_max{bounds(y_ref)}; const T y_min{std::get<0>(y_min_max)}; const T y_max{std::get<1>(y_min_max)}; + std::stringstream stream; + stream << std::setprecision(std::numeric_limits::digits10 + 1) + << "args: " << x_yref.x << "\n" + << "lower bound: " << y_min << "\n" + << "upper bound: " << y_max << "\n" + << "actual : " << y; + INFO(stream.str()); REQUIRE((y_min <= y and y <= y_max)); } @@ -222,6 +231,16 @@ TEMPLATE_TEST_CASE("special_functions", "", float, double, long double) { 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(16.5l), T(12.0l), T(1.2627126996732821737464036539300809019108e-01l)}, + tuple{T(16.5l), T(14.0l), T(2.8557520780884416274501697877174390654836e-01l)}, + tuple{T(16.5l), T(16.0l), T(4.8326039268967179332861124572559089314860e-01l)}, + tuple{T(16.5l), T(18.0l), T(6.7011079005720201592981130181526523422923e-01l)}, + tuple{T(16.5l), T(20.0l), T(8.1276165004754884951657099239430743008078e-01l)}, + tuple{T(28.5l), T(24.0l), T(2.0368352137506526380354602347413851063195e-01l)}, + tuple{T(28.5l), T(26.0l), T(3.3732537440409877898190669708355415940881e-01l)}, + tuple{T(28.5l), T(28.0l), T(4.8738370793105653978256547357820063511336e-01l)}, + tuple{T(28.5l), T(30.0l), T(6.3247470119696109362049958384311443851676e-01l)}, + tuple{T(28.5l), T(32.0l), T(7.5566381072794801096984673931665933212126e-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)}, @@ -254,6 +273,16 @@ TEMPLATE_TEST_CASE("special_functions", "", float, double, long double) { 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(16.5l), T(12.0l), T(8.7372873003267178262535963460699190980892e-01l)}, + tuple{T(16.5l), T(14.0l), T(7.1442479219115583725498302122825609345164e-01l)}, + tuple{T(16.5l), T(16.0l), T(5.1673960731032820667138875427440910685140e-01l)}, + tuple{T(16.5l), T(18.0l), T(3.2988920994279798407018869818473476577077e-01l)}, + tuple{T(16.5l), T(20.0l), T(1.8723834995245115048342900760569256991922e-01l)}, + tuple{T(28.5l), T(24.0l), T(7.9631647862493473619645397652586148936805e-01l)}, + tuple{T(28.5l), T(26.0l), T(6.6267462559590122101809330291644584059119e-01l)}, + tuple{T(28.5l), T(28.0l), T(5.1261629206894346021743452642179936488664e-01l)}, + tuple{T(28.5l), T(30.0l), T(3.6752529880303890637950041615688556148324e-01l)}, + tuple{T(28.5l), T(32.0l), T(2.4433618927205198903015326068334066787874e-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)}, diff --git a/trng/special_functions.hpp b/trng/special_functions.hpp index 2587f01..d9b4230 100644 --- a/trng/special_functions.hpp +++ b/trng/special_functions.hpp @@ -358,7 +358,7 @@ namespace trng { template<> class GammaPQ_asympt_coefficients { public: - static constexpr std::size_t n = 42; + static constexpr std::size_t n = 50; long double operator[](std::size_t index) const { // clang-format off @@ -398,13 +398,21 @@ namespace trng { -1.293256553803817501044325677449629120e-20l, // -6208770108287283939483943525987 / 480088045177548944685317694066691261071360000000000 -6.969230253185693380530546315855194520e-20l, // -6782242429223267933535073 / 97316951554628981439164518255878240000000000 2.835145432176936599923196187131687508e-20l, // 2355444393109967510921431436000087153 / 83080196394065199975683597593629056110920990720000000000 - -5.750982159007047500162666139520842346e-21, // -51748587106835353426330148693 / 8998217291595659510809468851493269705120000000000 + -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 + -2.522535663578433775881227291922899288e-24l, // -110855495796575034381969281033555329 / 43946056897097335807176660649672934064548596512000000000000 + 1.021727557887676825282765833070338486e-24l, // 34856851734234401648335623107688675640839679447003 / 34115602995281423626001889366894744876492284771185214084874240000000000000 + -2.065618928289515596162988977857018414e-25l, // -18447986573777204063499607563765439 / 89309728532810714704907407126754672453760050976000000000000 + 1.987728212387035132762438621481561682e-27l, // 909773124599542506852275229422593983242880452145053 / 457694929784695579366441347746259897263020492490220832162672803840000000000000 + 1.528011309299919423610824908396958132e-26l, // 38650132745379700438031566826935471987259957 / 2529440227971075696123684974426321512851113061316034836800000000000000 + -6.179660368053257853966883274340069798e-27l, // -1527335577854677023023224272800947125313629267269390501 / 247155262083735612857878327782980344522031065944719249367843314073600000000000000 + 1.247824052529354935357309838028423545e-27l, // 217784448556937372678947372805330071920344629 / 174531375730004223032534263235416184386726801230806403739200000000000000 + -1.099129014345020826298528140114925449e-29l, // -183856455668177802003316143799518064719008299958634826921 / 16727468137827226278221205224352109717251062543138598797215635496501248000000000000000 + -9.289074058313414187516698493744758188e-29l // -1167289109751840227800236733417523750884898531 / 12566259052560304058342466952949965275844329688618061069222400000000000000 }; // clang-format on return d[index]; @@ -475,8 +483,17 @@ 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 __cplusplus >= 201703L + if constexpr (numeric_limits::digits10 > 20) { +#else + if (numeric_limits::digits10 > 20) { +#endif + if (a > 28 and x > T(3) / T(10) * a and x < T(235) / T(100) * a) + return GammaP_asympt(a, x); + } else { + 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); @@ -497,8 +514,17 @@ 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 __cplusplus >= 201703L + if constexpr (numeric_limits::digits10 > 20) { +#else + if (numeric_limits::digits10 > 20) { +#endif + if (a > 28 and x > T(3) / T(10) * a and x < T(235) / T(100) * a) + return GammaQ_asympt(a, x); + } else { + 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);