sqrt(Float128) accuracy

[JeffreySarnoff]

There may be a slight flaw in the sqrt routine:

julia> using DoubleFloats, QD

julia> Float128(2) - sqrt(Float128(2))^2
9.86076131526264756764660706603503e-32

julia> Double64(2) - sqrt(Double64(2))^2
0.0

julia> Float64((sqrt(BigFloat(2)) - sqrt(Float128(2)))/sqrt(BigFloat(2)))
2.0358013360121176e-32

julia> Float64((sqrt(BigFloat(2)) - sqrt(Double64(2)))/sqrt(BigFloat(2)))
1.164224936801593e-32

[eschnett]

It seems that the qd library advertises an epsilon of 2^-104, or about 4.930380657631324e-32. Given this, the error is here acceptable (less than eps/2). That doesn't mean that a more accurate result wouldn't be desirable. In fact, I found other cases where the relative error is larger than eps.

It seems the algorithmic difference is that qd's sqrt uses one Newton iteration after the initial Float64 approximation, whereas DoubleFloats uses two Newton iterations.

[JeffreySarnoff]

Yes. There are many functions in qd that do not perform as advertised. Most of the time in developing DoubleFloats had been in providing algorithms that are more accurate, generally, than qd.

[eschnett]

@JeffreySarnoff That's good to know, thank you.