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