performance problems and imprecisions in sqrt_real

[cfbolz]

While working on support for the Sail Real type in Pydrofoil, I've hit on problems with sqrt_real in sail.c. I'm going to document them in this issue and would like opinions on what to do about them. The high level summary is that I discovered cases where the sqrt implementation is incredibly slow and cases where it returns really imprecise results. I am not a floating point/numerics expert, so please take all this with a grain of salt.

Background

Sail's Real type is implemented using gmp rationals. in the C backend. Almost the operations on Reals that Sail supports (addition, subtraction, ...) can be computed precisely on rationals and the implementation is performed by the really well tested gmp library.

The only exception to this is sqrt, which produces imprecise results on rational numbers (thanks Euclid). Sail's C and ML backend implement it by approximating the result up to an absolute precision of 1e-30, if the argument is not a perfect (integer) square. The approximation is performed using Heron's method. It's also notable that this implementation is really part of Sail, so a lot less tested than functionality that comes from GMP.

sqrt on Reals is used by at least the ARM model. As far as I understand it, the ARM model implements floating point operations in the following way: The operands of an fp instructions are are converted to a Sail Real, then the operation is performed on the Reals, then the result is converted back to the nearest representable floating point number taking rounding etc into account.

The problems with sqrt I'm describing below affect the ARM model directly and can quite likely lead to incorrect behaviour of the fsqrt instruction, but I didn't work on writing an ARM program that demonstrates this.

Problem 1: slow convergence due to bad initial guess for Heron's method

The first most obvious problem is that the sqrt implementation picks a bad initial guess for the result of sqrt for numbers that are smaller than 1. The guess is computed as

$$\frac{\mathrm{isqrt}\left\lceil \frac{a}{b}\right\rceil }{1}$$

whith is $= 1$ for numbers smaller than 1. This leads to extremely slow convergence in some cases, which can be seen by running the following Sail program:

$include <prelude.sail>
$include <arith.sail>
$include <real.sail>

val main : unit -> unit

function main() = {
    let a = 0.00000000015;
    print_real("a: ", a);
    let b = sqrt(a);
    print_real("sqrt: ", b);
}

Which takes more than two minutes on my (old) laptop.

This can be fixed using a better first guess

$$\frac{\mathrm{isqrt}a}{\mathrm{isqrt}b}$$

Which reduces the runtime of the program to less than a second.

Problem 2: Not enough precision

The next problem is that the convergence condition of Heron's method of sqrt uses an absolute difference of 1e-30, which leads to very imprecise results for reals that are very close to 0. Instead, a relative convergence condition should be used for small enough numbers. Otherwise, the results of taking the sqrt of a small number like 4.6270461409975375e-160 can yield completely wrong results.

Solution ideas

In my opinion the cleanest solution is to remove sqrt from Sail entirely, since the existence of the function promises a precise result that the implementation cannot deliver (unless you want to switch the implementation of Reals to some algebraic number package, but that is its own huge can of worms).

A more pragmatic option would be to have extra arguments to specify the required precision. Given that sqrt has to be imprecise on rationals it should be up to the model to know how much precision it needs.

Even more pragmatic would be to fix the bugs and pick a precision that is at least good enough for IEEE doubles (because that's what the ARM model supports). I am going to open a PR for that so you can see what it would look like.

Thoughts?

[Alasdair]

In practice we only have the real type to support the ARM specification. The RISC-V model doesn't use them (thankfully), although using softfloat as a library has it's own issues. I think your solution of having enough precision for doubles is probably good, but having the second option with a parameter for the required precision the model can specify seems slightly more principled, and we could patch the generated Sail from ASL to have that extra argument.

I vaguely remember talking with Alastair Reid about the sqrt function in the past, and I think he mentioned something about tracking the required precision along with the rational in the type, so their real type was something a (rational * required_precision) pair or a (precision -> rational) function, so you would always get the needed precision when you converted back to a bit vector.

[cfbolz]

Thanks Alasdair! Would you maybe approve the CI workflows on the PR so I can check that I didn't break anything? (though maybe I need to mark it as not draft first?)

Maybe I'll ping Alastair Reid and ask whether he still has opinions.

[bacam]

When we were speaking to some Arm spec people a while back they said that it should be possible to split out the floating point at a roughly IEEE-ish interface. Unfortunately we don't currently have a Sail floating point library to use instead.