Support complex numbers in eps (#21858)

This came up in
JuliaLinearAlgebra/IterativeSolvers.jl#113 (comment)
. JuliaDiffEq and IterativeSolvers.jl have to make sure that the
real-type is pulled out in order for `eps` to work:

```julia
eps(real(typeof(b)))
```

This detail can make many algorithms with tolerances that are written
generically that would otherwise work with complex numbers error. This
PR proposes to do just that trick, so that way `eps(1.0 + 1.0im)`
returns machine epsilon for a Float64 (and generally works for
`AbstractFloat` of course).

---------

Co-authored-by: Steven G. Johnson <[email protected]>

Author: ChrisRackauckas

base/complex.jl

@@ -1141,3 +1141,13 @@ function complex(A::AbstractArray{T}) where T
     end
     convert(AbstractArray{typeof(complex(zero(T)))}, A)
 end
+
+## Promotion to complex ##
+
+_default_type(T::Type{Complex}) = Complex{Int}
+
+## Machine epsilon for complex ##
+
+eps(z::Complex{<:AbstractFloat}) = hypot(eps(real(z)), eps(imag(z)))
+
+eps(::Type{Complex{T}}) where {T<:AbstractFloat} = sqrt(2*one(T))*eps(T)

base/float.jl

@@ -1064,6 +1064,13 @@ julia> 1.0 + eps()
 julia> 1.0 + eps()/2
 1.0
 ```
+
+More generally, for any floating-point numeric type, `eps` corresponds to an
+upper bound on the distance to the nearest floating-point complex value: if ``\text{fl}(x)`` is the closest
+floating-point value to a number ``x`` (e.g. an arbitrary real number), then ``\text{fl}(x)``
+satisfies ``|x - \text{fl}(x)| ≤ \text{eps}(x)/2``, not including overflow cases.
+This allows the definition of `eps` to be extended to complex numbers,
+for which ``\text{fl}(a + ib) = \text{fl}(a) + i \text{fl}(b)``.
 """
 eps(::Type{<:AbstractFloat})
 

test/complex.jl

@@ -932,6 +932,13 @@ end
     end
 end
 
+@testset "eps" begin
+    @test eps(1.0+1.0im) === 3.1401849173675503e-16
+    @test eps(Complex128) === eps(1.0+1.0im)
+    @test eps(Complex64) === 1.6858739f-7
+    @test eps(Float32(1.0)+Float32(1.0)im) === eps(Complex{Float32})
+end
+
 @testset "cis" begin
     @test cis(0.0+1.0im) ≈ 0.367879441171442321595523770161460867445811131031767834507836+0.0im
     @test cis(1.0+0.0im) ≈ 0.54030230586813971740093660744297660373231042061+0.84147098480789650665250232163029899962256306079im