1 /* $NetBSD: t_fe_round.c,v 1.20 2024/05/15 00:02:57 riastradh Exp $ */
2
3 /*
4 * Written by Maya Rashish <maya@NetBSD.org>
5 * Public domain.
6 *
7 * Testing IEEE-754 rounding modes (and lrint)
8 */
9
10 #include <sys/cdefs.h>
11 __RCSID("$NetBSD: t_fe_round.c,v 1.20 2024/05/15 00:02:57 riastradh Exp $");
12
13 #include <atf-c.h>
14 #include <fenv.h>
15 #include <math.h>
16 #include <stdint.h>
17 #include <stdio.h>
18 #include <stdlib.h>
19
20 #ifdef __HAVE_FENV
21
22 /*#pragma STDC FENV_ACCESS ON gcc?? */
23
24 #define INT 9223L
25
26 static const char *
rmname(int rm)27 rmname(int rm)
28 {
29 switch (rm) {
30 case FE_TOWARDZERO:
31 return "FE_TOWARDZERO";
32 case FE_DOWNWARD:
33 return "FE_DOWNWARD";
34 case FE_UPWARD:
35 return "FE_UPWARD";
36 case FE_TONEAREST:
37 return "FE_TONEAREST";
38 default:
39 return "unknown";
40 }
41 }
42
43 /*
44 * Examples are chosen to fit within the smallest single-precision
45 * format any NetBSD port uses, so that we can write the examples once
46 * in type double, and convert to single without raising inexact-result
47 * exceptions when we're trying to test whether the integer-rounding
48 * functions raise them.
49 */
50 static const struct {
51 int round_mode;
52 double input;
53 long int expected;
54 } values[] = {
55 { FE_DOWNWARD, 3.75, 3},
56 { FE_DOWNWARD, -3.75, -4},
57 { FE_DOWNWARD, +0., 0},
58 { FE_DOWNWARD, -0., 0},
59 { FE_DOWNWARD, -INT-0.0625, -INT-1},
60 { FE_DOWNWARD, +INT-0.0625, INT-1},
61 { FE_DOWNWARD, -INT+0.0625, -INT},
62 { FE_DOWNWARD, +INT+0.0625, INT},
63
64 { FE_UPWARD, +0., 0},
65 { FE_UPWARD, -0., 0},
66 { FE_UPWARD, -123.75, -123},
67 { FE_UPWARD, 123.75, 124},
68 { FE_UPWARD, -INT-0.0625, -INT},
69 { FE_UPWARD, +INT-0.0625, INT},
70 { FE_UPWARD, -INT+0.0625, -INT+1},
71 { FE_UPWARD, +INT+0.0625, INT+1},
72
73 { FE_TOWARDZERO, 1.9375, 1},
74 { FE_TOWARDZERO, -1.9375, -1},
75 { FE_TOWARDZERO, 0.25, 0},
76 { FE_TOWARDZERO, INT+0.0625, INT},
77 { FE_TOWARDZERO, INT-0.0625, INT - 1},
78 { FE_TOWARDZERO, -INT+0.0625, -INT + 1},
79 { FE_TOWARDZERO, +0., 0},
80 { FE_TOWARDZERO, -0., 0},
81
82 { FE_TONEAREST, -INT-0.0625, -INT},
83 { FE_TONEAREST, +INT-0.0625, INT},
84 { FE_TONEAREST, -INT+0.0625, -INT},
85 { FE_TONEAREST, +INT+0.0625, INT},
86 { FE_TONEAREST, -INT-0.53125, -INT-1},
87 { FE_TONEAREST, +INT-0.53125, INT-1},
88 { FE_TONEAREST, -INT+0.53125, -INT+1},
89 { FE_TONEAREST, +INT+0.53125, INT+1},
90 { FE_TONEAREST, +0., 0},
91 { FE_TONEAREST, -0., 0},
92 };
93
94 ATF_TC(fe_lrint);
ATF_TC_HEAD(fe_lrint,tc)95 ATF_TC_HEAD(fe_lrint, tc)
96 {
97 atf_tc_set_md_var(tc, "descr",
98 "Checking IEEE 754 rounding modes using lrint(3)");
99 }
100
ATF_TC_BODY(fe_lrint,tc)101 ATF_TC_BODY(fe_lrint, tc)
102 {
103 enum {
104 LLRINT,
105 LLRINTF,
106 LRINT,
107 LRINTF,
108 N_FN,
109 } fn;
110 static const char *const fnname[] = {
111 [LLRINT] = "llrint",
112 [LLRINTF] = "llrintf",
113 [LRINT] = "lrint",
114 [LRINTF] = "lrintf",
115 };
116 long int received;
117 unsigned i;
118
119 for (i = 0; i < __arraycount(values); i++) {
120 for (fn = 0; fn < N_FN; fn++) {
121 /*
122 * Set the requested rounding mode.
123 */
124 fesetround(values[i].round_mode);
125
126 /*
127 * Call the lrint(3)-family function.
128 */
129 switch (fn) {
130 case LLRINT:
131 received = llrint(values[i].input);
132 break;
133 case LLRINTF:
134 received = llrintf(values[i].input);
135 break;
136 case LRINT:
137 received = lrint(values[i].input);
138 break;
139 case LRINTF:
140 received = lrintf(values[i].input);
141 break;
142 default:
143 atf_tc_fail("impossible");
144 }
145
146 /*
147 * Assuming the result we got has zero
148 * fractional part, casting to long int should
149 * have no rounding. Verify it matches the
150 * integer we expect.
151 */
152 ATF_CHECK_MSG((long int)received == values[i].expected,
153 "[%u] %s %s(%f): got %ld, expected %ld",
154 i, rmname(values[i].round_mode), fnname[fn],
155 values[i].input,
156 (long int)received, values[i].expected);
157
158 /* Do we get the same rounding mode out? */
159 ATF_CHECK_MSG(fegetround() == values[i].round_mode,
160 "[%u] %s: set %d (%s), got %d (%s)",
161 i, fnname[fn],
162 values[i].round_mode, rmname(values[i].round_mode),
163 fegetround(), rmname(fegetround()));
164 }
165 }
166 }
167
168 ATF_TC(fe_nearbyint_rint);
ATF_TC_HEAD(fe_nearbyint_rint,tc)169 ATF_TC_HEAD(fe_nearbyint_rint, tc)
170 {
171 atf_tc_set_md_var(tc, "descr",
172 "Checking IEEE 754 rounding modes using nearbyint/rint");
173 }
174
ATF_TC_BODY(fe_nearbyint_rint,tc)175 ATF_TC_BODY(fe_nearbyint_rint, tc)
176 {
177 enum {
178 NEARBYINT,
179 NEARBYINTF,
180 NEARBYINTL,
181 RINT,
182 RINTF,
183 RINTL,
184 N_FN,
185 } fn;
186 static const char *const fnname[] = {
187 [NEARBYINT] = "nearbyint",
188 [NEARBYINTF] = "nearbyintf",
189 [NEARBYINTL] = "nearbyintl",
190 [RINT] = "rint",
191 [RINTF] = "rintf",
192 [RINTL] = "rintl",
193 };
194 double received, ipart, fpart;
195 unsigned i;
196
197 for (i = 0; i < __arraycount(values); i++) {
198 for (fn = 0; fn < N_FN; fn++) {
199 bool expect_except =
200 values[i].input != (double)values[i].expected;
201
202 /*
203 * Set the requested rounding mode.
204 */
205 fesetround(values[i].round_mode);
206
207 #ifdef __ia64__
208 /*
209 * Without this barrier, we get:
210 *
211 * /tmp//ccJayu9g.s:2793: Warning: Use of 'mov.m' violates RAW dependency 'AR[FPSR].sf0.flags' (impliedf)
212 * /tmp//ccJayu9g.s:2793: Warning: Only the first path encountering the conflict is reported
213 * /tmp//ccJayu9g.s:2757: Warning: This is the location of the conflicting usage
214 *
215 * (If you fix this, remove the entry from doc/HACKS.)
216 */
217 __insn_barrier();
218 #endif
219
220 /*
221 * Clear sticky floating-point exception bits
222 * so we can verify whether the FE_INEXACT
223 * exception is raised.
224 */
225 feclearexcept(FE_ALL_EXCEPT);
226
227 /*
228 * Call the rint(3)-family function.
229 */
230 switch (fn) {
231 case NEARBYINT:
232 received = nearbyint(values[i].input);
233 expect_except = false;
234 break;
235 case NEARBYINTF:
236 received = nearbyintf(values[i].input);
237 expect_except = false;
238 break;
239 case NEARBYINTL:
240 received = nearbyintl(values[i].input);
241 expect_except = false;
242 break;
243 case RINT:
244 received = rint(values[i].input);
245 break;
246 case RINTF:
247 received = rintf(values[i].input);
248 break;
249 case RINTL:
250 received = rintl(values[i].input);
251 break;
252 default:
253 atf_tc_fail("impossible");
254 }
255
256 /*
257 * Verify FE_INEXACT was raised or not,
258 * depending on whether there was rounding and
259 * whether the function is supposed to raise
260 * exceptions.
261 */
262 if (expect_except) {
263 ATF_CHECK_MSG(fetestexcept(FE_INEXACT) != 0,
264 "[%u] %s %s(%f)"
265 " failed to raise FE_INEXACT",
266 i, rmname(values[i].round_mode),
267 fnname[fn], values[i].input);
268 } else {
269 ATF_CHECK_MSG(fetestexcept(FE_INEXACT) == 0,
270 "[%u] %s %s(%f)"
271 " spuriously raised FE_INEXACT",
272 i, rmname(values[i].round_mode),
273 fnname[fn], values[i].input);
274 }
275
276 /*
277 * Verify the fractional part of the result is
278 * zero -- the result of rounding to an integer
279 * is supposed to be an integer.
280 */
281 fpart = modf(received, &ipart);
282 ATF_CHECK_MSG(fpart == 0,
283 "[%u] %s %s(%f)=%f has fractional part %f"
284 " (integer part %f)",
285 i, rmname(values[i].round_mode), fnname[fn],
286 values[i].input, received, fpart, ipart);
287
288 /*
289 * Assuming the result we got has zero
290 * fractional part, casting to long int should
291 * have no rounding. Verify it matches the
292 * integer we expect.
293 */
294 ATF_CHECK_MSG((long int)received == values[i].expected,
295 "[%u] %s %s(%f): got %f, expected %ld",
296 i, rmname(values[i].round_mode), fnname[fn],
297 values[i].input, received, values[i].expected);
298
299 /* Do we get the same rounding mode out? */
300 ATF_CHECK_MSG(fegetround() == values[i].round_mode,
301 "[%u] %s: set %d (%s), got %d (%s)",
302 i, fnname[fn],
303 values[i].round_mode, rmname(values[i].round_mode),
304 fegetround(), rmname(fegetround()));
305 }
306 }
307 }
308
309 #ifdef __HAVE_LONG_DOUBLE
310
311 /*
312 * Use one bit more than fits in IEEE 754 binary64.
313 */
314 static const struct {
315 int round_mode;
316 long double input;
317 int64_t expected;
318 } valuesl[] = {
319 { FE_TOWARDZERO, 0x2.00000000000008p+52L, 0x20000000000000 },
320 { FE_DOWNWARD, 0x2.00000000000008p+52L, 0x20000000000000 },
321 { FE_UPWARD, 0x2.00000000000008p+52L, 0x20000000000001 },
322 { FE_TONEAREST, 0x2.00000000000008p+52L, 0x20000000000000 },
323 { FE_TOWARDZERO, 0x2.00000000000018p+52L, 0x20000000000001 },
324 { FE_DOWNWARD, 0x2.00000000000018p+52L, 0x20000000000001 },
325 { FE_UPWARD, 0x2.00000000000018p+52L, 0x20000000000002 },
326 { FE_TONEAREST, 0x2.00000000000018p+52L, 0x20000000000002 },
327 };
328
329 ATF_TC(fe_nearbyintl_rintl);
ATF_TC_HEAD(fe_nearbyintl_rintl,tc)330 ATF_TC_HEAD(fe_nearbyintl_rintl, tc)
331 {
332 atf_tc_set_md_var(tc, "descr",
333 "Checking IEEE 754 rounding modes using nearbyintl/rintl");
334 }
335
ATF_TC_BODY(fe_nearbyintl_rintl,tc)336 ATF_TC_BODY(fe_nearbyintl_rintl, tc)
337 {
338 enum {
339 RINTL,
340 NEARBYINTL,
341 N_FN,
342 } fn;
343 static const char *const fnname[] = {
344 [RINTL] = "rintl",
345 [NEARBYINTL] = "nearbyintl",
346 };
347 long double received, ipart, fpart;
348 unsigned i;
349
350 for (i = 0; i < __arraycount(valuesl); i++) {
351 for (fn = 0; fn < N_FN; fn++) {
352 bool expect_except =
353 (valuesl[i].input !=
354 (long double)valuesl[i].expected);
355
356 /*
357 * Set the requested rounding mode.
358 */
359 fesetround(valuesl[i].round_mode);
360
361 /*
362 * Clear sticky floating-point exception bits
363 * so we can verify whether the FE_INEXACT
364 * exception is raised.
365 */
366 feclearexcept(FE_ALL_EXCEPT);
367
368 /*
369 * Call the rint(3)-family function.
370 */
371 switch (fn) {
372 case NEARBYINTL:
373 received = nearbyintl(valuesl[i].input);
374 expect_except = false;
375 break;
376 case RINTL:
377 received = rintl(valuesl[i].input);
378 break;
379 default:
380 atf_tc_fail("impossible");
381 }
382
383 /*
384 * Verify FE_INEXACT was raised or not,
385 * depending on whether there was rounding and
386 * whether the function is supposed to raise
387 * exceptions.
388 */
389 if (expect_except) {
390 ATF_CHECK_MSG(fetestexcept(FE_INEXACT) != 0,
391 "[%u] %s %s(%Lf)"
392 " failed to raise FE_INEXACT",
393 i, rmname(valuesl[i].round_mode),
394 fnname[fn], valuesl[i].input);
395 } else {
396 ATF_CHECK_MSG(fetestexcept(FE_INEXACT) == 0,
397 "[%u] %s %s(%Lf)"
398 " spuriously raised FE_INEXACT",
399 i, rmname(valuesl[i].round_mode),
400 fnname[fn], valuesl[i].input);
401 }
402
403 /*
404 * Verify the fractional part of the result is
405 * zero -- the result of rounding to an integer
406 * is supposed to be an integer.
407 */
408 fpart = modfl(received, &ipart);
409 ATF_CHECK_MSG(fpart == 0,
410 "[%u] %s %s(%Lf)=%Lf has fractional part %Lf"
411 " (integer part %Lf)",
412 i, rmname(valuesl[i].round_mode), fnname[fn],
413 valuesl[i].input, received, fpart, ipart);
414
415 /*
416 * Assuming the result we got has zero
417 * fractional part, casting to int64_t should
418 * have no rounding. Verify it matches the
419 * integer we expect.
420 */
421 ATF_CHECK_MSG(((int64_t)received ==
422 valuesl[i].expected),
423 "[%u] %s %s(%Lf): got %Lf, expected %jd",
424 i, rmname(valuesl[i].round_mode), fnname[fn],
425 valuesl[i].input, received, valuesl[i].expected);
426
427 /* Do we get the same rounding mode out? */
428 ATF_CHECK_MSG(fegetround() == valuesl[i].round_mode,
429 "[%u] %s: set %d (%s), got %d (%s)",
430 i, fnname[fn],
431 valuesl[i].round_mode,
432 rmname(valuesl[i].round_mode),
433 fegetround(), rmname(fegetround()));
434 }
435 }
436 }
437
438 #endif /* __HAVE_LONG_DOUBLE */
439
ATF_TP_ADD_TCS(tp)440 ATF_TP_ADD_TCS(tp)
441 {
442
443 ATF_TP_ADD_TC(tp, fe_lrint);
444 ATF_TP_ADD_TC(tp, fe_nearbyint_rint);
445 #ifdef __HAVE_LONG_DOUBLE
446 ATF_TP_ADD_TC(tp, fe_nearbyintl_rintl);
447 #endif
448
449 return atf_no_error();
450 }
451
452 #else /* !__HAVE_FENV */
453
ATF_TP_ADD_TCS(tp)454 ATF_TP_ADD_TCS(tp)
455 {
456
457 /*
458 * No fenv, no fesetround to test.
459 */
460 return atf_no_error();
461 }
462
463 #endif /* __HAVE_FENV */
464