From 11134a8cc33e02510258374befae638a3013e790 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Mon, 21 Jul 2025 00:13:01 +0200
Subject: [PATCH 01/33] a bit of code deduplication

---
 base/special/trig.jl | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/base/special/trig.jl b/base/special/trig.jl
index 2137f1486204d..eb13a1fc49301 100644
--- a/base/special/trig.jl
+++ b/base/special/trig.jl
@@ -1104,6 +1104,9 @@ Return a `T(NaN)` if `isnan(x)`.
 See also [`sinc`](@ref).
 """
 cosc(x::Number) = _cosc(float(x))
+function _cosc_generic(x)
+    ((pi*x)*cospi(x)-sinpi(x))/((pi*x)*x)
+end
 function _cosc(x::Number)
     # naive cosc formula is susceptible to catastrophic
     # cancellation error near x=0, so we use the Taylor series
@@ -1124,17 +1127,17 @@ function _cosc(x::Number)
         end
         return π*s
     else
-        return isinf_real(x) ? zero(x) : ((pi*x)*cospi(x)-sinpi(x))/((pi*x)*x)
+        return isinf_real(x) ? zero(x) : _cosc_generic(x)
     end
 end
 # hard-code Float64/Float32 Taylor series, with coefficients
 #  Float64.([(-1)^n*big(pi)^(2n)/((2n+1)*factorial(2n-1)) for n = 1:6])
 _cosc(x::Union{Float64,ComplexF64}) =
     fastabs(x) < 0.14 ? x*evalpoly(x^2, (-3.289868133696453, 3.2469697011334144, -1.1445109447325053, 0.2091827825412384, -0.023460810354558236, 0.001781145516372852)) :
-    isinf_real(x) ? zero(x) : ((pi*x)*cospi(x)-sinpi(x))/((pi*x)*x)
+    isinf_real(x) ? zero(x) : _cosc_generic(x)
 _cosc(x::Union{Float32,ComplexF32}) =
     fastabs(x) < 0.26f0 ? x*evalpoly(x^2, (-3.289868f0, 3.2469697f0, -1.144511f0, 0.20918278f0)) :
-    isinf_real(x) ? zero(x) : ((pi*x)*cospi(x)-sinpi(x))/((pi*x)*x)
+    isinf_real(x) ? zero(x) : _cosc_generic(x)
 _cosc(x::Float16) = Float16(_cosc(Float32(x)))
 _cosc(x::ComplexF16) = ComplexF16(_cosc(ComplexF32(x)))
 

From 2e17039a9781d6b8e91e29833e62f5eae912d40e Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Mon, 21 Jul 2025 18:06:24 +0200
Subject: [PATCH 02/33] eliminate common subexpression

---
 base/special/trig.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/base/special/trig.jl b/base/special/trig.jl
index eb13a1fc49301..d1488ff45d9cd 100644
--- a/base/special/trig.jl
+++ b/base/special/trig.jl
@@ -1105,7 +1105,8 @@ See also [`sinc`](@ref).
 """
 cosc(x::Number) = _cosc(float(x))
 function _cosc_generic(x)
-    ((pi*x)*cospi(x)-sinpi(x))/((pi*x)*x)
+    pi_x = pi * x
+    (pi_x*cospi(x)-sinpi(x))/(pi_x*x)
 end
 function _cosc(x::Number)
     # naive cosc formula is susceptible to catastrophic

From 0011a4f73ddbbb84cd114a68aaaf158c5dc227c5 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Mon, 21 Jul 2025 18:09:00 +0200
Subject: [PATCH 03/33] better accuracy for `cosc` for `Float32`, `Float64`
 around `x = 0`

Improve accuracy for `cosc(::Float32)` and `cosc(::Float64)`. This is
accomplished by:

* Using minimax instead of Taylor polynomials. They're more efficient.
  The minimax polynomials are found using Sollya.

* Using multiple polynomials for different subintervals of the domain.
  Known as "domain splitting". Enables good accuracy for a wider region
  around the origin than possible with only a single polynomial.

The polynomial degrees are kept as-is.

Fixes #59065
---
 base/special/trig.jl | 130 ++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 128 insertions(+), 2 deletions(-)

diff --git a/base/special/trig.jl b/base/special/trig.jl
index d1488ff45d9cd..688ea32fcc86c 100644
--- a/base/special/trig.jl
+++ b/base/special/trig.jl
@@ -1131,12 +1131,138 @@ function _cosc(x::Number)
         return isinf_real(x) ? zero(x) : _cosc_generic(x)
     end
 end
+
+#=
+
+## `cosc(x)` for `x` around the first zero, at `x = 0`
+
+`Float32`:
+
+```sollya
+prec = 500!;
+accurate = ((pi * x) * cos(pi * x) - sin(pi * x)) / (pi * x * x);
+b1 = 0.27001953125;
+b2 = 0.449951171875;
+domain_0 = [-b1/2, b1];
+domain_1 = [b1, b2];
+machinePrecision = 24;
+freeMonomials = [|1, 3, 5, 7|];
+freeMonomialPrecisions = [|machinePrecision, machinePrecision, machinePrecision, machinePrecision|];
+polynomial_0 = fpminimax(accurate, freeMonomials, freeMonomialPrecisions, domain_0);
+polynomial_1 = fpminimax(accurate, freeMonomials, freeMonomialPrecisions, domain_1);
+polynomial_0;
+polynomial_1;
+```
+
+`Float64`:
+
+```sollya
+prec = 500!;
+accurate = ((pi * x) * cos(pi * x) - sin(pi * x)) / (pi * x * x);
+b1 = 0.1700439453125;
+b2 = 0.27001953125;
+b3 = 0.340087890625;
+b4 = 0.39990234375;
+domain_0 = [-b1/2, b1];
+domain_1 = [b1, b2];
+domain_2 = [b2, b3];
+domain_3 = [b3, b4];
+machinePrecision = 53;
+freeMonomials = [|1, 3, 5, 7, 9, 11|];
+freeMonomialPrecisions = [|machinePrecision, machinePrecision, machinePrecision, machinePrecision, machinePrecision, machinePrecision|];
+polynomial_0 = fpminimax(accurate, freeMonomials, freeMonomialPrecisions, domain_0);
+polynomial_1 = fpminimax(accurate, freeMonomials, freeMonomialPrecisions, domain_1);
+polynomial_2 = fpminimax(accurate, freeMonomials, freeMonomialPrecisions, domain_2);
+polynomial_3 = fpminimax(accurate, freeMonomials, freeMonomialPrecisions, domain_3);
+polynomial_0;
+polynomial_1;
+polynomial_2;
+polynomial_3;
+```
+
+=#
+
+struct _CosCardinalEvaluationScheme{
+    T <: AbstractFloat,
+    PolynomialsCloseToZero <: (Tuple{Tuple{T, P}, Vararg{Tuple{T, P}}} where {P <: Tuple{T, Vararg{T}}}),
+}
+    polynomials_close_to_origin::PolynomialsCloseToZero
+end
+
+function (sch::_CosCardinalEvaluationScheme)(x::AbstractFloat)
+    a = abs(x)
+    polynomials_close_to_origin = sch.polynomials_close_to_origin
+    if (polynomials_close_to_origin !== ()) && (a ≤ polynomials_close_to_origin[end][1])
+        let
+            if length(polynomials_close_to_origin) == 1  # hardcode for each allowed tuple size
+                p = only(polynomials_close_to_origin)[2]
+            elseif length(polynomials_close_to_origin) == 2
+                p = let ((b1, p0), (_, p1)) = polynomials_close_to_origin
+                    if a ≤ b1  # hardcoded binary search
+                        p0
+                    else
+                        p1
+                    end
+                end
+            elseif length(polynomials_close_to_origin) == 4
+                p = let ((b1, p0), (b2, p1), (b3, p2), (_, p3)) = polynomials_close_to_origin
+                    if a ≤ b2  # hardcoded binary search
+                        if a ≤ b1
+                            p0
+                        else
+                            p1
+                        end
+                    else
+                        if a ≤ b3
+                            p2
+                        else
+                            p3
+                        end
+                    end
+                end
+            end
+            x * evalpoly(x * x, p)
+        end
+    elseif isinf(x)
+        typeof(x)(0)
+    else
+        _cosc_generic(x)
+    end
+end
+
+const _cosc_f32 = let b = Float32 ∘ Float16
+    _CosCardinalEvaluationScheme(
+        (
+            (b(0.27), (-3.289868f0, 3.246966f0, -1.1443111f0, 0.20542027f0)),
+            (b(0.45), (-3.2898617f0, 3.2467577f0, -1.1420113f0, 0.1965574f0)),
+        ),
+    )
+end
+
+const _cosc_f64 = let b = Float64 ∘ Float16
+    _CosCardinalEvaluationScheme(
+        (
+            (b(0.17), (-3.289868133696453, 3.2469697011333203, -1.1445109446992934, 0.20918277797812262, -0.023460519561502552, 0.001772485141534688)),
+            (b(0.27), (-3.289868133695205, 3.246969700970421, -1.1445109360543062, 0.20918254132488637, -0.023457115021035743, 0.0017515112964895303)),
+            (b(0.34), (-3.289868133634355, 3.246969697075094, -1.1445108347839286, 0.209181201609773, -0.023448079433318045, 0.001726628430505518)),
+            (b(0.4),  (-3.289868133074254, 3.2469696736659346, -1.1445104406286049, 0.20917785794416457, -0.02343378376047161, 0.0017019796223768677)),
+        ),
+    )
+end
+
+function _cosc(x::Float32)
+    _cosc_f32(x)
+end
+function _cosc(x::Float64)
+    _cosc_f64(x)
+end
+
 # hard-code Float64/Float32 Taylor series, with coefficients
 #  Float64.([(-1)^n*big(pi)^(2n)/((2n+1)*factorial(2n-1)) for n = 1:6])
-_cosc(x::Union{Float64,ComplexF64}) =
+_cosc(x::ComplexF64) =
     fastabs(x) < 0.14 ? x*evalpoly(x^2, (-3.289868133696453, 3.2469697011334144, -1.1445109447325053, 0.2091827825412384, -0.023460810354558236, 0.001781145516372852)) :
     isinf_real(x) ? zero(x) : _cosc_generic(x)
-_cosc(x::Union{Float32,ComplexF32}) =
+_cosc(x::ComplexF32) =
     fastabs(x) < 0.26f0 ? x*evalpoly(x^2, (-3.289868f0, 3.2469697f0, -1.144511f0, 0.20918278f0)) :
     isinf_real(x) ? zero(x) : _cosc_generic(x)
 _cosc(x::Float16) = Float16(_cosc(Float32(x)))

From a3c4a73504e293784b29d3006433e2fb4b70fc58 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Mon, 4 Aug 2025 17:35:25 +0200
Subject: [PATCH 04/33] test accuracy

---
 test/math.jl                 | 11 ++++++++
 test/testhelpers/ULPError.jl | 55 ++++++++++++++++++++++++++++++++++++
 2 files changed, 66 insertions(+)
 create mode 100644 test/testhelpers/ULPError.jl

diff --git a/test/math.jl b/test/math.jl
index 529f4e01405a9..e3fe9a5a6d47e 100644
--- a/test/math.jl
+++ b/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
@@ -685,6 +688,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
diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
new file mode 100644
index 0000000000000..abb44a83e121a
--- /dev/null
+++ b/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

From 0ab1f6b05e448a4e2603638a931b88d6dc9ca78a Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Mon, 4 Aug 2025 17:52:14 +0200
Subject: [PATCH 05/33] add copyright notice

---
 test/testhelpers/ULPError.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index abb44a83e121a..7b13bdb23f67c 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -1,3 +1,5 @@
+# This file is a part of Julia. License is MIT: https://julialang.org/license
+
 module ULPError
     export ulp_error, ulp_error_maximum
     @noinline function throw_invalid()

From d991f1ef60e31e422f7dd0ea924f6ab909d6607d Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Mon, 4 Aug 2025 22:20:06 +0200
Subject: [PATCH 06/33] ULPError simplification as per Oscar's suggestion

---
 test/testhelpers/ULPError.jl | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 7b13bdb23f67c..f682c0389712c 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -6,15 +6,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 4459c67a7f8699d2d5be7d6995c5347c4e4bf8af Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 13:45:03 +0200
Subject: [PATCH 07/33] ULPError common subexpression elimination

---
 test/testhelpers/ULPError.jl | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index f682c0389712c..8f33ba3010a30 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -17,18 +17,20 @@ module ULPError
                 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))
+        let accur_is_inf = isinf(accurate)
+            if accur_is_inf || iszero(accurate)  # handle floating-point edge cases
+                if accur_is_inf
+                    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
-            # `iszero(accurate)`
-            if iszero(approximate)
-                return zero_return
-            end
-            return inf_return
         end
         # assuming `precision(BigFloat)` is great enough
         acc = if accurate isa BigFloat

From 2592ee90d36a277f57ecaa314cbea6730dc4a541 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 13:55:35 +0200
Subject: [PATCH 08/33] ULPError: handle `isinf(approximate)`

---
 test/testhelpers/ULPError.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 8f33ba3010a30..94ce2524da812 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -17,10 +17,11 @@ module ULPError
                 return inf_return
             end
         end
+        approx_is_inf = isinf(approximate)
         let accur_is_inf = isinf(accurate)
             if accur_is_inf || iszero(accurate)  # handle floating-point edge cases
                 if accur_is_inf
-                    if isinf(approximate) && (signbit(accurate) == signbit(approximate))
+                    if approx_is_inf && (signbit(accurate) == signbit(approximate))
                         return zero_return
                     end
                     return inf_return
@@ -32,6 +33,9 @@ module ULPError
                 return inf_return
             end
         end
+        if approx_is_inf
+            return inf_return
+        end
         # assuming `precision(BigFloat)` is great enough
         acc = if accurate isa BigFloat
             accurate

From 48b809286d98ae5e417ac4d856b108c503ce125a Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 14:05:52 +0200
Subject: [PATCH 09/33] put comment in correct position

---
 test/testhelpers/ULPError.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 94ce2524da812..e685e98588173 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -9,6 +9,7 @@ module ULPError
         # the ULP error is usually not required to great accuracy, so `Float32` should be precise enough
         zero_return = Float32(0)
         inf_return = Float32(Inf)
+        # handle floating-point edge cases
         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
@@ -19,7 +20,7 @@ module ULPError
         end
         approx_is_inf = isinf(approximate)
         let accur_is_inf = isinf(accurate)
-            if accur_is_inf || iszero(accurate)  # handle floating-point edge cases
+            if accur_is_inf || iszero(accurate)
                 if accur_is_inf
                     if approx_is_inf && (signbit(accurate) == signbit(approximate))
                         return zero_return

From 7140f6055924641a27f625182c6ce9cfa0dc178f Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 14:09:58 +0200
Subject: [PATCH 10/33] style: get rid of the `let`s/nesting

---
 test/testhelpers/ULPError.jl | 33 ++++++++++++++++-----------------
 1 file changed, 16 insertions(+), 17 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index e685e98588173..5a6d48977ec0c 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -10,29 +10,28 @@ module ULPError
         zero_return = Float32(0)
         inf_return = Float32(Inf)
         # handle floating-point edge cases
-        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
+        accur_is_nan = isnan(accurate)
+        accur_is_inf = isinf(accurate)
+        approx_is_inf = isinf(approximate)
+        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
-        approx_is_inf = isinf(approximate)
-        let accur_is_inf = isinf(accurate)
-            if accur_is_inf || iszero(accurate)
-                if accur_is_inf
-                    if approx_is_inf && (signbit(accurate) == signbit(approximate))
-                        return zero_return
-                    end
-                    return inf_return
-                end
-                # `iszero(accurate)`
-                if iszero(approximate)
+        if accur_is_inf || iszero(accurate)
+            if accur_is_inf
+                if approx_is_inf && (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
         if approx_is_inf
             return inf_return

From 11c29f3d9ab23b99d3dcc60f369dd1bc3abe1ed1 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 14:16:33 +0200
Subject: [PATCH 11/33] style: simplify `return`s

---
 test/testhelpers/ULPError.jl | 21 ++++++++-------------
 1 file changed, 8 insertions(+), 13 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 5a6d48977ec0c..ea7d51a3d260c 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -15,23 +15,18 @@ module ULPError
         approx_is_inf = isinf(approximate)
         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
         if accur_is_inf || iszero(accurate)
-            if accur_is_inf
-                if approx_is_inf && (signbit(accurate) == signbit(approximate))
-                    return zero_return
-                end
-                return inf_return
-            end
-            # `iszero(accurate)`
-            if iszero(approximate)
-                return zero_return
+            return if (approx_is_inf && (signbit(accurate) == signbit(approximate))) || iszero(approximate)
+                zero_return
+            else
+                inf_return
             end
-            return inf_return
         end
         if approx_is_inf
             return inf_return

From fb2ac2f31c29a67c46cb978b34ddedff51b62ffc Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 14:41:55 +0200
Subject: [PATCH 12/33] fix `ULPError` edge case regarding zero/inf

---
 test/testhelpers/ULPError.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index ea7d51a3d260c..972f7c04023e4 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -10,10 +10,12 @@ module ULPError
         zero_return = Float32(0)
         inf_return = Float32(Inf)
         # handle floating-point edge cases
+        accur_is_zero = iszero(accurate)
         accur_is_nan = isnan(accurate)
         accur_is_inf = isinf(accurate)
         approx_is_inf = isinf(approximate)
         approx_is_nan = isnan(approximate)
+        approx_is_zero = iszero(approximate)
         if accur_is_nan || approx_is_nan
             return if accur_is_nan === approx_is_nan
                 zero_return
@@ -21,8 +23,8 @@ module ULPError
                 inf_return
             end
         end
-        if accur_is_inf || iszero(accurate)
-            return if (approx_is_inf && (signbit(accurate) == signbit(approximate))) || iszero(approximate)
+        if accur_is_inf || accur_is_zero
+            return if (accur_is_inf && approx_is_inf && (signbit(accurate) == signbit(approximate))) || (accur_is_zero && approx_is_zero)
                 zero_return
             else
                 inf_return

From c65115a2bdd445948efdaf4aed1e36b2cbd65ed3 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 14:31:33 +0200
Subject: [PATCH 13/33] test `ULPError` edge cases

---
 test/math.jl | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/test/math.jl b/test/math.jl
index e3fe9a5a6d47e..8007abe649ecf 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -28,6 +28,38 @@ has_fma = Dict(
     BigFloat => true,
 )
 
+@testset "meta test: ULPError" begin
+    @testset "edge cases" begin
+        @testset "zero error - two equivalent values" begin
+            @test 0 == @inferred ulp_error(NaN, NaN)
+            @test 0 == @inferred ulp_error(-NaN, -NaN)
+            @test 0 == @inferred ulp_error(-NaN, NaN)
+            @test 0 == @inferred ulp_error(NaN, -NaN)
+            @test 0 == @inferred ulp_error(Inf, Inf)
+            @test 0 == @inferred ulp_error(-Inf, -Inf)
+            @test 0 == @inferred ulp_error(0.0, 0.0)
+            @test 0 == @inferred ulp_error(-0.0, -0.0)
+            @test 0 == @inferred ulp_error(-0.0, 0.0)
+            @test 0 == @inferred ulp_error(0.0, -0.0)
+            @test 0 == @inferred ulp_error(3.0, 3.0)
+            @test 0 == @inferred ulp_error(-3.0, -3.0)
+        end
+        @testset "infinite error" begin
+            @test Inf == @inferred ulp_error(NaN, 0.0)
+            @test Inf == @inferred ulp_error(0.0, NaN)
+            @test Inf == @inferred ulp_error(NaN, 3.0)
+            @test Inf == @inferred ulp_error(3.0, NaN)
+            @test Inf == @inferred ulp_error(NaN, Inf)
+            @test Inf == @inferred ulp_error(Inf, -Inf)
+            @test Inf == @inferred ulp_error(-Inf, Inf)
+            @test Inf == @inferred ulp_error(0.0, Inf)
+            @test Inf == @inferred ulp_error(Inf, 0.0)
+            @test Inf == @inferred ulp_error(3.0, Inf)
+            @test Inf == @inferred ulp_error(Inf, 3.0)
+        end
+    end
+end
+
 @testset "clamp" begin
     let
         @test clamp(0, 1, 3) == 1

From 33717f57e5d960b7ef1f4561020a5fcec675b96f Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 14:53:17 +0200
Subject: [PATCH 14/33] `ULPError`: simplify

---
 test/testhelpers/ULPError.jl | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 972f7c04023e4..978a5ae042380 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -23,16 +23,13 @@ module ULPError
                 inf_return
             end
         end
-        if accur_is_inf || accur_is_zero
-            return if (accur_is_inf && approx_is_inf && (signbit(accurate) == signbit(approximate))) || (accur_is_zero && approx_is_zero)
+        if approx_is_inf
+            return if accur_is_inf && (signbit(accurate) == signbit(approximate))
                 zero_return
             else
                 inf_return
             end
         end
-        if approx_is_inf
-            return inf_return
-        end
         # assuming `precision(BigFloat)` is great enough
         acc = if accurate isa BigFloat
             accurate

From f67706265e8f2c94acb138343738dc5c773c025f Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 14:56:33 +0200
Subject: [PATCH 15/33] `ULPError`: delete dead code

---
 test/testhelpers/ULPError.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 978a5ae042380..3c3f77568f225 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -10,12 +10,10 @@ module ULPError
         zero_return = Float32(0)
         inf_return = Float32(Inf)
         # handle floating-point edge cases
-        accur_is_zero = iszero(accurate)
         accur_is_nan = isnan(accurate)
         accur_is_inf = isinf(accurate)
         approx_is_inf = isinf(approximate)
         approx_is_nan = isnan(approximate)
-        approx_is_zero = iszero(approximate)
         if accur_is_nan || approx_is_nan
             return if accur_is_nan === approx_is_nan
                 zero_return

From c0340a84885bac5de0309b0870900b59e2b5faa7 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:06:52 +0200
Subject: [PATCH 16/33] `ULPError`: eliminate variables only used once

---
 test/testhelpers/ULPError.jl | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 3c3f77568f225..535c6f9e678f3 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -11,8 +11,6 @@ module ULPError
         inf_return = Float32(Inf)
         # handle floating-point edge cases
         accur_is_nan = isnan(accurate)
-        accur_is_inf = isinf(accurate)
-        approx_is_inf = isinf(approximate)
         approx_is_nan = isnan(approximate)
         if accur_is_nan || approx_is_nan
             return if accur_is_nan === approx_is_nan
@@ -21,8 +19,8 @@ module ULPError
                 inf_return
             end
         end
-        if approx_is_inf
-            return if accur_is_inf && (signbit(accurate) == signbit(approximate))
+        if isinf(approximate)
+            return if isinf(accurate) && (signbit(accurate) == signbit(approximate))
                 zero_return
             else
                 inf_return

From 85ace6cebe664e8e501d79c7efa183e532f44868 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:08:07 +0200
Subject: [PATCH 17/33] `ULPError`: use `Float32` literals

---
 test/testhelpers/ULPError.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 535c6f9e678f3..6638ff73f729d 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -7,8 +7,8 @@ module ULPError
     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)
+        zero_return = 0f0
+        inf_return = Inf32
         # handle floating-point edge cases
         accur_is_nan = isnan(accurate)
         approx_is_nan = isnan(approximate)

From 0b05e22a3cdb4ea36797a631739b6b258de5465b Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:12:59 +0200
Subject: [PATCH 18/33] `ULPError`: guard the edge case handling behind a
 single `if`

---
 test/testhelpers/ULPError.jl | 28 +++++++++++++++-------------
 1 file changed, 15 insertions(+), 13 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 6638ff73f729d..d607089b9b463 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -10,20 +10,22 @@ module ULPError
         zero_return = 0f0
         inf_return = Inf32
         # handle floating-point edge cases
-        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
+        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
-        end
-        if isinf(approximate)
-            return if isinf(accurate) && (signbit(accurate) == signbit(approximate))
-                zero_return
-            else
-                inf_return
+            if isinf(approximate)
+                return if isinf(accurate) && (signbit(accurate) == signbit(approximate))
+                    zero_return
+                else
+                    inf_return
+                end
             end
         end
         # assuming `precision(BigFloat)` is great enough

From 71126d584e55a701cf6d4eb4de6c0cc50b6d1c89 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:19:01 +0200
Subject: [PATCH 19/33] `ULPError`: use `Float64` instead of `BigFloat`

---
 test/testhelpers/ULPError.jl | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index d607089b9b463..176cf205405fa 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -28,12 +28,7 @@ module ULPError
                 end
             end
         end
-        # assuming `precision(BigFloat)` is great enough
-        acc = if accurate isa BigFloat
-            accurate
-        else
-            BigFloat(accurate)::BigFloat
-        end
+        acc = Float64(accurate)::Float64
         err = abs(Float32((approximate - acc) / eps(approximate))::Float32)
         if isnan(err)
             @noinline throw_invalid()  # unexpected

From cc22fc97c5b38c2de0e33609557b35d9b5d6ee94 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:19:40 +0200
Subject: [PATCH 20/33] `ULPError`: avoid unnecessary conversion to `BigFloat`

---
 test/testhelpers/ULPError.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 176cf205405fa..7dec77f53c19b 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -41,7 +41,7 @@ module ULPError
         ulp_error(acc, app)
     end
     function ulp_error(func::Func, x::AbstractFloat) where {Func}
-        ulp_error(func ∘ BigFloat, func, x)
+        ulp_error(func, func, x)
     end
     function ulp_error_maximum(func::Func, iterator) where {Func}
         function f(x::AbstractFloat)

From ca57e60d20614a18c123da07ec46cdc25af875b6 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:21:14 +0200
Subject: [PATCH 21/33] delete debugging throw

Should not be necessary now we have some tests in the test suite.
---
 test/testhelpers/ULPError.jl | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 7dec77f53c19b..ee26ba556c0e1 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -2,9 +2,6 @@
 
 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 = 0f0
@@ -29,11 +26,7 @@ module ULPError
             end
         end
         acc = Float64(accurate)::Float64
-        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 ce8d25d684311f5baaf3b31b1563321a4b257676 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:25:29 +0200
Subject: [PATCH 22/33] test: add missing edge case

---
 test/math.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/math.jl b/test/math.jl
index 8007abe649ecf..fada323a70474 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -50,6 +50,7 @@ has_fma = Dict(
             @test Inf == @inferred ulp_error(NaN, 3.0)
             @test Inf == @inferred ulp_error(3.0, NaN)
             @test Inf == @inferred ulp_error(NaN, Inf)
+            @test Inf == @inferred ulp_error(Inf, NaN)
             @test Inf == @inferred ulp_error(Inf, -Inf)
             @test Inf == @inferred ulp_error(-Inf, Inf)
             @test Inf == @inferred ulp_error(0.0, Inf)

From 8a60bf08bb7d4358a4f8b0cc57d9773bc4e52fc2 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:32:04 +0200
Subject: [PATCH 23/33] test: `ULPError`: faitfhul approximation

---
 test/math.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/test/math.jl b/test/math.jl
index fada323a70474..55322b7c025a8 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -59,6 +59,12 @@ has_fma = Dict(
             @test Inf == @inferred ulp_error(Inf, 3.0)
         end
     end
+    @testset "faithful" begin
+        for x in (-2.0, -0.3, 0.1, 1.0)
+            @test 1 == @inferred ulp_error(x, nextfloat(x, 1))
+            @test 1 == @inferred ulp_error(x, nextfloat(x, -1))
+        end
+    end
 end
 
 @testset "clamp" begin

From a95af7732f9996c3f1d3d9ffd9dc6081687aa58e Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Wed, 6 Aug 2025 15:35:17 +0200
Subject: [PATCH 24/33] Revert "`ULPError`: avoid unnecessary conversion to
 `BigFloat`"

This reverts commit 33eca5fefbf48e94b6d26a533d03dc8d481ac247.
---
 test/testhelpers/ULPError.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index ee26ba556c0e1..9ce3592b57a32 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -34,7 +34,7 @@ module ULPError
         ulp_error(acc, app)
     end
     function ulp_error(func::Func, x::AbstractFloat) where {Func}
-        ulp_error(func, func, x)
+        ulp_error(func ∘ BigFloat, func, x)
     end
     function ulp_error_maximum(func::Func, iterator) where {Func}
         function f(x::AbstractFloat)

From 92a937840a9e07a99b619b4f4fa2ff973f5e91cf Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Thu, 7 Aug 2025 21:24:52 +0200
Subject: [PATCH 25/33] handle precision better

Trust the precision of `accurate` is great enough, unless it's one of
the known-imprecise types.
---
 test/testhelpers/ULPError.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index 9ce3592b57a32..b10c980812709 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -25,7 +25,12 @@ module ULPError
                 end
             end
         end
-        acc = Float64(accurate)::Float64
+        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}

From cee5a83a9c7a1c7a67dbd448ed5b5b80ddadf20f Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Thu, 7 Aug 2025 21:47:51 +0200
Subject: [PATCH 26/33] expand `ULPError` test suite

---
 test/math.jl | 12 +++++++++++-
 1 file changed, 11 insertions(+), 1 deletion(-)

diff --git a/test/math.jl b/test/math.jl
index 55322b7c025a8..759d8ccd1e535 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -29,6 +29,10 @@ has_fma = Dict(
 )
 
 @testset "meta test: ULPError" begin
+    examples_f64 = (-3e0, -2e0, -1e0, -1e-1, -1e-10, -0e0, 0e0, 1e-10, 1e-1, 1e0, 2e0, 3e0)::Tuple{Vararg{Float64}}
+    examples_f16 = Float16.(examples_f64)
+    examples_f32 = Float32.(examples_f64)
+    examples = (examples_f16..., examples_f32..., examples_f64...)
     @testset "edge cases" begin
         @testset "zero error - two equivalent values" begin
             @test 0 == @inferred ulp_error(NaN, NaN)
@@ -60,11 +64,17 @@ has_fma = Dict(
         end
     end
     @testset "faithful" begin
-        for x in (-2.0, -0.3, 0.1, 1.0)
+        for x in examples
             @test 1 == @inferred ulp_error(x, nextfloat(x, 1))
             @test 1 == @inferred ulp_error(x, nextfloat(x, -1))
         end
     end
+    @testset "midpoint" begin
+        for x in examples
+            a = abs(x)
+            @test 1 == 2 * @inferred ulp_error(copysign((widen(a) + nextfloat(a, 1))/2, x), x)
+        end
+    end
 end
 
 @testset "clamp" begin

From e8319acb56980c77875028ae2684559ff0e1f17c Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Sat, 9 Aug 2025 04:26:43 +0200
Subject: [PATCH 27/33] `ulp_error_maximum`: use `Fix1`

As suggested by devmotion on a PR to LogExpFunctions.jl:

* https://github.com/JuliaStats/LogExpFunctions.jl/pull/99
---
 test/testhelpers/ULPError.jl | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index b10c980812709..baad9593703bc 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -42,9 +42,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 1ca04f36caec90b3e98123fa7b5c1fb70d736703 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Sat, 9 Aug 2025 04:29:49 +0200
Subject: [PATCH 28/33] use `ulp_error_maximum` instead of a loop

As suggested by devmotion on a PR to LogExpFunctions.jl:

* https://github.com/JuliaStats/LogExpFunctions.jl/pull/99
---
 test/math.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/test/math.jl b/test/math.jl
index 759d8ccd1e535..6be53cf642bf3 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -740,9 +740,7 @@ 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
+            @test ulp_error_maximum(cosc, range(start = t(-1), stop = t(1), length = 5000)) < 4
         end
     end
 end

From 4f86cf31dd3d354355e941d20c5e1483909e36e2 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Sat, 9 Aug 2025 05:22:11 +0200
Subject: [PATCH 29/33] `ULPError`: delete unnecessary method static parameters

Specialization on the callable argument types should happen anyway.
---
 test/testhelpers/ULPError.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/testhelpers/ULPError.jl b/test/testhelpers/ULPError.jl
index baad9593703bc..6ca7d365ebaff 100644
--- a/test/testhelpers/ULPError.jl
+++ b/test/testhelpers/ULPError.jl
@@ -33,7 +33,7 @@ module ULPError
         end
         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)

From e9331157b8584104521e451667240e8fc6a5b0e8 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Mon, 11 Aug 2025 06:05:35 +0200
Subject: [PATCH 30/33] improve style

---
 base/special/trig.jl | 55 ++++++++++++++++++++------------------------
 1 file changed, 25 insertions(+), 30 deletions(-)

diff --git a/base/special/trig.jl b/base/special/trig.jl
index 688ea32fcc86c..4ea42e625bd62 100644
--- a/base/special/trig.jl
+++ b/base/special/trig.jl
@@ -1190,39 +1190,34 @@ struct _CosCardinalEvaluationScheme{
 end
 
 function (sch::_CosCardinalEvaluationScheme)(x::AbstractFloat)
+    function choose_poly(a::AbstractFloat, polynomials_close_to_origin::NTuple{2})
+        ((b1, p0), (_, p1)) = polynomials_close_to_origin
+        if a ≤ b1
+            p0
+        else
+            p1
+        end
+    end
+    function choose_poly(a::AbstractFloat, polynomials_close_to_origin::NTuple{4})
+        ((b1, p0), (b2, p1), (b3, p2), (_, p3)) = polynomials_close_to_origin
+        if a ≤ b2  # hardcoded binary search
+            if a ≤ b1
+                p0
+            else
+                p1
+            end
+        else
+            if a ≤ b3
+                p2
+            else
+                p3
+            end
+        end
+    end
     a = abs(x)
     polynomials_close_to_origin = sch.polynomials_close_to_origin
     if (polynomials_close_to_origin !== ()) && (a ≤ polynomials_close_to_origin[end][1])
-        let
-            if length(polynomials_close_to_origin) == 1  # hardcode for each allowed tuple size
-                p = only(polynomials_close_to_origin)[2]
-            elseif length(polynomials_close_to_origin) == 2
-                p = let ((b1, p0), (_, p1)) = polynomials_close_to_origin
-                    if a ≤ b1  # hardcoded binary search
-                        p0
-                    else
-                        p1
-                    end
-                end
-            elseif length(polynomials_close_to_origin) == 4
-                p = let ((b1, p0), (b2, p1), (b3, p2), (_, p3)) = polynomials_close_to_origin
-                    if a ≤ b2  # hardcoded binary search
-                        if a ≤ b1
-                            p0
-                        else
-                            p1
-                        end
-                    else
-                        if a ≤ b3
-                            p2
-                        else
-                            p3
-                        end
-                    end
-                end
-            end
-            x * evalpoly(x * x, p)
-        end
+        x * evalpoly(x * x, choose_poly(a, polynomials_close_to_origin))
     elseif isinf(x)
         typeof(x)(0)
     else

From 57fea2b813757c1a0871794032476849f76cbe2c Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Sat, 16 Aug 2025 04:33:24 +0200
Subject: [PATCH 31/33] get rid of the `struct`

---
 base/special/trig.jl | 33 +++++++++++----------------------
 1 file changed, 11 insertions(+), 22 deletions(-)

diff --git a/base/special/trig.jl b/base/special/trig.jl
index 4ea42e625bd62..b7f27b341d166 100644
--- a/base/special/trig.jl
+++ b/base/special/trig.jl
@@ -1182,14 +1182,7 @@ polynomial_3;
 
 =#
 
-struct _CosCardinalEvaluationScheme{
-    T <: AbstractFloat,
-    PolynomialsCloseToZero <: (Tuple{Tuple{T, P}, Vararg{Tuple{T, P}}} where {P <: Tuple{T, Vararg{T}}}),
-}
-    polynomials_close_to_origin::PolynomialsCloseToZero
-end
-
-function (sch::_CosCardinalEvaluationScheme)(x::AbstractFloat)
+function _cos_cardinal_eval(x::AbstractFloat, polynomials_close_to_origin::NTuple)
     function choose_poly(a::AbstractFloat, polynomials_close_to_origin::NTuple{2})
         ((b1, p0), (_, p1)) = polynomials_close_to_origin
         if a ≤ b1
@@ -1226,30 +1219,26 @@ function (sch::_CosCardinalEvaluationScheme)(x::AbstractFloat)
 end
 
 const _cosc_f32 = let b = Float32 ∘ Float16
-    _CosCardinalEvaluationScheme(
-        (
-            (b(0.27), (-3.289868f0, 3.246966f0, -1.1443111f0, 0.20542027f0)),
-            (b(0.45), (-3.2898617f0, 3.2467577f0, -1.1420113f0, 0.1965574f0)),
-        ),
+    (
+        (b(0.27), (-3.289868f0, 3.246966f0, -1.1443111f0, 0.20542027f0)),
+        (b(0.45), (-3.2898617f0, 3.2467577f0, -1.1420113f0, 0.1965574f0)),
     )
 end
 
 const _cosc_f64 = let b = Float64 ∘ Float16
-    _CosCardinalEvaluationScheme(
-        (
-            (b(0.17), (-3.289868133696453, 3.2469697011333203, -1.1445109446992934, 0.20918277797812262, -0.023460519561502552, 0.001772485141534688)),
-            (b(0.27), (-3.289868133695205, 3.246969700970421, -1.1445109360543062, 0.20918254132488637, -0.023457115021035743, 0.0017515112964895303)),
-            (b(0.34), (-3.289868133634355, 3.246969697075094, -1.1445108347839286, 0.209181201609773, -0.023448079433318045, 0.001726628430505518)),
-            (b(0.4),  (-3.289868133074254, 3.2469696736659346, -1.1445104406286049, 0.20917785794416457, -0.02343378376047161, 0.0017019796223768677)),
-        ),
+    (
+        (b(0.17), (-3.289868133696453, 3.2469697011333203, -1.1445109446992934, 0.20918277797812262, -0.023460519561502552, 0.001772485141534688)),
+        (b(0.27), (-3.289868133695205, 3.246969700970421, -1.1445109360543062, 0.20918254132488637, -0.023457115021035743, 0.0017515112964895303)),
+        (b(0.34), (-3.289868133634355, 3.246969697075094, -1.1445108347839286, 0.209181201609773, -0.023448079433318045, 0.001726628430505518)),
+        (b(0.4),  (-3.289868133074254, 3.2469696736659346, -1.1445104406286049, 0.20917785794416457, -0.02343378376047161, 0.0017019796223768677)),
     )
 end
 
 function _cosc(x::Float32)
-    _cosc_f32(x)
+    _cos_cardinal_eval(x, _cosc_f32)
 end
 function _cosc(x::Float64)
-    _cosc_f64(x)
+    _cos_cardinal_eval(x, _cosc_f64)
 end
 
 # hard-code Float64/Float32 Taylor series, with coefficients

From b7adfd3655c2c84abcf69b467b556a4591468f01 Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Sat, 16 Aug 2025 15:10:14 +0200
Subject: [PATCH 32/33] fix

---
 base/special/trig.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/base/special/trig.jl b/base/special/trig.jl
index b7f27b341d166..6858446837b1d 100644
--- a/base/special/trig.jl
+++ b/base/special/trig.jl
@@ -1208,7 +1208,6 @@ function _cos_cardinal_eval(x::AbstractFloat, polynomials_close_to_origin::NTupl
         end
     end
     a = abs(x)
-    polynomials_close_to_origin = sch.polynomials_close_to_origin
     if (polynomials_close_to_origin !== ()) && (a ≤ polynomials_close_to_origin[end][1])
         x * evalpoly(x * x, choose_poly(a, polynomials_close_to_origin))
     elseif isinf(x)

From c1e72b0ef56d4ca7801b68dd863edeb2b0a1540d Mon Sep 17 00:00:00 2001
From: Neven Sajko <s@purelymail.com>
Date: Sat, 16 Aug 2025 17:53:22 +0200
Subject: [PATCH 33/33] merge two methods into one

---
 base/special/trig.jl | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/base/special/trig.jl b/base/special/trig.jl
index 6858446837b1d..a972c762975bc 100644
--- a/base/special/trig.jl
+++ b/base/special/trig.jl
@@ -1233,11 +1233,13 @@ const _cosc_f64 = let b = Float64 ∘ Float16
     )
 end
 
-function _cosc(x::Float32)
-    _cos_cardinal_eval(x, _cosc_f32)
-end
-function _cosc(x::Float64)
-    _cos_cardinal_eval(x, _cosc_f64)
+function _cosc(x::Union{Float32, Float64})
+    if x isa Float32
+        pols = _cosc_f32
+    elseif x isa Float64
+        pols = _cosc_f64
+    end
+    _cos_cardinal_eval(x, pols)
 end
 
 # hard-code Float64/Float32 Taylor series, with coefficients