Commit | Line | Data |
---|---|---|
49114fd4 LC |
1 | /* Round towards zero. |
2 | Copyright (C) 2007, 2010-2011 Free Software Foundation, Inc. | |
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 | ||
19 | #include <config.h> | |
20 | ||
21 | /* Specification. */ | |
22 | #include <math.h> | |
23 | ||
24 | #include <float.h> | |
25 | ||
26 | #undef MIN | |
27 | ||
28 | #ifdef USE_LONG_DOUBLE | |
29 | # define FUNC truncl | |
30 | # define DOUBLE long double | |
31 | # define MANT_DIG LDBL_MANT_DIG | |
32 | # define MIN LDBL_MIN | |
33 | # define L_(literal) literal##L | |
34 | #elif ! defined USE_FLOAT | |
35 | # define FUNC trunc | |
36 | # define DOUBLE double | |
37 | # define MANT_DIG DBL_MANT_DIG | |
38 | # define MIN DBL_MIN | |
39 | # define L_(literal) literal | |
40 | #else /* defined USE_FLOAT */ | |
41 | # define FUNC truncf | |
42 | # define DOUBLE float | |
43 | # define MANT_DIG FLT_MANT_DIG | |
44 | # define MIN FLT_MIN | |
45 | # define L_(literal) literal##f | |
46 | #endif | |
47 | ||
48 | /* -0.0. See minus-zero.h. */ | |
49 | #if defined __hpux || defined __sgi || defined __ICC | |
50 | # define MINUS_ZERO (-MIN * MIN) | |
51 | #else | |
52 | # define MINUS_ZERO L_(-0.0) | |
53 | #endif | |
54 | ||
55 | /* 2^(MANT_DIG-1). */ | |
56 | static const DOUBLE TWO_MANT_DIG = | |
57 | /* Assume MANT_DIG <= 5 * 31. | |
58 | Use the identity | |
59 | n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5). */ | |
60 | (DOUBLE) (1U << ((MANT_DIG - 1) / 5)) | |
61 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5)) | |
62 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5)) | |
63 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5)) | |
64 | * (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5)); | |
65 | ||
66 | DOUBLE | |
67 | FUNC (DOUBLE x) | |
68 | { | |
69 | /* The use of 'volatile' guarantees that excess precision bits are dropped | |
70 | at each addition step and before the following comparison at the caller's | |
71 | site. It is necessary on x86 systems where double-floats are not IEEE | |
72 | compliant by default, to avoid that the results become platform and compiler | |
73 | option dependent. 'volatile' is a portable alternative to gcc's | |
74 | -ffloat-store option. */ | |
75 | volatile DOUBLE y = x; | |
76 | volatile DOUBLE z = y; | |
77 | ||
78 | if (z > L_(0.0)) | |
79 | { | |
80 | /* For 0 < x < 1, return +0.0 even if the current rounding mode is | |
81 | FE_DOWNWARD. */ | |
82 | if (z < L_(1.0)) | |
83 | z = L_(0.0); | |
84 | /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */ | |
85 | else if (z < TWO_MANT_DIG) | |
86 | { | |
87 | /* Round to the next integer (nearest or up or down, doesn't matter). */ | |
88 | z += TWO_MANT_DIG; | |
89 | z -= TWO_MANT_DIG; | |
90 | /* Enforce rounding down. */ | |
91 | if (z > y) | |
92 | z -= L_(1.0); | |
93 | } | |
94 | } | |
95 | else if (z < L_(0.0)) | |
96 | { | |
97 | /* For -1 < x < 0, return -0.0 regardless of the current rounding | |
98 | mode. */ | |
99 | if (z > L_(-1.0)) | |
100 | z = MINUS_ZERO; | |
101 | /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */ | |
102 | else if (z > - TWO_MANT_DIG) | |
103 | { | |
104 | /* Round to the next integer (nearest or up or down, doesn't matter). */ | |
105 | z -= TWO_MANT_DIG; | |
106 | z += TWO_MANT_DIG; | |
107 | /* Enforce rounding up. */ | |
108 | if (z < y) | |
109 | z += L_(1.0); | |
110 | } | |
111 | } | |
112 | return z; | |
113 | } |