Skip to content

Fix the iop=2 Temme expansion variable and the ind=2 wide-band truncation - #538

Merged
devmotion merged 5 commits into
JuliaMath:masterfrom
gaoflow:fix-temme-iop2-and-iop3-wide
Aug 5, 2026
Merged

Fix the iop=2 Temme expansion variable and the ind=2 wide-band truncation#538
devmotion merged 5 commits into
JuliaMath:masterfrom
gaoflow:fix-temme-iop2-and-iop3-wide

Conversation

@gaoflow

@gaoflow gaoflow commented Jul 31, 2026

Copy link
Copy Markdown
Contributor

gamma_inc(a, x, 1) enters the iop==2 Temme branch for a >= 14 within |1 - x/a| <=
0.025/sqrt(a). That branch expands c0/c1 in the wrong variable (@evalpoly(d00, ...)
puts d00 in the variable slot), so the z-dependence is dropped entirely (NSWC gratio
label 270 is C0 = D00 + D0(1)Z + D0(2)Z^2). Results miss the documented 6-digit tier
by up to ~9% and no existing test enters the branch.

The ind=2 wide band (0.14/sqrt(a) < |1 - x/a| <= 0.4) had the same class of porting
error: the 4-term expansion (label 230) was computed and then discarded, falling back
to the narrow-band 2-term form (label 280) — up to ~1.6e-3 error against the 3-digit
tier.

Repro: gamma_inc(14.0, 13.906552106765321, 1) -> (0.4806, 0.5194); true P =
0.525615807496999.

Tests: 34 reference points across both bands plus ind=0 controls (68 assertions,
test/gamma_inc.jl); full suite passes.

…tion

The iop=2 branch expanded c0/c1 in d00/d10 instead of z (NSWC gratio label
270), dropping the z-dependence of the expansion, and the ind=2 wide band
discarded its computed 4-term expansion (label 230) in favour of the
narrow-band 2-term form (label 280). Add reference-point regression tests
for both branches.
@codecov

codecov Bot commented Jul 31, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 94.49%. Comparing base (172df8d) to head (e148d24).

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #538      +/-   ##
==========================================
+ Coverage   94.34%   94.49%   +0.14%     
==========================================
  Files          14       14              
  Lines        3026     3016      -10     
==========================================
- Hits         2855     2850       -5     
+ Misses        171      166       -5     
Flag Coverage Δ
unittests 94.49% <100.00%> (+0.14%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Comment thread test/gamma_inc.jl
Comment on lines +103 to +104
@test gamma_inc(a, x, 1)[1] ≈ p rtol=5e-7
@test gamma_inc(a, x, 1)[2] ≈ q rtol=5e-7

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are these tolerances (and the ones below) as strict as possible, to increase the chances of detecting regressions?

Comment thread src/gamma_inc.jl Outdated
Comment on lines +922 to +928
if l < 1.0
p = c*(w - rt2pin*t/rta)
return (p, 1.0 - p)
else
q = c*(w + rt2pin*t/rta)
return (1.0 - q, q)
end

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems this is the same code as in gamma_inc_temme_1 above (and possibly other places). Maybe could be worth moving to a separate function?

…d tolerances

Extract the c/w/t -> (P,Q) combine into _gamma_inc_temme_combine, used by
both gamma_inc_temme_1 and the ind=2 wide-band branch in _gamma_inc that
previously duplicated it inline.

Measured actual error against the mpmath references: wide-band tops out
at 1.3e-4, narrow-band at 4.6e-5, so tighten rtol from 5e-4 to 2e-4/1e-4
respectively (~1.5-2x margin) to catch regressions sooner.
@gaoflow

gaoflow commented Aug 3, 2026

Copy link
Copy Markdown
Contributor Author

Both fair. Extracted the shared c/w/t -> (P,Q) combine into _gamma_inc_temme_combine, used by both gamma_inc_temme_1 and this branch (2b880b8).

On tolerances: measured actual error against the mpmath references — wide-band tops out at 1.3e-4, narrow-band at 4.6e-5 — so tightened rtol from 5e-4 to 2e-4 / 1e-4 (roughly 1.5-2x margin over what's observed). Full suite still green.

Comment thread src/gamma_inc.jl Outdated

External links: [DLMF 8.12.8](https://dlmf.nist.gov/8.12.8)
"""
# Combine the Temme series value `t` with the Gaussian-tail term `c*w` into (P, Q);

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This removes the docstring (above) from gamma_inc_temme_1 (below). Moreover, this function shows up in more places in this file.

Fixes the docstring being pushed off gamma_inc_temme_1 by the
helper's insertion point; also folds in gamma_inc_minimax and
gamma_inc_temme, which had the same duplicated block.
@gaoflow

gaoflow commented Aug 4, 2026

Copy link
Copy Markdown
Contributor Author

Good catch on both — the helper landed between the docstring and gamma_inc_temme_1, so @doc was attaching to the wrong thing. Moved _gamma_inc_temme_combine above gamma_inc_minimax (the earliest caller) and confirmed @doc gamma_inc_temme_1 resolves correctly now (1397469).

And yes — gamma_inc_minimax and gamma_inc_temme had the identical c/w/t→(P,Q) block too, just not the two I'd touched. All 5 sites now go through the one helper; full suite green (1,001,909 gamma_inc assertions).

@devmotion
devmotion merged commit 42ec04a into JuliaMath:master Aug 5, 2026
15 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants