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

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)

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).

Start by reproducing the listed large-argument cases and reading GammaQ in pytensor/scalar/c_code/gamma.c, including its MAXITER-based series and continued-fraction paths. Compare the C results with scipy.special.gammaincc, then implement and validate a robust treatment for large, similar a and x so the reported discrepancies and downstream Poisson.logcdf error are resolved.