Mercurial > hg > octave-lojdl > gnulib-hg
annotate lib/trunc.c @ 16511:2074f2bf7216
math code: Add comments.
* lib/acosl.c: Add comment about related glibc source files.
* lib/asinl.c: Likewise.
* lib/atanl.c: Likewise.
* lib/expl.c: Likewise.
* lib/logl.c: Likewise.
* lib/sincosl.c: Likewise.
* lib/sinl.c: Likewise.
* lib/tanl.c: Likewise.
* lib/trigl.c: Likewise.
* lib/cosl.c: Likewise. Fix comments.
| author | Bruno Haible <bruno@clisp.org> |
|---|---|
| date | Wed, 29 Feb 2012 12:30:07 +0100 |
| parents | fa4e9b981eb4 |
| children | e542fd46ad6f |
| rev | line source |
|---|---|
| 9282 | 1 /* Round towards zero. |
|
16201
8250f2777afc
maint: update all copyright year number ranges
Jim Meyering <meyering@redhat.com>
parents:
15926
diff
changeset
|
2 Copyright (C) 2007, 2010-2012 Free Software Foundation, Inc. |
| 9282 | 3 |
|
9309
bbbbbf4cd1c5
Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents:
9285
diff
changeset
|
4 This program is free software: you can redistribute it and/or modify |
| 9282 | 5 it under the terms of the GNU General Public License as published by |
|
9309
bbbbbf4cd1c5
Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents:
9285
diff
changeset
|
6 the Free Software Foundation; either version 3 of the License, or |
|
bbbbbf4cd1c5
Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents:
9285
diff
changeset
|
7 (at your option) any later version. |
| 9282 | 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 General Public License for more details. | |
| 13 | |
|
9309
bbbbbf4cd1c5
Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents:
9285
diff
changeset
|
14 You should have received a copy of the GNU General Public License |
|
bbbbbf4cd1c5
Change copyright notice from GPLv2+ to GPLv3+.
Bruno Haible <bruno@clisp.org>
parents:
9285
diff
changeset
|
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */ |
| 9282 | 16 |
| 17 /* Written by Bruno Haible <bruno@clisp.org>, 2007. */ | |
| 18 | |
|
15926
23d150bba2eb
truncl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents:
14079
diff
changeset
|
19 #if ! defined USE_LONG_DOUBLE |
|
23d150bba2eb
truncl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents:
14079
diff
changeset
|
20 # include <config.h> |
|
23d150bba2eb
truncl: Simplify for platforms where 'long double' == 'double'.
Bruno Haible <bruno@clisp.org>
parents:
14079
diff
changeset
|
21 #endif |
| 9282 | 22 |
| 23 /* Specification. */ | |
| 24 #include <math.h> | |
| 25 | |
| 26 #include <float.h> | |
| 27 | |
|
14018
ac9f59901721
ceil, trunc, round: Fix gcc warnings.
Bruno Haible <bruno@clisp.org>
parents:
13993
diff
changeset
|
28 #undef MIN |
|
ac9f59901721
ceil, trunc, round: Fix gcc warnings.
Bruno Haible <bruno@clisp.org>
parents:
13993
diff
changeset
|
29 |
| 9285 | 30 #ifdef USE_LONG_DOUBLE |
| 31 # define FUNC truncl | |
| 32 # define DOUBLE long double | |
| 33 # define MANT_DIG LDBL_MANT_DIG | |
|
13993
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
34 # define MIN LDBL_MIN |
| 9285 | 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 | |
|
13993
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
40 # define MIN DBL_MIN |
| 9285 | 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 | |
|
13993
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
46 # define MIN FLT_MIN |
| 9285 | 47 # define L_(literal) literal##f |
| 48 #endif | |
| 49 | |
|
13993
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
50 /* -0.0. See minus-zero.h. */ |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
51 #if defined __hpux || defined __sgi || defined __ICC |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
52 # define MINUS_ZERO (-MIN * MIN) |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
53 #else |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
54 # define MINUS_ZERO L_(-0.0) |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
55 #endif |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
56 |
|
16506
fa4e9b981eb4
Avoid compilation errors with MSVC option -fp:strict.
Bruno Haible <bruno@clisp.org>
parents:
16201
diff
changeset
|
57 /* MSVC with option -fp:strict refuses to compile constant initializers that |
|
fa4e9b981eb4
Avoid compilation errors with MSVC option -fp:strict.
Bruno Haible <bruno@clisp.org>
parents:
16201
diff
changeset
|
58 contain floating-point operations. Pacify this compiler. */ |
|
fa4e9b981eb4
Avoid compilation errors with MSVC option -fp:strict.
Bruno Haible <bruno@clisp.org>
parents:
16201
diff
changeset
|
59 #ifdef _MSC_VER |
|
fa4e9b981eb4
Avoid compilation errors with MSVC option -fp:strict.
Bruno Haible <bruno@clisp.org>
parents:
16201
diff
changeset
|
60 # pragma fenv_access (off) |
|
fa4e9b981eb4
Avoid compilation errors with MSVC option -fp:strict.
Bruno Haible <bruno@clisp.org>
parents:
16201
diff
changeset
|
61 #endif |
|
fa4e9b981eb4
Avoid compilation errors with MSVC option -fp:strict.
Bruno Haible <bruno@clisp.org>
parents:
16201
diff
changeset
|
62 |
| 9285 | 63 /* 2^(MANT_DIG-1). */ |
| 9313 | 64 static const DOUBLE TWO_MANT_DIG = |
| 9285 | 65 /* Assume MANT_DIG <= 5 * 31. |
| 9282 | 66 Use the identity |
| 9285 | 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)); | |
| 9282 | 73 |
| 9285 | 74 DOUBLE |
| 75 FUNC (DOUBLE x) | |
| 9282 | 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. */ | |
| 9285 | 83 volatile DOUBLE y = x; |
| 84 volatile DOUBLE z = y; | |
| 9282 | 85 |
| 9285 | 86 if (z > L_(0.0)) |
| 9282 | 87 { |
|
13993
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
88 /* For 0 < x < 1, return +0.0 even if the current rounding mode is |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
89 FE_DOWNWARD. */ |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
90 if (z < L_(1.0)) |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
91 z = L_(0.0); |
|
9327
5437ef3873b5
Fix incorrect rounding of trunc, truncf, truncl in some cases. Add a new test.
Bruno Haible <bruno@clisp.org>
parents:
9313
diff
changeset
|
92 /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */ |
|
13993
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
93 else if (z < TWO_MANT_DIG) |
|
12421
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
94 { |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
95 /* Round to the next integer (nearest or up or down, doesn't matter). */ |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
96 z += TWO_MANT_DIG; |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
97 z -= TWO_MANT_DIG; |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
98 /* Enforce rounding down. */ |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
99 if (z > y) |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
100 z -= L_(1.0); |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
101 } |
| 9282 | 102 } |
| 9285 | 103 else if (z < L_(0.0)) |
| 9282 | 104 { |
|
13993
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
105 /* For -1 < x < 0, return -0.0 regardless of the current rounding |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
106 mode. */ |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
107 if (z > L_(-1.0)) |
|
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
108 z = MINUS_ZERO; |
|
9327
5437ef3873b5
Fix incorrect rounding of trunc, truncf, truncl in some cases. Add a new test.
Bruno Haible <bruno@clisp.org>
parents:
9313
diff
changeset
|
109 /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */ |
|
13993
9e9aaf4c9464
trunc: Implement result sign according to IEEE 754.
Bruno Haible <bruno@clisp.org>
parents:
12559
diff
changeset
|
110 else if (z > - TWO_MANT_DIG) |
|
12421
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
111 { |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
112 /* Round to the next integer (nearest or up or down, doesn't matter). */ |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
113 z -= TWO_MANT_DIG; |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
114 z += TWO_MANT_DIG; |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
115 /* Enforce rounding up. */ |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
116 if (z < y) |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
117 z += L_(1.0); |
|
e8d2c6fc33ad
Use spaces for indentation, not tabs.
Bruno Haible <bruno@clisp.org>
parents:
9327
diff
changeset
|
118 } |
| 9282 | 119 } |
| 120 return z; | |
| 121 } |