JuliaMath / SpecialFunctions.jl

Special mathematical functions in Julia
https://specialfunctions.juliamath.org/stable/
Other
352 stars 98 forks source link

Rewrite `gamma_inc_taylor` and `gamma_inc_asym` more concisely #443

Open BioTurboNick opened 1 year ago

BioTurboNick commented 1 year ago

Many functions were translated directly from Fortran. This PR rewrites two of them using more concise Julia forms.

Ideally would like to eliminate the array allocation as well. (See #442) And now does!

codecov[bot] commented 1 year ago

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.01 :tada:

Comparison is base (ae35d10) 94.32% compared to head (8d37a34) 94.34%.

Additional details and impacted files ```diff @@ Coverage Diff @@ ## master #443 +/- ## ========================================== + Coverage 94.32% 94.34% +0.01% ========================================== Files 14 14 Lines 2909 2881 -28 ========================================== - Hits 2744 2718 -26 + Misses 165 163 -2 ``` | Flag | Coverage Δ | | |---|---|---| | unittests | `94.34% <100.00%> (+0.01%)` | :arrow_up: | Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMath#carryforward-flags-in-the-pull-request-comment) to find out more. | [Impacted Files](https://app.codecov.io/gh/JuliaMath/SpecialFunctions.jl/pull/443?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMath) | Coverage Δ | | |---|---|---| | [src/gamma\_inc.jl](https://app.codecov.io/gh/JuliaMath/SpecialFunctions.jl/pull/443?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMath#diff-c3JjL2dhbW1hX2luYy5qbA==) | `93.49% <100.00%> (+0.03%)` | :arrow_up: |

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

BioTurboNick commented 1 year ago

Oops, left an nmc off that commit. Apologies to that person.

heltonmc commented 1 year ago

I know I originally suggested it but it would still be nice to know why it was implemented in this way to begin with as it is quite unusual and we are really just relying on the strength of the original test set here.

My last comment is just about the convergence check. My original suggestion (https://github.com/JuliaMath/SpecialFunctions.jl/issues/442#issuecomment-1534610427) checked for relative tolerance where this just checks for absolute tolerance. It probably doesn't matter but if the final result is much smaller than 1 they could be different. (edit: relative check will be slower so not saying just to do it instead)

BioTurboNick commented 1 year ago

Alright, thanks! :-) I reverted to the original version, with your smaller suggestions. I kept the other one on another branch, gamma-fix2, just in case.

heltonmc commented 1 year ago

Sounds good 😊 I’ll let other people comment on what they think of the two versions.

This version will allocate though which might not solve your initial problem. If you want to sum in reverse order your probably best bet is to to use @nexprs to remove the array then can sum in any order you want. I’m not for sure it’s really needed here though but it would be nice since we are putting in the effort to make this allocation free !

BioTurboNick commented 1 year ago

Ah, found an allocation-free way! Before I commit it:

# compute first 21 terms
ts = cumprod(ntuple(i -> x / (a + i), Val(21)))
last_t = findfirst(<(1.0e-3), ts)
last_t = last_t === nothing ? 20 : last_t - 1
next_t = last_t + 1
BioTurboNick commented 1 year ago

Ah, cumprod for tuples doesn't exist in Julia 1.3... any chance of dropping support for Julia before 1.6? (It was introduced in 1.5)

EDIT: Or can just define it here. Nevermind, too much to bring in.

BioTurboNick commented 1 year ago

Bumping - any more work needed on this?

BioTurboNick commented 1 year ago

What do I need to do to get this in? I have a couple bits of code that would really benefit from this change, and I'd like not to have to disable precompilation.

ViralBShah commented 8 months ago

@devmotion Is this good to merge?