Clamp igamma/igammac output to [0,1] for numerical stability

libeigen/eigen!2229

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-02-28 08:52:44 -08:00
parent eddb470a09
commit 444ae9761d

View File

@@ -1103,7 +1103,9 @@ struct igammac_impl {
}
if ((x < one) || (x < a)) {
return (one - igamma_series_impl<Scalar, VALUE>::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<Scalar, VALUE>::run(a, x)));
}
return igammac_cf_impl<Scalar, VALUE>::run(a, x);
@@ -1155,7 +1157,9 @@ struct igamma_generic_impl {
if ((x > one) && (x > a)) {
Scalar ret = igammac_cf_impl<Scalar, mode>::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;
}