special case complex division for zero imaginary part

[KristofferC]

Fixes #48457

[mikmoore]

What about inv(complex(0.2,1e-320))? I assume this will still show the same issue on the real part.

But fixing that is likely more involved and maybe should be a different PR. Discrepancies like that are at least slightly more defensible than with literal zeros.

[oscardssmith]

The slightly more general case is

if absd < sqrt(eps(Float64))*absc
    return inv(c)*complex(1.0, -d/c)
end

This works because in this case c^2+d^2 == c^2 so we can simplify. This adds roughly 1 cycle to the condition cost but I could see that as worth it.

[mikmoore]

Sorry for not taking a deeper look at this the first time. This doesn't seem like the right place for a fix. This change seems worse than doing nothing, since inv(complex(0.2,0.0)) won't agree with inv(complex(0.2,0.0)*im)*im.

It appears that the fastpath conj(w) / muladd(cd, cd, dc*dc) computation is to blame here. For complex(c,0.0) this becomes c/c^2 != 1/c. It looks like the solution is simply to drop the fastpath. The slowpath robust_cinv will compute this correctly.

[KristofferC]

This change seems worse than doing nothing, since inv(complex(0.2,0.0)) won't agree with inv(complex(0.2,0.0)*im)*im

We can add a special case for zero real as well ;)

[oscardssmith]

IMO the fast path is worth it. It is occasionally slightly inaccurate, but it is 2x faster and is always pretty accurate.

[mikmoore]

Alright. I'm having trouble imagining other inputs where the fastpath result would be egregious. Just make sure that the extra checks aren't eroding all the benefit over a pure-slowpath version (don't need to post benchmarks, just know so).

[KristofferC]

Just to get a second data point I tried Rust (which is not really known for its numerics but still):

fn main() {
    println!("Complex float: {}", 1.0 / num::complex::Complex::new(0.2, 0.0));
    println!("Float: {}", 1.0 / 0.2);
}

// Printing:
Complex float: 4.999999999999999+0i
Float: 5

[oscardssmith]

for the record, this case is only off by 0.3 ULP (although there are probably cases where it's off by more, my guess is up to ~1.5 but I haven't done the math)

[mikmoore]

1.0 / num::complex::Complex::new(0.2, 0.0)

Does Rust have an inv operation? I'd check that, too. I checked MATLAB and it at least manages to give the 5.0 answer (even accounting for the dirty tricks they play with rounding printed numbers).

In Julia, I would accept 1.0/x and inv(x) giving these different results (or at least 2.0/x and 2.0*inv(x), were we to make a special case for 1.0). The reason is that I would want inv(complex(a,0)) == 1/a while I would want complex(a,b)/complex(a,b) == 1. Satisfying these both with a unified function requires branching to different computations.

Another thought: our automatic differentiation packages have explicit definitions of these computations, right? They don't just derive the result via these implementations? These sorts of special cases ruin derivatives. EDIT: the sufficiently-small-d version by oscardssmith shouldn't cause a problem, I think, only a pure-zero version.

[oscardssmith]

Good point. That would be another point in favor of the generalized version I gave.

[mikmoore]

I think the resulting derivatives would still be incorrect based on complex(inv(c),-d/c). Neither of those terms includes the negligible-value parts that may have non-negligible derivatives. The errors would be small in many cases, but definitely present and almost certainly enough to cause problems.

[oscardssmith]

I've done the math (i.e. typed it into wolfram alpha) and the derivatives work out correct to within epsilon also.

[KristofferC]

These sorts of special cases ruin derivatives.

I think the resulting derivatives would still be incorrect

These are imo bugs in the AD libraries. For example, ForwardDiff handles these type of problems on its master branch (by not considering Dual numbers with non-zero dual part to be equal to a real number).

[mikmoore]

I've done the math (i.e. typed it into wolfram alpha) and the derivatives work out correct to within epsilon also.

julia> using ForwardDiff

julia> ForwardDiff.jacobian(((c,d),)->[c,-d]/(c^2+d^2),[0.2,0.0])
2×2 Matrix{Float64}:
 -25.0    0.0
   0.0  -25.0

julia> ForwardDiff.jacobian(((c,d),)->[inv(c),-d/c],[0.2,0.0])
2×2 Matrix{Float64}:
 -25.0  -0.0
   0.0  -5.0

EDIT: The equation used here was wrong and is corrected in a later post. Also, all autodiff is irrelevant here because the input type is ComplexF64.
The latter version produces an incorrect result. Higher derivatives will be worse still cause problems.

If ForwardDiff etc refutes equality with nonzero derivatives that's a start, but I can't imagine it does the same for inequalities like we're proposing for "small" d. And if it does refute equality then it will change the results of computations.

Sorry, this might be an unnecessary diversion. I'll let the issue go unless someone wants to engage further.

[oscardssmith]

yeah that's because I typod. It was supposed to be inv(c)*Complex(1, -d/c)

[mikmoore]

Thanks. I should have caught that.

The point on higher derivatives stands

julia> ForwardDiff.hessian(((c,d),)->inv(c),[0.2,0.0])
2×2 Matrix{Float64}:
 250.0  0.0
   0.0  0.0

julia> ForwardDiff.hessian(((c,d),)->c/(c^2+d^2),[0.2,0.0])
2×2 Matrix{Float64}:
 250.0     0.0
   0.0  -250.0

Of course the nicer version of this (that doesn't fall into our aforementioned c/c^2 trap) is inv(c+d^2/c) as the real part, but then we're basically running the slow path.

EDIT: This method is only used for ComplexF64. My digression into autodiff is totally off-base (unless maybe this is subvertable with a Cassette.jl-based autodiff or something). Sorry for the noise.

[oscardssmith]

Is this ready to merge?

[KristofferC]

I fixed the test upstream so I vote for not changing anything until someone else complains about it.