1 /* Test of rounding towards zero.
2 Copyright (C) 2007-2022 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 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 General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <https://www.gnu.org/licenses/>. */
16
17 /* Written by Bruno Haible <bruno@clisp.org>, 2007. */
18
19 /* When this test fails on some platform, build it together with the gnulib
20 module 'fprintf-posix' for optimal debugging output. */
21
22 #include <config.h>
23
24 #include <math.h>
25
26 #include <float.h>
27 #include <stdbool.h>
28 #include <stdint.h>
29 #include <stdio.h>
30
31 #include "isnand-nolibm.h"
32 #include "minus-zero.h"
33 #include "macros.h"
34
35 /* MSVC with option -fp:strict refuses to compile constant initializers that
36 contain floating-point operations. Pacify this compiler. */
37 #if defined _MSC_VER && !defined __clang__
38 # pragma fenv_access (off)
39 #endif
40
41
42 /* The reference implementation, taken from lib/trunc.c. */
43
44 #define DOUBLE double
45 #define MANT_DIG DBL_MANT_DIG
46 #define L_(literal) literal
47
48 /* -0.0. See minus-zero.h. */
49 #define MINUS_ZERO minus_zerod
50
51 /* 2^(MANT_DIG-1). */
52 static const DOUBLE TWO_MANT_DIG =
53 /* Assume MANT_DIG <= 5 * 31.
54 Use the identity
55 n = floor(n/5) + floor((n+1)/5) + ... + floor((n+4)/5). */
56 (DOUBLE) (1U << ((MANT_DIG - 1) / 5))
57 * (DOUBLE) (1U << ((MANT_DIG - 1 + 1) / 5))
58 * (DOUBLE) (1U << ((MANT_DIG - 1 + 2) / 5))
59 * (DOUBLE) (1U << ((MANT_DIG - 1 + 3) / 5))
60 * (DOUBLE) (1U << ((MANT_DIG - 1 + 4) / 5));
61
62 DOUBLE
63 trunc_reference (DOUBLE x)
64 {
65 /* The use of 'volatile' guarantees that excess precision bits are dropped
66 at each addition step and before the following comparison at the caller's
67 site. It is necessary on x86 systems where double-floats are not IEEE
68 compliant by default, to avoid that the results become platform and compiler
69 option dependent. 'volatile' is a portable alternative to gcc's
70 -ffloat-store option. */
71 volatile DOUBLE y = x;
72 volatile DOUBLE z = y;
73
74 if (z > L_(0.0))
75 {
76 /* For 0 < x < 1, return +0.0 even if the current rounding mode is
77 FE_DOWNWARD. */
78 if (z < L_(1.0))
79 z = L_(0.0);
80 /* Avoid rounding errors for values near 2^k, where k >= MANT_DIG-1. */
81 else if (z < TWO_MANT_DIG)
82 {
83 /* Round to the next integer (nearest or up or down, doesn't matter). */
84 z += TWO_MANT_DIG;
85 z -= TWO_MANT_DIG;
86 /* Enforce rounding down. */
87 if (z > y)
88 z -= L_(1.0);
89 }
90 }
91 else if (z < L_(0.0))
92 {
93 /* For -1 < x < 0, return -0.0 regardless of the current rounding
94 mode. */
95 if (z > L_(-1.0))
96 z = MINUS_ZERO;
97 /* Avoid rounding errors for values near -2^k, where k >= MANT_DIG-1. */
98 else if (z > - TWO_MANT_DIG)
99 {
100 /* Round to the next integer (nearest or up or down, doesn't matter). */
101 z -= TWO_MANT_DIG;
102 z += TWO_MANT_DIG;
103 /* Enforce rounding up. */
104 if (z < y)
105 z += L_(1.0);
106 }
107 }
108 return z;
109 }
110
111
112 /* Test for equality. */
113 static int
114 equal (DOUBLE x, DOUBLE y)
115 {
116 return (isnand (x) ? isnand (y) : x == y);
117 }
118
119 /* Test whether the result for a given argument is correct. */
120 static bool
121 correct_result_p (DOUBLE x, DOUBLE result)
122 {
123 return
124 (x >= 0
125 ? (x < 1 ? result == L_(0.0) :
126 x - 1 < x ? result <= x && result >= x - 1 && x - result < 1 :
127 equal (result, x))
128 : (x > -1 ? result == L_(0.0) :
129 x + 1 > x ? result >= x && result <= x + 1 && result - x < 1 :
130 equal (result, x)));
131 }
132
133 /* Test the function for a given argument. */
134 static int
135 check (double x)
136 {
137 /* If the reference implementation is incorrect, bail out immediately. */
138 double reference = trunc_reference (x);
139 ASSERT (correct_result_p (x, reference));
140 /* If the actual implementation is wrong, return an error code. */
141 {
142 double result = trunc (x);
143 if (correct_result_p (x, result))
144 return 0;
145 else
146 {
147 #if GNULIB_TEST_FPRINTF_POSIX
148 fprintf (stderr, "trunc %g(%a) = %g(%a) or %g(%a)?\n",
149 x, x, reference, reference, result, result);
150 #endif
151 return 1;
152 }
153 }
154 }
155
156 #define NUM_HIGHBITS 13
157 #define NUM_LOWBITS 4
158
159 int
160 main ()
161 {
162 unsigned int highbits;
163 unsigned int lowbits;
164 int error = 0;
165 for (highbits = 0; highbits < (1 << NUM_HIGHBITS); highbits++)
166 for (lowbits = 0; lowbits < (1 << NUM_LOWBITS); lowbits++)
167 {
168 /* Combine highbits and lowbits into a floating-point number,
169 sign-extending the lowbits to 32-NUM_HIGHBITS bits. */
170 union { double f; uint64_t i; } janus;
171 janus.i = ((uint64_t) highbits << (64 - NUM_HIGHBITS))
172 | ((uint64_t) ((int64_t) ((uint64_t) lowbits << (64 - NUM_LOWBITS))
173 >> (64 - NUM_LOWBITS - NUM_HIGHBITS))
174 >> NUM_HIGHBITS);
175 error |= check (janus.f);
176 }
177 return (error ? 1 : 0);
178 }