Change in ComplexF64 division

[KristofferC]

#47255 caused the following:

julia> inv(0.2 + 0.0im)
4.999999999999999 - 0.0im

In 1.8 this was:

julia> inv(0.2 + 0.0im)
5.0 - 0.0im

I don't know what the guarantees are for this operation but it seems kind of silly to stop being able to do this.

[gbaraldi]

Yeah it seems odd. I would imagine we guarantee that a complex number with no imaginary part would behave the same as a normal Float and it here it doesn't.

[oscardssmith]

We could add a branch for zero imaginary part, but IMO it's not worth it (it would be a noticeable slowdown)

[KristofferC]

We could add a branch for zero imaginary part, but IMO it's not worth it (it would be a noticeable slowdown)

I am not sure about that:

julia> using GFlops

julia> @count_ops inv(1.0+2.0im) 
Flop Counter: 13 flop
┌────────┬─────────┐
│        │ Float64 │
├────────┼─────────┤
│ muladd │       1 │
│    add │       1 │
│    mul │       1 │
│    div │       4 │
│    abs │       2 │
│    neg │       1 │
│   sqrt │       2 │
└────────┴─────────┘

Feels like a well predicted branch should be not very noticeable among that?

(@count_ops probably has some extra ops there that are not visible with const prop but still)

[oscardssmith]

The actual count after const-prop is much lower since none of the sqrts and half of the divisions aren't real.

[KristofferC]

Tried first benchmark in #47255 (comment) and it went from 1.933 μs to 2.037 μs. That feels worth it to me.

[oscardssmith]

fair enough. One thing to be wary of is that that benchmark will look better for this change since with that benchmark the cpu will learn that the imaginary part is never 0, and if you care about this you are probably in the case where the branch is harder to predict.

[KristofferC]

But if you have a lot of zero imaginary numbers, shouldn't the new implementation be faster since it only has to do a single fp division?

[mikmoore]

I would imagine we guarantee that a complex number with no imaginary part would behave the same as a normal Float and it here it doesn't.

For a long time (and at least up to v1.8.0) we've also had this extremely annoying result:

julia> complex(49.0)/complex(49.0)
0.9999999999999999 + 0.0im

I've been working on a new division algorithm and am hoping to open a PR for it in the next month.

[oscardssmith]

@KristofferC it will depend a bit on the distribution. If you are more than 50% real then it will probably be faster, but if you are between 5% and 50% (roughly) it will be slower since it will mess up the branch prediction
ooh nice! I hate our current division algorithm.