Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ComplexF64 division: combine four if-statements into two if-elseif-statements #29042

Merged
merged 4 commits into from
Sep 18, 2018
Merged

ComplexF64 division: combine four if-statements into two if-elseif-statements #29042

Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8913022
ComplexF64 division: combine four if-statements into two if-elseif-st…
thchr Sep 4, 2018
bf5df6e
add @fastmath to magnitude check; NaN handling not needed
thchr Sep 4, 2018
32a0eb9
don't use @fastmath; use explicit route instead
thchr Sep 5, 2018
283528f
remove unnecessary two & half variables
thchr Sep 5, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 23 additions & 11 deletions base/complex.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -360,19 +360,31 @@ inv(z::Complex{<:Union{Float16,Float32}}) =
# c + i*d
function /(z::ComplexF64, w::ComplexF64)
a, b = reim(z); c, d = reim(w)
half = 0.5
two = 2.0
ab = max(abs(a), abs(b))
cd = max(abs(c), abs(d))
ov = floatmax(a)
un = floatmin(a)
ϵ = eps(Float64)
@fastmath ab = max(abs(a), abs(b))
Copy link
Contributor

Choose a reason for hiding this comment

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

Note @fastmath isn't really valid here: it implies that we can assume a or b are never Inf, NaN or subnormal.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ah, you're quick. I thought I could already put my lessons above to use.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Maybe I'm misunderstanding the nature of the @fastmath macro here: I thought it simply swapped the max function for max_fast (and same for abs, though abs and abs_fast appear identical)?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Regardless, I'll swap it to your initial suggestion instead, thanks.

Copy link
Contributor

Choose a reason for hiding this comment

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

@fastmath is probably okay in this case, but occasionally it can cause some problems when the compiler gets carried away. Probably better to be explicit.

@fastmath cd = max(abs(c), abs(d))

# constants
ov = floatmax(Float64)
un = floatmin(Float64)
ϵ = eps(Float64)
half = 0.5; halfov = half*ov
two = 2.0; twounϵ = un*two/ϵ
bs = two/(ϵ*ϵ)

# scaling
s = 1.0
ab >= half*ov && (a=half*a; b=half*b; s=two*s ) # scale down a,b
cd >= half*ov && (c=half*c; d=half*d; s=s*half) # scale down c,d
ab <= un*two/ϵ && (a=a*bs; b=b*bs; s=s/bs ) # scale up a,b
cd <= un*two/ϵ && (c=c*bs; d=d*bs; s=s*bs ) # scale up c,d
if ab >= halfov
a*=half; b*=half; s*=two # scale down a,b
elseif ab <= twounϵ
a*=bs; b*=bs; s/=bs # scale up a,b
end
if cd >= halfov
c*=half; d*=half; s*=half # scale down c,d
elseif cd <= twounϵ
c*=bs; d*=bs; s*=bs # scale up c,d
end

# division operations
abs(d)<=abs(c) ? ((p,q)=robust_cdiv1(a,b,c,d) ) : ((p,q)=robust_cdiv1(b,a,d,c); q=-q)
return ComplexF64(p*s,q*s) # undo scaling
end
Expand Down