From: Brian Inglis <Brian.Inglis@SystematicSw.ab.ca>
To: newlib@sourceware.org
Subject: Re: incorrectly rounded square root
Date: Sat, 5 Jun 2021 07:25:27 -0600 [thread overview]
Message-ID: <0a785f1d-4a6f-f880-a60a-05c68948f932@SystematicSw.ab.ca> (raw)
In-Reply-To: <CAOox84vMKMnqi1uGyO8tBO+5mjmnq_oAxqsqogCBz0c9Q--Rbw@mail.gmail.com>
Great catch, analysis, and fix!
Now sqrtf rounds correctly on Cygwin!
$ ./test-sqrtf-round
Direction CW MX Input Hex
Input Decimal Sqrt Hex
Sqrt Decimal
RNDN 0 0 37f 1f80: 0x1.ff07fe00p+127
339638501828070541185766401939693633536 0x1.ff83f000p+63
18429283829060468736
RNDD 1 1 77f 3f80: 0x1.ff07fe00p+127
339638501828070541185766401939693633536 0x1.ff83ee00p+63
18429282729548840960
RNDU 2 2 b7f 5f80: 0x1.ff07fe00p+127
339638501828070541185766401939693633536 0x1.ff83f000p+63
18429283829060468736
RNDZ 3 3 f7f 7f80: 0x1.ff07fe00p+127
339638501828070541185766401939693633536 0x1.ff83ee00p+63
18429282729548840960
--
Take care. Thanks, Brian Inglis, Calgary, Alberta, Canada
This email may be disturbing to some readers as it contains
too much technical detail. Reader discretion is advised.
On 2021-06-04 12:44, Jeff Johnston wrote:
> Ok, I now know exactly what is happening.
> The compiler is optimizing out the rounding check in ef_sqrt.c, probably
> due to the operation using two constants.
> 86 ix += (m <<23);
> (gdb) list
> 81 else
> 82 q += (q&1);
> When I debug, it always does the else at line 81 without performing the
> one-tiny operation. The difference in the mxcsr
> register is the PE bit which I believe gets set when you do the one-tiny
> operation. Since we aren't doing it, it never gets
> set on and the difference of 0x20 in the mxcsr register is explained.
> By making the constants volatile, I am able to get the code working
> as it should. I have pushed a patch for this.
> On Fri, Jun 4, 2021 at 3:14 AM Paul Zimmermann wrote:
>>> 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?
>> yes it works as expected, but it doesn't work with Newlib's fenv.h and
>> libm.a (see below).
>>> 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.
>> here is how to reproduce the issue:
>> tar xf newlib-4.1.0.tar.gz
>> cd newlib-4.1.0
>> mkdir build
>> cd build
>> ../configure --prefix=/tmp --disable-multilib --target=x86_64
>> make -j4
>> make install
>> $ cat test_sqrt_2.c
>> #include <stdio.h>
>> #include <math.h>
>> #include <fenv.h>
>> #ifdef NEWLIB
>> /* RedHat's libm claims:
>> undefined reference to `__errno' in j1f/y1f */
>> int errno;
>> int* __errno () { return &errno; }
>> #endif
>> int main()
>> {
>> int rnd[4] = { FE_TONEAREST, FE_TOWARDZERO, FE_UPWARD, FE_DOWNWARD };
>> char Rnd[4] = "NZUD";
>> float x = 0x1.ff07fep+127f;
>> float y;
>> for (int i = 0; i < 4; i++)
>> {
>> unsigned short cw;
>> unsigned int mxcsr = 0;
>> fesetround (rnd[i]);
>> __asm__ volatile ("fnstcw %0" : "=m" (cw) : );
>> __asm__ volatile ("stmxcsr %0" : "=m" (mxcsr) : );
>> y = sqrtf (x);
>> printf ("RND%c: %a cw=%u mxcsr=%u\n", Rnd[i], y, cw, mxcsr);
>> }
>> }
>> With GNU libc:
>> $ gcc -fno-builtin test_sqrt_2.c -lm
>> $ ./a.out
>> RNDN: 0x1.ff83fp+63 cw=895 mxcsr=8064
>> RNDZ: 0x1.ff83eep+63 cw=3967 mxcsr=32672
>> RNDU: 0x1.ff83fp+63 cw=2943 mxcsr=24480
>> RNDD: 0x1.ff83eep+63 cw=1919 mxcsr=16288
>> With Newlib:
>> $ gcc -I/tmp/x86_64/include -DNEWLIB -fno-builtin test_sqrt_2.c /tmp/libm.a
>> $ ./a.out
>> RNDN: 0x1.ff83fp+63 cw=895 mxcsr=8064
>> RNDZ: 0x1.ff83fp+63 cw=3967 mxcsr=32640
>> RNDU: 0x1.ff83fp+63 cw=2943 mxcsr=24448
>> RNDD: 0x1.ff83fp+63 cw=1919 mxcsr=16256
>> Can you reproduce that on x86_64 Linux?
next prev parent reply other threads:[~2021-06-05 13:25 UTC|newest]
Thread overview: 19+ messages / expand[flat|nested] mbox.gz Atom feed top
2021-05-04 8:08 Paul Zimmermann
2021-05-31 20:52 ` Jeff Johnston
2021-06-01 7:11 ` Paul Zimmermann
2021-06-01 16:28 ` Jeff Johnston
2021-06-02 7:51 ` Paul Zimmermann
2021-06-02 13:07 ` Joel Sherrill
2021-06-02 18:43 ` Jeff Johnston
2021-06-02 19:07 ` Marco Atzeri
2021-06-02 19:12 ` Joel Sherrill
2021-06-03 3:01 ` Jeff Johnston
2021-06-03 10:21 ` Paul Zimmermann
2021-06-03 15:50 ` Jeff Johnston
2021-06-04 7:14 ` Paul Zimmermann
2021-06-04 18:44 ` Jeff Johnston
2021-06-04 18:59 ` Joel Sherrill
2021-06-05 13:25 ` Brian Inglis [this message]
2021-06-07 9:51 ` Paul Zimmermann
2021-06-12 22:31 ` Maciej W. Rozycki
2021-06-03 11:25 ` Paul Zimmermann
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=0a785f1d-4a6f-f880-a60a-05c68948f932@SystematicSw.ab.ca \
--to=brian.inglis@systematicsw.ab.ca \
--cc=newlib@sourceware.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).