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