fix(math): regularize removable singularity in hyp2f1 when c-a-b is an integer

[Sumu004]

What

hyp2f1_small_argument returns nan for z near 1 whenever c - a - b is exactly an integer:

import tensorflow_probability as tfp
hyp2f1 = tfp.math.hypergeometric.hyp2f1_small_argument

H = 1
a = tf.constant(1.0)
b = tf.constant(0.5 - H)
c = tf.constant(H + 1.5)
x = tf.constant(0.9901961)
hyp2f1(a, b, c, x)  # nan, should be ~ 0.753603

Root cause

_hyp2f1_z_near_one computes d = c - a - b and evaluates the standard Gauss connection formula for z near 1, which needs both lgamma(d) and lgamma(-d). Whenever d is an integer of either sign, one of these two always lands on a nonpositive-integer pole of the Gamma function (lgamma(d) for d <= 0, lgamma(-d) for d >= 0), producing nan — even though 2F1 itself is finite and analytic at these points. It's a removable singularity in the formula, not in the function.

Why not the closed form?

The exact value at d integer is given by a distinct closed form involving digamma corrections and an extra log(1-z) term (DLMF 15.8.8 for d >= 0, 15.8.9 for d < 0). I initially hand-derived and implemented this directly, verified it exactly against mpmath at 50-digit precision for m=2 — and then found it silently diverged from mpmath (only 3 correct digits) for m=3 with a different parameter set, indicating a transcription error in my own derivation that a quick spot-check would have missed. Given how easy it is to introduce a subtle-but-wrong error in a delicate closed form like this, I moved to a different, independently-verifiable approach instead of debugging the series further:

Fix

Regularize by nudging c off the exact degenerate point before evaluating the (already-correct, well-tested) existing formula:

epsilon = numpy_dtype(np.sqrt(np.finfo(numpy_dtype).eps))
is_degenerate = tf.math.equal(d, tf.math.round(d))
c = tf.where(is_degenerate, c + epsilon, c)
d = tf.where(is_degenerate, d + epsilon, d)

epsilon = sqrt(machine epsilon) is the standard step size balancing truncation against rounding error in a one-sided numerical derivative. Since 2F1 is analytic in c away from its own poles, this recovers the true value at a controlled, first-order-in-epsilon rate — no new closed form to get subtly wrong, and it reuses the exact code path already validated by the surrounding test suite for non-degenerate d.

As a side effect, the same nudge transitively resolves the identical degeneracy in the two recursive _hyp2f1_internal(..., 1 - d, ...) / (..., d + 1, ...) calls a few lines down, whose own c parameters land on the exact same integer boundary (this recursion is what actually crashes today — dividing by an exact-zero series denominator — rather than surfacing as a "clean" nan from lgamma alone).

Verification

Checked against mpmath at 50-digit precision:

• The exact reported repro (m=2): matches to the reported expected value 0.753603006025111.
• m swept across {-4, ..., 4} (covering both DLMF 15.8.8 and 15.8.9 regimes) with hand-picked test cases: relative error consistently ~1e-8 to 1e-9, matching the expected O(epsilon) bound.
• 300 randomized (a, b, m, z) combinations at epsilon=1e-8: worst-case relative error 8.6e-7, still within the expected bound.

Checklist

• New testHyp2F1DegenerateCMinusAMinusB regression test: the exact reported repro plus 5 additional cases sweeping m = c-a-b through {-3, -1, 0, 1, 3}, checked against scipy.special.hyp2f1 at both float32 and float64
• Existing hyp2f1 tests (all currently-passing parameter ranges) untouched — the fix only activates when d is exactly an integer
• CLA signed

Fixes #2001

[maciejskorski]

Thanks for this patch!
Note that there is one more issue that could be possibly fixed
#1861