X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Fceil.c;h=b30902cbfd6dc059e2afbc7d1cf974da67fc29a9;hb=b344de996cd51f8a2f2558a3172016b64d99c622;hp=169c68fe5e09dac269e1575f0381ccf48c00a95a;hpb=d9cd580560f5ba3bd4c669d78321eb9947f4d568;p=gnulib.git diff --git a/lib/ceil.c b/lib/ceil.c index 169c68fe5..b30902cbf 100644 --- a/lib/ceil.c +++ b/lib/ceil.c @@ -1,5 +1,5 @@ /* Round towards positive infinity. - Copyright (C) 2007 Free Software Foundation, Inc. + Copyright (C) 2007, 2010-2011 Free Software Foundation, Inc. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -23,23 +23,35 @@ #include +#undef MIN + #ifdef USE_LONG_DOUBLE # define FUNC ceill # define DOUBLE long double # define MANT_DIG LDBL_MANT_DIG +# define MIN LDBL_MIN # define L_(literal) literal##L #elif ! defined USE_FLOAT # define FUNC ceil # define DOUBLE double # define MANT_DIG DBL_MANT_DIG +# define MIN DBL_MIN # define L_(literal) literal #else /* defined USE_FLOAT */ # define FUNC ceilf # define DOUBLE float # define MANT_DIG FLT_MANT_DIG +# define MIN FLT_MIN # define L_(literal) literal##f #endif +/* -0.0. See minus-zero.h. */ +#if defined __hpux || defined __sgi || defined __ICC +# define MINUS_ZERO (-MIN * MIN) +#else +# define MINUS_ZERO L_(-0.0) +#endif + /* 2^(MANT_DIG-1). */ static const DOUBLE TWO_MANT_DIG = /* Assume MANT_DIG <= 5 * 31. @@ -63,20 +75,35 @@ FUNC (DOUBLE x) volatile DOUBLE y = x; volatile DOUBLE z = y; - /* Round to the next integer (nearest or up or down, doesn't matter). */ if (z > L_(0.0)) { - z += TWO_MANT_DIG; - z -= TWO_MANT_DIG; + /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */ + if (z < TWO_MANT_DIG) + { + /* Round to the next integer (nearest or up or down, doesn't matter). */ + z += TWO_MANT_DIG; + z -= TWO_MANT_DIG; + /* Enforce rounding up. */ + if (z < y) + z += L_(1.0); + } } else if (z < L_(0.0)) { - z -= TWO_MANT_DIG; - z += TWO_MANT_DIG; + /* For -1 < x < 0, return -0.0 regardless of the current rounding + mode. */ + if (z > L_(-1.0)) + z = MINUS_ZERO; + /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */ + else if (z > - TWO_MANT_DIG) + { + /* Round to the next integer (nearest or up or down, doesn't matter). */ + z -= TWO_MANT_DIG; + z += TWO_MANT_DIG; + /* Enforce rounding up. */ + if (z < y) + z += L_(1.0); + } } - /* Enforce rounding up. */ - if (z < y) - z += L_(1.0); - return z; }