Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added erf(x) Float64 Julia implementation #491

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Added erf(x) Float64 Julia implementation #491

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
6f554ef
Added erf(x) Float64/Float32 Julia implementation
AhmedYKadah Mar 31, 2025
7f4fd2d
changed erf to _erf, got rid of unnecessary branch
AhmedYKadah Mar 31, 2025
e784c9f
fixed syntax error in ccall
AhmedYKadah Mar 31, 2025
da16cb1
fixed syntax error in ccall 2
AhmedYKadah Mar 31, 2025
0a755b6
NaN edge case for erf(x)
AhmedYKadah Mar 31, 2025
6efcec8
added test cases for erf(x)
AhmedYKadah Mar 31, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
254 changes: 249 additions & 5 deletions src/erf.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,6 @@ for f in (:erf, :erfc)
@eval begin
$f(x::Number) = $internalf(float(x))

$internalf(x::Float64) = ccall(($libopenlibmf, libopenlibm), Float64, (Float64,), x)
$internalf(x::Float32) = ccall(($libopenlibmf0, libopenlibm), Float32, (Float32,), x)
$internalf(x::Float16) = Float16($internalf(Float32(x)))

$internalf(z::Complex{Float64}) = Complex{Float64}(ccall(($openspecfunf, libopenspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64)))
$internalf(z::Complex{Float32}) = Complex{Float32}(ccall(($openspecfunf, libopenspecfun), Complex{Float64}, (Complex{Float64}, Float64), Complex{Float64}(z), Float64(eps(Float32))))
Expand All @@ -28,6 +25,10 @@ for f in (:erf, :erfc)
end
end

_erfc(x::Float64) = ccall((:erfc, libopenlibm), Float64, (Float64,), x)
_erfc(x::Float32) = ccall((:erfcf, libopenlibm), Float32, (Float32,), x)
_erfc(x::Float16) = Float16(_erfc(Float32(x)))

for f in (:erfcx, :erfi, :dawson, :faddeeva)
internalf = Symbol(:_, f)
openspecfunfsym = Symbol(:Faddeeva_, f === :dawson ? :Dawson : f === :faddeeva ? :w : f)
Expand Down Expand Up @@ -96,10 +97,253 @@ See also:
[`erfinv(x)`](@ref erfinv), [`erfcinv(x)`](@ref erfcinv).

# Implementation by
- `Float32`/`Float64`: C standard math library
[libm](https://en.wikipedia.org/wiki/C_mathematical_functions#libm).
- `Float32`/`Float64`: Julia implementation of https://github.com/ARM-software/optimized-routines/blob/master/math/erf.c
- `BigFloat`: C library for multiple-precision floating-point [MPFR](https://www.mpfr.org/)
"""

# Fast erf implementation using a mix of
# rational and polynomial approximations.
# Highest measured error is 1.01 ULPs at 0x1.39956ac43382fp+0.
function _erf(x::Float64)
# Minimax approximation of erf
PA=(0x1.06eba8214db68p-3, -0x1.812746b037948p-2, 0x1.ce2f21a03872p-4,-0x1.b82ce30e6548p-6, 0x1.565bcc360a2f2p-8, -0x1.c02d812bc979ap-11,0x1.f99bddfc1ebe9p-14, -0x1.f42c457cee912p-17, 0x1.b0e414ec20ee9p-20,-0x1.18c47fd143c5ep-23)
# Rational approximation on [0x1p-28, 0.84375]
NA=(0x1.06eba8214db68p-3, -0x1.4cd7d691cb913p-2, -0x1.d2a51dbd7194fp-6,-0x1.7a291236668e4p-8, -0x1.8ead6120016acp-16)
DA=(0x1.97779cddadc09p-2, 0x1.0a54c5536cebap-4, 0x1.4d022c4d36b0fp-8,0x1.15dc9221c1a1p-13, -0x1.09c4342a2612p-18)
# Rational approximation on [0.84375, 1.25]
NB=( -0x1.359b8bef77538p-9, 0x1.a8d00ad92b34dp-2, -0x1.7d240fbb8c3f1p-2, 0x1.45fca805120e4p-2, -0x1.c63983d3e28ecp-4, 0x1.22a36599795ebp-5, -0x1.1bf380a96073fp-9 )
DB=( 0x1.b3e6618eee323p-4, 0x1.14af092eb6f33p-1, 0x1.2635cd99fe9a7p-4, 0x1.02660e763351fp-3, 0x1.bedc26b51dd1cp-7, 0x1.88b545735151dp-7 )

# Generated using Sollya::remez(f(c*x+d), deg, [(a-d)/c;(b-d)/c], 1, 1e-16), [|D ...|] with deg=15 a=1.25 b=2 c=1 d=1.25
PC=( 0x1.3bcd133aa0ffcp-4, -0x1.e4652fadcb702p-3, 0x1.2ebf3dcca0446p-2, -0x1.571d01c62d66p-3, 0x1.93a9a8f5b3413p-8, 0x1.8281cbcc2cd52p-5, -0x1.5cffd86b4de16p-6, -0x1.db4ccf595053ep-9, 0x1.757cbf8684edap-8, -0x1.ce7dfd2a9e56ap-11, -0x1.99ee3bc5a3263p-11, 0x1.3c57cf9213f5fp-12, 0x1.60692996bf254p-14, -0x1.6e44cb7c1fa2ap-14, 0x1.9d4484ac482b2p-16, -0x1.578c9e375d37p-19)
# Generated using Sollya::remez(f(c*x+d), deg, [(a-d)/c;(b-d)/c], 1, 1e-16), [|D ...|] with deg=17 a=2 b=3.25 c=2 d=2
PD=( 0x1.328f5ec350e5p-8, -0x1.529b9e8cf8e99p-5, 0x1.529b9e8cd9e71p-3, -0x1.8b0ae3a023bf2p-2, 0x1.1a2c592599d82p-1, -0x1.ace732477e494p-2, -0x1.e1a06a27920ffp-6, 0x1.bae92a6d27af6p-2, -0x1.a15470fcf5ce7p-2, 0x1.bafe45d18e213p-6, 0x1.0d950680d199ap-2, -0x1.8c9481e8f22e3p-3, -0x1.158450ed5c899p-4, 0x1.c01f2973b44p-3, -0x1.73ed2827546a7p-3, 0x1.47733687d1ff7p-4, -0x1.2dec70d00b8e1p-6, 0x1.a947ab83cd4fp-10 )
# Generated using Sollya::remez(f(c*x+d), deg, [(a-d)/c;(b-d)/c], 1, 1e-16), [|D ...|] with deg=13 a=3.25 b=4 c=1 d=3.25
PE=( 0x1.20c13035539e4p-18, -0x1.e9b5e8d16df7ep-16, 0x1.8de3cd4733bf9p-14, -0x1.9aa48beb8382fp-13, 0x1.2c7d713370a9fp-12, -0x1.490b12110b9e2p-12, 0x1.1459c5d989d23p-12, -0x1.64b28e9f1269p-13, 0x1.57c76d9d05cf8p-14, -0x1.bf271d9951cf8p-16, 0x1.db7ea4d4535c9p-19, 0x1.91c2e102d5e49p-20, -0x1.e9f0826c2149ep-21, 0x1.60eebaea236e1p-23 )
# Generated using Sollya::remez(f(c*x+d), deg, [(a-d)/c;(b-d)/c], 1, 1e-16), [|D ...|] with deg=16 a=4 b=5.90625 c=2 d=4
PF=( 0x1.08ddd130d1fa6p-26, -0x1.10b146f59ff06p-22, 0x1.10b135328b7b2p-19, -0x1.6039988e7575fp-17, 0x1.497d365e19367p-15, -0x1.da48d9afac83ep-14, 0x1.1024c9b1fbb48p-12, -0x1.fc962e7066272p-12, 0x1.87297282d4651p-11, -0x1.f057b255f8c59p-11, 0x1.0228d0eee063p-10, -0x1.b1b21b84ec41cp-11, 0x1.1ead8ae9e1253p-11, -0x1.1e708fba37fccp-12, 0x1.9559363991edap-14, -0x1.68c827b783d9cp-16, 0x1.2ec4adeccf4a2p-19 )

C = 0x1.b0ac16p-1

TwoOverSqrtPiMinusOne=0x1.06eba8214db69p-3


# # top 32 bits
ix::UInt32=reinterpret(UInt64,x)>>32
# # top 32, without sign bit
ia::UInt32=ix & 0x7fffffff
# # sign
# sign::UInt32=ix>>31

sign::Bool=x<0



if (ia < 0x3feb0000)
# a = |x| < 0.84375.

x2 = x * x

if (ia < 0x3fe00000)
## a < 0.5 - Use polynomial approximation.
r1 = fma(x2, PA[2], PA[1])
r2 = fma(x2, PA[4], PA[3])
r3 = fma(x2, PA[6], PA[5])
r4 = fma(x2, PA[8], PA[7])
r5 = fma(x2, PA[10], PA[9])

x4 = x2 * x2
r = r5
r = fma(x4, r, r4)
r = fma(x4, r, r3)
r = fma(x4, r, r2)
r = fma(x4, r, r1)
return fma(r, x, x) ## This fma is crucial for accuracy.
else
## 0.5 <= a < 0.84375 - Use rational approximation.

r1n = fma(x2, NA[2], NA[1])
x4 = x2 * x2
r2n = fma(x2, NA[4], NA[3])
x8 = x4 * x4
r1d = fma(x2, DA[1], 1.0)
r2d = fma(x2, DA[3], DA[2])
r3d = fma(x2, DA[5], DA[4])
P = r1n + x4 * r2n + x8 * NA[5]

Q = r1d + x4 * r2d + x8 * r3d
return fma(P / Q, x, x)
end
elseif (ia < 0x3ff40000)
## 0.84375 <= |x| < 1.25.

a = abs(x) - 1.0
r1n = fma(a, NB[2], NB[1])
a2 = a * a
r1d = fma(a, DB[1], 1.0)
a4 = a2 * a2
r2n = fma(a, NB[4], NB[3])
a6 = a4 * a2
r2d = fma(a, DB[3], DB[2])
r3n = fma(a, NB[6], NB[5])
r3d = fma(a, DB[5], DB[4])
r4n = NB[7]
r4d = DB[6]

P = r1n + a2 * r2n + a4 * r3n + a6 * r4n
Q = r1d + a2 * r2d + a4 * r3d + a6 * r4d
if (sign)
return -C - P / Q
else
return C + P / Q
end
elseif (ia < 0x40000000)
## 1.25 <= |x| < 2.0.
a = abs(x)
a = a - 1.25

r1 = fma(a, PC[2], PC[1])
r2 = fma(a, PC[4], PC[3])
r3 = fma(a, PC[6], PC[5])
r4 = fma(a, PC[8], PC[7])
r5 = fma(a, PC[10], PC[9])
r6 = fma(a, PC[12], PC[11])
r7 = fma(a, PC[14], PC[13])
r8 = fma(a, PC[16], PC[15])


a2 = a * a

r = r8
r = fma(a2, r, r7)
r = fma(a2, r, r6)
r = fma(a2, r, r5)
r = fma(a2, r, r4)
r = fma(a2, r, r3)
r = fma(a2, r, r2)
r = fma(a2, r, r1)

if (sign)
return -1.0 + r
else
return 1.0 - r
end
elseif (ia < 0x400a0000)
## 2 <= |x| < 3.25.
a = abs(x)
a = fma(0.5, a, -1.0)

r1 = fma(a, PD[2], PD[1])
r2 = fma(a, PD[4], PD[3])
r3 = fma(a, PD[6], PD[5])
r4 = fma(a, PD[8], PD[7])
r5 = fma(a, PD[10], PD[9])
r6 = fma(a, PD[12], PD[11])
r7 = fma(a, PD[14], PD[13])
r8 = fma(a, PD[16], PD[15])
r9 = fma(a, PD[18], PD[17])

a2 = a * a

r = r9
r = fma(a2, r, r8)
r = fma(a2, r, r7)
r = fma(a2, r, r6)
r = fma(a2, r, r5)
r = fma(a2, r, r4)
r = fma(a2, r, r3)
r = fma(a2, r, r2)
r = fma(a2, r, r1)

if (sign)
return -1.0 + r
else
return 1.0 - r
end
elseif (ia < 0x40100000)
## 3.25 <= |x| < 4.0.
a = abs(x)
a = a - 3.25

r1 = fma(a, PE[2], PE[1])
r2 = fma(a, PE[4], PE[3])
r3 = fma(a, PE[6], PE[5])
r4 = fma(a, PE[8], PE[7])
r5 = fma(a, PE[10], PE[9])
r6 = fma(a, PE[12], PE[11])
r7 = fma(a, PE[14], PE[13])


a2 = a * a

r = r7
r = fma(a2, r, r6)
r = fma(a2, r, r5)
r = fma(a2, r, r4)
r = fma(a2, r, r3)
r = fma(a2, r, r2)
r = fma(a2, r, r1)

if (sign)
return -1.0 + r
else
return 1.0 - r
end
elseif (ia < 0x4017a000)
## 4 <= |x| < 5.90625.
a = abs(x)
a = fma(0.5, a, -2.0)

r1 = fma(a, PF[2], PF[1])
r2 = fma(a, PF[4], PF[3])
r3 = fma(a, PF[6], PF[5])
r4 = fma(a, PF[8], PF[7])
r5 = fma(a, PF[10], PF[9])
r6 = fma(a, PF[12], PF[11])
r7 = fma(a, PF[14], PF[13])
r8 = fma(a, PF[16], PF[15])

r9 = PF[17]

a2 = a * a

r = r9
r = fma(a2, r, r8)
r = fma(a2, r, r7)
r = fma(a2, r, r6)
r = fma(a2, r, r5)
r = fma(a2, r, r4)
r = fma(a2, r, r3)
r = fma(a2, r, r2)
r = fma(a2, r, r1)

if (sign)
return -1.0 + r
else
return 1.0 - r
end
else

if(isnan(x))
return NaN
end

if (sign)
return -1.0
else
return 1.0
end

end


end

_erf(x::Float32)=Float32(_erf(Float64(x)))

_erf(x::Float16)=Float16(_erf(Float64(x)))


function erf end
"""
erf(x, y)
Expand Down
16 changes: 16 additions & 0 deletions test/erf.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,23 @@
@testset "real argument" begin
for T in (Float16, Float32, Float64)
@test @inferred(erf(T(1))) isa T
@test erf(T(0.25)) ≈ T(0.27632639016823696) rtol=2*eps(T)
@test erf(T(0.75)) ≈ T(0.7111556336535151) rtol=2*eps(T)
@test erf(T(1)) ≈ T(0.84270079294971486934) rtol=2*eps(T)
@test erf(T(1.5)) ≈ T(0.9661051464753108) rtol=2*eps(T)
@test erf(T(2.5)) ≈ T(0.9995930479825551) rtol=2*eps(T)
@test erf(T(3.5)) ≈ T(0.9999992569016276) rtol=2*eps(T)
@test erf(T(4.5)) ≈ T(0.9999999998033839) rtol=2*eps(T)
@test erf(T(6)) ≈ T(1.0) rtol=2*eps(T)

@test erf(T(-0.25)) ≈ T(-0.27632639016823696) rtol=2*eps(T)
@test erf(T(-0.75)) ≈ T(-0.7111556336535151) rtol=2*eps(T)
@test erf(T(-1)) ≈ T(-0.84270079294971486934) rtol=2*eps(T)
@test erf(T(-1.5)) ≈ T(-0.9661051464753108) rtol=2*eps(T)
@test erf(T(-2.5)) ≈ T(-0.9995930479825551) rtol=2*eps(T)
@test erf(T(-3.5)) ≈ T(-0.9999992569016276) rtol=2*eps(T)
@test erf(T(-4.5)) ≈ T(-0.9999999998033839) rtol=2*eps(T)
@test erf(T(-6)) ≈ T(-1.0) rtol=2*eps(T)

@test @inferred(erfc(T(1))) isa T
@test erfc(T(1)) ≈ T(0.15729920705028513066) rtol=2*eps(T)
Expand Down
Loading