Use the NSWC GAMLN1 rational approximation in the incomplete beta routines - #543
Merged
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #543 +/- ##
==========================================
+ Coverage 94.49% 94.50% +0.01%
==========================================
Files 14 14
Lines 3016 3023 +7
==========================================
+ Hits 2850 2857 +7
Misses 166 166
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…tines The incomplete beta port substituted loggamma1p (from the incomplete gamma implementation) where the original NSWC routines call GAMLN1. loggamma1p is more accurate but roughly eight times more expensive: an 18-term Chebyshev recurrence plus a log1p, versus a single rational function evaluation. Within the incomplete beta routines the extra accuracy does not affect the result: the argument is always a0 = min(a, b) in (0, 1), where |log(gamma(1 + a))| <= 0.1215, so even a few ulps of relative error is an absolute error below 2e-16, and the value is only consumed through exp(z - u). Test suite results are unchanged. This makes beta_inc about 25% faster on the min(a, b) < 1 branches.
andreasnoack
force-pushed
the
an/beta-inc-gamln1
branch
from
August 17, 2026 12:24
71f707e to
7240ec0
Compare
devmotion
reviewed
Aug 17, 2026
devmotion
left a comment
devmotion
left a comment
Member
There was a problem hiding this comment.
LGTM apart from the comments below. I compared the coefficients and implementation with the Fortran source, and everything seems to match.
Review feedback: replace the @evalpoly macro with the evalpoly function (coefficients passed as tuples, since evalpoly has no vararg method), and add tests backing the docstring claims: gamln1 matches a 256-bit loggamma(1 + a) reference to within 2e-16 absolute error on [0, 1] (5e-16 on the full [-0.2, 1.25] domain), and substituting gamln1 for loggamma1p perturbs exp(z - u) by less than one ulp.
devmotion
reviewed
Aug 17, 2026
Docstring: the domain argument now rests on beta_inc's input validation rather than on "all call sites" (beta_integrand is also reachable from ncbeta with unvalidated parameters, harmlessly since beta_inc throws a DomainError for the same input on the next line), the perturbation claim is stated as "at most about one ulp", and the measured ulp errors and cost ratio are quoted. Tests: assert both the absolute bound (binding near the minimum of loggamma(1 + a)) and a 16/32 ulp bound (binding near the zeros at a = 0 and a = 1), add log-spaced grids for a -> 0 and a -> 1, compare directly against the BigFloat reference, loosen the full-domain absolute bound from 5e-16 to a bit-portable 8e-16 (FMA contraction in evalpoly is target-dependent), and drop the expm1 cross-check against loggamma1p, which asserted agreement with another approximation rather than with the true value.
devmotion
approved these changes
Aug 17, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Second follow-up from JuliaStats/StatsFuns.jl#228 (companion to the truncation PR; the two are independent and conflict only trivially in one hunk).
Where the original NSWC
BPSER/BRCOMProutines callGAMLN1— a pair of minimax rational approximations oflog(Γ(1+a))on-0.2 ≤ a ≤ 1.25— thebeta_inctranslation substitutedloggamma1pfrom the incomplete gamma implementation (Gil–Segura–Temme: an 18-term Chebyshev recurrence viaauxgamplus alog1p).loggamma1pis the more accurate function, but it is ~8× more expensive (30 ns vs 4 ns measured), and inside the incomplete beta routines the extra accuracy cannot reach the final result. This PR portsGAMLN1from the NSWC source (https://github.com/jacobwilliams/nswc) and uses it at the fourbeta_inc.jlcall sites.loggamma1pitself is untouched and remains in use bygamma_inc/gamma_inc_inv, where its wider domain and higher relative accuracy matter.Why it is safe here (also documented in the docstring):
a0 = min(a, b)inside ana0 < 1branch, with zero arguments handled earlier, so the argument is always in(0, 1)— well insideGAMLN1's designed domain[-0.2, 1.25]and away from the negative edge where its error grows.(0, 1),|log Γ(1+a)| ≤ 0.1215. A 200k-point sweep against a 256-bit reference gives max 8.0 ulp / mean 2.6 ulp forgamln1(vs 4.6 / 0.9 forloggamma1p) — but relative to a value bounded by 0.1215, 8 ulps is an absolute error below2e-16. The result is only ever consumed throughexp(z - u), where an absolute perturbation of that size changes the final incomplete beta value by less than one ulp.Consistent with that, all results in the test suite and in spot checks are bit-identical before and after.
Benchmarks (Apple Silicon,
@belapsed, this branch alone vs current master):beta_inc(0.5, 2.5, 0.3103)StatsFuns.tdistcdf(5, 1.5)StatsFuns.tdistinvcdf(5, 1e-8)Disclosure: this PR was prepared by Claude Code at my direction; I have reviewed the changes, the accuracy sweep, and the benchmark methodology.
🤖 Generated with Claude Code