Is there any way to get correct rounding with the i387 fsqrt instruction?

Is there any way to get correct rounding with the i387 fsqrt instruction?...

...**aside from changing the precision mode** in the x87 control word - I know that's possible, but it's not a reasonable solution because it has nasty reentrancy-type issues where the precision mode will be wrong if the sqrt operation is interrupted.

The issue I'm dealing with is as follows: the x87 `fsqrt` opcode performs a correctly-rounded (per IEEE 754) square root operation in the precision of the fpu registers, which I'll assume is extended (80-bit) precision. However, I want to use it to implement efficient single and double precision square root functions with the results correctly rounded (per the current rounding mode). Since the result has excess precision, the second step of converting the result to single or double precision rounds again, possibly leaving a not-correctly-rounded result.

With some operations it's possible to work around this with biases. For instance, I can avoid excess precision in the results of addition by adding a bias in the form of a power of two that forces the 52 significant bits of a double precision value into the last 52 bits of the 63-bit extended-precision mantissa. But I don't see any obvious way to do such a trick with square root.

Any clever ideas?

(Also tagged C because the intended application is implementation of the C `sqrt` and `sqrtf` functions.)

It may not be what you want, as it doesn't take advantage of the 387 `fsqrt` instruction, but there's a surprisingly efficient `sqrtf(float)` in [glibc][1] implemented with 32-bit integer arithmetic. It even handles NaNs, Infs, subnormals correctly - it might be possible to eliminate some of these checks with real x87 instructions / FP control word flags. see: `glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c`

The `dbl-64/e_sqrt.c` code is not so friendly. It's hard to tell what assumptions are being made at a glance. Curiously, the library's i386 `sqrt[f|l]` implementations just call `fsqrt`, but load the value differently. `flds` for SP, `fldl` for DP.


[1]:

[To see links please register here]

First, let's get the obvious out of the way: you should be using SSE instead of x87. The SSE `sqrtss` and `sqrtsd` instructions do exactly what you want, are supported on all modern x86 systems, and are significantly faster as well.

Now, if you insist on using x87, I'll start with the good news: you don't need to do anything for float. You need `2p + 2` bits to compute a correctly rounded square-root in a p-bit floating-point format. Because `80 > 2*24 + 2`, the additional rounding to single-precision will always round correctly, and you have a correctly rounded square root.

Now the bad news: `80 < 2*53 + 2`, so no such luck for double precision. I can suggest several workarounds; here's a nice easy one off the top of my head.

1. let `y = round_to_double(x87_square_root(x));`
2. use a Dekker (head-tail) product to compute `a` and `b` such that `y*y = a + b` exactly.
3. compute the residual `r = x - a - b`.
4. `if (r == 0) return y`
5. `if (r > 0)`, let `y1 = y + 1 ulp`, and compute `a1`, `b1` s.t. `y1*y1 = a1 + b1`. Compare `r1 = x - a1 - b1` to `r`, and return either `y` or `y1`, depending on which has the smaller residual (or the one with zero low-order bit, if the residuals are equal in magnitude).
6. `if (r < 0)`, do the same thing for `y1 = y - 1 ulp`.

This proceedure only handles the default rounding mode; however, in the directed rounding modes, simply rounding to the destination format does the right thing.

OK, I think I have a better solution:

1. Compute `y=sqrt(x)` in extended precision (`fsqrt`).
2. If last 11 bits are not `0x400`, simply convert to double precision and return.
3. Add `0x100-(fpu_status_word&0x200)` to the low word of the extended precision representation.
4. Convert to double precision and return.

Step 3 is based on the fact that the C1 bit (0x200) of the status word is 1 if and only if `fsqrt`'s result was rounded up. This is valid because, due to the test in step 2, `x` was not a perfect square; if it were a perfect square, `y` would have no bits beyond double precision.

It may be faster to perform step 3 with a conditional floating point operating rather than working on the bit representation and reloading.

Here's the code (seems to work in all cases):

sqrt:
fldl 4(%esp)
fsqrt
fstsw %ax
sub $12,%esp
fld %st(0)
fstpt (%esp)
mov (%esp),%ecx
and $0x7ff,%ecx
cmp $0x400,%ecx
jnz 1f
and $0x200,%eax
sub $0x100,%eax
sub %eax,(%esp)
fstp %st(0)
fldt (%esp)
1: add $12,%esp
fstpl 4(%esp)
fldl 4(%esp)
ret