Repository navigation
Conversation
The continued fraction is not specific to Float64, so it can be reused for the BigFloat gamma_inc. The accuracy options of ind still apply to Float64, while other types iterate to eps(T). Float64 results are bit-identical (checked on 400k random arguments), and the speed is unchanged. For this, rgammax gets a BigFloat method, x^a e^(-x) / Γ(a), which cannot overflow in BigFloat's exponent range. Like the Float64 method it goes through Γ(a+1) for a < 1, where MPFR's Γ(a) is up to ~3x slower. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
As for gamma_inc_cf, the accuracy options of ind apply to Float64 only, and other types sum to eps(T). Float64 results are bit-identical (checked on 400k random arguments), and its speed and zero allocations are unchanged. The stopping rule (the next term is below the tolerance) is kept. Stopping on a bound for the remaining terms instead would change most Float64 results and make some of them less accurate, since at the Float64 tolerance it ends the sum earlier. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Negative a or x and a = x = 0 now throw a DomainError like the Float64
method, instead of returning NaN or meaningless values.
mpfr_gamma_inc and gamma(a) are both infinite for a = Inf, so the
BigFloat method returned (NaN, NaN), while the Float64 and Float32
methods return (0, 1) for finite x and (1, 0) for x = Inf. Match them.
Special values are handled in the same order as in the Float64 method,
and the tests check that a, x ∈ {0, 2, Inf, NaN} give the same results.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The BigFloat method obtained Q(a,x) from mpfr_gamma_inc and returned P(a,x) = 1 - Q(a,x), which cancels catastrophically whenever P is small: at 256 bits gamma_inc(big(10), big(0.1)) had a relative error of 1e-61, and gamma_inc(big(1), big(1e-300)) and gamma_inc(big(1e5), big(9e4)) returned P = 0 instead of 1e-300 and 2e-235. If Q >= 0.5, P is now computed by gamma_inc_taylor, whose terms are all positive. Q >= 0.5 implies x < a (the median of the gamma distribution is below its mean), so the terms decrease. If Q < 0.5, P = 1 - Q loses at most about one bit and is kept. The computation runs with 32 guard bits and is rounded once at the end. This also removes the error of up to ~1.9 ulp in Q from dividing the separately rounded Γ(a,x) and Γ(a). Against MPFR's 1 - Q evaluated with enough extra precision to absorb the cancellation, 300 random (a, x) at 64, 256 and 1000 bits now give at most 0.5 ulp for both P and Q. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
mpfr_gamma_inc is slow for large arguments: at 256 bits, gamma_inc(big(10.5), big(1e4)) took 630 ms and gamma_inc(big(1e6), big(999000)) 360 ms. Both ratios are now computed with the generic gamma_inc_taylor and gamma_inc_cf: - For x < max(a + 1, prec/16), P(a,x) from gamma_inc_taylor, and Q = 1 - P. If P ≈ 1 (small a, or x > a + 1), the subtraction loses bits, and the computation is repeated with that many more bits. - Otherwise, Q(a,x) from gamma_inc_cf, and P = 1 - Q ≥ 1/2. The continued fraction needs roughly prec^2/x iterations, and x ≈ prec/16 is about where it becomes cheaper than the series. Against MPFR's Γ(a,x)/Γ(a) evaluated with enough extra precision, 1000 random (a, x) at 64, 256, 1000 and 4000 bits give at most 0.5 ulp for both P and Q. This includes a up to 1e7 with x ≈ a, where the 32 guard bits absorb that the stopping rule of the series does not bound the remaining terms. At 256 bits, (10.5, 1e4) now takes 0.06 ms, (1e5, 1.1e5) 0.14 ms (was 180 ms), and (1e6, 999000) 8.4 ms. Most small arguments are faster as well; the largest slowdown is (3.5, 20) with 0.37 ms instead of 0.16 ms. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #561 +/- ##
=======================================
Coverage 94.67% 94.67%
=======================================
Files 14 14
Lines 3023 3045 +22
=======================================
+ Hits 2862 2883 +21
- Misses 161 162 +1
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:
|
This branch has not been deployed
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
1 participant
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.
I was running into some catastrophic cancellation that could be solved with BigFloat only by going to very high precision (2048 bits). This made the code slow. This PR fixes gamma_inc for that.
While Claude was used to write the PR, I kept a very close eye on the code (I put less effort into the commit messages). Thus, this PR is best reviewed one commit at a time.