Skip to content

Commit

Permalink
faster inv for normal sized ComplexF64 (#47255)
Browse files Browse the repository at this point in the history
* faster `inv` for normal sized `ComplexF64`
  • Loading branch information
oscardssmith authored Nov 9, 2022
1 parent 86bb1fc commit 83592cf
Show file tree
Hide file tree
Showing 2 changed files with 11 additions and 6 deletions.
15 changes: 10 additions & 5 deletions base/complex.jl
Original file line number Diff line number Diff line change
Expand Up @@ -472,9 +472,13 @@ function inv(z::Complex{T}) where T<:Union{Float16,Float32}
end
function inv(w::ComplexF64)
c, d = reim(w)
(isinf(c) | isinf(d)) && return complex(copysign(0.0, c), flipsign(-0.0, d))
absc, absd = abs(c), abs(d)
cd = ifelse(absc>absd, absc, absd) # cheap `max`: don't need sign- and nan-checks here
cd, dc = ifelse(absc>absd, (absc, absd), (absd, absc))
# no overflow from abs2
if sqrt(floatmin(Float64)/2) <= cd <= sqrt(floatmax(Float64)/2)
return conj(w) / muladd(cd, cd, dc*dc)
end
(isinf(c) | isinf(d)) && return complex(copysign(0.0, c), flipsign(-0.0, d))

ϵ = eps(Float64)
bs = 2/*ϵ)
Expand All @@ -493,12 +497,13 @@ function inv(w::ComplexF64)
else
q, p = robust_cinv(-d, -c)
end
return ComplexF64(p*s, q*s) # undo scaling
return ComplexF64(p*s, q*s)
end
function robust_cinv(c::Float64, d::Float64)
r = d/c
p = inv(muladd(d, r, c))
q = -r*p
z = muladd(d, r, c)
p = 1.0/z
q = -r/z
return p, q
end

Expand Down
2 changes: 1 addition & 1 deletion doc/src/manual/complex-and-rational-numbers.md
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ julia> 3(2 - 5im)^2
-63 - 60im
julia> 3(2 - 5im)^-1.0
0.20689655172413796 + 0.5172413793103449im
0.20689655172413793 + 0.5172413793103449im
```

The promotion mechanism ensures that combinations of operands of different types just work:
Expand Down

2 comments on commit 83592cf

@maleadt
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@nanosoldier runtests(ALL, isdaily = true, configuration=(rr=true,))

@nanosoldier
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Your package evaluation job has completed - possible new issues were detected. A full report can be found here.

Please sign in to comment.