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