erfinv(::Float16) and erfcinv(::Float16) are wrong for abs(x) <= 0.75 (erfinv(Float16(0)) == 0.2201)
Nobody has claimed this yet.
Assessment
- Difficulty
- 2/5
- Estimated time
- 1-3 hours
- Newbie friendliness
- 82/100
Research direction
Start in src/erf.jl at the Float16 branch of _erfinv(x), then inspect the existing tests for erfinv and erfcinv. Compare Float16 results across the stated input range with the corresponding Float64 results; done means the reported wrong values are corrected and a regression test covers the range.
Written by the indexing model from the issue text.
Description
erfinv(::Float16) returns wrong values for every argument with abs(x) <= 0.75, x == 0 included. erfcinv(::Float16) calls erfinv(1 - y) for y > 0.0625, so it is wrong for 0.25 <= y <= 1.75 too.
julia> using SpecialFunctions
julia> erfinv(Float16(0))
Float16(0.2201)
julia> erfinv(Float16(0.5)), erfinv(0.5)
(Float16(12.91), 0.4769362762044699)
julia> erfcinv(Float16(0.5)), erfcinv(0.5)
(Float16(12.91), 0.4769362762044699)
Out of the 14849 Float16 values in [0, 0.75], none match Float16(erfinv(Float64(x))). The branch 0.75 < abs(x) < 1 is fine: all 511 values there are within 2 ulps.
Cause
This is the first branch of _erfinv(x::Float16) in src/erf.jl:
if a32 <= 0.75f0
# Simpler and more accurate alternative to Table 7 in Blair et al.
# Ref: https://github.com/JuliaMath/SpecialFunctions.jl/pull/372#discussion_r1592832735
t = muladd(-6.73815f1, x32, 1f0) / muladd(-4.18798f0, x32, 4.54263f0)
y = copysign(muladd(0.88622695f0, x32, t), x32)
It came from a review suggestion in #372. It isn't odd in x, and it isn't zero at x = 0. The Table 7 approximation it replaced is correct:
t = x32^2 - 0.5625f0
y = x32 * @evalpoly(t, -0.10976_672f1, 0.53062_1f0) / @evalpoly(t, -0.10123_953f1, 0.1f1)
With it, 9070 of the 14849 values equal Float16(erfinv(Float64(x))) and the others are off by 1 ulp. Restoring it, or computing Float16(erfinv(Float32(x))), would fix both functions. A test comparing erfinv(::Float16) with erfinv(Float64(x)) over the whole range would catch this kind of regression.
Versions
SpecialFunctions 2.9.0, Julia 1.13.3. The code is the same on master.
Found because blsimpv of FinancialToolbox.jl returned an implied volatility of 2.75 instead of 0.2 for at the money Float16 options.
- Dominant language
- Julia
- Stars
- 384
- Forks
- 115
- Avg merge
- 2d 9h
- Merged PRs (30d)
- 2
Getting set up
This project ships no dev container, Dockerfile or contributing guide, so setting up is up to you: start from its README, and see our first-contribution guide for the general steps.
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
More from JuliaMath/SpecialFunctions.jl
-
Difficulty 2/5 1-3 hours Newbie friendliness 68/100
JuliaMath/SpecialFunctions.jl#570 ·
-
Difficulty 5/5 Over a week Newbie friendliness 40/100
JuliaMath/SpecialFunctions.jl#569 · 2 comments ·
-
Difficulty 4/5 3-5 days Newbie friendliness 25/100
JuliaMath/SpecialFunctions.jl#567 ·
-
Difficulty 4/5 1-2 days Newbie friendliness 56/100
JuliaMath/SpecialFunctions.jl#563 ·
-
Difficulty 3/5 1-2 days Newbie friendliness 72/100
JuliaMath/SpecialFunctions.jl#560 ·
All issues in JuliaMath/SpecialFunctions.jl
Similar issues
-
Difficulty 2/5 1-3 hours Newbie friendliness 73/100
Maintainers usually reply within 1 day
-
lux
Difficulty 2/5 1-3 hours Newbie friendliness 76/100
JuliaGraphs/GraphNeuralNetworks.jl#732 ·
Maintainers usually reply within 2 days
-
Difficulty 2/5 1-3 hours Newbie friendliness 68/100
Maintainers usually reply within 1 day
-
Difficulty 2/5 1-3 hours Newbie friendliness 74/100
-
Difficulty 2/5 1-3 hours Newbie friendliness 67/100