libm functions sensitive to current rounding mode

Steve Kargl sgk at
Fri Mar 13 05:34:41 UTC 2015

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 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
> 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?)

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.


More information about the freebsd-numerics mailing list