Date: Thu, 12 Mar 2015 22:17:50 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: enh <enh@google.com> Cc: freebsd-numerics@freebsd.org Subject: Re: libm functions sensitive to current rounding mode Message-ID: <20150313051749.GA63174@troutmask.apl.washington.edu> In-Reply-To: <CAJgzZor=EgG2%2BE3L6MW-6DD7geUy=Ensa7G1B=viQBXLZKciLw@mail.gmail.com> References: <CAJgzZor=EgG2%2BE3L6MW-6DD7geUy=Ensa7G1B=viQBXLZKciLw@mail.gmail.com>
next in thread | previous in thread | raw e-mail | index | archive | help
On Thu, Mar 12, 2015 at 08:58:06PM -0700, enh via freebsd-numerics wrote:
> various SoC vendors have been sending us (Android) assembler
> implementations of various libm functions, and we noticed when looking
> at one change that it was slower on certain inputs than the FreeBSD C
> implementation we were already using. the response was that the cost
> was coping with other rounding modes, and they gave this example:
>
> #include <math.h>
> #include <stdio.h>
> #include <fenv.h>
>
> int main(){
> double x = -0xb.6365e22ee46fp4;
> fesetround(FE_UPWARD);
> printf("sin(%f) = %.16e\n", x, sin(x));
> fesetround(FE_DOWNWARD);
> printf("sin(%f) = %.16e\n", x, sin(x));
>
> x = 0x1.921fb54442d16p0;
> fesetround(FE_UPWARD);
> printf("sin(%f) = %.16e\n", x, sin(x));
> fesetround(FE_DOWNWARD);
> printf("sin(%f) = %.16e\n", x, sin(x));
> return 0;
> }
>
> see https://android-review.googlesource.com/112700 for the original
> change and related conversation.
>
> glibc 2.19:
>
> sin(-182.212373) = -2.4759225463534308e-18
> sin(-182.212374) = -2.4759225463534309e-18
> sin(1.570797) = 1.0000000000000000e+00
> sin(1.570796) = 1.0000000000000000e+00
>
> (glibc's fixed bugs here in the past. see
> https://sourceware.org/bugzilla/show_bug.cgi?id=3976 for example.)
>
> the FreeBSD libm:
>
> sin(-182.212374) = -2.4759225463534304e-18
> sin(-182.212374) = 2.2575413760606011e-11
> sin(1.570796) = 1.0000000000000000e+00
> sin(1.570796) = 1.0000000002522276e+00
>
> it looks like e_sqrtl.c, s_fmal.c, s_fma.c, and s_fmaf.c save/restore
> the rounding mode but nothing else does.
>
> let me know if you'd like me to send trivial patches for any of the
> affected functions. (do all the libm functions need this?)
>
http://wiki.freebsd.org/Numerics
documents the progress on implementing missing libm functions.
It's somewhat out of date as only powl and tgammal as required
by C99/11 are missing for the long double functions. I don't
know the status of complex.h functions. I don't do wiki, so
haven't even tried to update the page.
At the bottom of the page, you'll see a list of future directions.
Dealing with rounding modes is on the list. As of now, FreeBSD
libm should work with round-to-nearest; other modes probably
need to be tested.
e_sqrtl.c was implemented to handle the rounding mode for 2 reasons:
(1) IEEE-754 mandates that square root is correctly rounded in
all rounding modes; and (2), fdlibm's src/e_sqrt.c (might be
wrong filename) has a long comment that explains an algorithm that
can be used to implement software sqrt() in all rounding modes.
That's what I implemented and then Bruce Evans and David Schultz
fixed what I did.
I suspect all other functions, which are implemented in software,
need to be tested and patches developed. If you have trivial
patches, the best thing to do is submit a bug report via FreeBSD's
bugzilla so that the patch doesn't get lost. Someday after I get
powl and tgammal written, and I become independently wealthy, I may
have time to work on libm.
BTW, n1785.pdf is adding a number new functions to the C library
(e.g., sinpi(x), etc). I suspect I would work on those before
dealing with rounding modes; as all of my personal codes are
expecting round-to-nearest.
--
Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20150313051749.GA63174>