Skip to content

Commit

Permalink
improve numerical accuracy to incomplete gamma functions
Browse files Browse the repository at this point in the history 
to reach machine precision also for 128-bit floating point numbers, e.g., on ARM64
  • Loading branch information
@rabauke
rabauke committed Feb 25, 2024
1 parent 89a8060 commit 230a6f4
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 7 deletions.
29 changes: 29 additions & 0 deletions tests/test_special_functions.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,8 @@
#include <utility>
#include <limits>
#include <cmath>
#include <sstream>
#include <iomanip>
#include <ciso646>

#include <catch2/catch.hpp>
Expand Down Expand Up @@ -110,6 +112,13 @@ void check_function(const T_ref &x_yref, const T y) {
const std::tuple<T, T> 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<T>::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));
}

Expand Down Expand Up @@ -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)},
Expand Down Expand Up @@ -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)},
Expand Down
40 changes: 33 additions & 7 deletions trng/special_functions.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -358,7 +358,7 @@ namespace trng {
template<>
class GammaPQ_asympt_coefficients<long double> {
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
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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<T>::digits10 > 20) {
#else
if (numeric_limits<T>::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<T, true>(a, x);
return 1 - GammaQ_cf<T, true>(a, x);
Expand All @@ -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<T>::digits10 > 20) {
#else
if (numeric_limits<T>::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<T, true>(a, x);
return GammaQ_cf<T, true>(a, x);
Expand Down

0 comments on commit 230a6f4

Please sign in to comment.