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

more efficient muladd for complex/real operations #15980

Merged
merged 2 commits into from
Apr 22, 2016
Merged

more efficient muladd for complex/real operations #15980

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6445767
fix #15969: specialized muladd function for mixed complex/real operat…
stevengj Apr 21, 2016
c833e8e
more efficient use of muladd in evalpoly macro
stevengj Apr 21, 2016
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
13 changes: 13 additions & 0 deletions base/complex.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,10 @@ inv{T<:Integer}(z::Complex{T}) = inv(float(z))
*(z::Complex, w::Complex) = Complex(real(z) * real(w) - imag(z) * imag(w),
real(z) * imag(w) + imag(z) * real(w))

muladd(z::Complex, w::Complex, x::Complex) =
Complex(muladd(real(z), real(w), real(x)) - imag(z)*imag(w), # TODO: use mulsub given #15985
muladd(real(z), imag(w), muladd(imag(z), real(w), imag(x))))

# handle Bool and Complex{Bool}
# avoid type signature ambiguity warnings
+(x::Bool, z::Complex{Bool}) = Complex(x + real(z), imag(z))
Expand Down Expand Up @@ -163,6 +167,15 @@ end
*(x::Real, z::Complex) = Complex(x * real(z), x * imag(z))
*(z::Complex, x::Real) = Complex(x * real(z), x * imag(z))

muladd(x::Real, z::Complex, y::Number) = muladd(z, x, y)
muladd(z::Complex, x::Real, y::Real) = Complex(muladd(real(z),x,y), imag(z)*x)
muladd(z::Complex, x::Real, w::Complex) =
Complex(muladd(real(z),x,real(w)), muladd(imag(z),x,imag(w)))
muladd(x::Real, y::Real, z::Complex) = Complex(muladd(x,y,real(z)), imag(z))
Copy link
Contributor

Choose a reason for hiding this comment

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

This leaves out the case (Complex, Complex, Real). Is this intentional? If so, a comment might make sense.

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, it was intentional. All of the benefit from the muladd(Complex, Complex, Real) case could be gained by using muladd in Complex*Complex

Copy link
Contributor

Choose a reason for hiding this comment

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

Wouldn't this save one addition?

Copy link
Member Author

Choose a reason for hiding this comment

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

Added a comment.

Copy link
Member Author

@stevengj stevengj Apr 21, 2016 •
edited
Loading

Choose a reason for hiding this comment

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

@eschnett, no, because there is already an addition in both the real and imaginary parts of Complex*Complex that could be fused with one of the multiplications.

Oh wait, you're right; you could save an extra addition.

Copy link
Member Author

Choose a reason for hiding this comment

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

Okay, added the Complex*Complex cases

muladd(z::Complex, w::Complex, x::Real) =
Complex(muladd(real(z), real(w), x) - imag(z)*imag(w), # TODO: use mulsub given #15985
muladd(real(z), imag(w), imag(z) * real(w)))

/(a::Real, z::Complex) = a*inv(z)
/(z::Complex, x::Real) = Complex(real(z)/x, imag(z)/x)

Expand Down
4 changes: 2 additions & 2 deletions base/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,15 +74,15 @@ macro evalpoly(z, p...)
ai = symbol("a", i)
push!(as, :($ai = $a))
a = :(muladd(r, $ai, $b))
b = :(muladd(-s, $ai, $(esc(p[i]))))
b = :($(esc(p[i])) - s * $ai) # see issue #15985 on fused mul-subtract
end
ai = :a0
push!(as, :($ai = $a))
C = Expr(:block,
:(x = real(tt)),
:(y = imag(tt)),
:(r = x + x),
:(s = x*x + y*y),
:(s = muladd(x, x, y*y)),
Copy link
Contributor

Choose a reason for hiding this comment

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

This breaks the symmetry between x and y. I assume @evalpoly is only used in circumstances where something like @fastmath would make sense?

Copy link
Member Author

@stevengj stevengj Apr 21, 2016 •
edited
Loading

Choose a reason for hiding this comment

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

Yes. (Note that we are already rearranging the computation quite a bit compared to, say, Horner's rule. And r = x+x already breaks the real/imaginary symmetry.)

as...,
:(muladd($ai, tt, $b)))
R = Expr(:macrocall, symbol("@horner"), :tt, map(esc, p)...)
Expand Down
5 changes: 5 additions & 0 deletions test/complex.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -904,6 +904,11 @@ end
# issue #10926
@test typeof(π - 1im) == Complex{Float64}

# issue #15969: specialized muladd for complex types
for x in (3, 3+13im), y in (2, 2+7im), z in (5, 5+11im)
@test muladd(x,y,z) == x*y + z
Copy link
Contributor

Choose a reason for hiding this comment

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

should probably be ===

Copy link
Member Author

Choose a reason for hiding this comment

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

They are Complex{Int} values, so what does it matter?

Copy link
Contributor

Choose a reason for hiding this comment

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

to test that the return type stays that way

Copy link
Member Author

Choose a reason for hiding this comment

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

Ah, I see.

Copy link
Member Author

Choose a reason for hiding this comment

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

Fixed in de6468b

end

# issue #11839: type stability for Complex{Int64}
let x = 1+im
@inferred sin(x)
Expand Down