Dekker multiplication isn't exact

[timholy]

I was surprised to discover this:

x = Float16(0.992)
y = Float16(6.0e-8)   # subnormal

julia> d = Single(x)*Single(y)
Double(6.0e-8, 0.0)
 - value: 5.9604644775390625e-08

julia> bits(widen(d.hi) + widen(d.lo))
"00110011100000000000000000000000"

julia> bits(widen(x)*widen(y))
"00110011011111100000000000000000"

It can be fixed this way:

julia> ys, ye = frexp(y)
(Float16(0.5), -23)

julia> d = Single(x)*Single(ys)
Double(0.496, 0.0)
 - value: 0.49609375

julia> bits(ldexp(widen(d.hi) + widen(d.lo), ye))
"00110011011111100000000000000000"

But given that many treatises say "it's exact unless the splitting overflows," this was a surprising discovery to me.

[timholy]

Here's another interesting case, one that doesn't involve subnormals:

julia> x, y = Single(Float16(0.9775)), Single(Float16(0.5156))
(DoubleDouble.Single{Float16}(Float16(0.9775)), DoubleDouble.Single{Float16}(Float16(0.5156)))

julia> d = x*y
Double(0.5044, -0.0001068)
 - value: 0.5042877197265625

julia> bits(widen(d.hi) + widen(d.lo))
"00111111000000010001100100000000"

julia> bits(widen(x.hi) * widen(y.hi))
"00111111000000010000100100000000"

I thought this should satisfy the theorem that it should be exact?

[timholy]

I think here's the origin of the trouble (referencing https://github.com/JuliaMath/DoubleDouble.jl/blob/67683c79604b266781afb27ea357a78df61a1d0a/src/DoubleDouble.jl#L141-L145):

julia> bits(widen(hx)*widen(hy))
"00111111000000011111000000000000"

julia> bits(widen(hx*hy))
"00111111000000100000000000000000"

Doesn't the algorithm assume that multiplication of the two hi terms is exact?

[simonbyrne]

Ah, the problem is that half16 should be 2^6+1 == 65 not 2^5+1 == 33 (see §4.4.1 of the Handbook of Floating Point Arithmetic). I remember the Dekker paper being a bit confusing as to whether you should round 11/2 up or down.

In any case, we should move to using fma instead of Veltkamp splitting when it is available...