Expand the minimax Temme series in 1/a rather than z - #534
Conversation
gamma_inc_minimax evaluates the |1 - x/a| <= 1e-3 branch of
T(a,lambda) = sum_k c_k(z) a^-k
with @evalpoly(z, c0, ..., d80). The outer sum is in a^-k, as the
function's own docstring states and as the other three code paths
built from the same d tables already do: the rational branch ten lines
below uses @evalpoly(1.0/a, ...), gamma_inc_temme uses @evalpoly(1.0/a, ...)
and gamma_inc_temme_1 uses @evalpoly(u, ...) with u = 1.0/a.
In the affected band |z| is 2-5% of 1/a, so every correction term past
c0 is weighted 20-40x too small. Over 300 points with a in [25, 1.2e6]
the worst relative error drops from 1.16e-5 to 8.2e-14; afterwards it
sits at or below the 1-ulp conditioning limit of x. No point regresses,
and every other branch is bit-identical.
The band was covered by exactly one existing assertion, gamma_inc(1e7,
1e7 + 1), where |z| and 1/a happen to coincide to 1e-7 and the defect
cancels. The new testset pins 16 points across the band against mpmath.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #534 +/- ##
==========================================
+ Coverage 94.28% 94.34% +0.06%
==========================================
Files 14 14
Lines 3026 3026
==========================================
+ Hits 2853 2855 +2
+ Misses 173 171 -2
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:
|
| @testset "a=$a, x=$x" for (a, x, p, q) in ( | ||
| (25.0, 24.97525, 0.5246323715433058506793, 0.4753676284566941493207), | ||
| (25.0, 25.02475, 0.5285687428518152079350, 0.4714312571481847920650), | ||
| (25.0, 24.9925, 0.5260050200501651608654, 0.4739949799498348391346), | ||
| (50.0, 50.0495, 0.5215950013484762091395, 0.4784049986515237908605), | ||
| (50.0, 49.975, 0.5173998411006592979910, 0.4826001588993407020090), | ||
| (134.0, 134.05, 0.5132100266748297555471, 0.4867899733251702444529), | ||
| (134.0, 133.86734, 0.5069170257485132637813, 0.4930829742514867362187), | ||
| (400.0, 400.396, 0.5145421195364036819570, 0.4854578804635963180430), | ||
| (400.0, 399.92, 0.5050535378251085332799, 0.4949464621748914667201), | ||
| (2000.0, 2001.98, 0.5206211429087569237891, 0.4793788570912430762109), | ||
| (2000.0, 1998.5, 0.4895906656999191847208, 0.5104093343000808152792), | ||
| (10000.0, 9990.1, 0.4618797888026423050784, 0.5381202111973576949216), | ||
| (10000.0, 10005.0, 0.5212634681357695888112, 0.4787365318642304111888), | ||
| (100000.0, 100099.0, 0.6232456486261727122096, 0.3767543513738272877904), | ||
| (100000.0, 99901.0, 0.3774766852491783530590, 0.6225233147508216469410), | ||
| (100000.0, 99985.0, 0.4815027198019921823376, 0.5184972801980078176624), | ||
| ) | ||
| @test gamma_inc(a, x)[1] ≈ p rtol=1e-13 | ||
| @test gamma_inc(a, x)[2] ≈ q rtol=1e-13 |
There was a problem hiding this comment.
AFAIU these tests are failing on master?
There was a problem hiding this comment.
Yes — all 32 assertions fail on master (02a51af), worst relative error 1.16e-5 at a=25, x=24.97525, shrinking to ~7.8e-10 by a=1e5. With the one-line change they pass at rtol=1e-13.
The red CI isn't them, though: that's the pre-existing gamma_inc allocations testset on Julia pre. It fails 12 times on master too — https://github.com/JuliaMath/SpecialFunctions.jl/actions/runs/27935023667 on 02a51af (2026-06-22) has the same 12 failures at test/gamma_inc.jl:301/303/305/310, all three OSes, everything else green. Here they show as :328/:330/:332/:337 only because this PR inserts 27 lines above them.
There was a problem hiding this comment.
Done — bumped to 2.8.1 in 8a25843, ready to tag after merge.
devmotion
left a comment
devmotion
left a comment
There was a problem hiding this comment.
Can you update the version number, so we can tag a release with the fix once it is merged?
|
Thank you! |
2 participants
The
abs(s) <= 1e-3branch ofgamma_inc_minimaxevaluates the outer sum ofT(a,λ) = Σ cₖ(z)·a⁻ᵏwith@evalpoly(z, c0, ..., d80), but that sum runs ina⁻ᵏ— the docstring three lines above says so, and the other three paths built from the samedtables already do it: the rational branch twenty lines below uses@evalpoly(1.0/a, ...),gamma_inc_temmeuses1.0/a, andgamma_inc_temme_1usesu = 1.0/aover coefficients constructed identically to these. The existing suite already pins that convention in all three — applying the same1.0/a → zswap togamma_inc_temme_1reddenstest/gamma_inc.jl:18, togamma_inc_temmereddens:27, and to the minimax rational branch reddens both:19and:26; the one path of the four that no existing test pins is the one that gotz. It is the same one-token transcription shape as #418 → #419, which fixedind_→indingamma_inc_temme_1; the minimax branch wasn't looked at then.In that band
|z|is only 2–5% of1/a, so every term pastc0is weighted 20–40× too small, and the worst relative error over the branch drops from1.16e-5to8.2e-14— e.g.gamma_inc(134.0, 134.05)[2]returns0.48679041818186286where the true value is0.4867899733251702445. The change is contained: over 2331 evaluations spanninga ∈ [25, 1.2e6],|1 - x/a| ∈ [1e-7, 0.4], both signs and all threeindvalues, only the 301 points inside this branch move at all — everyminimax_rational,temme_1,temmeandind=1/ind=2result is bit-identical. The band was effectively untested, which is presumably how it survived:test/gamma_inc.jlenters it twice with a single unique input,gamma_inc(1e7, 1e7 + 1), and that point happens to sit where|z|/(1/a) = 1.0000, so it is bit-identical either way. The new testset pins 16 points across the band atrtol=1e-13(reference values from mpmath at 50+ digits, cross-checked against the DLMF 8.7.1 series and the 8.9.2 continued fraction); fullPkg.test()is green on 1.12, withgamma_incgoing from 1001809 to 1001841 assertions.