Improve numerical stability and domain coverage of gamma_approx and ln_gamma_approx

[jzeuzs]

We introduce domain checks, pole handling, and Euler's reflection formula to both gamma_approx and ln_gamma_approx. These changes prevent integer overflows and avoid NaN crashes on valid inputs.

Previously, ln_gamma_approx blindly evaluated the Lanczos approximation for all inputs. Lanczos breaks down for $z \lt 0.5$ and fails completely for negative numbers, causing severe precision loss or NaN crashes.

And, at exact negative integers (e.g., $z = -1.0$), (PI * z).sin() yields a tiny non-zero float rather than exactly 0.0. This caused gamma_approx to return massive incorrect finite numbers instead of failing mathematically, and ln_gamma_approx to return finite bounds instead of infinity.

Further, the fast-path for exact integers utilized the factorial function. Because factorial returns an integer, this caused a hard overflow/panic for any $z > 21$.

The following changes have been made:

• ln_gamma_approx

  • We added a check for exact integers $\leq 0$ and returned infinity.
  • Implemented the logarithmic Euler reflection formula for $z \lt 0.5$.
  • Applied .abs() to the sine calculation in the reflection formula.

• gamma_approx

  • Added a check for $z = 0.0$ (returns infinity) and exact negative integers (returns NaN, as standard limits diverge from left vs. right).
  • Removed the factorial(z_int - 1) shortcut. By falling back to ln_gamma_approx(z).exp(), the function can now safely evaluate exact integers up to $z \approx 171.6$ using floating-point math without panicking.

[jzeuzs]

I reimplemented the integer fast-path for gamma_approx, but calculating via f64 iteratively rather than falling back to an integer factorial function, since $\Gamma(22) = 51090942171657870000$ on the previous implementation but $\Gamma(22) = 51090942171709440000$ on Wolfram, possibly due to approximation error. This has now been rectified.

The fix prevents usize panics for $z \gt 21$, allows calculating massive factorials up to the float limit $z \approx 171$ and bypasses the Lanczos approximation error for exact integers.

[Axect]

Thanks for the improvement. The pole handling and the log-reflection formula look correct, and the added tests cover the tricky sign cases well.

One thing I'd like to fix before merging.

On the integer fast-path (src/special/lanczos.rs).

if z > 0.0 && z.fract() == 0.0 {
    let n = (z - 1.0) as u64;
    for i in 2..=n {
        result *= i as f64;
    }
    return result;
}

• Any f64 with magnitude above 2^52 is already integral, so z.fract() == 0.0 is true for huge values too, and they fall into this loop.
• result overflows to INFINITY around i = 171, but the loop keeps running up to n. So gamma(1e16) spins ~1e16 iterations, and 1e300 as u64 saturates to u64::MAX, which effectively never returns.
• Since public gamma() reaches this path with no magnitude guard, this is triggerable from outside as a hang.

Could you add a magnitude guard before the loop? Something like returning INFINITY for z > 171.0 (since gamma(171) is finite and gamma(172) already overflows) would keep the exact factorial for small integers and avoid the loop for large inputs.

Everything else looks good to me. Thanks again.