From e0347ccda4670806593c54c304bda6a1d15b101a Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Fri, 1 Aug 2025 20:06:20 +0200 Subject: [PATCH 1/9] fix accuracy of `logit` Fixes #98 --- src/basicfuns.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/basicfuns.jl b/src/basicfuns.jl index 5cc0df2..dbb46f9 100644 --- a/src/basicfuns.jl +++ b/src/basicfuns.jl @@ -120,7 +120,13 @@ for ``0 < x < 1``. Its inverse is the [`logistic`](@ref) function. """ -logit(x::Real) = log(x / (one(x) - x)) +function logit(x::Real) + if 4 * x < 1 + -log(inv(x) - 1) + else + 2 * atanh(2*x - 1) + end +end """ $(SIGNATURES) From 51ed3c3b2e1cbc159c7bff9f468da850774251fc Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 2 Aug 2025 12:28:18 +0200 Subject: [PATCH 2/9] test accuracy --- test/basicfuns.jl | 7 ++++++ test/common/ULPError.jl | 55 +++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 3 +++ 3 files changed, 65 insertions(+) create mode 100644 test/common/ULPError.jl diff --git a/test/basicfuns.jl b/test/basicfuns.jl index 6a84b54..0b76115 100644 --- a/test/basicfuns.jl +++ b/test/basicfuns.jl @@ -86,6 +86,13 @@ end @test logistic(+750.0) === 1.0 @test iszero(logit(0.5)) @test logit(logistic(2)) ≈ 2.0 + @testset "accuracy of `logit`" begin + for t in (Float16, Float32, Float64) + for x in range(start = t(0), stop = t(1), length = 500) + @test 2 * ulp_error(logit, x) < 3 + end + end + end end @testset "logcosh and logabssinh" begin diff --git a/test/common/ULPError.jl b/test/common/ULPError.jl new file mode 100644 index 0000000..abb44a8 --- /dev/null +++ b/test/common/ULPError.jl @@ -0,0 +1,55 @@ +module ULPError + export ulp_error, ulp_error_maximum + @noinline function throw_invalid() + throw(ArgumentError("invalid")) + end + function ulp_error(accurate::AbstractFloat, approximate::AbstractFloat) + if isnan(accurate) + @noinline throw_invalid() + end + if isnan(approximate) + @noinline throw_invalid() + end + # the ULP error is usually not required to great accuracy, so `Float32` should be precise enough + zero_return = Float32(0) + inf_return = Float32(Inf) + if isinf(accurate) || iszero(accurate) # handle floating-point edge cases + if isinf(accurate) + if isinf(approximate) && (signbit(accurate) == signbit(approximate)) + return zero_return + end + return inf_return + end + # `iszero(accurate)` + if iszero(approximate) + return zero_return + end + return inf_return + end + # assuming `precision(BigFloat)` is great enough + acc = if accurate isa BigFloat + accurate + else + BigFloat(accurate)::BigFloat + end + err = abs(Float32((approximate - acc) / eps(approximate))::Float32) + if isnan(err) + @noinline throw_invalid() # unexpected + end + err + end + function ulp_error(accurate::Acc, approximate::App, x::AbstractFloat) where {Acc, App} + acc = accurate(x) + app = approximate(x) + ulp_error(acc, app) + end + function ulp_error(func::Func, x::AbstractFloat) where {Func} + ulp_error(func ∘ BigFloat, func, x) + end + function ulp_error_maximum(func::Func, iterator) where {Func} + function f(x::AbstractFloat) + ulp_error(func, x) + end + maximum(f, iterator) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 27b247f..5fcc6ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,6 +10,9 @@ using OffsetArrays using Random using Test +include("common/ULPError.jl") +using .ULPError + Random.seed!(1234) include("basicfuns.jl") From 0a66f16ef8e7a5d8349638230e09de76fa57e52e Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Tue, 5 Aug 2025 14:38:32 +0200 Subject: [PATCH 3/9] update ULPError --- test/common/ULPError.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/test/common/ULPError.jl b/test/common/ULPError.jl index abb44a8..f788bd9 100644 --- a/test/common/ULPError.jl +++ b/test/common/ULPError.jl @@ -4,15 +4,17 @@ module ULPError throw(ArgumentError("invalid")) end function ulp_error(accurate::AbstractFloat, approximate::AbstractFloat) - if isnan(accurate) - @noinline throw_invalid() - end - if isnan(approximate) - @noinline throw_invalid() - end # the ULP error is usually not required to great accuracy, so `Float32` should be precise enough zero_return = Float32(0) inf_return = Float32(Inf) + let accur_is_nan = isnan(accurate), approx_is_nan = isnan(approximate) + if accur_is_nan || approx_is_nan + if accur_is_nan === approx_is_nan + return zero_return + end + return inf_return + end + end if isinf(accurate) || iszero(accurate) # handle floating-point edge cases if isinf(accurate) if isinf(approximate) && (signbit(accurate) == signbit(approximate)) From b407ae2aa0fcc787db0309917daf52ea1d8e8a33 Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Wed, 6 Aug 2025 16:03:09 +0200 Subject: [PATCH 4/9] update `ULPError` --- test/common/ULPError.jl | 41 ++++++++++++++++------------------------- 1 file changed, 16 insertions(+), 25 deletions(-) diff --git a/test/common/ULPError.jl b/test/common/ULPError.jl index f788bd9..d4d9280 100644 --- a/test/common/ULPError.jl +++ b/test/common/ULPError.jl @@ -1,32 +1,27 @@ module ULPError export ulp_error, ulp_error_maximum - @noinline function throw_invalid() - throw(ArgumentError("invalid")) - end function ulp_error(accurate::AbstractFloat, approximate::AbstractFloat) # the ULP error is usually not required to great accuracy, so `Float32` should be precise enough - zero_return = Float32(0) - inf_return = Float32(Inf) - let accur_is_nan = isnan(accurate), approx_is_nan = isnan(approximate) + zero_return = 0f0 + inf_return = Inf32 + # handle floating-point edge cases + if !(isfinite(accurate) && isfinite(approximate)) + accur_is_nan = isnan(accurate) + approx_is_nan = isnan(approximate) if accur_is_nan || approx_is_nan - if accur_is_nan === approx_is_nan - return zero_return + return if accur_is_nan === approx_is_nan + zero_return + else + inf_return end - return inf_return end - end - if isinf(accurate) || iszero(accurate) # handle floating-point edge cases - if isinf(accurate) - if isinf(approximate) && (signbit(accurate) == signbit(approximate)) - return zero_return + if isinf(approximate) + return if isinf(accurate) && (signbit(accurate) == signbit(approximate)) + zero_return + else + inf_return end - return inf_return end - # `iszero(accurate)` - if iszero(approximate) - return zero_return - end - return inf_return end # assuming `precision(BigFloat)` is great enough acc = if accurate isa BigFloat @@ -34,11 +29,7 @@ module ULPError else BigFloat(accurate)::BigFloat end - err = abs(Float32((approximate - acc) / eps(approximate))::Float32) - if isnan(err) - @noinline throw_invalid() # unexpected - end - err + abs(Float32((approximate - acc) / eps(approximate))::Float32) end function ulp_error(accurate::Acc, approximate::App, x::AbstractFloat) where {Acc, App} acc = accurate(x) From 200ee6cb1ba75f3a9a5eea88cde048a48c191c7a Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Thu, 7 Aug 2025 21:56:48 +0200 Subject: [PATCH 5/9] update `ULPError` oscardssmith has been reviewing the `ULPError` code on https://github.com/JuliaLang/julia/pull/59087 So sync the code with the changes there. --- test/common/ULPError.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/common/ULPError.jl b/test/common/ULPError.jl index d4d9280..438ccb9 100644 --- a/test/common/ULPError.jl +++ b/test/common/ULPError.jl @@ -23,11 +23,11 @@ module ULPError end end end - # assuming `precision(BigFloat)` is great enough - acc = if accurate isa BigFloat - accurate + acc = if accurate isa Union{Float16, Float32} + # widen for better accuracy when doing so does not impact performance too much + widen(accurate) else - BigFloat(accurate)::BigFloat + accurate end abs(Float32((approximate - acc) / eps(approximate))::Float32) end From c896b880727e48de2806da4aa1f1b9c7e416b02e Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 9 Aug 2025 04:58:29 +0200 Subject: [PATCH 6/9] `ulp_error_maximum`: use `Fix1` --- test/common/ULPError.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/test/common/ULPError.jl b/test/common/ULPError.jl index 438ccb9..e45eaa3 100644 --- a/test/common/ULPError.jl +++ b/test/common/ULPError.jl @@ -40,9 +40,6 @@ module ULPError ulp_error(func ∘ BigFloat, func, x) end function ulp_error_maximum(func::Func, iterator) where {Func} - function f(x::AbstractFloat) - ulp_error(func, x) - end - maximum(f, iterator) + maximum(Base.Fix1(ulp_error, func), iterator) end end From 0e33deb0f68a820854aa9b67ed9b56a23bbae860 Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 9 Aug 2025 05:00:12 +0200 Subject: [PATCH 7/9] test suite: use `ulp_error_maximum` instead of `for` loop --- test/basicfuns.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/basicfuns.jl b/test/basicfuns.jl index 0b76115..af0ee66 100644 --- a/test/basicfuns.jl +++ b/test/basicfuns.jl @@ -88,9 +88,7 @@ end @test logit(logistic(2)) ≈ 2.0 @testset "accuracy of `logit`" begin for t in (Float16, Float32, Float64) - for x in range(start = t(0), stop = t(1), length = 500) - @test 2 * ulp_error(logit, x) < 3 - end + @test 2 * ulp_error_maximum(logit, range(start = t(0), stop = t(1), length = 500)) < 3 end end end From e7e62267d02f1f4bca0a9343158ecfca10da291d Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 9 Aug 2025 05:01:46 +0200 Subject: [PATCH 8/9] `ULPError`: adjust whitespace style as per devmotion --- test/common/ULPError.jl | 82 ++++++++++++++++++++++------------------- 1 file changed, 44 insertions(+), 38 deletions(-) diff --git a/test/common/ULPError.jl b/test/common/ULPError.jl index e45eaa3..d9e64a2 100644 --- a/test/common/ULPError.jl +++ b/test/common/ULPError.jl @@ -1,45 +1,51 @@ module ULPError - export ulp_error, ulp_error_maximum - function ulp_error(accurate::AbstractFloat, approximate::AbstractFloat) - # the ULP error is usually not required to great accuracy, so `Float32` should be precise enough - zero_return = 0f0 - inf_return = Inf32 - # handle floating-point edge cases - if !(isfinite(accurate) && isfinite(approximate)) - accur_is_nan = isnan(accurate) - approx_is_nan = isnan(approximate) - if accur_is_nan || approx_is_nan - return if accur_is_nan === approx_is_nan - zero_return - else - inf_return - end - end - if isinf(approximate) - return if isinf(accurate) && (signbit(accurate) == signbit(approximate)) - zero_return - else - inf_return - end + +export ulp_error, ulp_error_maximum + +function ulp_error(accurate::AbstractFloat, approximate::AbstractFloat) + # the ULP error is usually not required to great accuracy, so `Float32` should be precise enough + zero_return = 0f0 + inf_return = Inf32 + # handle floating-point edge cases + if !(isfinite(accurate) && isfinite(approximate)) + accur_is_nan = isnan(accurate) + approx_is_nan = isnan(approximate) + if accur_is_nan || approx_is_nan + return if accur_is_nan === approx_is_nan + zero_return + else + inf_return end end - acc = if accurate isa Union{Float16, Float32} - # widen for better accuracy when doing so does not impact performance too much - widen(accurate) - else - accurate + if isinf(approximate) + return if isinf(accurate) && (signbit(accurate) == signbit(approximate)) + zero_return + else + inf_return + end end - abs(Float32((approximate - acc) / eps(approximate))::Float32) - end - function ulp_error(accurate::Acc, approximate::App, x::AbstractFloat) where {Acc, App} - acc = accurate(x) - app = approximate(x) - ulp_error(acc, app) - end - function ulp_error(func::Func, x::AbstractFloat) where {Func} - ulp_error(func ∘ BigFloat, func, x) end - function ulp_error_maximum(func::Func, iterator) where {Func} - maximum(Base.Fix1(ulp_error, func), iterator) + acc = if accurate isa Union{Float16, Float32} + # widen for better accuracy when doing so does not impact performance too much + widen(accurate) + else + accurate end + abs(Float32((approximate - acc) / eps(approximate))::Float32) +end + +function ulp_error(accurate::Acc, approximate::App, x::AbstractFloat) where {Acc, App} + acc = accurate(x) + app = approximate(x) + ulp_error(acc, app) +end + +function ulp_error(func::Func, x::AbstractFloat) where {Func} + ulp_error(func ∘ BigFloat, func, x) +end + +function ulp_error_maximum(func::Func, iterator) where {Func} + maximum(Base.Fix1(ulp_error, func), iterator) +end + end From a1eb57eaf29620ec51c7393b976329c18c3d1f68 Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 9 Aug 2025 05:25:18 +0200 Subject: [PATCH 9/9] `ULPError`: delete unnecessary method static parameters Specialization on the callable argument types should happen anyway. --- test/common/ULPError.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/common/ULPError.jl b/test/common/ULPError.jl index d9e64a2..ee62eca 100644 --- a/test/common/ULPError.jl +++ b/test/common/ULPError.jl @@ -34,7 +34,7 @@ function ulp_error(accurate::AbstractFloat, approximate::AbstractFloat) abs(Float32((approximate - acc) / eps(approximate))::Float32) end -function ulp_error(accurate::Acc, approximate::App, x::AbstractFloat) where {Acc, App} +function ulp_error(accurate, approximate, x::AbstractFloat) acc = accurate(x) app = approximate(x) ulp_error(acc, app)