Commit | Line | Data |
---|---|---|
49114fd4 | 1 | /* Round towards zero. |
5e69ceb7 | 2 | Copyright (C) 2007, 2010-2014 Free Software Foundation, Inc. |
49114fd4 LC |
3 | |
4 | This program is free software: you can redistribute it and/or modify | |
5 | it under the terms of the GNU Lesser General Public License as published by | |
6 | the Free Software Foundation; either version 3 of the License, or | |
7 | (at your option) any later version. | |
8 | ||
9 | This program is distributed in the hope that it will be useful, | |
10 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
11 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
12 | GNU Lesser General Public License for more details. | |
13 | ||
14 | You should have received a copy of the GNU Lesser General Public License | |
15 | along with this program. If not, see <http://www.gnu.org/licenses/>. */ | |
16 | ||
17 | /* Written by Bruno Haible <bruno@clisp.org>, 2007. */ | |
18 | ||
35428fb6 LC |
19 | #if ! defined USE_LONG_DOUBLE |
20 | # include <config.h> | |
21 | #endif | |
49114fd4 LC |
22 | |
23 | /* Specification. */ | |
24 | #include <math.h> | |
25 | ||
26 | #include <float.h> | |
27 | ||
28 | #undef MIN | |
29 | ||
30 | #ifdef USE_LONG_DOUBLE | |
31 | # define FUNC truncl | |
32 | # define DOUBLE long double | |
33 | # define MANT_DIG LDBL_MANT_DIG | |
34 | # define MIN LDBL_MIN | |
35 | # define L_(literal) literal##L | |
36 | #elif ! defined USE_FLOAT | |
37 | # define FUNC trunc | |
38 | # define DOUBLE double | |
39 | # define MANT_DIG DBL_MANT_DIG | |
40 | # define MIN DBL_MIN | |
41 | # define L_(literal) literal | |
42 | #else /* defined USE_FLOAT */ | |
43 | # define FUNC truncf | |
44 | # define DOUBLE float | |
45 | # define MANT_DIG FLT_MANT_DIG | |
46 | # define MIN FLT_MIN | |
47 | # define L_(literal) literal##f | |
48 | #endif | |
49 | ||
50 | /* -0.0. See minus-zero.h. */ | |
51 | #if defined __hpux || defined __sgi || defined __ICC | |
52 | # define MINUS_ZERO (-MIN * MIN) | |
53 | #else | |
54 | # define MINUS_ZERO L_(-0.0) | |
55 | #endif | |
56 | ||
005de2e8 LC |
57 | /* MSVC with option -fp:strict refuses to compile constant initializers that |
58 | contain floating-point operations. Pacify this compiler. */ | |
59 | #ifdef _MSC_VER | |
60 | # pragma fenv_access (off) | |
61 | #endif | |
62 | ||
49114fd4 LC |
63 | /* 2^(MANT_DIG-1). */ |
64 | static const DOUBLE TWO_MANT_DIG = | |
65 | /* Assume MANT_DIG <= 5 * 31. | |
66 | Use the identity | |
67 | n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5). */ | |
68 | (DOUBLE) (1U << ((MANT_DIG - 1) / 5)) | |
69 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5)) | |
70 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5)) | |
71 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5)) | |
72 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5)); | |
73 | ||
74 | DOUBLE | |
75 | FUNC (DOUBLE x) | |
76 | { | |
77 | /* The use of 'volatile' guarantees that excess precision bits are dropped | |
78 | at each addition step and before the following comparison at the caller's | |
79 | site. It is necessary on x86 systems where double-floats are not IEEE | |
80 | compliant by default, to avoid that the results become platform and compiler | |
81 | option dependent. 'volatile' is a portable alternative to gcc's | |
82 | -ffloat-store option. */ | |
83 | volatile DOUBLE y = x; | |
84 | volatile DOUBLE z = y; | |
85 | ||
86 | if (z > L_(0.0)) | |
87 | { | |
88 | /* For 0 < x < 1, return +0.0 even if the current rounding mode is | |
89 | FE_DOWNWARD. */ | |
90 | if (z < L_(1.0)) | |
91 | z = L_(0.0); | |
92 | /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */ | |
93 | else if (z < TWO_MANT_DIG) | |
94 | { | |
95 | /* Round to the next integer (nearest or up or down, doesn't matter). */ | |
96 | z += TWO_MANT_DIG; | |
97 | z -= TWO_MANT_DIG; | |
98 | /* Enforce rounding down. */ | |
99 | if (z > y) | |
100 | z -= L_(1.0); | |
101 | } | |
102 | } | |
103 | else if (z < L_(0.0)) | |
104 | { | |
105 | /* For -1 < x < 0, return -0.0 regardless of the current rounding | |
106 | mode. */ | |
107 | if (z > L_(-1.0)) | |
108 | z = MINUS_ZERO; | |
109 | /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */ | |
110 | else if (z > - TWO_MANT_DIG) | |
111 | { | |
112 | /* Round to the next integer (nearest or up or down, doesn't matter). */ | |
113 | z -= TWO_MANT_DIG; | |
114 | z += TWO_MANT_DIG; | |
115 | /* Enforce rounding up. */ | |
116 | if (z < y) | |
117 | z += L_(1.0); | |
118 | } | |
119 | } | |
120 | return z; | |
121 | } |