diff options
| author | Ulrich Drepper <drepper@redhat.com> | 2009-08-25 23:26:16 -0700 |
|---|---|---|
| committer | Ulrich Drepper <drepper@redhat.com> | 2009-08-25 23:26:16 -0700 |
| commit | 77a1e0873f98b57d9e50138b471be0d8db8da99f (patch) | |
| tree | 7fd21f70c118a34ea3264647b7790a84d579d516 | |
| parent | e6fd9b2490295b3c6901e06b94ea077ac1f9fc79 (diff) | |
| download | glibc-77a1e0873f98b57d9e50138b471be0d8db8da99f.tar.gz glibc-77a1e0873f98b57d9e50138b471be0d8db8da99f.tar.xz glibc-77a1e0873f98b57d9e50138b471be0d8db8da99f.zip | |
Optimized nearbyint for 64-bit.
| -rw-r--r-- | ChangeLog | 2 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c | 86 |
2 files changed, 88 insertions, 0 deletions
diff --git a/ChangeLog b/ChangeLog index af3d0578bd..a97b7d9604 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,7 @@ 2009-08-25 Ulrich Drepper <drepper@redhat.com> + * sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c: New file + * sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c: New file. * sysdeps/ieee754/dbl-64/wordsize-64/s_round.c: New file. diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c new file mode 100644 index 0000000000..cb49019ddb --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c @@ -0,0 +1,86 @@ +/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * rint(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rint(x). + */ + +#include <fenv.h> +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const double +#else +static double +#endif +TWO52[2]={ + 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ + -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ +}; + +#ifdef __STDC__ + double __nearbyint(double x) +#else + double __nearbyint(x) + double x; +#endif +{ + fenv_t env; + int64_t i0,sx; + int32_t j0; + EXTRACT_WORDS64(i0,x); + sx = (i0>>63)&1; + j0 = ((i0>>52)&0x7ff)-0x3ff; + if(j0<52) { + if(j0<0) { + if((i0&UINT64_C(0x7fffffffffffffff))==0) return x; + uint64_t i = i0 & UINT64_C(0xfffffffffffff); + i0 &= UINT64_C(0xfffe000000000000); + i0 |= (((i|-i) >> 12) & UINT64_C(0x8000000000000)); + INSERT_WORDS64(x,i0); + feholdexcept (&env); + double w = TWO52[sx]+x; + double t = w-TWO52[sx]; + fesetenv (&env); + EXTRACT_WORDS64(i0,t); + INSERT_WORDS64(t,(i0&UINT64_C(0x7fffffffffffffff))|(sx<<63)); + return t; + } else { + uint64_t i = UINT64_C(0x000fffffffffffff)>>j0; + if((i0&i)==0) return x; /* x is integral */ + i>>=1; + if((i0&i)!=0) + i0 = (i0&(~i))|(UINT64_C(0x4000000000000)>>j0); + } + } else { + if(j0==0x400) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } + INSERT_WORDS64(x,i0); + feholdexcept (&env); + double w = TWO52[sx]+x; + double t = w-TWO52[sx]; + fesetenv (&env); + return t; +} +weak_alias (__nearbyint, nearbyint) +#ifdef NO_LONG_DOUBLE +strong_alias (__nearbyint, __nearbyintl) +weak_alias (__nearbyint, nearbyintl) +#endif |