Skip to content
This repository has been archived by the owner on Dec 10, 2023. It is now read-only.

Dekker multiplication isn't exact #30

Open
timholy opened this issue Aug 14, 2017 · 3 comments · May be fixed by #31
Open

Dekker multiplication isn't exact #30

timholy opened this issue Aug 14, 2017 · 3 comments · May be fixed by #31

Comments

@timholy
Copy link
Member

timholy commented Aug 14, 2017

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 timholy changed the title Dekker multiplication doesn't handle subnormals Dekker multiplication isn't exact Aug 14, 2017
@timholy
Copy link
Member Author

timholy commented Aug 14, 2017

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
Copy link
Member Author

timholy commented Aug 14, 2017

I think here's the origin of the trouble (referencing

DoubleDouble.jl/src/DoubleDouble.jl

Lines 141 to 145 in 67683c7

function *{T}(x::Single{T},y::Single{T})
hx,lx = splitprec(x.hi)
hy,ly = splitprec(y.hi)
z = x.hi*y.hi
Double(z, ((hx*hy-z) + hx*ly + lx*hy) + lx*ly)
):

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
Copy link
Member

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...

simonbyrne added a commit that referenced this issue Aug 14, 2017
@simonbyrne
Use correct calculation for Valtkamp splitting
6a52562 
Should fix #30.
@simonbyrne simonbyrne linked a pull request Aug 14, 2017 that will close this issue
Open
simonbyrne added a commit that referenced this issue Aug 14, 2017
@simonbyrne
Use correct calculation for Veltkamp splitting
d8f392a 
Should fix #30.
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

Use correct calculation for Veltkamp splitting JuliaMath/DoubleDouble.jl
2 participants
@timholy @simonbyrne