Change the implementation of rotate! to avoid unary minus

[lanceXwq]

The current implementation of rotate!(x, y, c, s) overwrites x with c*x + s*y and y with -conj(s)*x + c*y. There exists this edge case of rotate!([NaN], [NaN], false, false)

julia> rotate!([NaN], [NaN], false, false)
([0.0], [NaN])

The new x becomes 0.0, because false acts as a "strong zero" in Julia (so false * NaN == 0.0). When calculating the new y, the code executes -conj(false)*NaN + false*NaN. However, the unary minus, -, turns the false::Bool to 0::Int, which is no longer a strong zero, and thus, we get NaN. This is inconsistent with

julia> reflect!([NaN], [NaN], false, false)
([0.0], [0.0])

This PR changes -conj(s)*x + c*y to c*y - conj(s)*x to avoid unary minus, so that

julia> rotate!([NaN], [NaN], false, false) # This PR
([0.0], [0.0])

I noted this from JuliaGPU/CUDA.jl#2607 (comment). More discussion in https://discourse.julialang.org/t/negative-false-no-longer-acts-as-a-strong-zero/128506/1

[jishnub]

This looks fine to me, but is this case really encountered in practice where s and c are both zero?

[lanceXwq]

Probably not something people would use for real computations. 😅 But I imagine there will occasionally be folks exploring these edge cases just to see how the "strong zeroness" of false is handled in a package like LinearAlgebra.jl, similar to the example I shared earlier from CUDA.jl. From that perspective, I think maintaining some level of consistency could be helpful for users deciding how to structure their own code. 😃

[mcabbott]

Some chance that more of these functions ought to branch on zero.

I don't see this documented, but mul!(fill(NaN,2,2), ones(2,2), ones(2,2), 3.0, 0.0) doesn't need false to understand that you mean to overwrite the first argument. Should rotate! follow that convention?

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.73%. Comparing base (c33045c) to head (67a065c).
Report is 1 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1323   +/-   ##
=======================================
  Coverage   93.73%   93.73%           
=======================================
  Files          34       34           
  Lines       15746    15746           
=======================================
  Hits        14760    14760           
  Misses        986      986           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[lanceXwq]

From what I can tell, mul!(fill(NaN,2,2), ones(2,2), ones(2,2), 3.0, 0.0) simply overwrites the NaNs without involving them in the computation. A fairer comparison might be

julia> mul!(ones(2,2), ones(2,2), fill(NaN,2,2), 3.0, false)
2×2 Matrix{Float64}:
 NaN  NaN
 NaN  NaN

which shows that "strong zeroness" doesn't propagate in this case. But then we also have

julia> rmul!(fill(NaN, 2,2), false)
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

julia> lmul!(false, fill(NaN, 2,2))
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

So yeah, behavior seems to vary a bit depending on the function. It might be too much to expect full consistency across all of LinearAlgebra.jl, but IMHO it would be nice if closely related functions, like rotate! and reflect!, rmul! and lmul!, were at least consistent with each other.

[mcabbott]

That mul! example is like false * 1.0 + 3.0 * NaN , so I think this issue of 0 being strong doesn't arise.

I agree rmul! doesn't follow mul! here. It could, but unlike mul! this is a slightly odd use... rmul!(A, 0) is better written fill!(A, 0). The point of the function is things like rmul!(fill(NaN,2,2), UpperTriangular(ones(2,2))), which is mul!(rand(2,2), fill(NaN,2,2), UpperTriangular(ones(2,2))) i.e. the NaN is unambiguously part of the calculation here.