Re: incorrectly rounded square root

[Jeff Johnston <jjohnstn@redhat.com>]

On Wed, Jun 2, 2021 at 9:08 AM Joel Sherrill <joel@rtems.org> wrote:

>
>
> On Wed, Jun 2, 2021 at 2:51 AM Paul Zimmermann <Paul.Zimmermann@inria.fr>
> wrote:
>
>>        Hi Jeff,
>>
>> thank you for your answer. After investigating more it seems fesetround()
>> is ineffective with newlib. Indeed if I print the result of fegetround()
>> just after the fesetround() calls I get:
>>
>> fegetround: 0
>> RNDN: 0x1.ff83fp+63
>> fegetround: 0
>> RNDZ: 0x1.ff83fp+63
>> fegetround: 0
>> RNDU: 0x1.ff83fp+63
>> fegetround: 0
>> RNDD: 0x1.ff83fp+63
>>
>> whereas with GNU libc I get:
>>
>> fegetround: 0
>> RNDN: 0x1.ff83fp+63
>> fegetround: 3072
>> RNDZ: 0x1.ff83eep+63
>> fegetround: 2048
>> RNDU: 0x1.ff83fp+63
>> fegetround: 1024
>> RNDD: 0x1.ff83eep+63
>>
>> Thus it seems this has nothing to do with the square root.
>>
>> According to gdb the fesetround() code used is the following:
>>
>> (gdb) break fesetround
>> Breakpoint 1 at 0x1650: file
>> ../../../../../../newlib/libm/machine/x86_64/fenv.c, line 371.
>>
>
> That is close to what I see in the source. That matches a variable
> declaration in fesetround().
>
> I was concerned that maybe the stub magic was resulting in empty
> methods getting in for one of these. But doing an objdump on i386-rtems,
> that doesn't appear to be the case.
>
> The implementation came from Cygwin and was lightly massaged to
> get here. A student and I migrated that code to newlib in Feb 2020.
> Corinna did some work in March 2021 to use this for Cygwin but I
> don't see any reason it is broken. The i386-rtems installed sys/fenv.h
> is the one for x86.
>
> fegetround() is quite simple and fesetround isn't much more than that.
>
>
> https://sourceware.org/git/?p=newlib-cygwin.git;a=blob;f=newlib/libm/machine/shared_x86/fenv.c;h=ccc08e2d8103ab974baf8d9591f71e5565d73ace;hb=HEAD#l350
>
>  I really expected to see a build issue and you using a stub. Could you
> confirm that the code we expect to be in those methods really is there
> for you? The line number indicates you have the expected code but....
>
> The rounding is set in the x87 control register and then if SSE is there
> (dynamic check), it is also set in the SSE control register. The rounding
> more is read from the x86 control register and the SSE control register is
> ignored.
>
> Could you step through get and set and see if it looks like the x87
> control register is actually changed?
>
> I'm confused. This looks like it should work.
>
> --joel
>
>
I'll second Joel's comment.  The code is extremely close to the glibc code
both in sqrtf and fesetround.  The only
thing I can think of is that the glibc code does the x87 stuff first and
does the set back into FPU state before doing the
SSE stuff.  The newlib code sets back the FPU state at the end after the
SSE stuff.  Don't know if this is relevant or not.

Any Cygwin users out there who can verify that the code is working/not
working for them?

-- Jeff J.


>> Paul
>>
>>
>>