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

Question

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.)

Answer 1

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.

Answer 2

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

Answer 3

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 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.