From 444ae9761d3055430e64ca61c63a42d4f14a65d7 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen <4643818-rmlarsen1@users.noreply.gitlab.com> Date: Sat, 28 Feb 2026 08:52:44 -0800 Subject: [PATCH] Clamp igamma/igammac output to [0,1] for numerical stability libeigen/eigen!2229 Co-authored-by: Rasmus Munk Larsen --- .../Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h index 5553bbccf..8e88b9dd9 100644 --- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h +++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h @@ -1103,7 +1103,9 @@ struct igammac_impl { } if ((x < one) || (x < a)) { - return (one - igamma_series_impl::run(a, x)); + // Clamp to [0,1] since 1-igamma can produce tiny negative values + // due to floating-point cancellation for extreme arguments. + return numext::mini(one, numext::maxi(zero, one - igamma_series_impl::run(a, x))); } return igammac_cf_impl::run(a, x); @@ -1155,7 +1157,9 @@ struct igamma_generic_impl { if ((x > one) && (x > a)) { Scalar ret = igammac_cf_impl::run(a, x); if (mode == VALUE) { - return one - ret; + // Clamp to [0,1] since 1-igammac can produce tiny negative values + // due to floating-point cancellation for extreme arguments. + return numext::mini(one, numext::maxi(zero, one - ret)); } else { return -ret; }