Re: incorrectly rounded square root
Jeff Johnston <[email protected]>
| Newsgroups | gmane.comp.lib.newlib |
|---|---|
| Message-ID | <CAOox84v54CvMzcmBY_Ax1aYV4rZDE7VQRKsug5-5qAcaFU+H0Q@mail.gmail.com> |
Ok, not sure what is happening.
You should be using the newlib/libm/machine/shared_x86/fenv.c which
includes fesetround which you can confirm.
Comparing to glibc's version of fesetround, it is very similar. The only
difference appears to be that glibc does the x87 fpu first, then the MSCSR
register for sse,
where newlib's version intermixes the code, but the same insn's are used.
Looking at the sqrtf implementation of newlib's libm/math/ef_sqrt.c vs
glibc's sysdeps/ieee754/flt-32/e_sqrtf.c
the rounding logic is identical.
/* use floating add to find out rounding direction */
if(ix!=0) {
z = one-tiny; /* trigger inexact flag */
if (z>=one) {
z = one+tiny;
if (z>one)
q += 2;
else
q += (q&1);
}
}
ix = (q>>1)+0x3f000000;
ix += (m <<23);
SET_FLOAT_WORD(z,ix);
return z;
In fact, other than the normalization code, the code is identical (not sure
that difference is applicable here).
Corinna is on vacation so I'm cc'ing Joel in case he can offer some input.
-- Jeff J.
On Tue, Jun 1, 2021 at 3:12 AM Paul Zimmermann <[email protected]>
wrote:
> Hi Jeff,
>
> > Not all platforms supply a machine implementation of fesetround. Which
> > platform are you using?
>
> this is on a Intel(R) Core(TM) i5-4590, with gcc 10.2.1, under
> Debian GNU/Linux 11 (bullseye).
>
> Best regards,
> Paul
>
>