`tf.math.igamma` and `tf.math.igammac` return values outside [0, 1] for large arguments

Issue type

Bug

Have you reproduced the bug with TensorFlow Nightly?

Yes

Source

source

TensorFlow version

TensorFlow 2.22.0-dev20260421 (nightly)

Custom code

Yes

OS platform and distribution

No response

Mobile device

No response

Python version

No response

Bazel version

No response

GCC/compiler version

No response

CUDA/cuDNN version

No response

GPU model and memory

No response

Current behavior?

tf.math.igamma and tf.math.igammac return values outside [0, 1] for large arguments

Summary

tf.math.igamma(a, x) computes the regularized lower incomplete gamma function P(a, x), which is mathematically bounded to [0, 1] for all valid inputs (a > 0, x ≥ 0). For large arguments (a ≈ x ≥ ~5×10⁴ in float32), TF returns values well outside this range — up to 1900× the correct value.

tf.math.igammac (the complement Q = 1 - P) is similarly affected. The identity P(a,x) + Q(a,x) = 1 breaks down: at a = x = 800000, P + Q ≈ 3.03 instead of 1.0.

Scale of the problem

a = x	igamma	igammac	sum (should be 1.0)
100	0.5133	0.4867	1.0000
1,000	0.5044	0.4960	1.0003
10,000	0.5020	0.4993	1.0013
100,000	0.5195	0.5187	1.0382
800,000	1.4981	1.5362	3.0343
1,000,000	1.3122	1.3741	2.6863
8,685,113	1900.3	—	—

Degradation starts around a ≈ 5×10⁴ and gets worse with increasing values. At a = 8.7×10⁶, the result is off by a factor of ~3800.

Affected devices

Both CPU and GPU eager kernels return identical wrong results. XLA (jit_compile=True) uses a different implementation that stays in range longer but also degrades at very large values.

Expected behavior

For a = x, the regularized incomplete gamma P(a, a) → 0.5 as a → ∞ (by the central limit theorem applied to the Gamma distribution). The function should always return values in [0, 1], regardless of input magnitude.

Root cause hypothesis

The Eigen igamma implementation likely uses a series expansion (for x < a+1) and a continued fraction (for x ≥ a+1). For a = x at the boundary of these two regimes, the series may require an extremely large number of terms to converge, leading to accumulated error and eventual overflow of intermediate sums. A numerically stable implementation would use the asymptotic uniform expansion for large a ≈ x (see DLMF §8.12 or Temme's method).

Environment

• TensorFlow 2.22.0-dev20260421 (nightly)
• Also confirmed on TensorFlow 2.21.0 (stable)
• CPU: x86_64
• GPU: NVIDIA T4 (CUDA 12)
• Both CPU and GPU affected identically

Standalone code to reproduce the issue

## Reproduction


import tensorflow as tf

a = tf.constant(800000.0, dtype=tf.float32)
x = tf.constant(800000.0, dtype=tf.float32)

ig = tf.math.igamma(a, x)
igc = tf.math.igammac(a, x)

print(f"igamma  = {ig.numpy():.6f}")   # 1.498142 (should be ~0.5)
print(f"igammac = {igc.numpy():.6f}")  # 1.536196 (should be ~0.5)
print(f"sum     = {(ig + igc).numpy():.6f}")  # 3.034338 (should be 1.0)

Relevant log output