X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Ftrunc.c;h=d28398b16190243ddce617ffa3a4e6e81199b373;hb=23eecb48e39afd0d267d64d40ba6bf97aa865e13;hp=385e2d139cdf03cc4b4cd94c50be9ea534f425e1;hpb=57fdfd3f8ec62b105c53bcdf6f127c35c7fe7391;p=gnulib.git diff --git a/lib/trunc.c b/lib/trunc.c index 385e2d139..d28398b16 100644 --- a/lib/trunc.c +++ b/lib/trunc.c @@ -1,5 +1,5 @@ /* Round towards zero. - Copyright (C) 2007 Free Software Foundation, Inc. + Copyright (C) 2007, 2010-2013 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 @@ -16,32 +16,52 @@ /* Written by Bruno Haible , 2007. */ -#include +#if ! defined USE_LONG_DOUBLE +# include +#endif /* Specification. */ #include #include +#undef MIN + #ifdef USE_LONG_DOUBLE # define FUNC truncl # 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 trunc # define DOUBLE double # define MANT_DIG DBL_MANT_DIG +# define MIN DBL_MIN # define L_(literal) literal #else /* defined USE_FLOAT */ # define FUNC truncf # 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 + +/* MSVC with option -fp:strict refuses to compile constant initializers that + contain floating-point operations. Pacify this compiler. */ +#ifdef _MSC_VER +# pragma fenv_access (off) +#endif + /* 2^(MANT_DIG-1). */ -static const double TWO_MANT_DIG = +static const DOUBLE TWO_MANT_DIG = /* Assume MANT_DIG <= 5 * 31. Use the identity n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5). */ @@ -65,21 +85,37 @@ FUNC (DOUBLE x) if (z > L_(0.0)) { - /* Round to the next integer (nearest or up or down, doesn't matter). */ - z += TWO_MANT_DIG; - z -= TWO_MANT_DIG; - /* Enforce rounding down. */ - if (z > y) - z -= L_(1.0); + /* For 0 < x < 1, return +0.0 even if the current rounding mode is + FE_DOWNWARD. */ + if (z < L_(1.0)) + z = L_(0.0); + /* 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 down. */ + if (z > y) + z -= L_(1.0); + } } else if (z < L_(0.0)) { - /* 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); + /* 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); + } } return z; }