incorrectly rounded square root

[Jeff Johnston]

Hi Paul,

I figured the values were off when I had to hard-code them in my own
test_sqrt.c but forgot to include that info in my note.

Now, that said, using the code I attached earlier, I am seeing the exact
values you are quoting above for glibc for the mxcsr register and the round
is working.  Have your
tried running that code?

The mxcsr values you are seeing that are different are not due to the
fesetround code.  The code is shifting the round value 13 bits
and for 3, that ends up being 0x6000.  It is masking mxcsr with 0xffff9fff
first so when you start with 0x1fxx and end up with 0x7fxx, the code is
doing what is supposed to do.
The difference in values above is 0x20 (e.g. 0x7fa0 vs 0x7f80) which is a
bit in the last 2 hex digits which isn't touched by the code logic.

-- Jeff J.



On Thu, Jun 3, 2021 at 6:22 AM Paul Zimmermann <Paul.Zimmermann@inria.fr>
wrote:

>        Hi Jeff,
>
> I've investigated and found (partially) the cause of the problem. When I
> compile my program, gcc uses the system fenv.h, which contains
> FE_TONEAREST=0,
> FE_DOWNWARD=0x400, FE_UPWARD=0x800, and FE_TOWARDZERO=0xc00 in
> /usr/include/bits/fenv.h.
>
> However, the values in x86_64/newlib/targ-include/sys/fenv.h are different:
> FE_TONEAREST=0, FE_DOWNWARD=1, FE_UPWARD=2, FE_TOWARDZERO=3. Thus the test
> at the beginning of fesetround in newlib/libm/machine/x86_64/fenv.c fails
> except for round=FE_TONEAREST:
>
>   /* Will succeed for any valid value of the input parameter.  */
>   if (round < FE_TONEAREST || round > FE_TOWARDZERO)
>     return EINVAL;
>
> enter fesetround, round=0 use_sse=1
> enter fesetround, round=3072 use_sse=1
> enter fesetround, round=2048 use_sse=1
> enter fesetround, round=1024 use_sse=1
>
> This explains why the other modes are not set.
>
> A first question is why Newlib and the system libm (here glibc) use
> different
> constants for the rounding modes.
>
> Now if I compile with -I/tmp/x86_64/include (where the Newlib fenv.h is
> installed) I get:
>
> enter fesetround, round=0 use_sse=1
> enter fesetround, round=3 use_sse=1
> enter fesetround, round=2 use_sse=1
> enter fesetround, round=1 use_sse=1
>
> but I still do get the same results whatever the rounding mode!
>
> If I look at the cw and mxcsr registers like in Newlib fenv.c, I get with
> glibc
> (after the fesetround call):
>
> fegetround: 0
> cw=895 mxcsr=8064
> RNDN: 0x1.ff83fp+63
> fegetround: 3072
> cw=3967 mxcsr=32672
> RNDZ: 0x1.ff83eep+63
> fegetround: 2048
> cw=2943 mxcsr=24480
> RNDU: 0x1.ff83fp+63
> fegetround: 1024
> cw=1919 mxcsr=16288
> RNDD: 0x1.ff83eep+63
>
> and with Newlib:
>
> fegetround: 0
> cw=895 mxcsr=8064
> RNDN: 0x1.ff83fp+63
> fegetround: 3
> cw=3967 mxcsr=32640
> RNDZ: 0x1.ff83fp+63
> fegetround: 2
> cw=2943 mxcsr=24448
> RNDU: 0x1.ff83fp+63
> fegetround: 1
> cw=1919 mxcsr=16256
> RNDD: 0x1.ff83fp+63
>
> The cw values do match but not the mxcsr values.
>
> Best regards,
> Paul
>
> > Something  weird is going on.  I condensed the newlib code into the
> > following files and when I compile and run the test,
> > it works as expected:
> >
> > gcc -std=c99 -g3 -O0 test_sqrt.c fesetround.c ef_sqrt.c
> > [jjohnstn@localhost shared_x86]$ ./a.out
> > RNDN: 0x1.ff83fp+63
> > RNDZ: 0x1.ff83eep+63
> > RNDU: 0x1.ff83fp+63
> > RNDD: 0x1.ff83eep+63
> >
> > Can you verify that you are calling the fesetround in fenv.c and what
> > configuration call did you use for building newlib?
>
>