Fuzzing test failure with TwicePrecision in test/ranges.jl

[mbauman]

This fuzzing loop seems to sometimes find failures.

Found in the wild in https://julia.iblis.cnmc.tw/#/builders/1/builds/1617/steps/6/logs/stdio, reproduced locally by increasing the number of iterations.

Here are the magic numbers:

x, y = Base.TwicePrecision{Float32}(0.5767872f0, 1.9906205f-8), Base.TwicePrecision{Float32}(0.2913009f0, 1.4812992f-8)
xw, yw = asww(x), asww(y)
T = Float32
Tw = widen(T)
slopbits = 5

julia> cmp_sn2(Tw(xw/yw), astuple(x/y)..., slopbits)
false

[timholy]

In this method everything is exact until the computation of lo, which is off by a surprising margin:

julia> bits(widen(lo))
"0011111000101010110111101101001111100000000000000000000000000000"

julia> bits(((((widen(x.hi) - uh) - ul) + x.lo) - hi*y.lo)/y.hi)
"0011111000101010110111101101011101000110001010101000110110100101"

Unfortunately, it's going to be a while before I have time to figure this out (debugging this high-precision arithmetic is surprisingly time consuming). If people are desperate, just bump slopbits by 1 and it should pass.

[timholy]

FYI the 5 comes from 53 - 2*24, where 53 = number of precision bits of Float64 and 24 = number of precision bits of Float32. So 5 is taking the precision very, very seriously. I thought one can get away with doing that, but it's possible that the answer is no. It's also possible that some intermediate is generating a subnormal and that's the real reason.

[simonbyrne]

I had a quick look, the problem is the line:

lo = ((((x.hi - uh) - ul) + x.lo) - hi*y.lo)/y.hi

In this case x.hi - uh == 0. The sum of the remaining part has cancellation issues which loses some precision.

[simonbyrne]

So I found this:
https://hal.archives-ouvertes.fr/hal-01351529v2/document

I'm putting together a PR, should have something tomorrow.