From d1a2ead878c27ac4ec600740320f8b76e1f961e9 Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Tue, 3 Sep 2013 03:27:02 +0000 Subject: math: rewrite rounding functions (ceil, floor, trunc, round, rint) * faster, smaller, cleaner implementation than the bit hacks of fdlibm * use arithmetics like y=(double)(x+0x1p52)-0x1p52, which is an integer neighbor of x in all rounding modes (0<=x<0x1p52) and only use bithacks when that's faster and smaller (for float it usually is) * the code assumes standard excess precision handling for casts * long double code supports both ld80 and ld128 * nearbyint is not changed (it is a wrapper around rint) --- src/math/ceil.c | 92 +++++++++------------------------------------- src/math/ceilf.c | 65 ++++++++++---------------------- src/math/ceill.c | 108 +++++++++++------------------------------------------- src/math/floor.c | 92 +++++++++------------------------------------- src/math/floorf.c | 75 ++++++++++--------------------------- src/math/floorl.c | 108 +++++++++++------------------------------------------- src/math/rint.c | 100 ++++++++------------------------------------------ src/math/rintf.c | 58 ++++++++--------------------- src/math/rintl.c | 86 ++++++++----------------------------------- src/math/round.c | 62 +++++++++++-------------------- src/math/roundf.c | 61 ++++++++++-------------------- src/math/roundl.c | 69 +++++++++++++--------------------- src/math/trunc.c | 70 +++++++---------------------------- src/math/truncf.c | 59 +++++++---------------------- src/math/truncl.c | 74 +++++++++++-------------------------- 15 files changed, 274 insertions(+), 905 deletions(-) diff --git a/src/math/ceil.c b/src/math/ceil.c index 19555180..22dc224c 100644 --- a/src/math/ceil.c +++ b/src/math/ceil.c @@ -1,82 +1,24 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_ceil.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* - * ceil(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to ceil(x). - */ - #include "libm.h" -static const double huge = 1.0e300; - double ceil(double x) { - int32_t i0,i1,j0; - uint32_t i,j; + union {double f; uint64_t i;} u = {x}; + int e = u.i >> 52 & 0x7ff; + double_t y; - EXTRACT_WORDS(i0, i1, x); - // FIXME signed shift - j0 = ((i0>>20)&0x7ff) - 0x3ff; - if (j0 < 20) { - if (j0 < 0) { - /* raise inexact if x != 0 */ - if (huge+x > 0.0) { - if (i0 < 0) { - i0 = 0x80000000; - i1=0; - } else if ((i0|i1) != 0) { - i0=0x3ff00000; - i1=0; - } - } - } else { - i = 0x000fffff>>j0; - if (((i0&i)|i1) == 0) /* x is integral */ - return x; - /* raise inexact flag */ - if (huge+x > 0.0) { - if (i0 > 0) - i0 += 0x00100000>>j0; - i0 &= ~i; - i1 = 0; - } - } - } else if (j0 > 51) { - if (j0 == 0x400) /* inf or NaN */ - return x+x; - return x; /* x is integral */ - } else { - i = (uint32_t)0xffffffff>>(j0-20); - if ((i1&i) == 0) - return x; /* x is integral */ - /* raise inexact flag */ - if (huge+x > 0.0) { - if (i0 > 0) { - if (j0 == 20) - i0 += 1; - else { - j = i1 + (1<<(52-j0)); - if (j < i1) /* got a carry */ - i0 += 1; - i1 = j; - } - } - i1 &= ~i; - } + if (e >= 0x3ff+52 || x == 0) + return x; + /* y = int(x) - x, where int(x) is an integer neighbor of x */ + if (u.i >> 63) + y = (double)(x - 0x1p52) + 0x1p52 - x; + else + y = (double)(x + 0x1p52) - 0x1p52 - x; + /* special case because of non-nearest rounding modes */ + if (e <= 0x3ff-1) { + FORCE_EVAL(y); + return u.i >> 63 ? -0.0 : 1; } - INSERT_WORDS(x, i0, i1); - return x; + if (y < 0) + return x + y + 1; + return x + y; } diff --git a/src/math/ceilf.c b/src/math/ceilf.c index fec945b6..869835f3 100644 --- a/src/math/ceilf.c +++ b/src/math/ceilf.c @@ -1,54 +1,27 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_ceilf.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@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. - * ==================================================== - */ - #include "libm.h" -static const float huge = 1.0e30; - float ceilf(float x) { - int32_t i0,j0; - uint32_t i; + union {float f; uint32_t i;} u = {x}; + int e = (int)(u.i >> 23 & 0xff) - 0x7f; + uint32_t m; - GET_FLOAT_WORD(i0, x); - j0 = ((i0>>23)&0xff) - 0x7f; - if (j0 < 23) { - if (j0 < 0) { - /* raise inexact if x != 0 */ - if (huge+x > 0.0f) { - if (i0 < 0) - i0 = 0x80000000; - else if(i0 != 0) - i0 = 0x3f800000; - } - } else { - i = 0x007fffff>>j0; - if ((i0&i) == 0) - return x; /* x is integral */ - /* raise inexact flag */ - if (huge+x > 0.0f) { - if (i0 > 0) - i0 += 0x00800000>>j0; - i0 &= ~i; - } - } + if (e >= 23) + return x; + if (e >= 0) { + m = 0x007fffff >> e; + if ((u.i & m) == 0) + return x; + FORCE_EVAL(x + 0x1p120f); + if (u.i >> 31 == 0) + u.i += m; + u.i &= ~m; } else { - if (j0 == 0x80) /* inf or NaN */ - return x+x; - return x; /* x is integral */ + FORCE_EVAL(x + 0x1p120f); + if (u.i >> 31) + u.f = -0.0; + else if (u.i << 1) + u.f = 1.0; } - SET_FLOAT_WORD(x, i0); - return x; + return u.f; } diff --git a/src/math/ceill.c b/src/math/ceill.c index a3523f9d..a2cb0a7f 100644 --- a/src/math/ceill.c +++ b/src/math/ceill.c @@ -1,23 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_ceill.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* - * ceill(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to ceill(x). - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -26,77 +6,31 @@ long double ceill(long double x) return ceil(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -#ifdef LDBL_IMPLICIT_NBIT -#define MANH_SIZE (LDBL_MANH_SIZE + 1) -#define INC_MANH(u, c) do { \ - uint64_t o = u.bits.manh; \ - u.bits.manh += (c); \ - if (u.bits.manh < o) \ - u.bits.exp++; \ -} while (0) -#else -#define MANH_SIZE LDBL_MANH_SIZE -#define INC_MANH(u, c) do { \ - uint64_t o = u.bits.manh; \ - u.bits.manh += (c); \ - if (u.bits.manh < o) { \ - u.bits.exp++; \ - u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \ - } \ -} while (0) +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 #endif - -static const long double huge = 1.0e300; - long double ceill(long double x) { - union IEEEl2bits u = { .e = x }; - int e = u.bits.exp - LDBL_MAX_EXP + 1; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; + long double y; - if (e < MANH_SIZE - 1) { - if (e < 0) { - /* raise inexact if x != 0 */ - if (huge + x > 0.0) - if (u.bits.exp > 0 || - (u.bits.manh | u.bits.manl) != 0) - u.e = u.bits.sign ? -0.0 : 1.0; - } else { - uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); - if (((u.bits.manh & m) | u.bits.manl) == 0) - return x; /* x is integral */ - if (!u.bits.sign) { -#ifdef LDBL_IMPLICIT_NBIT - if (e == 0) - u.bits.exp++; - else -#endif - INC_MANH(u, 1llu << (MANH_SIZE - e - 1)); - } - /* raise inexact flag */ - if (huge + x > 0.0) { - u.bits.manh &= ~m; - u.bits.manl = 0; - } - } - } else if (e < LDBL_MANT_DIG - 1) { - uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); - if ((u.bits.manl & m) == 0) - return x; /* x is integral */ - if (!u.bits.sign) { - if (e == MANH_SIZE - 1) - INC_MANH(u, 1); - else { - uint64_t o = u.bits.manl; - u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1); - if (u.bits.manl < o) /* got a carry */ - INC_MANH(u, 1); - } - } - /* raise inexact flag */ - if (huge + x > 0.0) - u.bits.manl &= ~m; + if (e >= 0x3fff+LDBL_MANT_DIG-1 || x == 0) + return x; + /* y = int(x) - x, where int(x) is an integer neighbor of x */ + if (u.i.se >> 15) + y = x - TOINT + TOINT - x; + else + y = x + TOINT - TOINT - x; + /* special case because of non-nearest rounding modes */ + if (e <= 0x3fff-1) { + FORCE_EVAL(y); + return u.i.se >> 15 ? -0.0 : 1; } - return u.e; + if (y < 0) + return x + y + 1; + return x + y; } #endif diff --git a/src/math/floor.c b/src/math/floor.c index ecb9dde8..ebc9fabe 100644 --- a/src/math/floor.c +++ b/src/math/floor.c @@ -1,82 +1,24 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_floor.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* - * floor(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floor(x). - */ - #include "libm.h" -static const double huge = 1.0e300; - double floor(double x) { - int32_t i0,i1,j0; - uint32_t i,j; + union {double f; uint64_t i;} u = {x}; + int e = u.i >> 52 & 0x7ff; + double_t y; - EXTRACT_WORDS(i0, i1, x); - // FIXME: signed shift - j0 = ((i0>>20)&0x7ff) - 0x3ff; - if (j0 < 20) { - if (j0 < 0) { /* |x| < 1 */ - /* raise inexact if x != 0 */ - if (huge+x > 0.0) { - if (i0 >= 0) { /* x >= 0 */ - i0 = i1 = 0; - } else if (((i0&0x7fffffff)|i1) != 0) { - i0 = 0xbff00000; - i1 = 0; - } - } - } else { - i = 0x000fffff>>j0; - if (((i0&i)|i1) == 0) - return x; /* x is integral */ - /* raise inexact flag */ - if (huge+x > 0.0) { - if (i0 < 0) - i0 += 0x00100000>>j0; - i0 &= ~i; - i1 = 0; - } - } - } else if (j0 > 51) { - if (j0 == 0x400) - return x+x; /* inf or NaN */ - else - return x; /* x is integral */ - } else { - i = (uint32_t)0xffffffff>>(j0-20); - if ((i1&i) == 0) - return x; /* x is integral */ - /* raise inexact flag */ - if (huge+x > 0.0) { - if (i0 < 0) { - if (j0 == 20) - i0++; - else { - j = i1+(1<<(52-j0)); - if (j < i1) - i0++; /* got a carry */ - i1 = j; - } - } - i1 &= ~i; - } + if (e >= 0x3ff+52 || x == 0) + return x; + /* y = int(x) - x, where int(x) is an integer neighbor of x */ + if (u.i >> 63) + y = (double)(x - 0x1p52) + 0x1p52 - x; + else + y = (double)(x + 0x1p52) - 0x1p52 - x; + /* special case because of non-nearest rounding modes */ + if (e <= 0x3ff-1) { + FORCE_EVAL(y); + return u.i >> 63 ? -1 : 0; } - INSERT_WORDS(x, i0, i1); - return x; + if (y > 0) + return x + y - 1; + return x + y; } diff --git a/src/math/floorf.c b/src/math/floorf.c index 23fa5e7d..dceec739 100644 --- a/src/math/floorf.c +++ b/src/math/floorf.c @@ -1,64 +1,27 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_floorf.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@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. - * ==================================================== - */ -/* - * floorf(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floorf(x). - */ - #include "libm.h" -static const float huge = 1.0e30; - float floorf(float x) { - int32_t i0,j0; - uint32_t i; + union {float f; uint32_t i;} u = {x}; + int e = (int)(u.i >> 23 & 0xff) - 0x7f; + uint32_t m; - GET_FLOAT_WORD(i0, x); - // FIXME: signed shift - j0 = ((i0>>23)&0xff) - 0x7f; - if (j0 < 23) { - if (j0 < 0) { /* |x| < 1 */ - /* raise inexact if x != 0 */ - if (huge+x > 0.0f) { - if (i0 >= 0) /* x >= 0 */ - i0 = 0; - else if ((i0&0x7fffffff) != 0) - i0 = 0xbf800000; - } - } else { - i = 0x007fffff>>j0; - if ((i0&i) == 0) - return x; /* x is integral */ - /* raise inexact flag */ - if (huge+x > 0.0f) { - if (i0 < 0) - i0 += 0x00800000>>j0; - i0 &= ~i; - } - } + if (e >= 23) + return x; + if (e >= 0) { + m = 0x007fffff >> e; + if ((u.i & m) == 0) + return x; + FORCE_EVAL(x + 0x1p120f); + if (u.i >> 31) + u.i += m; + u.i &= ~m; } else { - if (j0 == 0x80) /* inf or NaN */ - return x+x; - else - return x; /* x is integral */ + FORCE_EVAL(x + 0x1p120f); + if (u.i >> 31 == 0) + u.i = 0; + else if (u.i << 1) + u.f = -1.0; } - SET_FLOAT_WORD(x, i0); - return x; + return u.f; } diff --git a/src/math/floorl.c b/src/math/floorl.c index 3901b060..961f9e89 100644 --- a/src/math/floorl.c +++ b/src/math/floorl.c @@ -1,23 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_floorl.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* - * floorl(x) - * Return x rounded toward -inf to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floorl(x). - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -26,77 +6,31 @@ long double floorl(long double x) return floor(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -#ifdef LDBL_IMPLICIT_NBIT -#define MANH_SIZE (LDBL_MANH_SIZE + 1) -#define INC_MANH(u, c) do { \ - uint64_t o = u.bits.manh; \ - u.bits.manh += (c); \ - if (u.bits.manh < o) \ - u.bits.exp++; \ -} while (0) -#else -#define MANH_SIZE LDBL_MANH_SIZE -#define INC_MANH(u, c) do { \ - uint64_t o = u.bits.manh; \ - u.bits.manh += (c); \ - if (u.bits.manh < o) { \ - u.bits.exp++; \ - u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \ - } \ -} while (0) +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 #endif - -static const long double huge = 1.0e300; - long double floorl(long double x) { - union IEEEl2bits u = { .e = x }; - int e = u.bits.exp - LDBL_MAX_EXP + 1; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; + long double y; - if (e < MANH_SIZE - 1) { - if (e < 0) { - /* raise inexact if x != 0 */ - if (huge + x > 0.0) - if (u.bits.exp > 0 || - (u.bits.manh | u.bits.manl) != 0) - u.e = u.bits.sign ? -1.0 : 0.0; - } else { - uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); - if (((u.bits.manh & m) | u.bits.manl) == 0) - return x; /* x is integral */ - if (u.bits.sign) { -#ifdef LDBL_IMPLICIT_NBIT - if (e == 0) - u.bits.exp++; - else -#endif - INC_MANH(u, 1llu << (MANH_SIZE - e - 1)); - } - /* raise inexact flag */ - if (huge + x > 0.0) { - u.bits.manh &= ~m; - u.bits.manl = 0; - } - } - } else if (e < LDBL_MANT_DIG - 1) { - uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); - if ((u.bits.manl & m) == 0) - return x; /* x is integral */ - if (u.bits.sign) { - if (e == MANH_SIZE - 1) - INC_MANH(u, 1); - else { - uint64_t o = u.bits.manl; - u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1); - if (u.bits.manl < o) /* got a carry */ - INC_MANH(u, 1); - } - } - /* raise inexact flag */ - if (huge + x > 0.0) - u.bits.manl &= ~m; + if (e >= 0x3fff+LDBL_MANT_DIG-1 || x == 0) + return x; + /* y = int(x) - x, where int(x) is an integer neighbor of x */ + if (u.i.se >> 15) + y = x - TOINT + TOINT - x; + else + y = x + TOINT - TOINT - x; + /* special case because of non-nearest rounding modes */ + if (e <= 0x3fff-1) { + FORCE_EVAL(y); + return u.i.se >> 15 ? -1 : 0; } - return u.e; + if (y > 0) + return x + y - 1; + return x + y; } #endif diff --git a/src/math/rint.c b/src/math/rint.c index 775c7b8d..81f4e622 100644 --- a/src/math/rint.c +++ b/src/math/rint.c @@ -1,90 +1,20 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_rint.c */ -/* - * ==================================================== - * 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 "libm.h" - -static const double -TWO52[2] = { - 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ - -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ -}; +#include +#include double rint(double x) { - int32_t i0,j0,sx; - uint32_t i,i1; - double w,t; + union {double f; uint64_t i;} u = {x}; + int e = u.i>>52 & 0x7ff; + int s = u.i>>63; + double_t y; - EXTRACT_WORDS(i0, i1, x); - // FIXME: signed shift - sx = (i0>>31) & 1; - j0 = ((i0>>20)&0x7ff) - 0x3ff; - if (j0 < 20) { - if (j0 < 0) { - if (((i0&0x7fffffff)|i1) == 0) - return x; - i1 |= i0 & 0x0fffff; - i0 &= 0xfffe0000; - i0 |= ((i1|-i1)>>12) & 0x80000; - SET_HIGH_WORD(x, i0); - STRICT_ASSIGN(double, w, TWO52[sx] + x); - t = w - TWO52[sx]; - GET_HIGH_WORD(i0, t); - SET_HIGH_WORD(t, (i0&0x7fffffff)|(sx<<31)); - return t; - } else { - i = 0x000fffff>>j0; - if (((i0&i)|i1) == 0) - return x; /* x is integral */ - i >>= 1; - if (((i0&i)|i1) != 0) { - /* - * Some bit is set after the 0.5 bit. To avoid the - * possibility of errors from double rounding in - * w = TWO52[sx]+x, adjust the 0.25 bit to a lower - * guard bit. We do this for all j0<=51. The - * adjustment is trickiest for j0==18 and j0==19 - * since then it spans the word boundary. - */ - if (j0 == 19) - i1 = 0x40000000; - else if (j0 == 18) - i1 = 0x80000000; - else - i0 = (i0 & ~i)|(0x20000>>j0); - } - } - } else if (j0 > 51) { - if (j0 == 0x400) - return x+x; /* inf or NaN */ - return x; /* x is integral */ - } else { - i = (uint32_t)0xffffffff>>(j0-20); - if ((i1&i) == 0) - return x; /* x is integral */ - i >>= 1; - if ((i1&i) != 0) - i1 = (i1 & ~i)|(0x40000000>>(j0-20)); - } - INSERT_WORDS(x, i0, i1); - STRICT_ASSIGN(double, w, TWO52[sx] + x); - return w - TWO52[sx]; + if (e >= 0x3ff+52) + return x; + if (s) + y = (double)(x - 0x1p52) + 0x1p52; + else + y = (double)(x + 0x1p52) - 0x1p52; + if (y == 0) + return s ? -0.0 : 0; + return y; } diff --git a/src/math/rintf.c b/src/math/rintf.c index e8d44969..9cfc2a26 100644 --- a/src/math/rintf.c +++ b/src/math/rintf.c @@ -1,48 +1,20 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_rintf.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@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. - * ==================================================== - */ - -#include "libm.h" - -static const float -TWO23[2] = { - 8.3886080000e+06, /* 0x4b000000 */ - -8.3886080000e+06, /* 0xcb000000 */ -}; +#include +#include float rintf(float x) { - int32_t i0,j0,sx; - float w,t; + union {float f; uint32_t i;} u = {x}; + int e = u.i>>23 & 0xff; + int s = u.i>>31; + float_t y; - GET_FLOAT_WORD(i0, x); - sx = (i0>>31) & 1; - j0 = ((i0>>23)&0xff) - 0x7f; - if (j0 < 23) { - if (j0 < 0) { - if ((i0&0x7fffffff) == 0) - return x; - STRICT_ASSIGN(float, w, TWO23[sx] + x); - t = w - TWO23[sx]; - GET_FLOAT_WORD(i0, t); - SET_FLOAT_WORD(t, (i0&0x7fffffff)|(sx<<31)); - return t; - } - STRICT_ASSIGN(float, w, TWO23[sx] + x); - return w - TWO23[sx]; - } - if (j0 == 0x80) - return x+x; /* inf or NaN */ - return x; /* x is integral */ + if (e >= 0x7f+23) + return x; + if (s) + y = (float)(x - 0x1p23f) + 0x1p23f; + else + y = (float)(x + 0x1p23f) - 0x1p23f; + if (y == 0) + return s ? -0.0f : 0.0f; + return y; } diff --git a/src/math/rintl.c b/src/math/rintl.c index b13cfeb3..26725073 100644 --- a/src/math/rintl.c +++ b/src/math/rintl.c @@ -1,30 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_rintl.c */ -/*- - * Copyright (c) 2008 David Schultz - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice, this list of conditions and the following disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND - * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE - * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE - * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE - * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL - * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS - * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) - * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT - * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY - * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF - * SUCH DAMAGE. - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -33,53 +6,26 @@ long double rintl(long double x) return rint(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 - -#define BIAS (LDBL_MAX_EXP - 1) - -static const float -shift[2] = { #if LDBL_MANT_DIG == 64 - 0x1.0p63, -0x1.0p63 +#define TOINT 0x1p63 #elif LDBL_MANT_DIG == 113 - 0x1.0p112, -0x1.0p112 +#define TOINT 0x1p112 #endif -}; -static const float zero[2] = { 0.0, -0.0 }; - long double rintl(long double x) { - union IEEEl2bits u; - uint32_t expsign; - int ex, sign; - - u.e = x; - expsign = u.xbits.expsign; - ex = expsign & 0x7fff; - - if (ex >= BIAS + LDBL_MANT_DIG - 1) { - if (ex == BIAS + LDBL_MAX_EXP) - return x + x; /* Inf, NaN, or unsupported format */ - return x; /* finite and already an integer */ - } - sign = expsign >> 15; - - /* - * The following code assumes that intermediate results are - * evaluated in long double precision. If they are evaluated in - * greater precision, double rounding may occur, and if they are - * evaluated in less precision (as on i386), results will be - * wildly incorrect. - */ - x += shift[sign]; - x -= shift[sign]; - - /* - * If the result is +-0, then it must have the same sign as x, but - * the above calculation doesn't always give this. Fix up the sign. - */ - if (ex < BIAS && x == 0.0) - return zero[sign]; - - return x; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; + int s = u.i.se >> 15; + long double y; + + if (e >= 0x3fff+LDBL_MANT_DIG-1) + return x; + if (s) + y = x - TOINT + TOINT; + else + y = x + TOINT - TOINT; + if (y == 0) + return 0*x; + return y; } #endif diff --git a/src/math/round.c b/src/math/round.c index bf0b453f..4b38d1fd 100644 --- a/src/math/round.c +++ b/src/math/round.c @@ -1,48 +1,28 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_round.c */ -/*- - * Copyright (c) 2003, Steven G. Kargl - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice unmodified, this list of conditions, and the following - * disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR - * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. - * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, - * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT - * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, - * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY - * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF - * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - -#include +#include "libm.h" double round(double x) { - double t; + union {double f; uint64_t i;} u = {x}; + int e = u.i >> 52 & 0x7ff; + double_t y; - if (!isfinite(x)) + if (e >= 0x3ff+52) return x; - - if (x >= 0.0) { - t = floor(x); - if (t - x <= -0.5) - t += 1.0; - return t; - } else { - t = floor(-x); - if (t + x <= -0.5) - t += 1.0; - return -t; + if (u.i >> 63) + x = -x; + if (e < 0x3ff-1) { + /* raise inexact if x!=0 */ + FORCE_EVAL(x + 0x1p52); + return 0*u.f; } + y = (double)(x + 0x1p52) - 0x1p52 - x; + if (y > 0.5) + y = y + x - 1; + else if (y <= -0.5) + y = y + x + 1; + else + y = y + x; + if (u.i >> 63) + y = -y; + return y; } diff --git a/src/math/roundf.c b/src/math/roundf.c index 662a454b..c6b27797 100644 --- a/src/math/roundf.c +++ b/src/math/roundf.c @@ -1,48 +1,27 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_roundf.c */ -/*- - * Copyright (c) 2003, Steven G. Kargl - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice unmodified, this list of conditions, and the following - * disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR - * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. - * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, - * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT - * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, - * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY - * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF - * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - -#include +#include "libm.h" float roundf(float x) { - float t; + union {float f; uint32_t i;} u = {x}; + int e = u.i >> 23 & 0xff; + float_t y; - if (!isfinite(x)) + if (e >= 0x7f+23) return x; - - if (x >= 0.0) { - t = floorf(x); - if (t - x <= -0.5) - t += 1.0; - return t; - } else { - t = floorf(-x); - if (t + x <= -0.5) - t += 1.0; - return -t; + if (u.i >> 31) + x = -x; + if (e < 0x7f-1) { + FORCE_EVAL(x + 0x1p23f); + return 0*u.f; } + y = (float)(x + 0x1p23f) - 0x1p23f - x; + if (y > 0.5f) + y = y + x - 1; + else if (y <= -0.5f) + y = y + x + 1; + else + y = y + x; + if (u.i >> 31) + y = -y; + return y; } diff --git a/src/math/roundl.c b/src/math/roundl.c index 99146f07..8f3f28b3 100644 --- a/src/math/roundl.c +++ b/src/math/roundl.c @@ -1,56 +1,39 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_roundl.c */ -/*- - * Copyright (c) 2003, Steven G. Kargl - * All rights reserved. - * - * Redistribution and use in source and binary forms, with or without - * modification, are permitted provided that the following conditions - * are met: - * 1. Redistributions of source code must retain the above copyright - * notice unmodified, this list of conditions, and the following - * disclaimer. - * 2. Redistributions in binary form must reproduce the above copyright - * notice, this list of conditions and the following disclaimer in the - * documentation and/or other materials provided with the distribution. - * - * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR - * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES - * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. - * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, - * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT - * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, - * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY - * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF - * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - -#include -#include +#include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 long double roundl(long double x) { return round(x); } -#else +#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 +#endif long double roundl(long double x) { - long double t; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; + long double y; - if (!isfinite(x)) + if (e >= 0x3fff+LDBL_MANT_DIG-1) return x; - - if (x >= 0.0) { - t = floorl(x); - if (t - x <= -0.5) - t += 1.0; - return t; - } else { - t = floorl(-x); - if (t + x <= -0.5) - t += 1.0; - return -t; + if (u.i.se >> 15) + x = -x; + if (e < 0x3fff-1) { + FORCE_EVAL(x + TOINT); + return 0*u.f; } + y = x + TOINT - TOINT - x; + if (y > 0.5) + y = y + x - 1; + else if (y <= -0.5) + y = y + x + 1; + else + y = y + x; + if (u.i.se >> 15) + y = -y; + return y; } #endif diff --git a/src/math/trunc.c b/src/math/trunc.c index 44b04ecc..d13711b5 100644 --- a/src/math/trunc.c +++ b/src/math/trunc.c @@ -1,63 +1,19 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_trunc.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* - * trunc(x) - * Return x rounded toward 0 to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to trunc(x). - */ - #include "libm.h" -static const double huge = 1.0e300; - double trunc(double x) { - int32_t i0,i1,j0; - uint32_t i; + union {double f; uint64_t i;} u = {x}; + int e = (int)(u.i >> 52 & 0x7ff) - 0x3ff + 12; + uint64_t m; - EXTRACT_WORDS(i0, i1, x); - j0 = ((i0>>20)&0x7ff) - 0x3ff; - if (j0 < 20) { - if (j0 < 0) { /* |x|<1, return 0*sign(x) */ - /* raise inexact if x != 0 */ - if (huge+x > 0.0) { - i0 &= 0x80000000U; - i1 = 0; - } - } else { - i = 0x000fffff>>j0; - if (((i0&i)|i1) == 0) - return x; /* x is integral */ - /* raise inexact */ - if (huge+x > 0.0) { - i0 &= ~i; - i1 = 0; - } - } - } else if (j0 > 51) { - if (j0 == 0x400) - return x + x; /* inf or NaN */ - return x; /* x is integral */ - } else { - i = (uint32_t)0xffffffff>>(j0-20); - if ((i1&i) == 0) - return x; /* x is integral */ - /* raise inexact */ - if (huge+x > 0.0) - i1 &= ~i; - } - INSERT_WORDS(x, i0, i1); - return x; + if (e >= 52 + 12) + return x; + if (e < 12) + e = 1; + m = -1ULL >> e; + if ((u.i & m) == 0) + return x; + FORCE_EVAL(x + 0x1p120f); + u.i &= ~m; + return u.f; } diff --git a/src/math/truncf.c b/src/math/truncf.c index 0afcdfbf..1a7d03c3 100644 --- a/src/math/truncf.c +++ b/src/math/truncf.c @@ -1,52 +1,19 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_truncf.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* - * truncf(x) - * Return x rounded toward 0 to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to truncf(x). - */ - #include "libm.h" -static const float huge = 1.0e30f; - float truncf(float x) { - int32_t i0,j0; - uint32_t i; + union {float f; uint32_t i;} u = {x}; + int e = (int)(u.i >> 23 & 0xff) - 0x7f + 9; + uint32_t m; - GET_FLOAT_WORD(i0, x); - j0 = ((i0>>23)&0xff) - 0x7f; - if (j0 < 23) { - if (j0 < 0) { /* |x|<1, return 0*sign(x) */ - /* raise inexact if x != 0 */ - if (huge+x > 0.0f) - i0 &= 0x80000000; - } else { - i = 0x007fffff>>j0; - if ((i0&i) == 0) - return x; /* x is integral */ - /* raise inexact */ - if (huge+x > 0.0f) - i0 &= ~i; - } - } else { - if (j0 == 0x80) - return x + x; /* inf or NaN */ - return x; /* x is integral */ - } - SET_FLOAT_WORD(x, i0); - return x; + if (e >= 23 + 9) + return x; + if (e < 9) + e = 1; + m = -1U >> e; + if ((u.i & m) == 0) + return x; + FORCE_EVAL(x + 0x1p120f); + u.i &= ~m; + return u.f; } diff --git a/src/math/truncl.c b/src/math/truncl.c index d817e9ad..3eedb083 100644 --- a/src/math/truncl.c +++ b/src/math/truncl.c @@ -1,23 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_truncl.c */ -/* - * ==================================================== - * 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. - * ==================================================== - */ -/* - * truncl(x) - * Return x rounded toward 0 to integral value - * Method: - * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to truncl(x). - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -26,43 +6,31 @@ long double truncl(long double x) return trunc(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#ifdef LDBL_IMPLICIT_NBIT -#define MANH_SIZE (LDBL_MANH_SIZE + 1) -#else -#define MANH_SIZE LDBL_MANH_SIZE +#if LDBL_MANT_DIG == 64 +#define TOINT 0x1p63 +#elif LDBL_MANT_DIG == 113 +#define TOINT 0x1p112 #endif - -static const long double huge = 1.0e300; -static const float zero[] = { 0.0, -0.0 }; - long double truncl(long double x) { - union IEEEl2bits u = { .e = x }; - int e = u.bits.exp - LDBL_MAX_EXP + 1; + union ldshape u = {x}; + int e = u.i.se & 0x7fff; + int s = u.i.se >> 15; + long double y; - if (e < MANH_SIZE - 1) { - if (e < 0) { - /* raise inexact if x != 0 */ - if (huge + x > 0.0) - u.e = zero[u.bits.sign]; - } else { - uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); - if (((u.bits.manh & m) | u.bits.manl) == 0) - return x; /* x is integral */ - /* raise inexact */ - if (huge + x > 0.0) { - u.bits.manh &= ~m; - u.bits.manl = 0; - } - } - } else if (e < LDBL_MANT_DIG - 1) { - uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); - if ((u.bits.manl & m) == 0) - return x; /* x is integral */ - /* raise inexact */ - if (huge + x > 0.0) - u.bits.manl &= ~m; + if (e >= 0x3fff+LDBL_MANT_DIG-1) + return x; + if (e <= 0x3fff-1) { + FORCE_EVAL(x + 0x1p120f); + return x*0; } - return u.e; + /* y = int(|x|) - |x|, where int(|x|) is an integer neighbor of |x| */ + if (s) + x = -x; + y = x + TOINT - TOINT - x; + if (y > 0) + y -= 1; + x += y; + return s ? -x : x; } #endif -- cgit v1.2.1