[libc][math] Implement single precision lgamma function (#205490)
The implementation reduces to `|x|` and selects a path:
- **Tiny** ($|x| < 2^{-23}$): truncated Laurent series $\text{lgamma}(x)
= -\log|x| - \gamma\,|x|$.
- **Small** ($|x| < 0.66$): a centered-monomial minimax fit (degree 17
for $x>0$, degree 25 for $x<0$) of the smooth quotient $g(x) =
\frac{\text{lgamma}(x) + \log|x|}{x}$, recovering $\text{lgamma}(x) =
x\cdot g(x) - \log|x|$ with a single FMA rounding.
- **Medium** ($|x| \in [0.66, 3.37]$): centered-monomial minimax fits
split into three sub-ranges that factor out the function's roots —
- $[0.66, 1)$: $\text{lgamma}(t) = (t-1)\,P(t)$ (degree 14)
- $[1, 2)$: $\text{lgamma}(t) = (t-1)(t-2)\,P(t)$ (degree 18)
- $[2, 3.37)$: $\text{lgamma}(t) = (t-2)\,P(t)$ (degree 15)
- **Stirling** ($|x| > 3.37$): $(x-0.5)(\log x - 1) +
\frac{\log(2\pi)}{2} - \tfrac12$ plus a $1/x$ Bernoulli correction whose
term count is selected by sub-range (degree-10 residual poly on $(3.37,
73]$; 4-term, then 2-term; dropped entirely for $|x| > 2^{20}$).
Evaluated in plain `double` (≈$2^{-48}$ relative error).
For **negative** $x$:
- $|x| \le 3.37$ uses the Γ-recurrence — e.g. $\text{lgamma}(x) =
x(x+1)\,P_{M2}(x+0.5) - \log|x(x+1)|$ on $(-2, -0.66]$ (the recurrence's
polynomial prefactor is exactly the root factor, so relative accuracy is
preserved through lgamma's zeros, and the whole $\sin$ chain is
avoided). Narrow Taylor windows handle the two regular zeros in $(-3.37,
-2)$.
- $|x| > 3.37$ uses the reflection formula $\text{lgamma}(x) = \log\pi -
\log|\sin(\pi x)| - \log|x| - \text{lgamma}(|x|)$, reusing the
already-computed $\log|x|$.
Exhaustive tests pass in all 4 rounding modes.
Benchmark
Per-range throughput and latency, cycles/call (rdtsc, min of 20 trials;
core-math `cr_lgammaf` vs this implementation, both `-O3
-march=native`):
| range | interval | thru cr → llvm | lat cr → llvm |
| ------------- | --------------------------------------- |
--------------- | --------------- |
| tiny_pos | $[0, 2^{-23}]$ | 22.8 → **16.7** | 43.0 → **36.8** |
| small_pos | $[2^{-23}, 0.66]$ | 23.0 → **21.1** | 42.8 → **39.7** |
| m1_pos | $[0.66, 1)$ | 14.7 → **12.8** | 34.4 → **26.2** |
| m2_pos | $[1, 2)$ | **14.9** → 14.9 | 34.5 → **32.0** |
| m3_pos | $[2, 3.37)$ | 14.9 → **13.7** | 34.7 → **26.4** |
| stir_cheb_pos | $[3.37, 73]$ | **19.9** → 20.4 | **38.2** → 41.2 |
| bern4_pos | $[73, 2000]$ | 25.2 → **18.5** | **39.7** → 39.7 |
| bern2_pos | $[2000, 2\!\times\!10^7]$ | 26.0 → **25.0** | 44.4 →
**45.5** |
| huge_pos | $[2\!\times\!10^7, 4\!\times\!10^{36}]$ | **16.7** → 18.7 |
**35.5** → 38.7 |
| tiny_neg | $[-2^{-23}, 0]$ | 23.7 → **16.9** | 43.2 → **36.9** |
| small_neg | $[-0.66, -2^{-23}]$ | 23.6 → **22.2** | 43.6 → **41.1** |
| m1_neg | $[-1, -0.66]$ | 29.2 → **21.7** | 67.9 → **41.5** |
| m2_neg | $[-2, -1]$ | 29.3 → **23.3** | 67.9 → **45.6** |
| m3_neg | $[-3.37, -2]$ | 28.5 → **24.3** | 67.3 → **50.6** |
| stir_cheb_neg | $[-73, -3.37]$ | **32.5** → 39.6 | **63.5** → 73.0 |
| bern4_neg | $[-2000, -73]$ | **35.9** → 37.5 | **68.0** → 71.9 |
| bern2_neg | $[-10^5, -2000]$ | **28.5** → 35.8 | **61.3** → 69.8 |
GitOrigin-RevId: 665bc50f894248b4f676ff23381456683f2338c7
28 files changed