`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

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