Commit | Line | Data |
---|---|---|
005de2e8 | 1 | /* Round toward nearest, breaking ties away from zero. |
5e69ceb7 | 2 | Copyright (C) 2007, 2010-2014 Free Software Foundation, Inc. |
005de2e8 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 2, or (at your option) | |
7 | 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 along | |
15 | with this program; if not, see <http://www.gnu.org/licenses/>. */ | |
16 | ||
17 | /* Written by Ben Pfaff <blp@gnu.org>, 2007. | |
18 | Based heavily on code by Bruno Haible. */ | |
19 | ||
20 | #if ! defined USE_LONG_DOUBLE | |
21 | # include <config.h> | |
22 | #endif | |
23 | ||
24 | /* Specification. */ | |
25 | #include <math.h> | |
26 | ||
27 | #include <float.h> | |
28 | ||
29 | #undef MIN | |
30 | ||
31 | #ifdef USE_LONG_DOUBLE | |
32 | # define ROUND roundl | |
33 | # define FLOOR floorl | |
34 | # define CEIL ceill | |
35 | # define DOUBLE long double | |
36 | # define MANT_DIG LDBL_MANT_DIG | |
37 | # define MIN LDBL_MIN | |
38 | # define L_(literal) literal##L | |
39 | # define HAVE_FLOOR_AND_CEIL HAVE_FLOORL_AND_CEILL | |
40 | #elif ! defined USE_FLOAT | |
41 | # define ROUND round | |
42 | # define FLOOR floor | |
43 | # define CEIL ceil | |
44 | # define DOUBLE double | |
45 | # define MANT_DIG DBL_MANT_DIG | |
46 | # define MIN DBL_MIN | |
47 | # define L_(literal) literal | |
48 | # define HAVE_FLOOR_AND_CEIL 1 | |
49 | #else /* defined USE_FLOAT */ | |
50 | # define ROUND roundf | |
51 | # define FLOOR floorf | |
52 | # define CEIL ceilf | |
53 | # define DOUBLE float | |
54 | # define MANT_DIG FLT_MANT_DIG | |
55 | # define MIN FLT_MIN | |
56 | # define L_(literal) literal##f | |
57 | # define HAVE_FLOOR_AND_CEIL HAVE_FLOORF_AND_CEILF | |
58 | #endif | |
59 | ||
60 | /* -0.0. See minus-zero.h. */ | |
61 | #if defined __hpux || defined __sgi || defined __ICC | |
62 | # define MINUS_ZERO (-MIN * MIN) | |
63 | #else | |
64 | # define MINUS_ZERO L_(-0.0) | |
65 | #endif | |
66 | ||
67 | /* MSVC with option -fp:strict refuses to compile constant initializers that | |
68 | contain floating-point operations. Pacify this compiler. */ | |
69 | #ifdef _MSC_VER | |
70 | # pragma fenv_access (off) | |
71 | #endif | |
72 | ||
73 | /* If we're being included from test-round2[f].c, it already defined names for | |
74 | our round implementations. Otherwise, pick the preferred implementation for | |
75 | this machine. */ | |
76 | #if !defined FLOOR_BASED_ROUND && !defined FLOOR_FREE_ROUND | |
77 | # if HAVE_FLOOR_AND_CEIL | |
78 | # define FLOOR_BASED_ROUND ROUND | |
79 | # else | |
80 | # define FLOOR_FREE_ROUND ROUND | |
81 | # endif | |
82 | #endif | |
83 | ||
84 | #ifdef FLOOR_BASED_ROUND | |
85 | /* An implementation of the C99 round function based on floor and ceil. We use | |
86 | this when floor and ceil are available, on the assumption that they are | |
87 | faster than the open-coded versions below. */ | |
88 | DOUBLE | |
89 | FLOOR_BASED_ROUND (DOUBLE x) | |
90 | { | |
91 | if (x >= L_(0.0)) | |
92 | { | |
93 | DOUBLE y = FLOOR (x); | |
94 | if (x - y >= L_(0.5)) | |
95 | y += L_(1.0); | |
96 | return y; | |
97 | } | |
98 | else | |
99 | { | |
100 | DOUBLE y = CEIL (x); | |
101 | if (y - x >= L_(0.5)) | |
102 | y -= L_(1.0); | |
103 | return y; | |
104 | } | |
105 | } | |
106 | #endif /* FLOOR_BASED_ROUND */ | |
107 | ||
108 | #ifdef FLOOR_FREE_ROUND | |
109 | /* An implementation of the C99 round function without floor or ceil. | |
110 | We use this when floor or ceil is missing. */ | |
111 | DOUBLE | |
112 | FLOOR_FREE_ROUND (DOUBLE x) | |
113 | { | |
114 | /* 2^(MANT_DIG-1). */ | |
115 | static const DOUBLE TWO_MANT_DIG = | |
116 | /* Assume MANT_DIG <= 5 * 31. | |
117 | Use the identity | |
118 | n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5). */ | |
119 | (DOUBLE) (1U << ((MANT_DIG - 1) / 5)) | |
120 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5)) | |
121 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5)) | |
122 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5)) | |
123 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5)); | |
124 | ||
125 | /* The use of 'volatile' guarantees that excess precision bits are dropped at | |
126 | each addition step and before the following comparison at the caller's | |
127 | site. It is necessary on x86 systems where double-floats are not IEEE | |
128 | compliant by default, to avoid that the results become platform and | |
129 | compiler option dependent. 'volatile' is a portable alternative to gcc's | |
130 | -ffloat-store option. */ | |
131 | volatile DOUBLE y = x; | |
132 | volatile DOUBLE z = y; | |
133 | ||
134 | if (z > L_(0.0)) | |
135 | { | |
136 | /* Avoid rounding error for x = 0.5 - 2^(-MANT_DIG-1). */ | |
137 | if (z < L_(0.5)) | |
138 | z = L_(0.0); | |
139 | /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */ | |
140 | else if (z < TWO_MANT_DIG) | |
141 | { | |
142 | /* Add 0.5 to the absolute value. */ | |
143 | y = z += L_(0.5); | |
144 | /* Round to the next integer (nearest or up or down, doesn't | |
145 | matter). */ | |
146 | z += TWO_MANT_DIG; | |
147 | z -= TWO_MANT_DIG; | |
148 | /* Enforce rounding down. */ | |
149 | if (z > y) | |
150 | z -= L_(1.0); | |
151 | } | |
152 | } | |
153 | else if (z < L_(0.0)) | |
154 | { | |
155 | /* Avoid rounding error for x = -(0.5 - 2^(-MANT_DIG-1)). */ | |
156 | if (z > - L_(0.5)) | |
157 | z = MINUS_ZERO; | |
158 | /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */ | |
159 | else if (z > -TWO_MANT_DIG) | |
160 | { | |
161 | /* Add 0.5 to the absolute value. */ | |
162 | y = z -= L_(0.5); | |
163 | /* Round to the next integer (nearest or up or down, doesn't | |
164 | matter). */ | |
165 | z -= TWO_MANT_DIG; | |
166 | z += TWO_MANT_DIG; | |
167 | /* Enforce rounding up. */ | |
168 | if (z < y) | |
169 | z += L_(1.0); | |
170 | } | |
171 | } | |
172 | return z; | |
173 | } | |
174 | #endif /* FLOOR_FREE_ROUND */ |