From 96b5f3d9603fb47ebde47a9546f69daf626c2d95 Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Fri, 1 Aug 2025 22:46:38 +0200 Subject: [PATCH 01/10] fix accuracy of `logcosh(::Union{Float32, Float64})` Fixes #100 --- src/basicfuns.jl | 54 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/src/basicfuns.jl b/src/basicfuns.jl index 5cc0df2..a856077 100644 --- a/src/basicfuns.jl +++ b/src/basicfuns.jl @@ -131,9 +131,63 @@ The implementation ensures `logcosh(-x) = logcosh(x)`. """ function logcosh(x::Real) abs_x = abs(x) + if (x isa Union{Float32, Float64}) && (abs_x < oftype(x, 0.7373046875)) + return logcosh_ker(x) + end return abs_x + log1pexp(- 2 * abs_x) - IrrationalConstants.logtwo end +""" + logcosh_ker(x::Union{Float32, Float64}) + +The kernel of `logcosh`. + +The polynomial coefficients were found using Sollya: + +```sollya +prec = 500!; +points = 50001!; +accurate = log(cosh(x)); +domain = [-0.125, 0.7373046875]; +constrained_part = (x^2) / 2; +free_monomials_32 = [|4, 6, 8, 10, 12|]; +free_monomials_64 = [|4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24|]; +polynomial_32 = fpminimax(accurate, free_monomials_32, [|single...|], domain, constrained_part); +polynomial_64 = fpminimax(accurate, free_monomials_64, [|double...|], domain, constrained_part); +polynomial_32; +polynomial_64; +``` +""" +function logcosh_ker(x::Union{Float32, Float64}) + x² = x * x + if x isa Float32 + p = ( + 5f-1, + -0.083333164f0, + 0.022217678f0, + -0.0067060017f0, + 0.0020296266f0, + -0.00044135848f0, + ) + elseif x isa Float64 + p = ( + 5e-1, + -0.08333333333332801, + 0.02222222222164912, + -0.0067460317245250445, + 0.0021869484500251714, + -0.0007385985435311435, + 0.0002565500026777061, + -9.084985367586575e-5, + 3.2348259905568986e-5, + -1.1058814347469105e-5, + 3.16293199955507e-6, + -5.312230207322749e-7, + ) + end + evalpoly(x², p) * x² +end + """ $(SIGNATURES) From 596ec9579aae6fcca9109f7fa08ccc0274dd94fe Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 2 Aug 2025 14:56:52 +0200 Subject: [PATCH 02/10] add ULPError.jl --- test/common/ULPError.jl | 55 +++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 3 +++ 2 files changed, 58 insertions(+) create mode 100644 test/common/ULPError.jl 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 3656499c4a82a40abae655bfa03dadbd2394909a Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 2 Aug 2025 15:00:19 +0200 Subject: [PATCH 03/10] test accuracy --- test/basicfuns.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/basicfuns.jl b/test/basicfuns.jl index 6a84b54..cbc6157 100644 --- a/test/basicfuns.jl +++ b/test/basicfuns.jl @@ -107,6 +107,14 @@ end @test @inferred(logcosh(x)) === x @test @inferred(logabssinh(x)) === x end + + @testset "accuracy of `logcosh`" begin + for t in (Float32, Float64) + for x in range(start = t(-3), stop = t(3), length = 1000) + @test ulp_error(logcosh, x) < 3 + end + end + end end @testset "log1psq" begin From 30430c4f28f31ef02479d1dcf1992f1e893578fa Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Tue, 5 Aug 2025 14:35:24 +0200 Subject: [PATCH 04/10] 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 03037baa6ad4e4fe8e1f28501c02604a7ae6178c Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Wed, 6 Aug 2025 16:04:37 +0200 Subject: [PATCH 05/10] 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 4fc1b9ddfb31c7caf1ef79b05e7f037b300cfb1a Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Thu, 7 Aug 2025 21:59:22 +0200 Subject: [PATCH 06/10] 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 12ed7ce6c9c8408ea5de8561662861191f95995e Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 9 Aug 2025 04:38:13 +0200 Subject: [PATCH 07/10] `ulp_error_maximum`: use `Fix1` As suggested by devmotion on PR #99. --- 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 714e08a2a471965956656479caeaeb9f834c30e7 Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 9 Aug 2025 04:40:03 +0200 Subject: [PATCH 08/10] test suite: use `ulp_error_maximum` instead of `for` loop As suggested by devmotion on PR #99. --- test/basicfuns.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/basicfuns.jl b/test/basicfuns.jl index cbc6157..c4601d2 100644 --- a/test/basicfuns.jl +++ b/test/basicfuns.jl @@ -110,9 +110,7 @@ end @testset "accuracy of `logcosh`" begin for t in (Float32, Float64) - for x in range(start = t(-3), stop = t(3), length = 1000) - @test ulp_error(logcosh, x) < 3 - end + @test ulp_error_maximum(logcosh, range(start = t(-3), stop = t(3), length = 1000)) < 3 end end end From 297586b3b35edd3005d32b749d961dfed02ee1c2 Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 9 Aug 2025 04:51:05 +0200 Subject: [PATCH 09/10] `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 b754b402c75b80b0e7feeea3995b0fc4797a20b5 Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sat, 9 Aug 2025 05:31:17 +0200 Subject: [PATCH 10/10] `ulp_error`: 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)