faster inv and div for Complex{Union{Float16, Float32}}

[oscardssmith]

As discovered in https://discourse.julialang.org/t/efficient-calculation-of-many-complex-lorentzians-on-a-large-array/76150/21, Complex division for Float16 and Float32 are currently doing a lot of extra work. This still does the math in Float64, but ignores all the extra checks.

[oscardssmith]

Benchmarks look pretty good to me

#generates random numbers with widely varying values
myrand() = ComplexF32(reinterpret(Float32, rand(Int32)), reinterpret(Float32, rand(Int32)))
@btime myrand()
  3.528 ns (0 allocations: 0 bytes)

#old
@btime inv(myrand())  
  10.503 ns (0 allocations: 0 bytes)

@btime myrand()/myrand()
  17.030 ns (0 allocations: 0 bytes)
  
#new
@btime inv(myrand())  
  3.983 ns (0 allocations: 0 bytes)

@btime myrand()/myrand()
  6.893 ns (0 allocations: 0 bytes)

[LaurentPlagne]

Maybe I have missed something but isn't it a problem that CFloat32/16 semantic differs from CFloat64 ? (not the same checks).

[oscardssmith]

This isn't quite right yet, but a lot of the checks have to do with numeric issues that don't exist when you do the math in wider precision.

[giordano]

Can we please have correct functions rather than fast ones?

julia> function Base.inv(z::Complex{<:Union{Float16,Float32}})
           c, d = reim(widen(z))
           (isinf(c) | isinf(d)) && return complex(copysign(zero(z), c), flipsign(-zero(z), d))
           mag = inv(muladd(c, c, d^2))
           return oftype(z, Complex(c*mag, -d*mag))
       end

julia> inv(complex(Inf32))
ERROR: MethodError: no method matching copysign(::ComplexF32, ::Float64)
Closest candidates are:
  copysign(::Signed, ::Float64) at /usr/share/julia/base/int.jl:149
  copysign(::Signed, ::Real) at /usr/share/julia/base/int.jl:150
  copysign(::Rational, ::Real) at /usr/share/julia/base/rational.jl:256
  ...
Stacktrace:
 [1] inv(z::ComplexF32)
   @ Main ./REPL[1]:3
 [2] top-level scope
   @ REPL[3]:1

[oscardssmith]

spelling is hard.

[giordano]

Relatedly, division of complex numbers has been plagued by several issues, some of them have been fixed over time, but for example #32992 still stands. Might be a good occasion to fix it?

[oscardssmith]

Can you write me a test script for what the behavior should be? I think I want to merge this without fixing the Float64 edge cases (1.8 feature freeze is coming up), but I'd be happy to fix these in a follow-up PR.

[giordano]

@testset "Complex division by infinity - $T" for T in (Float16, Float32, Float64)
       @test isequal(complex(one(T)) / complex(T(Inf), T(Inf)), complex(zero(T)))
end

Note that with latest version of this PR inv of infinity is still broken:

julia> function Base.:/(z::Complex{T}, w::Complex{T}) where {T<:Union{Float16,Float32}}
           c, d = reim(widen(w))
           a, b = reim(widen(z))
           (isinf(c) | isinf(d)) && return complex(copysign(zero(T), c), flipsign(zero(T), d)) * z
           mag = inv(muladd(c, c, d^2))
           re_part = muladd(a, c, b*d)
           im_part = muladd(b, c, -a*d)
           return oftype(z, Complex(re_part*mag, im_part*mag))
       end

julia> function Base.inv(z::T) where T<:Complex{<:Union{Float16,Float32}}
           c, d = reim(widen(z))
           (isinf(c) | isinf(d)) && return complex(copysign(zero(T), c), flipsign(-zero(T), d))
           mag = inv(muladd(c, c, d^2))
           return oftype(z, Complex(c*mag, -d*mag))
       end

julia> inv(complex(Inf32))
ERROR: MethodError: no method matching copysign(::ComplexF32, ::Float64)
Closest candidates are:
  copysign(::Signed, ::Float64) at /usr/share/julia/base/int.jl:149
  copysign(::Signed, ::Real) at /usr/share/julia/base/int.jl:150
  copysign(::Rational, ::Real) at /usr/share/julia/base/rational.jl:256
  ...
Stacktrace:
 [1] inv(z::ComplexF32)
   @ Main ./REPL[5]:3
 [2] top-level scope
   @ REPL[6]:1

Same error as before

[oscardssmith]

@giordano how would you feel about changing the sign of the zeros for a/b to match those of a*inv(b). Currently, big(1.0 + 0.0im)/big(-Inf + Inf*im) != inv(big(-Inf + Inf*im)) which seems wrong.

[giordano]

What you suggest makes sense, but in v1.7.2 I get:

julia> big(1.0 + 0.0im)/big(-Inf + Inf*im) == inv(big(-Inf + Inf*im))
true

so that's already the case? If in master you get something different that looks like a regression

[oscardssmith]

huh, I clearly need to experiment more...

[giordano]

Why the sign change?

[oscardssmith]

Because I was having trouble with math. Complex numbers in IEEE don't have a multiplicative identity. Complex(1.,0.)*Complex(-0.,-0.)==Complex(0.,-0.) and Complex(1.,-0.)*Complex(-0.,-0.)==Complex(-0.,0.), so my idea for getting the signs wasn't working...

[oscardssmith]

do you know a good algorithm for computing the sign of the zeros when we know that the divisor is some type of infinity? I've looked for a bit, but can't come up with anything nice.

[giordano]

I think it should be the following: in (a + b * im) / (c + d * im), assuming a and b finite and non-zero, and c and d infinite:

• real part: sign(a) * sign(c)
• imaginary part: sign(b) * -sign(d)

[giordano]

For completeness, you should likely refer to the ISO/IEC 10967-3:2006 standard, section "5.2
Imaginary and complex floating point datatypes and operations"

[KristofferC]

It's a bit worrying to me how e.g. signs in tests are being changed with an off-hand "fix tests" commit. It's to me a bit of a Chesterton's fence situation. A significant amount of work has gone in to making sure that these things are working the way they are and one needs to be quite careful when changing these pretty low-level operations.

[oscardssmith]

@giordano thanks for finding the thing that actually works.

[oscardssmith]

Is this ready to merge (assuming CI finally works?)

[giordano]

Failures like

  return type ComplexF16 does not match inferred return type Union{ComplexF16, ComplexF32}

and

  Test threw exception
  Expression: (D \ S)::Tridiagonal{elty} == Tridiagonal(Matrix(D) \ Matrix(S))
  TypeError: in typeassert, expected LinearAlgebra.Tridiagonal{ComplexF32, V} where V<:AbstractVector{ComplexF32}, got a value of type LinearAlgebra.Tridiagonal{Union{ComplexF32, ComplexF64}, Vector{Union{ComplexF32, ComplexF64}}}

look relevant

[giordano]

Ok, performance improvements look good and we're reconciling with math for division by infinity across all types. Failure on Linux looks unrelated as far as I could tell, but FreeBSD and macOS tests are stalling, I don't know what's their status nowadays.

[oscardssmith]

Ok. I think this is finally ready.