bug: "clean" mode has worse accuracy than regular Float64

[nsajko]

Here's a case where MultiFloats "clean" mode is worse than "standard" mode. To make matters worse, "clean" mode is worse than Float64 in this case:

julia> using MultiFloats

julia> MultiFloats.use_clean_multifloat_arithmetic(20)

julia> const F = MultiFloat{Float64, n} where {n}
Float64x (alias for MultiFloat{Float64})

julia> f(x) = isone(x*inv(x))
f (generic function with 1 method)

julia> f(3.0^100)
true

julia> f(F{13}(3.0^100))
false

julia> MultiFloats.use_standard_multifloat_arithmetic(20)

julia> f(F{13}(3.0^100))
true

julia> MultiFloats.use_clean_multifloat_arithmetic(20)

julia> f(F{13}(3.0^100))
false

So there seems to be a bug somewhere in "clean" mode.

[dzhang314]

Hey @nsajko, thanks for discovering this surprising behavior! It seems that, in this case, the extra limb computed in clean mode is not quite cancelled out in the multiplication x * inv(x), producing a number that is very slightly different from 1.0.

This behavior is counterintuitive but actually correct. The number (1/3)^100 cannot be exactly represented as a terminating decimal in base 2, so you should not expect inv(3.0^100) to be the exact multiplicative inverse of 3^100. Indeed, if we compare using BigFloat:

julia> MultiFloats.use_clean_multifloat_arithmetic(20)

julia> Float64x{13}(3.0^100) * inv(Float64x{13}(3.0^100))
1.00000000000000000 [ ... lots of zeroes ... ] 000009233803766632279

julia> setprecision(BigFloat, 1000)

julia> BigFloat(3.0^100) * BigFloat(inv(Float64x{13}(3.0^100)))
1.00000000000000000 [ ... lots of zeroes ... ] 000009233803766632279027870659735403696888168679729487436063148964080335784965191272233867971

You can see that clean mode correctly produces a more accurate result. Therefore, there is no bug here. The only issue is the false expectation that isone(x * inv(x)) hold for any x whose multiplicative inverse cannot be exactly represented.