Dear GSL maintainer, In our team we found some errors in the GSL incomplete gamma function. These are huge and manifest when you call gsl_sf_gamma_inc_Q_e() and gsl_sf_gamma_inc_P_e() with values for (a, x) around 950000 each.
Specifically when the variable "a" values slightly higher than "(x-a)^2" and where "10e5 < x < 10e6". I attach a patch that fixes this problem, as well as a patch that adds tests to exercise the problem (and acknowledge the fix): 0001-added-tests-for-incomplete-gamma-near-x-a-2-a.txt and 0002-fixed-errors-in-the-incomplete-gamma-function-at-x-a.txt After that there are a few more patches. Patch 0003 tests for other edge cases similar to the ones that had significant errors. Patch 0004 tests for errors near P=0 as there seems to be some floating point error happening there and patch 0005 fixes the error bounds of the current implementation in this area without actually increasing its accuracy. The test data was generated using SciPy's implementation of the incomplete gamma function and Fortran code from this paper: https://arxiv.org/pdf/1306.1754 . If there are any questions please let me know but I really do suggest adding at least the first 2 patches as the error is truly huge. Thanks, Noah Unger-Schulz
From db4b66874b78c0c91f7a9a842278908f5819d10d Mon Sep 17 00:00:00 2001 From: Noah Unger-Schulz <[email protected]> Date: Mon, 21 Jul 2025 13:10:06 -0600 Subject: [PATCH 3/5] added tests for incomplete gamma function near other edge cases --- specfunc/test_gamma.c | 45 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/specfunc/test_gamma.c b/specfunc/test_gamma.c index 50173d44..9dbfa7c2 100644 --- a/specfunc/test_gamma.c +++ b/specfunc/test_gamma.c @@ -297,6 +297,51 @@ int test_gamma(void) TEST_SF(s, gsl_sf_gamma_inc_Q_e, (847000.0, 846160.0, &r), 0.8192892296443739, TEST_TOL0, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_Q_e, (978740.0, 977760.0, &r), 0.8390547330533766, TEST_TOL0, GSL_SUCCESS); + +/*The following tests test in the region right below (x-a)^2=a which currently have error around 10e-14. + If implementations change these values should be checked. */ + + TEST_SF(s, gsl_sf_gamma_inc_P_e, (858353.4, 857416.2, &r), 0.1558679723188377, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (903565.3, 902606.8, &r), 0.1566411958884375, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (862237.3, 861236.2, &r), 0.1404787341181084, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (915812.8, 914811.7, &r), 0.1477492968528649, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (915184.7, 914162.3, &r), 0.1425853107051554, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (884894.3, 883850.6, &r), 0.1335880944683347, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (965419.9, 964376.2, &r), 0.14405627117667602, TEST_TOL1, GSL_SUCCESS); + +/* The following tests test at x=0.2*a as this is another branch edge that + has slightly higher error than the surrounding area. */ + TEST_SF(s, gsl_sf_gamma_inc_P_e, (104.4, 20.00, &r), 2.5824266399145001e-040, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (226.0, 50.00, &r), 8.0563758802602748e-074, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (250.0, 50.00, &r), 4.1168220544357860e-090, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (477.4, 110.0, &r), 3.9534827781417115e-147, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (525.8, 110.0, &r), 4.7703402616217646e-179, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (572.4, 120.0, &r), 2.5991048520232467e-194, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (633.6, 120.0, &r), 3.0460087195879340e-237, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (659.4, 140.0, &r), 1.1947328166286074e-220, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (783.0, 150.0, &r), 1.6541772234483730e-289, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (696.0, 160.0, &r), 4.8537373853443579e-214, TEST_TOL1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (774.4, 160.0, &r), 5.5334098010985248e-266, TEST_TOL1, GSL_SUCCESS); + +/* The following tests test at x=0.5*a with similar reasons to the tests above. */ + + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7125, 19.6625, &r), 9.4942649542195976e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7125, 19.6500, &r), 9.3745107453207659e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (37.3375, 18.9125, &r), 1.1816285005746677e-4, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (37.0250, 18.7250, &r), 1.2203993541570113e-4, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7750, 19.5375, &r), 7.9891762227154276e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7125, 19.3500, &r), 6.8804750433475805e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (39.3375, 19.6625, &r), 6.0533918681359682e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7750, 19.3500, &r), 6.5740992245717479e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (37.9625, 18.9250, &r), 7.6144060398413831e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7125, 19.2875, &r), 6.4437279980465698e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.0000, 18.9125, &r), 7.3112985566458930e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (39.4000, 19.5375, &r), 5.0695437496483094e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7125, 19.1625, &r), 5.6448572778536202e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (37.9625, 18.7250, &r), 6.1565788962559479e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7750, 19.1000, &r), 5.0411017292971622e-5, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (39.2000, 19.2875, &r), 4.4988168191115288e-5, TEST_TOL0, GSL_SUCCESS); + /* non-normalized "Q" function */ TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0/1048576.0, 1.0/1048576.0, &r), 13.285819596290624271, TEST_TOL0, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_e, (-0.001, 1.0/1048576.0, &r), 13.381275128625328858, TEST_TOL0, GSL_SUCCESS); -- 2.50.1
From cce89ca01495dc334d10da0b6ae91d22e4966d5f Mon Sep 17 00:00:00 2001 From: Noah Unger-Schulz <[email protected]> Date: Mon, 21 Jul 2025 13:06:17 -0600 Subject: [PATCH 2/5] fixed errors in the incomplete gamma function at (x-a)^2>a --- specfunc/gamma_inc.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/specfunc/gamma_inc.c b/specfunc/gamma_inc.c index 99027516..18a359fc 100644 --- a/specfunc/gamma_inc.c +++ b/specfunc/gamma_inc.c @@ -524,7 +524,7 @@ gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result) result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); return stat_P; } - else if(a >= 1.0e+06 && (x-a)*(x-a) < a) { + else if(a >= 1.0e+04 && (x-a)*(x-a) < a) { /* Then try the difficult asymptotic regime. * This is the only way to do this region. */ @@ -593,7 +593,7 @@ gsl_sf_gamma_inc_P_e(const double a, const double x, gsl_sf_result * result) */ return gamma_inc_P_series(a, x, result); } - else if(a > 1.0e+06 && (x-a)*(x-a) < a) { + else if(a > 1.0e+04 && (x-a)*(x-a) < a) { /* Crossover region. Note that Q and P are * roughly the same order of magnitude here, * so the subtraction is stable. -- 2.50.1
From 5c2af8f247dfc3a4b6771160590a132851c0b2ec Mon Sep 17 00:00:00 2001 From: Noah Unger-Schulz <[email protected]> Date: Mon, 21 Jul 2025 13:15:45 -0600 Subject: [PATCH 5/5] fixed errors in the incomplete gamma function error estimate near P=0 --- specfunc/gamma_inc.c | 1 + 1 file changed, 1 insertion(+) diff --git a/specfunc/gamma_inc.c b/specfunc/gamma_inc.c index 18a359fc..69f21469 100644 --- a/specfunc/gamma_inc.c +++ b/specfunc/gamma_inc.c @@ -140,6 +140,7 @@ gamma_inc_P_series(const double a, const double x, gsl_sf_result * result) result->val = D.val * sum; result->err = D.err * fabs(sum) + fabs(D.val * remainder); + result->err += fabs(GSL_DBL_MIN * remainder); result->err += (1.0 + n) * GSL_DBL_EPSILON * fabs(result->val); if(n == nmax && fabs(remainder/sum) > GSL_SQRT_DBL_EPSILON) -- 2.50.1
From ea81a4922db1746b472906b5670632bc70f5a5ec Mon Sep 17 00:00:00 2001 From: Noah Unger-Schulz <[email protected]> Date: Mon, 21 Jul 2025 13:12:50 -0600 Subject: [PATCH 4/5] added tests for incomplete gamma function for spots near P=0 --- specfunc/test_gamma.c | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/specfunc/test_gamma.c b/specfunc/test_gamma.c index 9dbfa7c2..cb13afce 100644 --- a/specfunc/test_gamma.c +++ b/specfunc/test_gamma.c @@ -298,7 +298,7 @@ int test_gamma(void) TEST_SF(s, gsl_sf_gamma_inc_Q_e, (978740.0, 977760.0, &r), 0.8390547330533766, TEST_TOL0, GSL_SUCCESS); -/*The following tests test in the region right below (x-a)^2=a which currently have error around 10e-14. + /*The following tests test in the region right below (x-a)^2=a which currently have error around 10e-14. If implementations change these values should be checked. */ TEST_SF(s, gsl_sf_gamma_inc_P_e, (858353.4, 857416.2, &r), 0.1558679723188377, TEST_TOL1, GSL_SUCCESS); @@ -309,7 +309,7 @@ int test_gamma(void) TEST_SF(s, gsl_sf_gamma_inc_P_e, (884894.3, 883850.6, &r), 0.1335880944683347, TEST_TOL1, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_P_e, (965419.9, 964376.2, &r), 0.14405627117667602, TEST_TOL1, GSL_SUCCESS); -/* The following tests test at x=0.2*a as this is another branch edge that + /* The following tests test at x=0.2*a as this is another branch edge that has slightly higher error than the surrounding area. */ TEST_SF(s, gsl_sf_gamma_inc_P_e, (104.4, 20.00, &r), 2.5824266399145001e-040, TEST_TOL1, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_P_e, (226.0, 50.00, &r), 8.0563758802602748e-074, TEST_TOL1, GSL_SUCCESS); @@ -323,7 +323,7 @@ int test_gamma(void) TEST_SF(s, gsl_sf_gamma_inc_P_e, (696.0, 160.0, &r), 4.8537373853443579e-214, TEST_TOL1, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_P_e, (774.4, 160.0, &r), 5.5334098010985248e-266, TEST_TOL1, GSL_SUCCESS); -/* The following tests test at x=0.5*a with similar reasons to the tests above. */ + /* The following tests test at x=0.5*a with similar reasons to the tests above. */ TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7125, 19.6625, &r), 9.4942649542195976e-5, TEST_TOL0, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7125, 19.6500, &r), 9.3745107453207659e-5, TEST_TOL0, GSL_SUCCESS); @@ -342,6 +342,25 @@ int test_gamma(void) TEST_SF(s, gsl_sf_gamma_inc_P_e, (38.7750, 19.1000, &r), 5.0411017292971622e-5, TEST_TOL0, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_P_e, (39.2000, 19.2875, &r), 4.4988168191115288e-5, TEST_TOL0, GSL_SUCCESS); + + /*The following tests test numerical imprecision near floating point round off. */ + + TEST_SF(s, gsl_sf_gamma_inc_P_e, (27965.40, 22020.00, &r), 1.4821969375237396e-323, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (102852.6, 91020.00, &r), 3.9525251667299724e-323, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (118742.4, 106020.0, &r), 1.0869444208507424e-321, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (139882.2, 126020.0, &r), 2.6679544875427313e-322, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (167222.0, 152020.0, &r), 2.5691413583744820e-322, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (203851.8, 187020.0, &r), 3.3596463917204765e-322, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (254901.6, 236020.0, &r), 2.7173610521268560e-322, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (329581.4, 308020.0, &r), 4.9406564584124654e-323, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (440981.2, 416020.0, &r), 1.3389179002297781e-321, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (442041.2, 417020.0, &r), 2.3221085354538588e-322, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (627921.0, 598020.0, &r), 4.4959973771553436e-322, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (628971.0, 599020.0, &r), 1.2845706791872410e-322, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (967220.8, 930020.0, &r), 1.3290365873129532e-321, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (968260.8, 931020.0, &r), 6.0770074438473325e-322, 1, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (969300.8, 932020.0, &r), 2.7667676167109806e-322, 1, GSL_SUCCESS); + /* non-normalized "Q" function */ TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0/1048576.0, 1.0/1048576.0, &r), 13.285819596290624271, TEST_TOL0, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_e, (-0.001, 1.0/1048576.0, &r), 13.381275128625328858, TEST_TOL0, GSL_SUCCESS); -- 2.50.1
From dc91bef68238c59c98c50c6ddad1e7cdd713b61d Mon Sep 17 00:00:00 2001 From: Noah Unger-Schulz <[email protected]> Date: Mon, 21 Jul 2025 12:56:13 -0600 Subject: [PATCH 1/5] added tests for incomplete gamma near (x-a)^2<~a --- specfunc/test_gamma.c | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/specfunc/test_gamma.c b/specfunc/test_gamma.c index 97993553..50173d44 100644 --- a/specfunc/test_gamma.c +++ b/specfunc/test_gamma.c @@ -266,6 +266,37 @@ int test_gamma(void) TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0e+06, 1.0e+06-2.0, &r), 0.50066490399940144811, TEST_TOL2, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_Q_e, (1.0e+07, 1.0e+07-2.0, &r), 0.50021026104978614908, TEST_TOL2, GSL_SUCCESS); + + /*The following tests test in the region (x-a)^2<~a because the continued fraction method + * doesn't always converge for high values when too close to this curve*/ + + + TEST_SF(s, gsl_sf_gamma_inc_P_e, (970100.0, 969120.0, &r), 0.15987184308570648, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (967420.0, 966440.0, &r), 0.15953692653917345, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (872880.0, 872060.0, &r), 0.19008070116271375, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (930480.0, 930020.0, &r), 0.316820727586329, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (832460.0, 832780.0, &r), 0.6372233296338622, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (825420.0, 824640.0, &r), 0.19532563380564685, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (969780.0, 969000.0, &r), 0.21419964517274928, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (906220.0, 905420.0, &r), 0.20037785977318115, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (860920.0, 860000.0, &r), 0.16071454049177544, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (964180.0, 963360.0, &r), 0.20186159428235603, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (847000.0, 846160.0, &r), 0.18071077035562616, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_P_e, (978740.0, 977760.0, &r), 0.16094526694662334, TEST_TOL0, GSL_SUCCESS); + + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (970100.0, 969120.0, &r), 0.8401281569142935, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (967420.0, 966440.0, &r), 0.8404630734608265, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (872880.0, 872060.0, &r), 0.8099192988372863, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (930480.0, 930020.0, &r), 0.6831792724136709, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (832460.0, 832780.0, &r), 0.3627766703661378, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (825420.0, 824640.0, &r), 0.8046743661943532, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (969780.0, 969000.0, &r), 0.7858003548272507, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (906220.0, 905420.0, &r), 0.7996221402268189, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (860920.0, 860000.0, &r), 0.8392854595082245, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (964180.0, 963360.0, &r), 0.7981384057176439, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (847000.0, 846160.0, &r), 0.8192892296443739, TEST_TOL0, GSL_SUCCESS); + TEST_SF(s, gsl_sf_gamma_inc_Q_e, (978740.0, 977760.0, &r), 0.8390547330533766, TEST_TOL0, GSL_SUCCESS); + /* non-normalized "Q" function */ TEST_SF(s, gsl_sf_gamma_inc_e, (-1.0/1048576.0, 1.0/1048576.0, &r), 13.285819596290624271, TEST_TOL0, GSL_SUCCESS); TEST_SF(s, gsl_sf_gamma_inc_e, (-0.001, 1.0/1048576.0, &r), 13.381275128625328858, TEST_TOL0, GSL_SUCCESS); -- 2.50.1
