C implementation of `gammaincc` silently loses precision for large `a ≈ x`

[jacobpstein]

Description

The C kernel for gammaincc (GammaQ in pytensor/scalar/c_code/gamma.c) returns badly wrong values when both arguments are large and of similar magnitude, without any warning:

import pytensor
import pytensor.tensor as pt
from pytensor.compile.mode import Mode
import scipy.special as sp

a, x = pt.dscalars("a", "x")
f = pytensor.function([a, x], pt.gammaincc(a, x), mode=Mode(linker="c", optimizer=None))

f(1_000_001.0, 1_000_000.0)             # 0.6528316434
sp.gammaincc(1_000_001.0, 1_000_000.0)  # 0.5002659615 (correct)

a = x+1	C kernel	scipy	abs diff
1e4	0.5026595812	0.5026595812	2.5e-12
1e5	0.5014430600	0.5008410431	6.0e-4
3e5	0.5311209048	0.5004855769	3.1e-2
1e6	0.6528316434	0.5002659615	1.5e-1

Cause: gamma.c uses the Numerical Recipes approach (series for x < a+1, Lentz continued fraction otherwise) with MAXITER = 1024. Near a ≈ x both methods need O(√x) iterations to converge, so from roughly x ≳ 1e5 the loop hits the cap and the partial result is returned as-is. scipy switches to a uniform asymptotic (Temme) expansion in this regime.

Impact: any graph that evaluates gammaincc through the C backend gets wrong values in this regime, e.g. fused composites inside scan. Concretely, Poisson.logcdf in PyMC is off by ~1e-4 (relative, in probability space) at mu = 1e6, which shifts quantile searches by several integers. Confirmed on 3.0.7 and 3.1.3; the code is unchanged on main.

A minimal fix could be raising MAXITER for the affected range, but the more robust fix is an asymptotic expansion for the large-a ≈ x regime (as in scipy/cephes igam.c).