test accuracy

Author: nsajko

test/math.jl

@@ -3,6 +3,9 @@
 include("testhelpers/EvenIntegers.jl")
 using .EvenIntegers
 
+include("testhelpers/ULPError.jl")
+using .ULPError
+
 using Random
 using LinearAlgebra
 using Base.Experimental: @force_compile
@@ -652,6 +655,14 @@ end
         @test cosc(big"0.5") ≈ big"-1.273239544735162686151070106980114896275677165923651589981338752471174381073817" rtol=1e-76
         @test cosc(big"0.499") ≈ big"-1.272045747741181369948389133250213864178198918667041860771078493955590574971317" rtol=1e-76
     end
+
+    @testset "accuracy of `cosc` around the origin" begin
+        for t in (Float32, Float64)
+            for x in range(start = t(-1), stop = t(1), length = 5000)
+                @test ulp_error(cosc, x) < 4
+            end
+        end
+    end
 end
 
 @testset "Irrational args to sinpi/cospi/tanpi/sinc/cosc" begin

test/testhelpers/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