Home | History | Annotate | Line # | Download | only in libm
      1 /*	$NetBSD: t_fe_round.c,v 1.21 2025/04/17 13:45:22 riastradh Exp $	*/
      2 
      3 /*
      4  * Written by Maya Rashish <maya (at) 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.21 2025/04/17 13:45:22 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 #ifdef __sparc64__
    198 	atf_tc_expect_fail("PR port-sparc64/59310:"
    199 	    " t_fe_round:fe_nearbyint_rint tests are failing");
    200 #endif
    201 
    202 	for (i = 0; i < __arraycount(values); i++) {
    203 		for (fn = 0; fn < N_FN; fn++) {
    204 			bool expect_except =
    205 			    values[i].input != (double)values[i].expected;
    206 
    207 			/*
    208 			 * Set the requested rounding mode.
    209 			 */
    210 			fesetround(values[i].round_mode);
    211 
    212 #ifdef __ia64__
    213 /*
    214  * Without this barrier, we get:
    215  *
    216  * /tmp//ccJayu9g.s:2793: Warning: Use of 'mov.m' violates RAW dependency 'AR[FPSR].sf0.flags' (impliedf)
    217  * /tmp//ccJayu9g.s:2793: Warning: Only the first path encountering the conflict is reported
    218  * /tmp//ccJayu9g.s:2757: Warning: This is the location of the conflicting usage
    219  *
    220  * (If you fix this, remove the entry from doc/HACKS.)
    221  */
    222 			__insn_barrier();
    223 #endif
    224 
    225 			/*
    226 			 * Clear sticky floating-point exception bits
    227 			 * so we can verify whether the FE_INEXACT
    228 			 * exception is raised.
    229 			 */
    230 			feclearexcept(FE_ALL_EXCEPT);
    231 
    232 			/*
    233 			 * Call the rint(3)-family function.
    234 			 */
    235 			switch (fn) {
    236 			case NEARBYINT:
    237 				received = nearbyint(values[i].input);
    238 				expect_except = false;
    239 				break;
    240 			case NEARBYINTF:
    241 				received = nearbyintf(values[i].input);
    242 				expect_except = false;
    243 				break;
    244 			case NEARBYINTL:
    245 				received = nearbyintl(values[i].input);
    246 				expect_except = false;
    247 				break;
    248 			case RINT:
    249 				received = rint(values[i].input);
    250 				break;
    251 			case RINTF:
    252 				received = rintf(values[i].input);
    253 				break;
    254 			case RINTL:
    255 				received = rintl(values[i].input);
    256 				break;
    257 			default:
    258 				atf_tc_fail("impossible");
    259 			}
    260 
    261 			/*
    262 			 * Verify FE_INEXACT was raised or not,
    263 			 * depending on whether there was rounding and
    264 			 * whether the function is supposed to raise
    265 			 * exceptions.
    266 			 */
    267 			if (expect_except) {
    268 				ATF_CHECK_MSG(fetestexcept(FE_INEXACT) != 0,
    269 				    "[%u] %s %s(%f)"
    270 				    " failed to raise FE_INEXACT",
    271 				    i, rmname(values[i].round_mode),
    272 				    fnname[fn], values[i].input);
    273 			} else {
    274 				ATF_CHECK_MSG(fetestexcept(FE_INEXACT) == 0,
    275 				    "[%u] %s %s(%f)"
    276 				    " spuriously raised FE_INEXACT",
    277 				    i, rmname(values[i].round_mode),
    278 				    fnname[fn], values[i].input);
    279 			}
    280 
    281 			/*
    282 			 * Verify the fractional part of the result is
    283 			 * zero -- the result of rounding to an integer
    284 			 * is supposed to be an integer.
    285 			 */
    286 			fpart = modf(received, &ipart);
    287 			ATF_CHECK_MSG(fpart == 0,
    288 			    "[%u] %s %s(%f)=%f has fractional part %f"
    289 			    " (integer part %f)",
    290 			    i, rmname(values[i].round_mode), fnname[fn],
    291 			    values[i].input, received, fpart, ipart);
    292 
    293 			/*
    294 			 * Assuming the result we got has zero
    295 			 * fractional part, casting to long int should
    296 			 * have no rounding.  Verify it matches the
    297 			 * integer we expect.
    298 			 */
    299 			ATF_CHECK_MSG((long int)received == values[i].expected,
    300 			    "[%u] %s %s(%f): got %f, expected %ld",
    301 			    i, rmname(values[i].round_mode), fnname[fn],
    302 			    values[i].input, received, values[i].expected);
    303 
    304 			/* Do we get the same rounding mode out? */
    305 			ATF_CHECK_MSG(fegetround() == values[i].round_mode,
    306 			    "[%u] %s: set %d (%s), got %d (%s)",
    307 			    i, fnname[fn],
    308 			    values[i].round_mode, rmname(values[i].round_mode),
    309 			    fegetround(), rmname(fegetround()));
    310 		}
    311 	}
    312 }
    313 
    314 #ifdef __HAVE_LONG_DOUBLE
    315 
    316 /*
    317  * Use one bit more than fits in IEEE 754 binary64.
    318  */
    319 static const struct {
    320 	int round_mode;
    321 	long double input;
    322 	int64_t expected;
    323 } valuesl[] = {
    324 	{ FE_TOWARDZERO,	0x2.00000000000008p+52L, 0x20000000000000 },
    325 	{ FE_DOWNWARD,		0x2.00000000000008p+52L, 0x20000000000000 },
    326 	{ FE_UPWARD,		0x2.00000000000008p+52L, 0x20000000000001 },
    327 	{ FE_TONEAREST,		0x2.00000000000008p+52L, 0x20000000000000 },
    328 	{ FE_TOWARDZERO,	0x2.00000000000018p+52L, 0x20000000000001 },
    329 	{ FE_DOWNWARD,		0x2.00000000000018p+52L, 0x20000000000001 },
    330 	{ FE_UPWARD,		0x2.00000000000018p+52L, 0x20000000000002 },
    331 	{ FE_TONEAREST,		0x2.00000000000018p+52L, 0x20000000000002 },
    332 };
    333 
    334 ATF_TC(fe_nearbyintl_rintl);
    335 ATF_TC_HEAD(fe_nearbyintl_rintl, tc)
    336 {
    337 	atf_tc_set_md_var(tc, "descr",
    338 	    "Checking IEEE 754 rounding modes using nearbyintl/rintl");
    339 }
    340 
    341 ATF_TC_BODY(fe_nearbyintl_rintl, tc)
    342 {
    343 	enum {
    344 		RINTL,
    345 		NEARBYINTL,
    346 		N_FN,
    347 	} fn;
    348 	static const char *const fnname[] = {
    349 		[RINTL] = "rintl",
    350 		[NEARBYINTL] = "nearbyintl",
    351 	};
    352 	long double received, ipart, fpart;
    353 	unsigned i;
    354 
    355 #ifdef __sparc64__
    356 	atf_tc_expect_fail("PR port-sparc64/59310:"
    357 	    " t_fe_round:fe_nearbyint_rint tests are failing");
    358 #endif
    359 
    360 	for (i = 0; i < __arraycount(valuesl); i++) {
    361 		for (fn = 0; fn < N_FN; fn++) {
    362 			bool expect_except =
    363 			    (valuesl[i].input !=
    364 				(long double)valuesl[i].expected);
    365 
    366 			/*
    367 			 * Set the requested rounding mode.
    368 			 */
    369 			fesetround(valuesl[i].round_mode);
    370 
    371 			/*
    372 			 * Clear sticky floating-point exception bits
    373 			 * so we can verify whether the FE_INEXACT
    374 			 * exception is raised.
    375 			 */
    376 			feclearexcept(FE_ALL_EXCEPT);
    377 
    378 			/*
    379 			 * Call the rint(3)-family function.
    380 			 */
    381 			switch (fn) {
    382 			case NEARBYINTL:
    383 				received = nearbyintl(valuesl[i].input);
    384 				expect_except = false;
    385 				break;
    386 			case RINTL:
    387 				received = rintl(valuesl[i].input);
    388 				break;
    389 			default:
    390 				atf_tc_fail("impossible");
    391 			}
    392 
    393 			/*
    394 			 * Verify FE_INEXACT was raised or not,
    395 			 * depending on whether there was rounding and
    396 			 * whether the function is supposed to raise
    397 			 * exceptions.
    398 			 */
    399 			if (expect_except) {
    400 				ATF_CHECK_MSG(fetestexcept(FE_INEXACT) != 0,
    401 				    "[%u] %s %s(%Lf)"
    402 				    " failed to raise FE_INEXACT",
    403 				    i, rmname(valuesl[i].round_mode),
    404 				    fnname[fn], valuesl[i].input);
    405 			} else {
    406 				ATF_CHECK_MSG(fetestexcept(FE_INEXACT) == 0,
    407 				    "[%u] %s %s(%Lf)"
    408 				    " spuriously raised FE_INEXACT",
    409 				    i, rmname(valuesl[i].round_mode),
    410 				    fnname[fn], valuesl[i].input);
    411 			}
    412 
    413 			/*
    414 			 * Verify the fractional part of the result is
    415 			 * zero -- the result of rounding to an integer
    416 			 * is supposed to be an integer.
    417 			 */
    418 			fpart = modfl(received, &ipart);
    419 			ATF_CHECK_MSG(fpart == 0,
    420 			    "[%u] %s %s(%Lf)=%Lf has fractional part %Lf"
    421 			    " (integer part %Lf)",
    422 			    i, rmname(valuesl[i].round_mode), fnname[fn],
    423 			    valuesl[i].input, received, fpart, ipart);
    424 
    425 			/*
    426 			 * Assuming the result we got has zero
    427 			 * fractional part, casting to int64_t should
    428 			 * have no rounding.  Verify it matches the
    429 			 * integer we expect.
    430 			 */
    431 			ATF_CHECK_MSG(((int64_t)received ==
    432 				valuesl[i].expected),
    433 			    "[%u] %s %s(%Lf): got %Lf, expected %jd",
    434 			    i, rmname(valuesl[i].round_mode), fnname[fn],
    435 			    valuesl[i].input, received, valuesl[i].expected);
    436 
    437 			/* Do we get the same rounding mode out? */
    438 			ATF_CHECK_MSG(fegetround() == valuesl[i].round_mode,
    439 			    "[%u] %s: set %d (%s), got %d (%s)",
    440 			    i, fnname[fn],
    441 			    valuesl[i].round_mode,
    442 			    rmname(valuesl[i].round_mode),
    443 			    fegetround(), rmname(fegetround()));
    444 		}
    445 	}
    446 }
    447 
    448 #endif	/* __HAVE_LONG_DOUBLE */
    449 
    450 ATF_TP_ADD_TCS(tp)
    451 {
    452 
    453 	ATF_TP_ADD_TC(tp, fe_lrint);
    454 	ATF_TP_ADD_TC(tp, fe_nearbyint_rint);
    455 #ifdef __HAVE_LONG_DOUBLE
    456 	ATF_TP_ADD_TC(tp, fe_nearbyintl_rintl);
    457 #endif
    458 
    459 	return atf_no_error();
    460 }
    461 
    462 #else  /* !__HAVE_FENV */
    463 
    464 ATF_TP_ADD_TCS(tp)
    465 {
    466 
    467 	/*
    468 	 * No fenv, no fesetround to test.
    469 	 */
    470 	return atf_no_error();
    471 }
    472 
    473 #endif	/* __HAVE_FENV */
    474