search

dzhang314 / MultiFloats.jl

Fast, SIMD-accelerated extended-precision arithmetic for Julia
MIT License
77 stars 10 forks source link

Catastrophic loss of precision in `hypot` #30

Closed mwallerb closed 1 year ago

mwallerb commented 2 years ago

hypot is just not correctly implemented. it loses all digits of precision for small as well as large values.

using MultiFloats
small = 1e-300
@assert hypot(small, small) / sqrt(2) ≈ 1e-300
xsmall = Float64x2(small)
@assert hypot(xsmall, xsmall) / sqrt(2) ≈ 1e-300

One should probably have some really robust precision audit as part of the unit tests.

SamuelBadr commented 2 years ago

The implementation in Base is generic, so funnily enough, if you just remove the stopgap implementation, everything works. (This is what DoubleFloats.jl does.) The computation is roughly 3.3x faster than with DoubleFloats.jl on my machine.

lrnv commented 2 years ago

It is amazing how a stop-gap needed 3 years ago to obtain compatibility with GenericLinearAlgebra.jl is now not needed anymore. The generics have improved a lot!

dzhang314 commented 1 year ago

Thanks for pointing this out @mwallerb! As @lrnv said, this stopgap was added years ago to make MultiFloats.jl work with GenericLinearAlgebra.jl, and was never intended to actually do the job hypot(x, y) (i.e., avoiding undue overflow and underflow). I'm very glad to see that the generics have improved to the point that these sorts of hacks are no longer necessary.