Fix NaN arithmetic

[PatrickHaecker]

Builds on the infinity_direction branch, so we can rebase the current PR on master once #80 has landed to reduce the overall diff size (we both have the permissions). I recommend to review again commit-by-commit. The most recent 4 commits are new.

1. Base promotes every Real to BigFloat, so Base._promote converted the infinity away and big(2.0) + ∞, big(2.0) * ∞ and (+∞)^big(2.0) were all MethodErrors, which also made the NaN work below untestable at that precision.
2. NaN + ∞ gave ∞, NaN * ∞ gave +∞, div(NaN, ∞) gave 0.0 and (+∞)^NaN gave 0.0, so a NaN marking a failed computation silently became a plausible infinity; every entry point now returns its NaN argument unchanged, precision included.
3. signbit(ComplexInfinity(0.5)) returned the field itself, 0.5, so any caller branching on it hit a TypeError. The three methods are replaced by one correct method.
4. Add tests for such cases to avoid these kind of errors in the future.

This PR changes inference, but it was necessary to have a correct result:

• x::Float64 + ∞ now infers Union{Float64, Infinity} rather than Infinity (but no allocations)
• x::Int + ∞ is unchanged because isnan folds away for types that have no NaN.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 100.00%. Comparing base (13f9e1f) to head (cfba27d).

Additional details and impacted files

@@            Coverage Diff            @@
##            master       #81   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files            6         6           
  Lines          266       272    +6     
=========================================
+ Hits           266       272    +6     

☔ View full report in Codecov by Harness.
📢 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.

[dlfivefifty]

Hahaha this is exactly what I just suggested!

[dlfivefifty]

Btw For every θ that is "natural" in base 2 (1, -1, 1/2, 1/4, 3/4, etc.) it turns out angle(exp(im*θ)) == θ. This seems like a a very lucky coincidence that rounding goes in the right direction!

[PatrickHaecker]

Hahaha this is exactly what I just suggested!

Yes, but it only covers part of it, because we still have no defined interval the value is in and we still have a very strange format so store an angle when we have an integer. So this needs more thought.

[PatrickHaecker]

Btw For every θ that is "natural" in base 2 (1, -1, 1/2, 1/4, 3/4, etc.) it turns out angle(exp(im*θ)) == θ. This seems like a a very lucky coincidence that rounding goes in the right direction!

You are definitely onto something here. While it is true for all your numbers, it interestingly is not true, e.g., for 1/8:

julia> angle(exp(im*π/8))/π
0.12500000000000003

So although 1/8 is exact in floating points, the calculation is not exact. That not only means that our tests are incomplete here (I already added some, but have not yet pushed). It also means that we need to check whether this loss of precision is a problem anywhere. And there are quite some locations where we do the exact comparison.

But I need to do this with a fresh mind tomorrow.

Especially as this lead me to toinf and I currently think we have a bug in toinf:

julia> angle(toinf(complex(0,Inf)))
4.934802200544679

julia> angle(complex(0,Inf))
1.5707963267948966

And I don't understand how we can have such a large bug here without the test suite detecting it. So I probably need to have a fresh look tomorrow to see why this is not a bug and only my current understanding is incomplete.

[PatrickHaecker]

It was indeed a bug and there is an explanation why the tests haven't caught it. Fix, details and off-axis tests in #85. So for all I know this conversation is resolved.

[dlfivefifty]

The really important case is that for r > 0 where r is any real we always have _sb(r+0im) = zero(r), _sb(-r±0im) = one(r), _sb(r*im) = one(r)/2, _sb(-r*im) = -one(r)/2

so we want this to be correct for BigFloat and other types too. Actually, I'm not sure what the types should be (I guess for <: Real we want them to be Bool)

[PatrickHaecker]

I assume you mean == as comparison. What do you mean with "types"? The storage type or the return type?

I am currently trying a combination of UInt64 as a storage type and using sign with a Float64 as a possible user interface in addition to other methods for Base functions.
This seems to work well up to now (unique values, some precision everywhere, wrapping for free, usage without rounding possible), but I can't yet see whether this will work in the whole package, which is the relevant question in the end.

[dlfivefifty]

This is a bit of a fringe question so I don't think any choice makes a big practical choice. But what I mean is if _sb(x::Real) = signbit(x) then the return type is a Bool, but for _sb(::Complex{<:Union{Int,Float64}) it should probably be a Float64. But then what about _sb(::Complex32) and _sb(::Complex{BigFloat})? Probably Float32 and BigFloat. (which I believe your implementation does)

[dlfivefifty]

I guess the concern is: do we actually know that rounding works correctly for all BigFloat? Eg, we have

julia> setprecision(10) do
       angle(one(BigFloat)im)
       end
1.5703

julia> setprecision(20) do
       angle(one(BigFloat)im)
       end
1.570797

julia> setprecision(40) do
       angle(one(BigFloat)im)
       end
1.5707963267942

There is no obvious guarantee that when we divide this by π we get exactly 1/2.

Except.....in tests we do up to precision 10_000:

julia> for p in 1:10_000
       setprecision(p) do
       @test angle(one(BigFloat)im)/π == 1/2
       end
       end

Which is to say: I think by luck we are fine.

[PatrickHaecker]

I assume you talking about #81 and not about the possible future solution.

I think in general we can't be exact with floats. The user can always have any arbitrary transformations between two points in time and when they compare then according to what should be identical with two mathematical reals it can be different in our implementation.

That being said, I think in this case it is probably more than luck. I am speculating, because I haven't followed into MPFR's implementation, but your example boils down to:

atan(one(BigFloat), zero(BigFloat))/π == 1/2

which is

atan(one(BigFloat), zero(BigFloat)) == π/2

and I strongly assume that this special cases such that atan will return the best rounding of π/2 according to the current precision if the first argument is positive and the second is 0. This triggers reliably because 1 and 0 are exact. π and π/2 will then have the same number in the significant and only differ in the mantissa by exactly 1 leading to 0.5 exactly in the division.

So while it would be possible to implement that in a different way, I think an implementation which cares about precision will always return the exact value here, despite being a floating point operation. That's just a special case which is exact when using binary floating points. It wouldn't for some other base in general, but it is with base 2.

[dlfivefifty]

Once the conflict is fixed I can merge this

[PatrickHaecker]

Once the conflict is fixed I can merge this

Let me please have a look at your rounding finding first. Because this might indeed not be the best way to do it. I am not sure whether this is guaranteed in all cases and whether we are missing a "not on one of the complex axes" case.

[PatrickHaecker]

With the fix and the missing tests which were identified during review being handled in #85, the current PR (#81) is ready again.