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 * 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); 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 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); 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 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); 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 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 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 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