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