1 /* 2 * dotest.c - actually generate mathlib test cases 3 * 4 * Copyright (c) 1999-2024, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include <stdio.h> 9 #include <string.h> 10 #include <stdlib.h> 11 #include <stdint.h> 12 #include <assert.h> 13 #include <limits.h> 14 15 #include "semi.h" 16 #include "intern.h" 17 #include "random.h" 18 19 #define MPFR_PREC 96 /* good enough for float or double + a few extra bits */ 20 21 #if MPFR_VERSION < MPFR_VERSION_NUM(4, 2, 0) 22 int 23 mpfr_tanpi (mpfr_t ret, const mpfr_t arg, mpfr_rnd_t rnd) 24 { 25 MPFR_DECL_INIT (frd, MPFR_PREC); 26 mpfr_const_pi (frd, GMP_RNDN); 27 mpfr_mul (frd, frd, arg, GMP_RNDN); 28 return mpfr_tan (ret, frd, GMP_RNDN); 29 } 30 31 int 32 mpfr_sinpi (mpfr_t ret, const mpfr_t arg, mpfr_rnd_t rnd) 33 { 34 MPFR_DECL_INIT (frd, MPFR_PREC); 35 mpfr_const_pi (frd, GMP_RNDN); 36 mpfr_mul (frd, frd, arg, GMP_RNDN); 37 return mpfr_sin (ret, frd, GMP_RNDN); 38 } 39 40 int 41 mpfr_cospi (mpfr_t ret, const mpfr_t arg, mpfr_rnd_t rnd) 42 { 43 MPFR_DECL_INIT (frd, MPFR_PREC); 44 mpfr_const_pi (frd, GMP_RNDN); 45 mpfr_mul (frd, frd, arg, GMP_RNDN); 46 return mpfr_cos (ret, frd, GMP_RNDN); 47 } 48 #endif 49 50 extern int lib_fo, lib_no_arith, ntests; 51 52 /* 53 * Prototypes. 54 */ 55 static void cases_biased(uint32 *, uint32, uint32); 56 static void cases_biased_positive(uint32 *, uint32, uint32); 57 static void cases_biased_float(uint32 *, uint32, uint32); 58 static void cases_uniform(uint32 *, uint32, uint32); 59 static void cases_uniform_positive(uint32 *, uint32, uint32); 60 static void cases_uniform_float(uint32 *, uint32, uint32); 61 static void cases_uniform_float_positive(uint32 *, uint32, uint32); 62 static void log_cases(uint32 *, uint32, uint32); 63 static void log_cases_float(uint32 *, uint32, uint32); 64 static void log1p_cases(uint32 *, uint32, uint32); 65 static void log1p_cases_float(uint32 *, uint32, uint32); 66 static void minmax_cases(uint32 *, uint32, uint32); 67 static void minmax_cases_float(uint32 *, uint32, uint32); 68 static void atan2_cases(uint32 *, uint32, uint32); 69 static void atan2_cases_float(uint32 *, uint32, uint32); 70 static void pow_cases(uint32 *, uint32, uint32); 71 static void pow_cases_float(uint32 *, uint32, uint32); 72 static void rred_cases(uint32 *, uint32, uint32); 73 static void rred_cases_float(uint32 *, uint32, uint32); 74 static void cases_semi1(uint32 *, uint32, uint32); 75 static void cases_semi1_float(uint32 *, uint32, uint32); 76 static void cases_semi2(uint32 *, uint32, uint32); 77 static void cases_semi2_float(uint32 *, uint32, uint32); 78 static void cases_ldexp(uint32 *, uint32, uint32); 79 static void cases_ldexp_float(uint32 *, uint32, uint32); 80 81 static void complex_cases_uniform(uint32 *, uint32, uint32); 82 static void complex_cases_uniform_float(uint32 *, uint32, uint32); 83 static void complex_cases_biased(uint32 *, uint32, uint32); 84 static void complex_cases_biased_float(uint32 *, uint32, uint32); 85 static void complex_log_cases(uint32 *, uint32, uint32); 86 static void complex_log_cases_float(uint32 *, uint32, uint32); 87 static void complex_pow_cases(uint32 *, uint32, uint32); 88 static void complex_pow_cases_float(uint32 *, uint32, uint32); 89 static void complex_arithmetic_cases(uint32 *, uint32, uint32); 90 static void complex_arithmetic_cases_float(uint32 *, uint32, uint32); 91 92 static uint32 doubletop(int x, int scale); 93 static uint32 floatval(int x, int scale); 94 95 /* 96 * Convert back and forth between IEEE bit patterns and the 97 * mpfr_t/mpc_t types. 98 */ 99 static void set_mpfr_d(mpfr_t x, uint32 h, uint32 l) 100 { 101 uint64_t hl = ((uint64_t)h << 32) | l; 102 uint32 exp = (hl >> 52) & 0x7ff; 103 int64_t mantissa = hl & (((uint64_t)1 << 52) - 1); 104 int sign = (hl >> 63) ? -1 : +1; 105 if (exp == 0x7ff) { 106 if (mantissa == 0) 107 mpfr_set_inf(x, sign); 108 else 109 mpfr_set_nan(x); 110 } else if (exp == 0 && mantissa == 0) { 111 mpfr_set_ui(x, 0, GMP_RNDN); 112 mpfr_setsign(x, x, sign < 0, GMP_RNDN); 113 } else { 114 if (exp != 0) 115 mantissa |= ((uint64_t)1 << 52); 116 else 117 exp++; 118 mpfr_set_sj_2exp(x, mantissa * sign, (int)exp - 0x3ff - 52, GMP_RNDN); 119 } 120 } 121 static void set_mpfr_f(mpfr_t x, uint32 f) 122 { 123 uint32 exp = (f >> 23) & 0xff; 124 int32 mantissa = f & ((1 << 23) - 1); 125 int sign = (f >> 31) ? -1 : +1; 126 if (exp == 0xff) { 127 if (mantissa == 0) 128 mpfr_set_inf(x, sign); 129 else 130 mpfr_set_nan(x); 131 } else if (exp == 0 && mantissa == 0) { 132 mpfr_set_ui(x, 0, GMP_RNDN); 133 mpfr_setsign(x, x, sign < 0, GMP_RNDN); 134 } else { 135 if (exp != 0) 136 mantissa |= (1 << 23); 137 else 138 exp++; 139 mpfr_set_sj_2exp(x, mantissa * sign, (int)exp - 0x7f - 23, GMP_RNDN); 140 } 141 } 142 static void set_mpc_d(mpc_t z, uint32 rh, uint32 rl, uint32 ih, uint32 il) 143 { 144 mpfr_t x, y; 145 mpfr_init2(x, MPFR_PREC); 146 mpfr_init2(y, MPFR_PREC); 147 set_mpfr_d(x, rh, rl); 148 set_mpfr_d(y, ih, il); 149 mpc_set_fr_fr(z, x, y, MPC_RNDNN); 150 mpfr_clear(x); 151 mpfr_clear(y); 152 } 153 static void set_mpc_f(mpc_t z, uint32 r, uint32 i) 154 { 155 mpfr_t x, y; 156 mpfr_init2(x, MPFR_PREC); 157 mpfr_init2(y, MPFR_PREC); 158 set_mpfr_f(x, r); 159 set_mpfr_f(y, i); 160 mpc_set_fr_fr(z, x, y, MPC_RNDNN); 161 mpfr_clear(x); 162 mpfr_clear(y); 163 } 164 static void get_mpfr_d(const mpfr_t x, uint32 *h, uint32 *l, uint32 *extra) 165 { 166 uint32_t sign, expfield, mantfield; 167 mpfr_t significand; 168 int exp; 169 170 if (mpfr_nan_p(x)) { 171 *h = 0x7ff80000; 172 *l = 0; 173 *extra = 0; 174 return; 175 } 176 177 sign = mpfr_signbit(x) ? 0x80000000U : 0; 178 179 if (mpfr_inf_p(x)) { 180 *h = 0x7ff00000 | sign; 181 *l = 0; 182 *extra = 0; 183 return; 184 } 185 186 if (mpfr_zero_p(x)) { 187 *h = 0x00000000 | sign; 188 *l = 0; 189 *extra = 0; 190 return; 191 } 192 193 mpfr_init2(significand, MPFR_PREC); 194 mpfr_set(significand, x, GMP_RNDN); 195 exp = mpfr_get_exp(significand); 196 mpfr_set_exp(significand, 0); 197 198 /* Now significand is in [1/2,1), and significand * 2^exp == x. 199 * So the IEEE exponent corresponding to exp==0 is 0x3fe. */ 200 if (exp > 0x400) { 201 /* overflow to infinity anyway */ 202 *h = 0x7ff00000 | sign; 203 *l = 0; 204 *extra = 0; 205 mpfr_clear(significand); 206 return; 207 } 208 209 if (exp <= -0x3fe || mpfr_zero_p(x)) 210 exp = -0x3fd; /* denormalise */ 211 expfield = exp + 0x3fd; /* offset to cancel leading mantissa bit */ 212 213 mpfr_div_2si(significand, x, exp - 21, GMP_RNDN); 214 mpfr_abs(significand, significand, GMP_RNDN); 215 mantfield = mpfr_get_ui(significand, GMP_RNDZ); 216 *h = sign + ((uint64_t)expfield << 20) + mantfield; 217 mpfr_sub_ui(significand, significand, mantfield, GMP_RNDN); 218 mpfr_mul_2ui(significand, significand, 32, GMP_RNDN); 219 mantfield = mpfr_get_ui(significand, GMP_RNDZ); 220 *l = mantfield; 221 mpfr_sub_ui(significand, significand, mantfield, GMP_RNDN); 222 mpfr_mul_2ui(significand, significand, 32, GMP_RNDN); 223 mantfield = mpfr_get_ui(significand, GMP_RNDZ); 224 *extra = mantfield; 225 226 mpfr_clear(significand); 227 } 228 static void get_mpfr_f(const mpfr_t x, uint32 *f, uint32 *extra) 229 { 230 uint32_t sign, expfield, mantfield; 231 mpfr_t significand; 232 int exp; 233 234 if (mpfr_nan_p(x)) { 235 *f = 0x7fc00000; 236 *extra = 0; 237 return; 238 } 239 240 sign = mpfr_signbit(x) ? 0x80000000U : 0; 241 242 if (mpfr_inf_p(x)) { 243 *f = 0x7f800000 | sign; 244 *extra = 0; 245 return; 246 } 247 248 if (mpfr_zero_p(x)) { 249 *f = 0x00000000 | sign; 250 *extra = 0; 251 return; 252 } 253 254 mpfr_init2(significand, MPFR_PREC); 255 mpfr_set(significand, x, GMP_RNDN); 256 exp = mpfr_get_exp(significand); 257 mpfr_set_exp(significand, 0); 258 259 /* Now significand is in [1/2,1), and significand * 2^exp == x. 260 * So the IEEE exponent corresponding to exp==0 is 0x7e. */ 261 if (exp > 0x80) { 262 /* overflow to infinity anyway */ 263 *f = 0x7f800000 | sign; 264 *extra = 0; 265 mpfr_clear(significand); 266 return; 267 } 268 269 if (exp <= -0x7e || mpfr_zero_p(x)) 270 exp = -0x7d; /* denormalise */ 271 expfield = exp + 0x7d; /* offset to cancel leading mantissa bit */ 272 273 mpfr_div_2si(significand, x, exp - 24, GMP_RNDN); 274 mpfr_abs(significand, significand, GMP_RNDN); 275 mantfield = mpfr_get_ui(significand, GMP_RNDZ); 276 *f = sign + ((uint64_t)expfield << 23) + mantfield; 277 mpfr_sub_ui(significand, significand, mantfield, GMP_RNDN); 278 mpfr_mul_2ui(significand, significand, 32, GMP_RNDN); 279 mantfield = mpfr_get_ui(significand, GMP_RNDZ); 280 *extra = mantfield; 281 282 mpfr_clear(significand); 283 } 284 static void get_mpc_d(const mpc_t z, 285 uint32 *rh, uint32 *rl, uint32 *rextra, 286 uint32 *ih, uint32 *il, uint32 *iextra) 287 { 288 mpfr_t x, y; 289 mpfr_init2(x, MPFR_PREC); 290 mpfr_init2(y, MPFR_PREC); 291 mpc_real(x, z, GMP_RNDN); 292 mpc_imag(y, z, GMP_RNDN); 293 get_mpfr_d(x, rh, rl, rextra); 294 get_mpfr_d(y, ih, il, iextra); 295 mpfr_clear(x); 296 mpfr_clear(y); 297 } 298 static void get_mpc_f(const mpc_t z, 299 uint32 *r, uint32 *rextra, 300 uint32 *i, uint32 *iextra) 301 { 302 mpfr_t x, y; 303 mpfr_init2(x, MPFR_PREC); 304 mpfr_init2(y, MPFR_PREC); 305 mpc_real(x, z, GMP_RNDN); 306 mpc_imag(y, z, GMP_RNDN); 307 get_mpfr_f(x, r, rextra); 308 get_mpfr_f(y, i, iextra); 309 mpfr_clear(x); 310 mpfr_clear(y); 311 } 312 313 /* 314 * Implementation of mathlib functions that aren't trivially 315 * implementable using an existing mpfr or mpc function. 316 */ 317 int test_rred(mpfr_t ret, const mpfr_t x, int *quadrant) 318 { 319 mpfr_t halfpi; 320 long quo; 321 int status; 322 323 /* 324 * In the worst case of range reduction, we get an input of size 325 * around 2^1024, and must find its remainder mod pi, which means 326 * we need 1024 bits of pi at least. Plus, the remainder might 327 * happen to come out very very small if we're unlucky. How 328 * unlucky can we be? Well, conveniently, I once went through and 329 * actually worked that out using Paxson's modular minimisation 330 * algorithm, and it turns out that the smallest exponent you can 331 * get out of a nontrivial[1] double precision range reduction is 332 * 0x3c2, i.e. of the order of 2^-61. So we need 1024 bits of pi 333 * to get us down to the units digit, another 61 or so bits (say 334 * 64) to get down to the highest set bit of the output, and then 335 * some bits to make the actual mantissa big enough. 336 * 337 * [1] of course the output of range reduction can have an 338 * arbitrarily small exponent in the trivial case, where the 339 * input is so small that it's the identity function. That 340 * doesn't count. 341 */ 342 mpfr_init2(halfpi, MPFR_PREC + 1024 + 64); 343 mpfr_const_pi(halfpi, GMP_RNDN); 344 mpfr_div_ui(halfpi, halfpi, 2, GMP_RNDN); 345 346 status = mpfr_remquo(ret, &quo, x, halfpi, GMP_RNDN); 347 *quadrant = quo & 3; 348 349 mpfr_clear(halfpi); 350 351 return status; 352 } 353 int test_lgamma(mpfr_t ret, const mpfr_t x, mpfr_rnd_t rnd) 354 { 355 /* 356 * mpfr_lgamma takes an extra int * parameter to hold the output 357 * sign. We don't bother testing that, so this wrapper throws away 358 * the sign and hence fits into the same function prototype as all 359 * the other real->real mpfr functions. 360 * 361 * There is also mpfr_lngamma which has no sign output and hence 362 * has the right prototype already, but unfortunately it returns 363 * NaN in cases where gamma(x) < 0, so it's no use to us. 364 */ 365 int sign; 366 return mpfr_lgamma(ret, &sign, x, rnd); 367 } 368 int test_cpow(mpc_t ret, const mpc_t x, const mpc_t y, mpc_rnd_t rnd) 369 { 370 /* 371 * For complex pow, we must bump up the precision by a huge amount 372 * if we want it to get the really difficult cases right. (Not 373 * that we expect the library under test to be getting those cases 374 * right itself, but we'd at least like the test suite to report 375 * them as wrong for the _right reason_.) 376 * 377 * This works around a bug in mpc_pow(), fixed by r1455 in the MPC 378 * svn repository (2014-10-14) and expected to be in any MPC 379 * release after 1.0.2 (which was the latest release already made 380 * at the time of the fix). So as and when we update to an MPC 381 * with the fix in it, we could remove this workaround. 382 * 383 * For the reasons for choosing this amount of extra precision, 384 * see analysis in complex/cpownotes.txt for the rationale for the 385 * amount. 386 */ 387 mpc_t xbig, ybig, retbig; 388 int status; 389 390 mpc_init2(xbig, 1034 + 53 + 60 + MPFR_PREC); 391 mpc_init2(ybig, 1034 + 53 + 60 + MPFR_PREC); 392 mpc_init2(retbig, 1034 + 53 + 60 + MPFR_PREC); 393 394 mpc_set(xbig, x, MPC_RNDNN); 395 mpc_set(ybig, y, MPC_RNDNN); 396 status = mpc_pow(retbig, xbig, ybig, rnd); 397 mpc_set(ret, retbig, rnd); 398 399 mpc_clear(xbig); 400 mpc_clear(ybig); 401 mpc_clear(retbig); 402 403 return status; 404 } 405 406 /* 407 * Identify 'hard' values (NaN, Inf, nonzero denormal) for deciding 408 * whether microlib will decline to run a test. 409 */ 410 #define is_shard(in) ( \ 411 (((in)[0] & 0x7F800000) == 0x7F800000 || \ 412 (((in)[0] & 0x7F800000) == 0 && ((in)[0]&0x7FFFFFFF) != 0))) 413 414 #define is_dhard(in) ( \ 415 (((in)[0] & 0x7FF00000) == 0x7FF00000 || \ 416 (((in)[0] & 0x7FF00000) == 0 && (((in)[0] & 0xFFFFF) | (in)[1]) != 0))) 417 418 /* 419 * Identify integers. 420 */ 421 int is_dinteger(uint32 *in) 422 { 423 uint32 out[3]; 424 if ((0x7FF00000 & ~in[0]) == 0) 425 return 0; /* not finite, hence not integer */ 426 test_ceil(in, out); 427 return in[0] == out[0] && in[1] == out[1]; 428 } 429 int is_sinteger(uint32 *in) 430 { 431 uint32 out[3]; 432 if ((0x7F800000 & ~in[0]) == 0) 433 return 0; /* not finite, hence not integer */ 434 test_ceilf(in, out); 435 return in[0] == out[0]; 436 } 437 438 /* 439 * Identify signalling NaNs. 440 */ 441 int is_dsnan(const uint32 *in) 442 { 443 if ((in[0] & 0x7FF00000) != 0x7FF00000) 444 return 0; /* not the inf/nan exponent */ 445 if ((in[0] << 12) == 0 && in[1] == 0) 446 return 0; /* inf */ 447 if (in[0] & 0x00080000) 448 return 0; /* qnan */ 449 return 1; 450 } 451 int is_ssnan(const uint32 *in) 452 { 453 if ((in[0] & 0x7F800000) != 0x7F800000) 454 return 0; /* not the inf/nan exponent */ 455 if ((in[0] << 9) == 0) 456 return 0; /* inf */ 457 if (in[0] & 0x00400000) 458 return 0; /* qnan */ 459 return 1; 460 } 461 int is_snan(const uint32 *in, int size) 462 { 463 return size == 2 ? is_dsnan(in) : is_ssnan(in); 464 } 465 466 /* 467 * Wrapper functions called to fix up unusual results after the main 468 * test function has run. 469 */ 470 void universal_wrapper(wrapperctx *ctx) 471 { 472 /* 473 * Any SNaN input gives rise to a QNaN output. 474 */ 475 int op; 476 for (op = 0; op < wrapper_get_nops(ctx); op++) { 477 int size = wrapper_get_size(ctx, op); 478 479 if (!wrapper_is_complex(ctx, op) && 480 is_snan(wrapper_get_ieee(ctx, op), size)) { 481 wrapper_set_nan(ctx); 482 } 483 } 484 } 485 486 /* clang-format off */ 487 Testable functions[] = { 488 /* 489 * Trig functions: sin, cos, tan. We test the core function 490 * between -16 and +16: we assume that range reduction exists 491 * and will be used for larger arguments, and we'll test that 492 * separately. Also we only go down to 2^-27 in magnitude, 493 * because below that sin(x)=tan(x)=x and cos(x)=1 as far as 494 * double precision can tell, which is boring. 495 */ 496 {"sin", (funcptr)mpfr_sin, args1, {NULL}, 497 cases_uniform, 0x3e400000, 0x40300000}, 498 {"sinf", (funcptr)mpfr_sin, args1f, {NULL}, 499 cases_uniform_float, 0x39800000, 0x41800000}, 500 {"cos", (funcptr)mpfr_cos, args1, {NULL}, 501 cases_uniform, 0x3e400000, 0x40300000}, 502 {"cosf", (funcptr)mpfr_cos, args1f, {NULL}, 503 cases_uniform_float, 0x39800000, 0x41800000}, 504 {"tan", (funcptr)mpfr_tan, args1, {NULL}, 505 cases_uniform, 0x3e400000, 0x40300000}, 506 {"tanf", (funcptr)mpfr_tan, args1f, {NULL}, 507 cases_uniform_float, 0x39800000, 0x41800000}, 508 {"sincosf_sinf", (funcptr)mpfr_sin, args1f, {NULL}, 509 cases_uniform_float, 0x39800000, 0x41800000}, 510 {"sincosf_cosf", (funcptr)mpfr_cos, args1f, {NULL}, 511 cases_uniform_float, 0x39800000, 0x41800000}, 512 {"sinpi", (funcptr)mpfr_sinpi, args1, {NULL}, 513 cases_uniform, 0x3e400000, 0x40300000}, 514 {"sinpif", (funcptr)mpfr_sinpi, args1f, {NULL}, 515 cases_uniform_float, 0x39800000, 0x41800000}, 516 {"cospi", (funcptr)mpfr_cospi, args1, {NULL}, 517 cases_uniform, 0x3e400000, 0x40300000}, 518 {"cospif", (funcptr)mpfr_cospi, args1f, {NULL}, 519 cases_uniform_float, 0x39800000, 0x41800000}, 520 {"tanpi", (funcptr)mpfr_tanpi, args1, {NULL}, 521 cases_uniform, 0x3e400000, 0x40300000}, 522 {"tanpif", (funcptr)mpfr_tanpi, args1f, {NULL}, 523 cases_uniform_float, 0x39800000, 0x41800000}, 524 /* 525 * Inverse trig: asin, acos. Between 1 and -1, of course. acos 526 * goes down to 2^-54, asin to 2^-27. 527 */ 528 {"asin", (funcptr)mpfr_asin, args1, {NULL}, 529 cases_uniform, 0x3e400000, 0x3fefffff}, 530 {"asinf", (funcptr)mpfr_asin, args1f, {NULL}, 531 cases_uniform_float, 0x39800000, 0x3f7fffff}, 532 {"acos", (funcptr)mpfr_acos, args1, {NULL}, 533 cases_uniform, 0x3c900000, 0x3fefffff}, 534 {"acosf", (funcptr)mpfr_acos, args1f, {NULL}, 535 cases_uniform_float, 0x33800000, 0x3f7fffff}, 536 /* 537 * Inverse trig: atan. atan is stable (in double prec) with 538 * argument magnitude past 2^53, so we'll test up to there. 539 * atan(x) is boringly just x below 2^-27. 540 */ 541 {"atan", (funcptr)mpfr_atan, args1, {NULL}, 542 cases_uniform, 0x3e400000, 0x43400000}, 543 {"atanf", (funcptr)mpfr_atan, args1f, {NULL}, 544 cases_uniform_float, 0x39800000, 0x4b800000}, 545 /* 546 * atan2. Interesting cases arise when the exponents of the 547 * arguments differ by at most about 50. 548 */ 549 {"atan2", (funcptr)mpfr_atan2, args2, {NULL}, 550 atan2_cases, 0}, 551 {"atan2f", (funcptr)mpfr_atan2, args2f, {NULL}, 552 atan2_cases_float, 0}, 553 /* 554 * The exponentials: exp, sinh, cosh. They overflow at around 555 * 710. exp and sinh are boring below 2^-54, cosh below 2^-27. 556 */ 557 {"exp", (funcptr)mpfr_exp, args1, {NULL}, 558 cases_uniform, 0x3c900000, 0x40878000}, 559 {"expf", (funcptr)mpfr_exp, args1f, {NULL}, 560 cases_uniform_float, 0x33800000, 0x42dc0000}, 561 {"sinh", (funcptr)mpfr_sinh, args1, {NULL}, 562 cases_uniform, 0x3c900000, 0x40878000}, 563 {"sinhf", (funcptr)mpfr_sinh, args1f, {NULL}, 564 cases_uniform_float, 0x33800000, 0x42dc0000}, 565 {"cosh", (funcptr)mpfr_cosh, args1, {NULL}, 566 cases_uniform, 0x3e400000, 0x40878000}, 567 {"coshf", (funcptr)mpfr_cosh, args1f, {NULL}, 568 cases_uniform_float, 0x39800000, 0x42dc0000}, 569 /* 570 * tanh is stable past around 20. It's boring below 2^-27. 571 */ 572 {"tanh", (funcptr)mpfr_tanh, args1, {NULL}, 573 cases_uniform, 0x3e400000, 0x40340000}, 574 {"tanhf", (funcptr)mpfr_tanh, args1f, {NULL}, 575 cases_uniform, 0x39800000, 0x41100000}, 576 /* 577 * log must be tested only on positive numbers, but can cover 578 * the whole range of positive nonzero finite numbers. It never 579 * gets boring. 580 */ 581 {"log", (funcptr)mpfr_log, args1, {NULL}, log_cases, 0}, 582 {"logf", (funcptr)mpfr_log, args1f, {NULL}, log_cases_float, 0}, 583 {"log10", (funcptr)mpfr_log10, args1, {NULL}, log_cases, 0}, 584 {"log10f", (funcptr)mpfr_log10, args1f, {NULL}, log_cases_float, 0}, 585 /* 586 * pow. 587 */ 588 {"pow", (funcptr)mpfr_pow, args2, {NULL}, pow_cases, 0}, 589 {"powf", (funcptr)mpfr_pow, args2f, {NULL}, pow_cases_float, 0}, 590 /* 591 * Trig range reduction. We are able to test this for all 592 * finite values, but will only bother for things between 2^-3 593 * and 2^+52. 594 */ 595 {"rred", (funcptr)test_rred, rred, {NULL}, rred_cases, 0}, 596 {"rredf", (funcptr)test_rred, rredf, {NULL}, rred_cases_float, 0}, 597 /* 598 * Square and cube root. 599 */ 600 {"sqrt", (funcptr)mpfr_sqrt, args1, {NULL}, log_cases, 0}, 601 {"sqrtf", (funcptr)mpfr_sqrt, args1f, {NULL}, log_cases_float, 0}, 602 {"cbrt", (funcptr)mpfr_cbrt, args1, {NULL}, log_cases, 0}, 603 {"cbrtf", (funcptr)mpfr_cbrt, args1f, {NULL}, log_cases_float, 0}, 604 {"hypot", (funcptr)mpfr_hypot, args2, {NULL}, atan2_cases, 0}, 605 {"hypotf", (funcptr)mpfr_hypot, args2f, {NULL}, atan2_cases_float, 0}, 606 /* 607 * Seminumerical functions. 608 */ 609 {"ceil", (funcptr)test_ceil, semi1, {NULL}, cases_semi1}, 610 {"ceilf", (funcptr)test_ceilf, semi1f, {NULL}, cases_semi1_float}, 611 {"floor", (funcptr)test_floor, semi1, {NULL}, cases_semi1}, 612 {"floorf", (funcptr)test_floorf, semi1f, {NULL}, cases_semi1_float}, 613 {"fmod", (funcptr)test_fmod, semi2, {NULL}, cases_semi2}, 614 {"fmodf", (funcptr)test_fmodf, semi2f, {NULL}, cases_semi2_float}, 615 {"ldexp", (funcptr)test_ldexp, t_ldexp, {NULL}, cases_ldexp}, 616 {"ldexpf", (funcptr)test_ldexpf, t_ldexpf, {NULL}, cases_ldexp_float}, 617 {"frexp", (funcptr)test_frexp, t_frexp, {NULL}, cases_semi1}, 618 {"frexpf", (funcptr)test_frexpf, t_frexpf, {NULL}, cases_semi1_float}, 619 {"modf", (funcptr)test_modf, t_modf, {NULL}, cases_semi1}, 620 {"modff", (funcptr)test_modff, t_modff, {NULL}, cases_semi1_float}, 621 622 /* 623 * Classification and more semi-numericals 624 */ 625 {"copysign", (funcptr)test_copysign, semi2, {NULL}, cases_semi2}, 626 {"copysignf", (funcptr)test_copysignf, semi2f, {NULL}, cases_semi2_float}, 627 {"isfinite", (funcptr)test_isfinite, classify, {NULL}, cases_uniform, 0, 0x7fffffff}, 628 {"isfinitef", (funcptr)test_isfinitef, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 629 {"isinf", (funcptr)test_isinf, classify, {NULL}, cases_uniform, 0, 0x7fffffff}, 630 {"isinff", (funcptr)test_isinff, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 631 {"isnan", (funcptr)test_isnan, classify, {NULL}, cases_uniform, 0, 0x7fffffff}, 632 {"isnanf", (funcptr)test_isnanf, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 633 {"isnormal", (funcptr)test_isnormal, classify, {NULL}, cases_uniform, 0, 0x7fffffff}, 634 {"isnormalf", (funcptr)test_isnormalf, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 635 {"signbit", (funcptr)test_signbit, classify, {NULL}, cases_uniform, 0, 0x7fffffff}, 636 {"signbitf", (funcptr)test_signbitf, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 637 {"fpclassify", (funcptr)test_fpclassify, classify, {NULL}, cases_uniform, 0, 0x7fffffff}, 638 {"fpclassifyf", (funcptr)test_fpclassifyf, classifyf, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 639 /* 640 * Comparisons 641 */ 642 {"isgreater", (funcptr)test_isgreater, compare, {NULL}, cases_uniform, 0, 0x7fffffff}, 643 {"isgreaterequal", (funcptr)test_isgreaterequal, compare, {NULL}, cases_uniform, 0, 0x7fffffff}, 644 {"isless", (funcptr)test_isless, compare, {NULL}, cases_uniform, 0, 0x7fffffff}, 645 {"islessequal", (funcptr)test_islessequal, compare, {NULL}, cases_uniform, 0, 0x7fffffff}, 646 {"islessgreater", (funcptr)test_islessgreater, compare, {NULL}, cases_uniform, 0, 0x7fffffff}, 647 {"isunordered", (funcptr)test_isunordered, compare, {NULL}, cases_uniform, 0, 0x7fffffff}, 648 649 {"isgreaterf", (funcptr)test_isgreaterf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 650 {"isgreaterequalf", (funcptr)test_isgreaterequalf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 651 {"islessf", (funcptr)test_islessf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 652 {"islessequalf", (funcptr)test_islessequalf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 653 {"islessgreaterf", (funcptr)test_islessgreaterf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 654 {"isunorderedf", (funcptr)test_isunorderedf, comparef, {NULL}, cases_uniform_float, 0, 0x7fffffff}, 655 656 /* 657 * Inverse Hyperbolic functions 658 */ 659 {"atanh", (funcptr)mpfr_atanh, args1, {NULL}, cases_uniform, 0x3e400000, 0x3fefffff}, 660 {"asinh", (funcptr)mpfr_asinh, args1, {NULL}, cases_uniform, 0x3e400000, 0x3fefffff}, 661 {"acosh", (funcptr)mpfr_acosh, args1, {NULL}, cases_uniform_positive, 0x3ff00000, 0x7fefffff}, 662 663 {"atanhf", (funcptr)mpfr_atanh, args1f, {NULL}, cases_uniform_float, 0x32000000, 0x3f7fffff}, 664 {"asinhf", (funcptr)mpfr_asinh, args1f, {NULL}, cases_uniform_float, 0x32000000, 0x3f7fffff}, 665 {"acoshf", (funcptr)mpfr_acosh, args1f, {NULL}, cases_uniform_float_positive, 0x3f800000, 0x7f800000}, 666 667 /* 668 * Everything else (sitting in a section down here at the bottom 669 * because historically they were not tested because we didn't 670 * have reference implementations for them) 671 */ 672 {"csin", (funcptr)mpc_sin, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 673 {"csinf", (funcptr)mpc_sin, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 674 {"ccos", (funcptr)mpc_cos, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 675 {"ccosf", (funcptr)mpc_cos, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 676 {"ctan", (funcptr)mpc_tan, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 677 {"ctanf", (funcptr)mpc_tan, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 678 679 {"casin", (funcptr)mpc_asin, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 680 {"casinf", (funcptr)mpc_asin, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 681 {"cacos", (funcptr)mpc_acos, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 682 {"cacosf", (funcptr)mpc_acos, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 683 {"catan", (funcptr)mpc_atan, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 684 {"catanf", (funcptr)mpc_atan, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 685 686 {"csinh", (funcptr)mpc_sinh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 687 {"csinhf", (funcptr)mpc_sinh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 688 {"ccosh", (funcptr)mpc_cosh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 689 {"ccoshf", (funcptr)mpc_cosh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 690 {"ctanh", (funcptr)mpc_tanh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 691 {"ctanhf", (funcptr)mpc_tanh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 692 693 {"casinh", (funcptr)mpc_asinh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 694 {"casinhf", (funcptr)mpc_asinh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 695 {"cacosh", (funcptr)mpc_acosh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 696 {"cacoshf", (funcptr)mpc_acosh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 697 {"catanh", (funcptr)mpc_atanh, args1c, {NULL}, complex_cases_uniform, 0x3f000000, 0x40300000}, 698 {"catanhf", (funcptr)mpc_atanh, args1fc, {NULL}, complex_cases_uniform_float, 0x38000000, 0x41800000}, 699 700 {"cexp", (funcptr)mpc_exp, args1c, {NULL}, complex_cases_uniform, 0x3c900000, 0x40862000}, 701 {"cpow", (funcptr)test_cpow, args2c, {NULL}, complex_pow_cases, 0x3fc00000, 0x40000000}, 702 {"clog", (funcptr)mpc_log, args1c, {NULL}, complex_log_cases, 0, 0}, 703 {"csqrt", (funcptr)mpc_sqrt, args1c, {NULL}, complex_log_cases, 0, 0}, 704 705 {"cexpf", (funcptr)mpc_exp, args1fc, {NULL}, complex_cases_uniform_float, 0x24800000, 0x42b00000}, 706 {"cpowf", (funcptr)test_cpow, args2fc, {NULL}, complex_pow_cases_float, 0x3e000000, 0x41000000}, 707 {"clogf", (funcptr)mpc_log, args1fc, {NULL}, complex_log_cases_float, 0, 0}, 708 {"csqrtf", (funcptr)mpc_sqrt, args1fc, {NULL}, complex_log_cases_float, 0, 0}, 709 710 {"cdiv", (funcptr)mpc_div, args2c, {NULL}, complex_arithmetic_cases, 0, 0}, 711 {"cmul", (funcptr)mpc_mul, args2c, {NULL}, complex_arithmetic_cases, 0, 0}, 712 {"cadd", (funcptr)mpc_add, args2c, {NULL}, complex_arithmetic_cases, 0, 0}, 713 {"csub", (funcptr)mpc_sub, args2c, {NULL}, complex_arithmetic_cases, 0, 0}, 714 715 {"cdivf", (funcptr)mpc_div, args2fc, {NULL}, complex_arithmetic_cases_float, 0, 0}, 716 {"cmulf", (funcptr)mpc_mul, args2fc, {NULL}, complex_arithmetic_cases_float, 0, 0}, 717 {"caddf", (funcptr)mpc_add, args2fc, {NULL}, complex_arithmetic_cases_float, 0, 0}, 718 {"csubf", (funcptr)mpc_sub, args2fc, {NULL}, complex_arithmetic_cases_float, 0, 0}, 719 720 {"cabsf", (funcptr)mpc_abs, args1fcr, {NULL}, complex_arithmetic_cases_float, 0, 0}, 721 {"cabs", (funcptr)mpc_abs, args1cr, {NULL}, complex_arithmetic_cases, 0, 0}, 722 {"cargf", (funcptr)mpc_arg, args1fcr, {NULL}, complex_arithmetic_cases_float, 0, 0}, 723 {"carg", (funcptr)mpc_arg, args1cr, {NULL}, complex_arithmetic_cases, 0, 0}, 724 {"cimagf", (funcptr)mpc_imag, args1fcr, {NULL}, complex_arithmetic_cases_float, 0, 0}, 725 {"cimag", (funcptr)mpc_imag, args1cr, {NULL}, complex_arithmetic_cases, 0, 0}, 726 {"conjf", (funcptr)mpc_conj, args1fc, {NULL}, complex_arithmetic_cases_float, 0, 0}, 727 {"conj", (funcptr)mpc_conj, args1c, {NULL}, complex_arithmetic_cases, 0, 0}, 728 {"cprojf", (funcptr)mpc_proj, args1fc, {NULL}, complex_arithmetic_cases_float, 0, 0}, 729 {"cproj", (funcptr)mpc_proj, args1c, {NULL}, complex_arithmetic_cases, 0, 0}, 730 {"crealf", (funcptr)mpc_real, args1fcr, {NULL}, complex_arithmetic_cases_float, 0, 0}, 731 {"creal", (funcptr)mpc_real, args1cr, {NULL}, complex_arithmetic_cases, 0, 0}, 732 {"erfcf", (funcptr)mpfr_erfc, args1f, {NULL}, cases_biased_float, 0x1e800000, 0x41000000}, 733 {"erfc", (funcptr)mpfr_erfc, args1, {NULL}, cases_biased, 0x3bd00000, 0x403c0000}, 734 {"erff", (funcptr)mpfr_erf, args1f, {NULL}, cases_biased_float, 0x03800000, 0x40700000}, 735 {"erf", (funcptr)mpfr_erf, args1, {NULL}, cases_biased, 0x00800000, 0x40200000}, 736 {"exp2f", (funcptr)mpfr_exp2, args1f, {NULL}, cases_uniform_float, 0x33800000, 0x43c00000}, 737 {"exp2", (funcptr)mpfr_exp2, args1, {NULL}, cases_uniform, 0x3ca00000, 0x40a00000}, 738 {"expm1f", (funcptr)mpfr_expm1, args1f, {NULL}, cases_uniform_float, 0x33000000, 0x43800000}, 739 {"expm1", (funcptr)mpfr_expm1, args1, {NULL}, cases_uniform, 0x3c900000, 0x409c0000}, 740 {"fmaxf", (funcptr)mpfr_max, args2f, {NULL}, minmax_cases_float, 0, 0x7f7fffff}, 741 {"fmax", (funcptr)mpfr_max, args2, {NULL}, minmax_cases, 0, 0x7fefffff}, 742 {"fminf", (funcptr)mpfr_min, args2f, {NULL}, minmax_cases_float, 0, 0x7f7fffff}, 743 {"fmin", (funcptr)mpfr_min, args2, {NULL}, minmax_cases, 0, 0x7fefffff}, 744 {"lgammaf", (funcptr)test_lgamma, args1f, {NULL}, cases_uniform_float, 0x01800000, 0x7f800000}, 745 {"lgamma", (funcptr)test_lgamma, args1, {NULL}, cases_uniform, 0x00100000, 0x7ff00000}, 746 {"log1pf", (funcptr)mpfr_log1p, args1f, {NULL}, log1p_cases_float, 0, 0}, 747 {"log1p", (funcptr)mpfr_log1p, args1, {NULL}, log1p_cases, 0, 0}, 748 {"log2f", (funcptr)mpfr_log2, args1f, {NULL}, log_cases_float, 0, 0}, 749 {"log2", (funcptr)mpfr_log2, args1, {NULL}, log_cases, 0, 0}, 750 {"tgammaf", (funcptr)mpfr_gamma, args1f, {NULL}, cases_uniform_float, 0x2f800000, 0x43000000}, 751 {"tgamma", (funcptr)mpfr_gamma, args1, {NULL}, cases_uniform, 0x3c000000, 0x40800000}, 752 }; 753 /* clang-format on */ 754 755 const int nfunctions = ( sizeof(functions)/sizeof(*functions) ); 756 757 #define random_sign ( random_upto(1) ? 0x80000000 : 0 ) 758 759 static int iszero(uint32 *x) { 760 return !((x[0] & 0x7FFFFFFF) || x[1]); 761 } 762 763 764 static void complex_log_cases(uint32 *out, uint32 param1, 765 uint32 param2) { 766 cases_uniform(out,0x00100000,0x7fefffff); 767 cases_uniform(out+2,0x00100000,0x7fefffff); 768 } 769 770 771 static void complex_log_cases_float(uint32 *out, uint32 param1, 772 uint32 param2) { 773 cases_uniform_float(out,0x00800000,0x7f7fffff); 774 cases_uniform_float(out+2,0x00800000,0x7f7fffff); 775 } 776 777 static void complex_cases_biased(uint32 *out, uint32 lowbound, 778 uint32 highbound) { 779 cases_biased(out,lowbound,highbound); 780 cases_biased(out+2,lowbound,highbound); 781 } 782 783 static void complex_cases_biased_float(uint32 *out, uint32 lowbound, 784 uint32 highbound) { 785 cases_biased_float(out,lowbound,highbound); 786 cases_biased_float(out+2,lowbound,highbound); 787 } 788 789 static void complex_cases_uniform(uint32 *out, uint32 lowbound, 790 uint32 highbound) { 791 cases_uniform(out,lowbound,highbound); 792 cases_uniform(out+2,lowbound,highbound); 793 } 794 795 static void complex_cases_uniform_float(uint32 *out, uint32 lowbound, 796 uint32 highbound) { 797 cases_uniform_float(out,lowbound,highbound); 798 cases_uniform(out+2,lowbound,highbound); 799 } 800 801 static void complex_pow_cases(uint32 *out, uint32 lowbound, 802 uint32 highbound) { 803 /* 804 * Generating non-overflowing cases for complex pow: 805 * 806 * Our base has both parts within the range [1/2,2], and hence 807 * its magnitude is within [1/2,2*sqrt(2)]. The magnitude of its 808 * logarithm in base 2 is therefore at most the magnitude of 809 * (log2(2*sqrt(2)) + i*pi/log(2)), or in other words 810 * hypot(3/2,pi/log(2)) = 4.77. So the magnitude of the exponent 811 * input must be at most our output magnitude limit (as a power 812 * of two) divided by that. 813 * 814 * I also set the output magnitude limit a bit low, because we 815 * don't guarantee (and neither does glibc) to prevent internal 816 * overflow in cases where the output _magnitude_ overflows but 817 * scaling it back down by cos and sin of the argument brings it 818 * back in range. 819 */ 820 cases_uniform(out,0x3fe00000, 0x40000000); 821 cases_uniform(out+2,0x3fe00000, 0x40000000); 822 cases_uniform(out+4,0x3f800000, 0x40600000); 823 cases_uniform(out+6,0x3f800000, 0x40600000); 824 } 825 826 static void complex_pow_cases_float(uint32 *out, uint32 lowbound, 827 uint32 highbound) { 828 /* 829 * Reasoning as above, though of course the detailed numbers are 830 * all different. 831 */ 832 cases_uniform_float(out,0x3f000000, 0x40000000); 833 cases_uniform_float(out+2,0x3f000000, 0x40000000); 834 cases_uniform_float(out+4,0x3d600000, 0x41900000); 835 cases_uniform_float(out+6,0x3d600000, 0x41900000); 836 } 837 838 static void complex_arithmetic_cases(uint32 *out, uint32 lowbound, 839 uint32 highbound) { 840 cases_uniform(out,0,0x7fefffff); 841 cases_uniform(out+2,0,0x7fefffff); 842 cases_uniform(out+4,0,0x7fefffff); 843 cases_uniform(out+6,0,0x7fefffff); 844 } 845 846 static void complex_arithmetic_cases_float(uint32 *out, uint32 lowbound, 847 uint32 highbound) { 848 cases_uniform_float(out,0,0x7f7fffff); 849 cases_uniform_float(out+2,0,0x7f7fffff); 850 cases_uniform_float(out+4,0,0x7f7fffff); 851 cases_uniform_float(out+6,0,0x7f7fffff); 852 } 853 854 /* 855 * Included from fplib test suite, in a compact self-contained 856 * form. 857 */ 858 859 void float32_case(uint32 *ret) { 860 int n, bits; 861 uint32 f; 862 static int premax, preptr; 863 static uint32 *specifics = NULL; 864 865 if (!ret) { 866 if (specifics) 867 free(specifics); 868 specifics = NULL; 869 premax = preptr = 0; 870 return; 871 } 872 873 if (!specifics) { 874 int exps[] = { 875 -127, -126, -125, -24, -4, -3, -2, -1, 0, 1, 2, 3, 4, 876 24, 29, 30, 31, 32, 61, 62, 63, 64, 126, 127, 128 877 }; 878 int sign, eptr; 879 uint32 se, j; 880 /* 881 * We want a cross product of: 882 * - each of two sign bits (2) 883 * - each of the above (unbiased) exponents (25) 884 * - the following list of fraction parts: 885 * * zero (1) 886 * * all bits (1) 887 * * one-bit-set (23) 888 * * one-bit-clear (23) 889 * * one-bit-and-above (20: 3 are duplicates) 890 * * one-bit-and-below (20: 3 are duplicates) 891 * (total 88) 892 * (total 4400) 893 */ 894 specifics = malloc(4400 * sizeof(*specifics)); 895 preptr = 0; 896 for (sign = 0; sign <= 1; sign++) { 897 for (eptr = 0; eptr < sizeof(exps)/sizeof(*exps); eptr++) { 898 se = (sign ? 0x80000000 : 0) | ((exps[eptr]+127) << 23); 899 /* 900 * Zero. 901 */ 902 specifics[preptr++] = se | 0; 903 /* 904 * All bits. 905 */ 906 specifics[preptr++] = se | 0x7FFFFF; 907 /* 908 * One-bit-set. 909 */ 910 for (j = 1; j && j <= 0x400000; j <<= 1) 911 specifics[preptr++] = se | j; 912 /* 913 * One-bit-clear. 914 */ 915 for (j = 1; j && j <= 0x400000; j <<= 1) 916 specifics[preptr++] = se | (0x7FFFFF ^ j); 917 /* 918 * One-bit-and-everything-below. 919 */ 920 for (j = 2; j && j <= 0x100000; j <<= 1) 921 specifics[preptr++] = se | (2*j-1); 922 /* 923 * One-bit-and-everything-above. 924 */ 925 for (j = 4; j && j <= 0x200000; j <<= 1) 926 specifics[preptr++] = se | (0x7FFFFF ^ (j-1)); 927 /* 928 * Done. 929 */ 930 } 931 } 932 assert(preptr == 4400); 933 premax = preptr; 934 } 935 936 /* 937 * Decide whether to return a pre or a random case. 938 */ 939 n = random32() % (premax+1); 940 if (n < preptr) { 941 /* 942 * Return pre[n]. 943 */ 944 uint32 t; 945 t = specifics[n]; 946 specifics[n] = specifics[preptr-1]; 947 specifics[preptr-1] = t; /* (not really needed) */ 948 preptr--; 949 *ret = t; 950 } else { 951 /* 952 * Random case. 953 * Sign and exponent: 954 * - FIXME 955 * Significand: 956 * - with prob 1/5, a totally random bit pattern 957 * - with prob 1/5, all 1s down to some point and then random 958 * - with prob 1/5, all 1s up to some point and then random 959 * - with prob 1/5, all 0s down to some point and then random 960 * - with prob 1/5, all 0s up to some point and then random 961 */ 962 n = random32() % 5; 963 f = random32(); /* some random bits */ 964 bits = random32() % 22 + 1; /* 1-22 */ 965 switch (n) { 966 case 0: 967 break; /* leave f alone */ 968 case 1: 969 f |= (1<<bits)-1; 970 break; 971 case 2: 972 f &= ~((1<<bits)-1); 973 break; 974 case 3: 975 f |= ~((1<<bits)-1); 976 break; 977 case 4: 978 f &= (1<<bits)-1; 979 break; 980 } 981 f &= 0x7FFFFF; 982 f |= (random32() & 0xFF800000);/* FIXME - do better */ 983 *ret = f; 984 } 985 } 986 static void float64_case(uint32 *ret) { 987 int n, bits; 988 uint32 f, g; 989 static int premax, preptr; 990 static uint32 (*specifics)[2] = NULL; 991 992 if (!ret) { 993 if (specifics) 994 free(specifics); 995 specifics = NULL; 996 premax = preptr = 0; 997 return; 998 } 999 1000 if (!specifics) { 1001 int exps[] = { 1002 -1023, -1022, -1021, -129, -128, -127, -126, -53, -4, -3, -2, 1003 -1, 0, 1, 2, 3, 4, 29, 30, 31, 32, 53, 61, 62, 63, 64, 127, 1004 128, 129, 1022, 1023, 1024 1005 }; 1006 int sign, eptr; 1007 uint32 se, j; 1008 /* 1009 * We want a cross product of: 1010 * - each of two sign bits (2) 1011 * - each of the above (unbiased) exponents (32) 1012 * - the following list of fraction parts: 1013 * * zero (1) 1014 * * all bits (1) 1015 * * one-bit-set (52) 1016 * * one-bit-clear (52) 1017 * * one-bit-and-above (49: 3 are duplicates) 1018 * * one-bit-and-below (49: 3 are duplicates) 1019 * (total 204) 1020 * (total 13056) 1021 */ 1022 specifics = malloc(13056 * sizeof(*specifics)); 1023 preptr = 0; 1024 for (sign = 0; sign <= 1; sign++) { 1025 for (eptr = 0; eptr < sizeof(exps)/sizeof(*exps); eptr++) { 1026 se = (sign ? 0x80000000 : 0) | ((exps[eptr]+1023) << 20); 1027 /* 1028 * Zero. 1029 */ 1030 specifics[preptr][0] = 0; 1031 specifics[preptr][1] = 0; 1032 specifics[preptr++][0] |= se; 1033 /* 1034 * All bits. 1035 */ 1036 specifics[preptr][0] = 0xFFFFF; 1037 specifics[preptr][1] = ~0; 1038 specifics[preptr++][0] |= se; 1039 /* 1040 * One-bit-set. 1041 */ 1042 for (j = 1; j && j <= 0x80000000; j <<= 1) { 1043 specifics[preptr][0] = 0; 1044 specifics[preptr][1] = j; 1045 specifics[preptr++][0] |= se; 1046 if (j & 0xFFFFF) { 1047 specifics[preptr][0] = j; 1048 specifics[preptr][1] = 0; 1049 specifics[preptr++][0] |= se; 1050 } 1051 } 1052 /* 1053 * One-bit-clear. 1054 */ 1055 for (j = 1; j && j <= 0x80000000; j <<= 1) { 1056 specifics[preptr][0] = 0xFFFFF; 1057 specifics[preptr][1] = ~j; 1058 specifics[preptr++][0] |= se; 1059 if (j & 0xFFFFF) { 1060 specifics[preptr][0] = 0xFFFFF ^ j; 1061 specifics[preptr][1] = ~0; 1062 specifics[preptr++][0] |= se; 1063 } 1064 } 1065 /* 1066 * One-bit-and-everything-below. 1067 */ 1068 for (j = 2; j && j <= 0x80000000; j <<= 1) { 1069 specifics[preptr][0] = 0; 1070 specifics[preptr][1] = 2*j-1; 1071 specifics[preptr++][0] |= se; 1072 } 1073 for (j = 1; j && j <= 0x20000; j <<= 1) { 1074 specifics[preptr][0] = 2*j-1; 1075 specifics[preptr][1] = ~0; 1076 specifics[preptr++][0] |= se; 1077 } 1078 /* 1079 * One-bit-and-everything-above. 1080 */ 1081 for (j = 4; j && j <= 0x80000000; j <<= 1) { 1082 specifics[preptr][0] = 0xFFFFF; 1083 specifics[preptr][1] = ~(j-1); 1084 specifics[preptr++][0] |= se; 1085 } 1086 for (j = 1; j && j <= 0x40000; j <<= 1) { 1087 specifics[preptr][0] = 0xFFFFF ^ (j-1); 1088 specifics[preptr][1] = 0; 1089 specifics[preptr++][0] |= se; 1090 } 1091 /* 1092 * Done. 1093 */ 1094 } 1095 } 1096 assert(preptr == 13056); 1097 premax = preptr; 1098 } 1099 1100 /* 1101 * Decide whether to return a pre or a random case. 1102 */ 1103 n = (uint32) random32() % (uint32) (premax+1); 1104 if (n < preptr) { 1105 /* 1106 * Return pre[n]. 1107 */ 1108 uint32 t; 1109 t = specifics[n][0]; 1110 specifics[n][0] = specifics[preptr-1][0]; 1111 specifics[preptr-1][0] = t; /* (not really needed) */ 1112 ret[0] = t; 1113 t = specifics[n][1]; 1114 specifics[n][1] = specifics[preptr-1][1]; 1115 specifics[preptr-1][1] = t; /* (not really needed) */ 1116 ret[1] = t; 1117 preptr--; 1118 } else { 1119 /* 1120 * Random case. 1121 * Sign and exponent: 1122 * - FIXME 1123 * Significand: 1124 * - with prob 1/5, a totally random bit pattern 1125 * - with prob 1/5, all 1s down to some point and then random 1126 * - with prob 1/5, all 1s up to some point and then random 1127 * - with prob 1/5, all 0s down to some point and then random 1128 * - with prob 1/5, all 0s up to some point and then random 1129 */ 1130 n = random32() % 5; 1131 f = random32(); /* some random bits */ 1132 g = random32(); /* some random bits */ 1133 bits = random32() % 51 + 1; /* 1-51 */ 1134 switch (n) { 1135 case 0: 1136 break; /* leave f alone */ 1137 case 1: 1138 if (bits <= 32) 1139 f |= (1<<bits)-1; 1140 else { 1141 bits -= 32; 1142 g |= (1<<bits)-1; 1143 f = ~0; 1144 } 1145 break; 1146 case 2: 1147 if (bits <= 32) 1148 f &= ~((1<<bits)-1); 1149 else { 1150 bits -= 32; 1151 g &= ~((1<<bits)-1); 1152 f = 0; 1153 } 1154 break; 1155 case 3: 1156 if (bits <= 32) 1157 g &= (1<<bits)-1; 1158 else { 1159 bits -= 32; 1160 f &= (1<<bits)-1; 1161 g = 0; 1162 } 1163 break; 1164 case 4: 1165 if (bits <= 32) 1166 g |= ~((1<<bits)-1); 1167 else { 1168 bits -= 32; 1169 f |= ~((1<<bits)-1); 1170 g = ~0; 1171 } 1172 break; 1173 } 1174 g &= 0xFFFFF; 1175 g |= (random32() & 0xFFF00000);/* FIXME - do better */ 1176 ret[0] = g; 1177 ret[1] = f; 1178 } 1179 } 1180 1181 static void cases_biased(uint32 *out, uint32 lowbound, 1182 uint32 highbound) { 1183 do { 1184 out[0] = highbound - random_upto_biased(highbound-lowbound, 8); 1185 out[1] = random_upto(0xFFFFFFFF); 1186 out[0] |= random_sign; 1187 } while (iszero(out)); /* rule out zero */ 1188 } 1189 1190 static void cases_biased_positive(uint32 *out, uint32 lowbound, 1191 uint32 highbound) { 1192 do { 1193 out[0] = highbound - random_upto_biased(highbound-lowbound, 8); 1194 out[1] = random_upto(0xFFFFFFFF); 1195 } while (iszero(out)); /* rule out zero */ 1196 } 1197 1198 static void cases_biased_float(uint32 *out, uint32 lowbound, 1199 uint32 highbound) { 1200 do { 1201 out[0] = highbound - random_upto_biased(highbound-lowbound, 8); 1202 out[1] = 0; 1203 out[0] |= random_sign; 1204 } while (iszero(out)); /* rule out zero */ 1205 } 1206 1207 static void cases_semi1(uint32 *out, uint32 param1, 1208 uint32 param2) { 1209 float64_case(out); 1210 } 1211 1212 static void cases_semi1_float(uint32 *out, uint32 param1, 1213 uint32 param2) { 1214 float32_case(out); 1215 } 1216 1217 static void cases_semi2(uint32 *out, uint32 param1, 1218 uint32 param2) { 1219 float64_case(out); 1220 float64_case(out+2); 1221 } 1222 1223 static void cases_semi2_float(uint32 *out, uint32 param1, 1224 uint32 param2) { 1225 float32_case(out); 1226 float32_case(out+2); 1227 } 1228 1229 static void cases_ldexp(uint32 *out, uint32 param1, 1230 uint32 param2) { 1231 float64_case(out); 1232 out[2] = random_upto(2048)-1024; 1233 } 1234 1235 static void cases_ldexp_float(uint32 *out, uint32 param1, 1236 uint32 param2) { 1237 float32_case(out); 1238 out[2] = random_upto(256)-128; 1239 } 1240 1241 static void cases_uniform(uint32 *out, uint32 lowbound, 1242 uint32 highbound) { 1243 do { 1244 out[0] = highbound - random_upto(highbound-lowbound); 1245 out[1] = random_upto(0xFFFFFFFF); 1246 out[0] |= random_sign; 1247 } while (iszero(out)); /* rule out zero */ 1248 } 1249 static void cases_uniform_float(uint32 *out, uint32 lowbound, 1250 uint32 highbound) { 1251 do { 1252 out[0] = highbound - random_upto(highbound-lowbound); 1253 out[1] = 0; 1254 out[0] |= random_sign; 1255 } while (iszero(out)); /* rule out zero */ 1256 } 1257 1258 static void cases_uniform_positive(uint32 *out, uint32 lowbound, 1259 uint32 highbound) { 1260 do { 1261 out[0] = highbound - random_upto(highbound-lowbound); 1262 out[1] = random_upto(0xFFFFFFFF); 1263 } while (iszero(out)); /* rule out zero */ 1264 } 1265 static void cases_uniform_float_positive(uint32 *out, uint32 lowbound, 1266 uint32 highbound) { 1267 do { 1268 out[0] = highbound - random_upto(highbound-lowbound); 1269 out[1] = 0; 1270 } while (iszero(out)); /* rule out zero */ 1271 } 1272 1273 1274 static void log_cases(uint32 *out, uint32 param1, 1275 uint32 param2) { 1276 do { 1277 out[0] = random_upto(0x7FEFFFFF); 1278 out[1] = random_upto(0xFFFFFFFF); 1279 } while (iszero(out)); /* rule out zero */ 1280 } 1281 1282 static void log_cases_float(uint32 *out, uint32 param1, 1283 uint32 param2) { 1284 do { 1285 out[0] = random_upto(0x7F7FFFFF); 1286 out[1] = 0; 1287 } while (iszero(out)); /* rule out zero */ 1288 } 1289 1290 static void log1p_cases(uint32 *out, uint32 param1, uint32 param2) 1291 { 1292 uint32 sign = random_sign; 1293 if (sign == 0) { 1294 cases_uniform_positive(out, 0x3c700000, 0x43400000); 1295 } else { 1296 cases_uniform_positive(out, 0x3c000000, 0x3ff00000); 1297 } 1298 out[0] |= sign; 1299 } 1300 1301 static void log1p_cases_float(uint32 *out, uint32 param1, uint32 param2) 1302 { 1303 uint32 sign = random_sign; 1304 if (sign == 0) { 1305 cases_uniform_float_positive(out, 0x32000000, 0x4c000000); 1306 } else { 1307 cases_uniform_float_positive(out, 0x30000000, 0x3f800000); 1308 } 1309 out[0] |= sign; 1310 } 1311 1312 static void minmax_cases(uint32 *out, uint32 param1, uint32 param2) 1313 { 1314 do { 1315 out[0] = random_upto(0x7FEFFFFF); 1316 out[1] = random_upto(0xFFFFFFFF); 1317 out[0] |= random_sign; 1318 out[2] = random_upto(0x7FEFFFFF); 1319 out[3] = random_upto(0xFFFFFFFF); 1320 out[2] |= random_sign; 1321 } while (iszero(out)); /* rule out zero */ 1322 } 1323 1324 static void minmax_cases_float(uint32 *out, uint32 param1, uint32 param2) 1325 { 1326 do { 1327 out[0] = random_upto(0x7F7FFFFF); 1328 out[1] = 0; 1329 out[0] |= random_sign; 1330 out[2] = random_upto(0x7F7FFFFF); 1331 out[3] = 0; 1332 out[2] |= random_sign; 1333 } while (iszero(out)); /* rule out zero */ 1334 } 1335 1336 static void rred_cases(uint32 *out, uint32 param1, 1337 uint32 param2) { 1338 do { 1339 out[0] = ((0x3fc00000 + random_upto(0x036fffff)) | 1340 (random_upto(1) << 31)); 1341 out[1] = random_upto(0xFFFFFFFF); 1342 } while (iszero(out)); /* rule out zero */ 1343 } 1344 1345 static void rred_cases_float(uint32 *out, uint32 param1, 1346 uint32 param2) { 1347 do { 1348 out[0] = ((0x3e000000 + random_upto(0x0cffffff)) | 1349 (random_upto(1) << 31)); 1350 out[1] = 0; /* for iszero */ 1351 } while (iszero(out)); /* rule out zero */ 1352 } 1353 1354 static void atan2_cases(uint32 *out, uint32 param1, 1355 uint32 param2) { 1356 do { 1357 int expdiff = random_upto(101)-51; 1358 int swap; 1359 if (expdiff < 0) { 1360 expdiff = -expdiff; 1361 swap = 2; 1362 } else 1363 swap = 0; 1364 out[swap ^ 0] = random_upto(0x7FEFFFFF-((expdiff+1)<<20)); 1365 out[swap ^ 2] = random_upto(((expdiff+1)<<20)-1) + out[swap ^ 0]; 1366 out[1] = random_upto(0xFFFFFFFF); 1367 out[3] = random_upto(0xFFFFFFFF); 1368 out[0] |= random_sign; 1369 out[2] |= random_sign; 1370 } while (iszero(out) || iszero(out+2));/* rule out zero */ 1371 } 1372 1373 static void atan2_cases_float(uint32 *out, uint32 param1, 1374 uint32 param2) { 1375 do { 1376 int expdiff = random_upto(44)-22; 1377 int swap; 1378 if (expdiff < 0) { 1379 expdiff = -expdiff; 1380 swap = 2; 1381 } else 1382 swap = 0; 1383 out[swap ^ 0] = random_upto(0x7F7FFFFF-((expdiff+1)<<23)); 1384 out[swap ^ 2] = random_upto(((expdiff+1)<<23)-1) + out[swap ^ 0]; 1385 out[0] |= random_sign; 1386 out[2] |= random_sign; 1387 out[1] = out[3] = 0; /* for iszero */ 1388 } while (iszero(out) || iszero(out+2));/* rule out zero */ 1389 } 1390 1391 static void pow_cases(uint32 *out, uint32 param1, 1392 uint32 param2) { 1393 /* 1394 * Pick an exponent e (-0x33 to +0x7FE) for x, and here's the 1395 * range of numbers we can use as y: 1396 * 1397 * For e < 0x3FE, the range is [-0x400/(0x3FE-e),+0x432/(0x3FE-e)] 1398 * For e > 0x3FF, the range is [-0x432/(e-0x3FF),+0x400/(e-0x3FF)] 1399 * 1400 * For e == 0x3FE or e == 0x3FF, the range gets infinite at one 1401 * end or the other, so we have to be cleverer: pick a number n 1402 * of useful bits in the mantissa (1 thru 52, so 1 must imply 1403 * 0x3ff00000.00000001 whereas 52 is anything at least as big 1404 * as 0x3ff80000.00000000; for e == 0x3fe, 1 necessarily means 1405 * 0x3fefffff.ffffffff and 52 is anything at most as big as 1406 * 0x3fe80000.00000000). Then, as it happens, a sensible 1407 * maximum power is 2^(63-n) for e == 0x3fe, and 2^(62-n) for 1408 * e == 0x3ff. 1409 * 1410 * We inevitably get some overflows in approximating the log 1411 * curves by these nasty step functions, but that's all right - 1412 * we do want _some_ overflows to be tested. 1413 * 1414 * Having got that, then, it's just a matter of inventing a 1415 * probability distribution for all of this. 1416 */ 1417 int e, n; 1418 uint32 dmin, dmax; 1419 const uint32 pmin = 0x3e100000; 1420 1421 /* 1422 * Generate exponents in a slightly biased fashion. 1423 */ 1424 e = (random_upto(1) ? /* is exponent small or big? */ 1425 0x3FE - random_upto_biased(0x431,2) : /* small */ 1426 0x3FF + random_upto_biased(0x3FF,2)); /* big */ 1427 1428 /* 1429 * Now split into cases. 1430 */ 1431 if (e < 0x3FE || e > 0x3FF) { 1432 uint32 imin, imax; 1433 if (e < 0x3FE) 1434 imin = 0x40000 / (0x3FE - e), imax = 0x43200 / (0x3FE - e); 1435 else 1436 imin = 0x43200 / (e - 0x3FF), imax = 0x40000 / (e - 0x3FF); 1437 /* Power range runs from -imin to imax. Now convert to doubles */ 1438 dmin = doubletop(imin, -8); 1439 dmax = doubletop(imax, -8); 1440 /* Compute the number of mantissa bits. */ 1441 n = (e > 0 ? 53 : 52+e); 1442 } else { 1443 /* Critical exponents. Generate a top bit index. */ 1444 n = 52 - random_upto_biased(51, 4); 1445 if (e == 0x3FE) 1446 dmax = 63 - n; 1447 else 1448 dmax = 62 - n; 1449 dmax = (dmax << 20) + 0x3FF00000; 1450 dmin = dmax; 1451 } 1452 /* Generate a mantissa. */ 1453 if (n <= 32) { 1454 out[0] = 0; 1455 out[1] = random_upto((1 << (n-1)) - 1) + (1 << (n-1)); 1456 } else if (n == 33) { 1457 out[0] = 1; 1458 out[1] = random_upto(0xFFFFFFFF); 1459 } else if (n > 33) { 1460 out[0] = random_upto((1 << (n-33)) - 1) + (1 << (n-33)); 1461 out[1] = random_upto(0xFFFFFFFF); 1462 } 1463 /* Negate the mantissa if e == 0x3FE. */ 1464 if (e == 0x3FE) { 1465 out[1] = -out[1]; 1466 out[0] = -out[0]; 1467 if (out[1]) out[0]--; 1468 } 1469 /* Put the exponent on. */ 1470 out[0] &= 0xFFFFF; 1471 out[0] |= ((e > 0 ? e : 0) << 20); 1472 /* Generate a power. Powers don't go below 2^-30. */ 1473 if (random_upto(1)) { 1474 /* Positive power */ 1475 out[2] = dmax - random_upto_biased(dmax-pmin, 10); 1476 } else { 1477 /* Negative power */ 1478 out[2] = (dmin - random_upto_biased(dmin-pmin, 10)) | 0x80000000; 1479 } 1480 out[3] = random_upto(0xFFFFFFFF); 1481 } 1482 static void pow_cases_float(uint32 *out, uint32 param1, 1483 uint32 param2) { 1484 /* 1485 * Pick an exponent e (-0x16 to +0xFE) for x, and here's the 1486 * range of numbers we can use as y: 1487 * 1488 * For e < 0x7E, the range is [-0x80/(0x7E-e),+0x95/(0x7E-e)] 1489 * For e > 0x7F, the range is [-0x95/(e-0x7F),+0x80/(e-0x7F)] 1490 * 1491 * For e == 0x7E or e == 0x7F, the range gets infinite at one 1492 * end or the other, so we have to be cleverer: pick a number n 1493 * of useful bits in the mantissa (1 thru 23, so 1 must imply 1494 * 0x3f800001 whereas 23 is anything at least as big as 1495 * 0x3fc00000; for e == 0x7e, 1 necessarily means 0x3f7fffff 1496 * and 23 is anything at most as big as 0x3f400000). Then, as 1497 * it happens, a sensible maximum power is 2^(31-n) for e == 1498 * 0x7e, and 2^(30-n) for e == 0x7f. 1499 * 1500 * We inevitably get some overflows in approximating the log 1501 * curves by these nasty step functions, but that's all right - 1502 * we do want _some_ overflows to be tested. 1503 * 1504 * Having got that, then, it's just a matter of inventing a 1505 * probability distribution for all of this. 1506 */ 1507 int e, n; 1508 uint32 dmin, dmax; 1509 const uint32 pmin = 0x38000000; 1510 1511 /* 1512 * Generate exponents in a slightly biased fashion. 1513 */ 1514 e = (random_upto(1) ? /* is exponent small or big? */ 1515 0x7E - random_upto_biased(0x94,2) : /* small */ 1516 0x7F + random_upto_biased(0x7f,2)); /* big */ 1517 1518 /* 1519 * Now split into cases. 1520 */ 1521 if (e < 0x7E || e > 0x7F) { 1522 uint32 imin, imax; 1523 if (e < 0x7E) 1524 imin = 0x8000 / (0x7e - e), imax = 0x9500 / (0x7e - e); 1525 else 1526 imin = 0x9500 / (e - 0x7f), imax = 0x8000 / (e - 0x7f); 1527 /* Power range runs from -imin to imax. Now convert to doubles */ 1528 dmin = floatval(imin, -8); 1529 dmax = floatval(imax, -8); 1530 /* Compute the number of mantissa bits. */ 1531 n = (e > 0 ? 24 : 23+e); 1532 } else { 1533 /* Critical exponents. Generate a top bit index. */ 1534 n = 23 - random_upto_biased(22, 4); 1535 if (e == 0x7E) 1536 dmax = 31 - n; 1537 else 1538 dmax = 30 - n; 1539 dmax = (dmax << 23) + 0x3F800000; 1540 dmin = dmax; 1541 } 1542 /* Generate a mantissa. */ 1543 out[0] = random_upto((1 << (n-1)) - 1) + (1 << (n-1)); 1544 out[1] = 0; 1545 /* Negate the mantissa if e == 0x7E. */ 1546 if (e == 0x7E) { 1547 out[0] = -out[0]; 1548 } 1549 /* Put the exponent on. */ 1550 out[0] &= 0x7FFFFF; 1551 out[0] |= ((e > 0 ? e : 0) << 23); 1552 /* Generate a power. Powers don't go below 2^-15. */ 1553 if (random_upto(1)) { 1554 /* Positive power */ 1555 out[2] = dmax - random_upto_biased(dmax-pmin, 10); 1556 } else { 1557 /* Negative power */ 1558 out[2] = (dmin - random_upto_biased(dmin-pmin, 10)) | 0x80000000; 1559 } 1560 out[3] = 0; 1561 } 1562 1563 void vet_for_decline(Testable *fn, uint32 *args, uint32 *result, int got_errno_in) { 1564 int declined = 0; 1565 1566 switch (fn->type) { 1567 case args1: 1568 case rred: 1569 case semi1: 1570 case t_frexp: 1571 case t_modf: 1572 case classify: 1573 case t_ldexp: 1574 declined |= lib_fo && is_dhard(args+0); 1575 break; 1576 case args1f: 1577 case rredf: 1578 case semi1f: 1579 case t_frexpf: 1580 case t_modff: 1581 case classifyf: 1582 declined |= lib_fo && is_shard(args+0); 1583 break; 1584 case args2: 1585 case semi2: 1586 case args1c: 1587 case args1cr: 1588 case compare: 1589 declined |= lib_fo && is_dhard(args+0); 1590 declined |= lib_fo && is_dhard(args+2); 1591 break; 1592 case args2f: 1593 case semi2f: 1594 case t_ldexpf: 1595 case comparef: 1596 case args1fc: 1597 case args1fcr: 1598 declined |= lib_fo && is_shard(args+0); 1599 declined |= lib_fo && is_shard(args+2); 1600 break; 1601 case args2c: 1602 declined |= lib_fo && is_dhard(args+0); 1603 declined |= lib_fo && is_dhard(args+2); 1604 declined |= lib_fo && is_dhard(args+4); 1605 declined |= lib_fo && is_dhard(args+6); 1606 break; 1607 case args2fc: 1608 declined |= lib_fo && is_shard(args+0); 1609 declined |= lib_fo && is_shard(args+2); 1610 declined |= lib_fo && is_shard(args+4); 1611 declined |= lib_fo && is_shard(args+6); 1612 break; 1613 } 1614 1615 switch (fn->type) { 1616 case args1: /* return an extra-precise result */ 1617 case args2: 1618 case rred: 1619 case semi1: /* return a double result */ 1620 case semi2: 1621 case t_ldexp: 1622 case t_frexp: /* return double * int */ 1623 case args1cr: 1624 declined |= lib_fo && is_dhard(result); 1625 break; 1626 case args1f: 1627 case args2f: 1628 case rredf: 1629 case semi1f: 1630 case semi2f: 1631 case t_ldexpf: 1632 case args1fcr: 1633 declined |= lib_fo && is_shard(result); 1634 break; 1635 case t_modf: /* return double * double */ 1636 declined |= lib_fo && is_dhard(result+0); 1637 declined |= lib_fo && is_dhard(result+2); 1638 break; 1639 case t_modff: /* return float * float */ 1640 declined |= lib_fo && is_shard(result+2); 1641 /* fall through */ 1642 case t_frexpf: /* return float * int */ 1643 declined |= lib_fo && is_shard(result+0); 1644 break; 1645 case args1c: 1646 case args2c: 1647 declined |= lib_fo && is_dhard(result+0); 1648 declined |= lib_fo && is_dhard(result+4); 1649 break; 1650 case args1fc: 1651 case args2fc: 1652 declined |= lib_fo && is_shard(result+0); 1653 declined |= lib_fo && is_shard(result+4); 1654 break; 1655 } 1656 1657 /* Expect basic arithmetic tests to be declined if the command 1658 * line said that would happen */ 1659 declined |= (lib_no_arith && (fn->func == (funcptr)mpc_add || 1660 fn->func == (funcptr)mpc_sub || 1661 fn->func == (funcptr)mpc_mul || 1662 fn->func == (funcptr)mpc_div)); 1663 1664 if (!declined) { 1665 if (got_errno_in) 1666 ntests++; 1667 else 1668 ntests += 3; 1669 } 1670 } 1671 1672 void docase(Testable *fn, uint32 *args) { 1673 uint32 result[8]; /* real part in first 4, imaginary part in last 4 */ 1674 char *errstr = NULL; 1675 mpfr_t a, b, r; 1676 mpc_t ac, bc, rc; 1677 int rejected, printextra; 1678 wrapperctx ctx; 1679 1680 mpfr_init2(a, MPFR_PREC); 1681 mpfr_init2(b, MPFR_PREC); 1682 mpfr_init2(r, MPFR_PREC); 1683 mpc_init2(ac, MPFR_PREC); 1684 mpc_init2(bc, MPFR_PREC); 1685 mpc_init2(rc, MPFR_PREC); 1686 1687 printf("func=%s", fn->name); 1688 1689 rejected = 0; /* FIXME */ 1690 1691 switch (fn->type) { 1692 case args1: 1693 case rred: 1694 case semi1: 1695 case t_frexp: 1696 case t_modf: 1697 case classify: 1698 printf(" op1=%08x.%08x", args[0], args[1]); 1699 break; 1700 case args1f: 1701 case rredf: 1702 case semi1f: 1703 case t_frexpf: 1704 case t_modff: 1705 case classifyf: 1706 printf(" op1=%08x", args[0]); 1707 break; 1708 case args2: 1709 case semi2: 1710 case compare: 1711 printf(" op1=%08x.%08x", args[0], args[1]); 1712 printf(" op2=%08x.%08x", args[2], args[3]); 1713 break; 1714 case args2f: 1715 case semi2f: 1716 case t_ldexpf: 1717 case comparef: 1718 printf(" op1=%08x", args[0]); 1719 printf(" op2=%08x", args[2]); 1720 break; 1721 case t_ldexp: 1722 printf(" op1=%08x.%08x", args[0], args[1]); 1723 printf(" op2=%08x", args[2]); 1724 break; 1725 case args1c: 1726 case args1cr: 1727 printf(" op1r=%08x.%08x", args[0], args[1]); 1728 printf(" op1i=%08x.%08x", args[2], args[3]); 1729 break; 1730 case args2c: 1731 printf(" op1r=%08x.%08x", args[0], args[1]); 1732 printf(" op1i=%08x.%08x", args[2], args[3]); 1733 printf(" op2r=%08x.%08x", args[4], args[5]); 1734 printf(" op2i=%08x.%08x", args[6], args[7]); 1735 break; 1736 case args1fc: 1737 case args1fcr: 1738 printf(" op1r=%08x", args[0]); 1739 printf(" op1i=%08x", args[2]); 1740 break; 1741 case args2fc: 1742 printf(" op1r=%08x", args[0]); 1743 printf(" op1i=%08x", args[2]); 1744 printf(" op2r=%08x", args[4]); 1745 printf(" op2i=%08x", args[6]); 1746 break; 1747 default: 1748 fprintf(stderr, "internal inconsistency?!\n"); 1749 abort(); 1750 } 1751 1752 if (rejected == 2) { 1753 printf(" - test case rejected\n"); 1754 goto cleanup; 1755 } 1756 1757 wrapper_init(&ctx); 1758 1759 if (rejected == 0) { 1760 switch (fn->type) { 1761 case args1: 1762 set_mpfr_d(a, args[0], args[1]); 1763 wrapper_op_real(&ctx, a, 2, args); 1764 ((testfunc1)(fn->func))(r, a, GMP_RNDN); 1765 get_mpfr_d(r, &result[0], &result[1], &result[2]); 1766 wrapper_result_real(&ctx, r, 2, result); 1767 if (wrapper_run(&ctx, fn->wrappers)) 1768 get_mpfr_d(r, &result[0], &result[1], &result[2]); 1769 break; 1770 case args1cr: 1771 set_mpc_d(ac, args[0], args[1], args[2], args[3]); 1772 wrapper_op_complex(&ctx, ac, 2, args); 1773 ((testfunc1cr)(fn->func))(r, ac, GMP_RNDN); 1774 get_mpfr_d(r, &result[0], &result[1], &result[2]); 1775 wrapper_result_real(&ctx, r, 2, result); 1776 if (wrapper_run(&ctx, fn->wrappers)) 1777 get_mpfr_d(r, &result[0], &result[1], &result[2]); 1778 break; 1779 case args1f: 1780 set_mpfr_f(a, args[0]); 1781 wrapper_op_real(&ctx, a, 1, args); 1782 ((testfunc1)(fn->func))(r, a, GMP_RNDN); 1783 get_mpfr_f(r, &result[0], &result[1]); 1784 wrapper_result_real(&ctx, r, 1, result); 1785 if (wrapper_run(&ctx, fn->wrappers)) 1786 get_mpfr_f(r, &result[0], &result[1]); 1787 break; 1788 case args1fcr: 1789 set_mpc_f(ac, args[0], args[2]); 1790 wrapper_op_complex(&ctx, ac, 1, args); 1791 ((testfunc1cr)(fn->func))(r, ac, GMP_RNDN); 1792 get_mpfr_f(r, &result[0], &result[1]); 1793 wrapper_result_real(&ctx, r, 1, result); 1794 if (wrapper_run(&ctx, fn->wrappers)) 1795 get_mpfr_f(r, &result[0], &result[1]); 1796 break; 1797 case args2: 1798 set_mpfr_d(a, args[0], args[1]); 1799 wrapper_op_real(&ctx, a, 2, args); 1800 set_mpfr_d(b, args[2], args[3]); 1801 wrapper_op_real(&ctx, b, 2, args+2); 1802 ((testfunc2)(fn->func))(r, a, b, GMP_RNDN); 1803 get_mpfr_d(r, &result[0], &result[1], &result[2]); 1804 wrapper_result_real(&ctx, r, 2, result); 1805 if (wrapper_run(&ctx, fn->wrappers)) 1806 get_mpfr_d(r, &result[0], &result[1], &result[2]); 1807 break; 1808 case args2f: 1809 set_mpfr_f(a, args[0]); 1810 wrapper_op_real(&ctx, a, 1, args); 1811 set_mpfr_f(b, args[2]); 1812 wrapper_op_real(&ctx, b, 1, args+2); 1813 ((testfunc2)(fn->func))(r, a, b, GMP_RNDN); 1814 get_mpfr_f(r, &result[0], &result[1]); 1815 wrapper_result_real(&ctx, r, 1, result); 1816 if (wrapper_run(&ctx, fn->wrappers)) 1817 get_mpfr_f(r, &result[0], &result[1]); 1818 break; 1819 case rred: 1820 set_mpfr_d(a, args[0], args[1]); 1821 wrapper_op_real(&ctx, a, 2, args); 1822 ((testrred)(fn->func))(r, a, (int *)&result[3]); 1823 get_mpfr_d(r, &result[0], &result[1], &result[2]); 1824 wrapper_result_real(&ctx, r, 2, result); 1825 /* We never need to mess about with the integer auxiliary 1826 * output. */ 1827 if (wrapper_run(&ctx, fn->wrappers)) 1828 get_mpfr_d(r, &result[0], &result[1], &result[2]); 1829 break; 1830 case rredf: 1831 set_mpfr_f(a, args[0]); 1832 wrapper_op_real(&ctx, a, 1, args); 1833 ((testrred)(fn->func))(r, a, (int *)&result[3]); 1834 get_mpfr_f(r, &result[0], &result[1]); 1835 wrapper_result_real(&ctx, r, 1, result); 1836 /* We never need to mess about with the integer auxiliary 1837 * output. */ 1838 if (wrapper_run(&ctx, fn->wrappers)) 1839 get_mpfr_f(r, &result[0], &result[1]); 1840 break; 1841 case semi1: 1842 case semi1f: 1843 errstr = ((testsemi1)(fn->func))(args, result); 1844 break; 1845 case semi2: 1846 case compare: 1847 errstr = ((testsemi2)(fn->func))(args, args+2, result); 1848 break; 1849 case semi2f: 1850 case comparef: 1851 case t_ldexpf: 1852 errstr = ((testsemi2f)(fn->func))(args, args+2, result); 1853 break; 1854 case t_ldexp: 1855 errstr = ((testldexp)(fn->func))(args, args+2, result); 1856 break; 1857 case t_frexp: 1858 errstr = ((testfrexp)(fn->func))(args, result, result+2); 1859 break; 1860 case t_frexpf: 1861 errstr = ((testfrexp)(fn->func))(args, result, result+2); 1862 break; 1863 case t_modf: 1864 errstr = ((testmodf)(fn->func))(args, result, result+2); 1865 break; 1866 case t_modff: 1867 errstr = ((testmodf)(fn->func))(args, result, result+2); 1868 break; 1869 case classify: 1870 errstr = ((testclassify)(fn->func))(args, &result[0]); 1871 break; 1872 case classifyf: 1873 errstr = ((testclassifyf)(fn->func))(args, &result[0]); 1874 break; 1875 case args1c: 1876 set_mpc_d(ac, args[0], args[1], args[2], args[3]); 1877 wrapper_op_complex(&ctx, ac, 2, args); 1878 ((testfunc1c)(fn->func))(rc, ac, MPC_RNDNN); 1879 get_mpc_d(rc, &result[0], &result[1], &result[2], &result[4], &result[5], &result[6]); 1880 wrapper_result_complex(&ctx, rc, 2, result); 1881 if (wrapper_run(&ctx, fn->wrappers)) 1882 get_mpc_d(rc, &result[0], &result[1], &result[2], &result[4], &result[5], &result[6]); 1883 break; 1884 case args2c: 1885 set_mpc_d(ac, args[0], args[1], args[2], args[3]); 1886 wrapper_op_complex(&ctx, ac, 2, args); 1887 set_mpc_d(bc, args[4], args[5], args[6], args[7]); 1888 wrapper_op_complex(&ctx, bc, 2, args+4); 1889 ((testfunc2c)(fn->func))(rc, ac, bc, MPC_RNDNN); 1890 get_mpc_d(rc, &result[0], &result[1], &result[2], &result[4], &result[5], &result[6]); 1891 wrapper_result_complex(&ctx, rc, 2, result); 1892 if (wrapper_run(&ctx, fn->wrappers)) 1893 get_mpc_d(rc, &result[0], &result[1], &result[2], &result[4], &result[5], &result[6]); 1894 break; 1895 case args1fc: 1896 set_mpc_f(ac, args[0], args[2]); 1897 wrapper_op_complex(&ctx, ac, 1, args); 1898 ((testfunc1c)(fn->func))(rc, ac, MPC_RNDNN); 1899 get_mpc_f(rc, &result[0], &result[1], &result[4], &result[5]); 1900 wrapper_result_complex(&ctx, rc, 1, result); 1901 if (wrapper_run(&ctx, fn->wrappers)) 1902 get_mpc_f(rc, &result[0], &result[1], &result[4], &result[5]); 1903 break; 1904 case args2fc: 1905 set_mpc_f(ac, args[0], args[2]); 1906 wrapper_op_complex(&ctx, ac, 1, args); 1907 set_mpc_f(bc, args[4], args[6]); 1908 wrapper_op_complex(&ctx, bc, 1, args+4); 1909 ((testfunc2c)(fn->func))(rc, ac, bc, MPC_RNDNN); 1910 get_mpc_f(rc, &result[0], &result[1], &result[4], &result[5]); 1911 wrapper_result_complex(&ctx, rc, 1, result); 1912 if (wrapper_run(&ctx, fn->wrappers)) 1913 get_mpc_f(rc, &result[0], &result[1], &result[4], &result[5]); 1914 break; 1915 default: 1916 fprintf(stderr, "internal inconsistency?!\n"); 1917 abort(); 1918 } 1919 } 1920 1921 switch (fn->type) { 1922 case args1: /* return an extra-precise result */ 1923 case args2: 1924 case args1cr: 1925 case rred: 1926 printextra = 1; 1927 if (rejected == 0) { 1928 errstr = NULL; 1929 if (!mpfr_zero_p(a)) { 1930 if ((result[0] & 0x7FFFFFFF) == 0 && result[1] == 0) { 1931 /* 1932 * If the output is +0 or -0 apart from the extra 1933 * precision in result[2], then there's a tricky 1934 * judgment call about what we require in the 1935 * output. If we output the extra bits and set 1936 * errstr="?underflow" then mathtest will tolerate 1937 * the function under test rounding down to zero 1938 * _or_ up to the minimum denormal; whereas if we 1939 * suppress the extra bits and set 1940 * errstr="underflow", then mathtest will enforce 1941 * that the function really does underflow to zero. 1942 * 1943 * But where to draw the line? It seems clear to 1944 * me that numbers along the lines of 1945 * 00000000.00000000.7ff should be treated 1946 * similarly to 00000000.00000000.801, but on the 1947 * other hand, we must surely be prepared to 1948 * enforce a genuine underflow-to-zero in _some_ 1949 * case where the true mathematical output is 1950 * nonzero but absurdly tiny. 1951 * 1952 * I think a reasonable place to draw the 1953 * distinction is at 00000000.00000000.400, i.e. 1954 * one quarter of the minimum positive denormal. 1955 * If a value less than that rounds up to the 1956 * minimum denormal, that must mean the function 1957 * under test has managed to make an error of an 1958 * entire factor of two, and that's something we 1959 * should fix. Above that, you can misround within 1960 * the limits of your accuracy bound if you have 1961 * to. 1962 */ 1963 if (result[2] < 0x40000000) { 1964 /* Total underflow (ERANGE + UFL) is required, 1965 * and we suppress the extra bits to make 1966 * mathtest enforce that the output is really 1967 * zero. */ 1968 errstr = "underflow"; 1969 printextra = 0; 1970 } else { 1971 /* Total underflow is not required, but if the 1972 * function rounds down to zero anyway, then 1973 * we should be prepared to tolerate it. */ 1974 errstr = "?underflow"; 1975 } 1976 } else if (!(result[0] & 0x7ff00000)) { 1977 /* 1978 * If the output is denormal, we usually expect a 1979 * UFL exception, warning the user of partial 1980 * underflow. The exception is if the denormal 1981 * being returned is just one of the input values, 1982 * unchanged even in principle. I bodgily handle 1983 * this by just special-casing the functions in 1984 * question below. 1985 */ 1986 if (!strcmp(fn->name, "fmax") || 1987 !strcmp(fn->name, "fmin") || 1988 !strcmp(fn->name, "creal") || 1989 !strcmp(fn->name, "cimag")) { 1990 /* no error expected */ 1991 } else { 1992 errstr = "u"; 1993 } 1994 } else if ((result[0] & 0x7FFFFFFF) > 0x7FEFFFFF) { 1995 /* 1996 * Infinite results are usually due to overflow, 1997 * but one exception is lgamma of a negative 1998 * integer. 1999 */ 2000 if (!strcmp(fn->name, "lgamma") && 2001 (args[0] & 0x80000000) != 0 && /* negative */ 2002 is_dinteger(args)) { 2003 errstr = "ERANGE status=z"; 2004 } else { 2005 errstr = "overflow"; 2006 } 2007 printextra = 0; 2008 } 2009 } else { 2010 /* lgamma(0) is also a pole. */ 2011 if (!strcmp(fn->name, "lgamma")) { 2012 errstr = "ERANGE status=z"; 2013 printextra = 0; 2014 } 2015 } 2016 } 2017 2018 if (!printextra || (rejected && !(rejected==1 && result[2]!=0))) { 2019 printf(" result=%08x.%08x", 2020 result[0], result[1]); 2021 } else { 2022 printf(" result=%08x.%08x.%03x", 2023 result[0], result[1], (result[2] >> 20) & 0xFFF); 2024 } 2025 if (fn->type == rred) { 2026 printf(" res2=%08x", result[3]); 2027 } 2028 break; 2029 case args1f: 2030 case args2f: 2031 case args1fcr: 2032 case rredf: 2033 printextra = 1; 2034 if (rejected == 0) { 2035 errstr = NULL; 2036 if (!mpfr_zero_p(a)) { 2037 if ((result[0] & 0x7FFFFFFF) == 0) { 2038 /* 2039 * Decide whether to print the extra bits based on 2040 * just how close to zero the number is. See the 2041 * big comment in the double-precision case for 2042 * discussion. 2043 */ 2044 if (result[1] < 0x40000000) { 2045 errstr = "underflow"; 2046 printextra = 0; 2047 } else { 2048 errstr = "?underflow"; 2049 } 2050 } else if (!(result[0] & 0x7f800000)) { 2051 /* 2052 * Functions which do not report partial overflow 2053 * are listed here as special cases. (See the 2054 * corresponding double case above for a fuller 2055 * comment.) 2056 */ 2057 if (!strcmp(fn->name, "fmaxf") || 2058 !strcmp(fn->name, "fminf") || 2059 !strcmp(fn->name, "crealf") || 2060 !strcmp(fn->name, "cimagf")) { 2061 /* no error expected */ 2062 } else { 2063 errstr = "u"; 2064 } 2065 } else if ((result[0] & 0x7FFFFFFF) > 0x7F7FFFFF) { 2066 /* 2067 * Infinite results are usually due to overflow, 2068 * but one exception is lgamma of a negative 2069 * integer. 2070 */ 2071 if (!strcmp(fn->name, "lgammaf") && 2072 (args[0] & 0x80000000) != 0 && /* negative */ 2073 is_sinteger(args)) { 2074 errstr = "ERANGE status=z"; 2075 } else { 2076 errstr = "overflow"; 2077 } 2078 printextra = 0; 2079 } 2080 } else { 2081 /* lgamma(0) is also a pole. */ 2082 if (!strcmp(fn->name, "lgammaf")) { 2083 errstr = "ERANGE status=z"; 2084 printextra = 0; 2085 } 2086 } 2087 } 2088 2089 if (!printextra || (rejected && !(rejected==1 && result[1]!=0))) { 2090 printf(" result=%08x", 2091 result[0]); 2092 } else { 2093 printf(" result=%08x.%03x", 2094 result[0], (result[1] >> 20) & 0xFFF); 2095 } 2096 if (fn->type == rredf) { 2097 printf(" res2=%08x", result[3]); 2098 } 2099 break; 2100 case semi1: /* return a double result */ 2101 case semi2: 2102 case t_ldexp: 2103 printf(" result=%08x.%08x", result[0], result[1]); 2104 break; 2105 case semi1f: 2106 case semi2f: 2107 case t_ldexpf: 2108 printf(" result=%08x", result[0]); 2109 break; 2110 case t_frexp: /* return double * int */ 2111 printf(" result=%08x.%08x res2=%08x", result[0], result[1], 2112 result[2]); 2113 break; 2114 case t_modf: /* return double * double */ 2115 printf(" result=%08x.%08x res2=%08x.%08x", 2116 result[0], result[1], result[2], result[3]); 2117 break; 2118 case t_modff: /* return float * float */ 2119 /* fall through */ 2120 case t_frexpf: /* return float * int */ 2121 printf(" result=%08x res2=%08x", result[0], result[2]); 2122 break; 2123 case classify: 2124 case classifyf: 2125 case compare: 2126 case comparef: 2127 printf(" result=%x", result[0]); 2128 break; 2129 case args1c: 2130 case args2c: 2131 if (0/* errstr */) { 2132 printf(" resultr=%08x.%08x", result[0], result[1]); 2133 printf(" resulti=%08x.%08x", result[4], result[5]); 2134 } else { 2135 printf(" resultr=%08x.%08x.%03x", 2136 result[0], result[1], (result[2] >> 20) & 0xFFF); 2137 printf(" resulti=%08x.%08x.%03x", 2138 result[4], result[5], (result[6] >> 20) & 0xFFF); 2139 } 2140 /* Underflow behaviour doesn't seem to be specified for complex arithmetic */ 2141 errstr = "?underflow"; 2142 break; 2143 case args1fc: 2144 case args2fc: 2145 if (0/* errstr */) { 2146 printf(" resultr=%08x", result[0]); 2147 printf(" resulti=%08x", result[4]); 2148 } else { 2149 printf(" resultr=%08x.%03x", 2150 result[0], (result[1] >> 20) & 0xFFF); 2151 printf(" resulti=%08x.%03x", 2152 result[4], (result[5] >> 20) & 0xFFF); 2153 } 2154 /* Underflow behaviour doesn't seem to be specified for complex arithmetic */ 2155 errstr = "?underflow"; 2156 break; 2157 } 2158 2159 if (errstr && *(errstr+1) == '\0') { 2160 printf(" errno=0 status=%c",*errstr); 2161 } else if (errstr && *errstr == '?') { 2162 printf(" maybeerror=%s", errstr+1); 2163 } else if (errstr && errstr[0] == 'E') { 2164 printf(" errno=%s", errstr); 2165 } else { 2166 printf(" error=%s", errstr && *errstr ? errstr : "0"); 2167 } 2168 2169 printf("\n"); 2170 2171 vet_for_decline(fn, args, result, 0); 2172 2173 cleanup: 2174 mpfr_clear(a); 2175 mpfr_clear(b); 2176 mpfr_clear(r); 2177 mpc_clear(ac); 2178 mpc_clear(bc); 2179 mpc_clear(rc); 2180 } 2181 2182 void gencases(Testable *fn, int number) { 2183 int i; 2184 uint32 args[8]; 2185 2186 float32_case(NULL); 2187 float64_case(NULL); 2188 2189 printf("random=on\n"); /* signal to runtests.pl that the following tests are randomly generated */ 2190 for (i = 0; i < number; i++) { 2191 /* generate test point */ 2192 fn->cases(args, fn->caseparam1, fn->caseparam2); 2193 docase(fn, args); 2194 } 2195 printf("random=off\n"); 2196 } 2197 2198 static uint32 doubletop(int x, int scale) { 2199 int e = 0x412 + scale; 2200 while (!(x & 0x100000)) 2201 x <<= 1, e--; 2202 return (e << 20) + x; 2203 } 2204 2205 static uint32 floatval(int x, int scale) { 2206 int e = 0x95 + scale; 2207 while (!(x & 0x800000)) 2208 x <<= 1, e--; 2209 return (e << 23) + x; 2210 } 2211