1 /*
2 * dotest.c - actually generate mathlib test cases
3 *
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-exception
7 */
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
22 extern int lib_fo, lib_no_arith, ntests;
23
24 /*
25 * Prototypes.
26 */
27 static void cases_biased(uint32 *, uint32, uint32);
28 static void cases_biased_positive(uint32 *, uint32, uint32);
29 static void cases_biased_float(uint32 *, uint32, uint32);
30 static void cases_uniform(uint32 *, uint32, uint32);
31 static void cases_uniform_positive(uint32 *, uint32, uint32);
32 static void cases_uniform_float(uint32 *, uint32, uint32);
33 static void cases_uniform_float_positive(uint32 *, uint32, uint32);
34 static void log_cases(uint32 *, uint32, uint32);
35 static void log_cases_float(uint32 *, uint32, uint32);
36 static void log1p_cases(uint32 *, uint32, uint32);
37 static void log1p_cases_float(uint32 *, uint32, uint32);
38 static void minmax_cases(uint32 *, uint32, uint32);
39 static void minmax_cases_float(uint32 *, uint32, uint32);
40 static void atan2_cases(uint32 *, uint32, uint32);
41 static void atan2_cases_float(uint32 *, uint32, uint32);
42 static void pow_cases(uint32 *, uint32, uint32);
43 static void pow_cases_float(uint32 *, uint32, uint32);
44 static void rred_cases(uint32 *, uint32, uint32);
45 static void rred_cases_float(uint32 *, uint32, uint32);
46 static void cases_semi1(uint32 *, uint32, uint32);
47 static void cases_semi1_float(uint32 *, uint32, uint32);
48 static void cases_semi2(uint32 *, uint32, uint32);
49 static void cases_semi2_float(uint32 *, uint32, uint32);
50 static void cases_ldexp(uint32 *, uint32, uint32);
51 static void cases_ldexp_float(uint32 *, uint32, uint32);
52
53 static void complex_cases_uniform(uint32 *, uint32, uint32);
54 static void complex_cases_uniform_float(uint32 *, uint32, uint32);
55 static void complex_cases_biased(uint32 *, uint32, uint32);
56 static void complex_cases_biased_float(uint32 *, uint32, uint32);
57 static void complex_log_cases(uint32 *, uint32, uint32);
58 static void complex_log_cases_float(uint32 *, uint32, uint32);
59 static void complex_pow_cases(uint32 *, uint32, uint32);
60 static void complex_pow_cases_float(uint32 *, uint32, uint32);
61 static void complex_arithmetic_cases(uint32 *, uint32, uint32);
62 static void complex_arithmetic_cases_float(uint32 *, uint32, uint32);
63
64 static uint32 doubletop(int x, int scale);
65 static uint32 floatval(int x, int scale);
66
67 /*
68 * Convert back and forth between IEEE bit patterns and the
69 * mpfr_t/mpc_t types.
70 */
set_mpfr_d(mpfr_t x,uint32 h,uint32 l)71 static 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 else
81 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 else
89 exp++;
90 mpfr_set_sj_2exp(x, mantissa * sign, (int)exp - 0x3ff - 52, GMP_RNDN);
91 }
92 }
set_mpfr_f(mpfr_t x,uint32 f)93 static 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 else
102 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 else
110 exp++;
111 mpfr_set_sj_2exp(x, mantissa * sign, (int)exp - 0x7f - 23, GMP_RNDN);
112 }
113 }
set_mpc_d(mpc_t z,uint32 rh,uint32 rl,uint32 ih,uint32 il)114 static 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 }
set_mpc_f(mpc_t z,uint32 r,uint32 i)125 static 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 }
get_mpfr_d(const mpfr_t x,uint32 * h,uint32 * l,uint32 * extra)136 static 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 }
get_mpfr_f(const mpfr_t x,uint32 * f,uint32 * extra)200 static 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 }
get_mpc_d(const mpc_t z,uint32 * rh,uint32 * rl,uint32 * rextra,uint32 * ih,uint32 * il,uint32 * iextra)256 static 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 }
get_mpc_f(const mpc_t z,uint32 * r,uint32 * rextra,uint32 * i,uint32 * iextra)270 static 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 trivially
287 * implementable using an existing mpfr or mpc function.
288 */
test_rred(mpfr_t ret,const mpfr_t x,int * quadrant)289 int 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 size
297 * around 2^1024, and must find its remainder mod pi, which means
298 * we need 1024 bits of pi at least. Plus, the remainder might
299 * happen to come out very very small if we're unlucky. How
300 * unlucky can we be? Well, conveniently, I once went through and
301 * actually worked that out using Paxson's modular minimisation
302 * algorithm, and it turns out that the smallest exponent you can
303 * get out of a nontrivial[1] double precision range reduction is
304 * 0x3c2, i.e. of the order of 2^-61. So we need 1024 bits of pi
305 * to get us down to the units digit, another 61 or so bits (say
306 * 64) to get down to the highest set bit of the output, and then
307 * some bits to make the actual mantissa big enough.
308 *
309 * [1] of course the output of range reduction can have an
310 * arbitrarily small exponent in the trivial case, where the
311 * input is so small that it's the identity function. That
312 * 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 }
test_lgamma(mpfr_t ret,const mpfr_t x,mpfr_rnd_t rnd)325 int 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 output
329 * sign. We don't bother testing that, so this wrapper throws away
330 * the sign and hence fits into the same function prototype as all
331 * the other real->real mpfr functions.
332 *
333 * There is also mpfr_lngamma which has no sign output and hence
334 * has the right prototype already, but unfortunately it returns
335 * 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 }
test_cpow(mpc_t ret,const mpc_t x,const mpc_t y,mpc_rnd_t rnd)340 int 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 amount
344 * if we want it to get the really difficult cases right. (Not
345 * that we expect the library under test to be getting those cases
346 * right itself, but we'd at least like the test suite to report
347 * them as wrong for the _right reason_.)
348 *
349 * This works around a bug in mpc_pow(), fixed by r1455 in the MPC
350 * svn repository (2014-10-14) and expected to be in any MPC
351 * release after 1.0.2 (which was the latest release already made
352 * at the time of the fix). So as and when we update to an MPC
353 * 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 the
357 * 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 deciding
380 * 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 */
is_dinteger(uint32 * in)393 int 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 }
is_sinteger(uint32 * in)401 int 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 */
is_dsnan(const uint32 * in)413 int 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 }
is_ssnan(const uint32 * in)423 int 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 }
is_snan(const uint32 * in,int size)433 int 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 main
440 * test function has run.
441 */
universal_wrapper(wrapperctx * ctx)442 void 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
458 Testable functions[] = {
459 /*
460 * Trig functions: sin, cos, tan. We test the core function
461 * between -16 and +16: we assume that range reduction exists
462 * and will be used for larger arguments, and we'll test that
463 * 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 as
465 * 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. acos
485 * 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) with
497 * 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 the
506 * 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 around
514 * 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 cover
537 * the whole range of positive nonzero finite numbers. It never
538 * 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 all
551 * finite values, but will only bother for things between 2^-3
552 * 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-numericals
583 */
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 * Comparisons
600 */
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 functions
617 */
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 bottom
628 * because historically they were not tested because we didn't
629 * 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
713 const int nfunctions = ( sizeof(functions)/sizeof(*functions) );
714
715 #define random_sign ( random_upto(1) ? 0x80000000 : 0 )
716
iszero(uint32 * x)717 static int iszero(uint32 *x) {
718 return !((x[0] & 0x7FFFFFFF) || x[1]);
719 }
720
721
complex_log_cases(uint32 * out,uint32 param1,uint32 param2)722 static 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
complex_log_cases_float(uint32 * out,uint32 param1,uint32 param2)729 static 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
complex_cases_biased(uint32 * out,uint32 lowbound,uint32 highbound)735 static 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
complex_cases_biased_float(uint32 * out,uint32 lowbound,uint32 highbound)741 static 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
complex_cases_uniform(uint32 * out,uint32 lowbound,uint32 highbound)747 static 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
complex_cases_uniform_float(uint32 * out,uint32 lowbound,uint32 highbound)753 static 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
complex_pow_cases(uint32 * out,uint32 lowbound,uint32 highbound)759 static 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 hence
765 * its magnitude is within [1/2,2*sqrt(2)]. The magnitude of its
766 * logarithm in base 2 is therefore at most the magnitude of
767 * (log2(2*sqrt(2)) + i*pi/log(2)), or in other words
768 * hypot(3/2,pi/log(2)) = 4.77. So the magnitude of the exponent
769 * input must be at most our output magnitude limit (as a power
770 * of two) divided by that.
771 *
772 * I also set the output magnitude limit a bit low, because we
773 * don't guarantee (and neither does glibc) to prevent internal
774 * overflow in cases where the output _magnitude_ overflows but
775 * scaling it back down by cos and sin of the argument brings it
776 * 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
complex_pow_cases_float(uint32 * out,uint32 lowbound,uint32 highbound)784 static void complex_pow_cases_float(uint32 *out, uint32 lowbound,
785 uint32 highbound) {
786 /*
787 * Reasoning as above, though of course the detailed numbers are
788 * 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
complex_arithmetic_cases(uint32 * out,uint32 lowbound,uint32 highbound)796 static 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
complex_arithmetic_cases_float(uint32 * out,uint32 lowbound,uint32 highbound)804 static 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-contained
814 * form.
815 */
816
float32_case(uint32 * ret)817 void 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, 128
835 };
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 * - FIXME
913 * Significand:
914 * - with prob 1/5, a totally random bit pattern
915 * - with prob 1/5, all 1s down to some point and then random
916 * - with prob 1/5, all 1s up to some point and then random
917 * - with prob 1/5, all 0s down to some point and then random
918 * - with prob 1/5, all 0s up to some point and then random
919 */
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 }
float64_case(uint32 * ret)944 static 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, 1024
963 };
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 * - FIXME
1081 * Significand:
1082 * - with prob 1/5, a totally random bit pattern
1083 * - with prob 1/5, all 1s down to some point and then random
1084 * - with prob 1/5, all 1s up to some point and then random
1085 * - with prob 1/5, all 0s down to some point and then random
1086 * - with prob 1/5, all 0s up to some point and then random
1087 */
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
cases_biased(uint32 * out,uint32 lowbound,uint32 highbound)1139 static 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
cases_biased_positive(uint32 * out,uint32 lowbound,uint32 highbound)1148 static 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
cases_biased_float(uint32 * out,uint32 lowbound,uint32 highbound)1156 static 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
cases_semi1(uint32 * out,uint32 param1,uint32 param2)1165 static void cases_semi1(uint32 *out, uint32 param1,
1166 uint32 param2) {
1167 float64_case(out);
1168 }
1169
cases_semi1_float(uint32 * out,uint32 param1,uint32 param2)1170 static void cases_semi1_float(uint32 *out, uint32 param1,
1171 uint32 param2) {
1172 float32_case(out);
1173 }
1174
cases_semi2(uint32 * out,uint32 param1,uint32 param2)1175 static void cases_semi2(uint32 *out, uint32 param1,
1176 uint32 param2) {
1177 float64_case(out);
1178 float64_case(out+2);
1179 }
1180
cases_semi2_float(uint32 * out,uint32 param1,uint32 param2)1181 static void cases_semi2_float(uint32 *out, uint32 param1,
1182 uint32 param2) {
1183 float32_case(out);
1184 float32_case(out+2);
1185 }
1186
cases_ldexp(uint32 * out,uint32 param1,uint32 param2)1187 static void cases_ldexp(uint32 *out, uint32 param1,
1188 uint32 param2) {
1189 float64_case(out);
1190 out[2] = random_upto(2048)-1024;
1191 }
1192
cases_ldexp_float(uint32 * out,uint32 param1,uint32 param2)1193 static void cases_ldexp_float(uint32 *out, uint32 param1,
1194 uint32 param2) {
1195 float32_case(out);
1196 out[2] = random_upto(256)-128;
1197 }
1198
cases_uniform(uint32 * out,uint32 lowbound,uint32 highbound)1199 static 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 }
cases_uniform_float(uint32 * out,uint32 lowbound,uint32 highbound)1207 static 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
cases_uniform_positive(uint32 * out,uint32 lowbound,uint32 highbound)1216 static 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 }
cases_uniform_float_positive(uint32 * out,uint32 lowbound,uint32 highbound)1223 static 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
log_cases(uint32 * out,uint32 param1,uint32 param2)1232 static 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
log_cases_float(uint32 * out,uint32 param1,uint32 param2)1240 static 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
log1p_cases(uint32 * out,uint32 param1,uint32 param2)1248 static 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
log1p_cases_float(uint32 * out,uint32 param1,uint32 param2)1259 static 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
minmax_cases(uint32 * out,uint32 param1,uint32 param2)1270 static 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
minmax_cases_float(uint32 * out,uint32 param1,uint32 param2)1282 static 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
rred_cases(uint32 * out,uint32 param1,uint32 param2)1294 static 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
rred_cases_float(uint32 * out,uint32 param1,uint32 param2)1303 static 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
atan2_cases(uint32 * out,uint32 param1,uint32 param2)1312 static 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 } else
1321 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
atan2_cases_float(uint32 * out,uint32 param1,uint32 param2)1331 static 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 } else
1340 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
pow_cases(uint32 * out,uint32 param1,uint32 param2)1349 static void pow_cases(uint32 *out, uint32 param1,
1350 uint32 param2) {
1351 /*
1352 * Pick an exponent e (-0x33 to +0x7FE) for x, and here's the
1353 * 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 one
1359 * end or the other, so we have to be cleverer: pick a number n
1360 * of useful bits in the mantissa (1 thru 52, so 1 must imply
1361 * 0x3ff00000.00000001 whereas 52 is anything at least as big
1362 * as 0x3ff80000.00000000; for e == 0x3fe, 1 necessarily means
1363 * 0x3fefffff.ffffffff and 52 is anything at most as big as
1364 * 0x3fe80000.00000000). Then, as it happens, a sensible
1365 * maximum power is 2^(63-n) for e == 0x3fe, and 2^(62-n) for
1366 * e == 0x3ff.
1367 *
1368 * We inevitably get some overflows in approximating the log
1369 * 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 a
1373 * 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 else
1394 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 else
1406 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 }
pow_cases_float(uint32 * out,uint32 param1,uint32 param2)1440 static 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 the
1444 * 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 one
1450 * end or the other, so we have to be cleverer: pick a number n
1451 * of useful bits in the mantissa (1 thru 23, so 1 must imply
1452 * 0x3f800001 whereas 23 is anything at least as big as
1453 * 0x3fc00000; for e == 0x7e, 1 necessarily means 0x3f7fffff
1454 * and 23 is anything at most as big as 0x3f400000). Then, as
1455 * 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 log
1459 * 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 a
1463 * 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 else
1484 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 else
1496 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
vet_for_decline(Testable * fn,uint32 * args,uint32 * result,int got_errno_in)1521 void 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 command
1616 * 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 else
1626 ntests += 3;
1627 }
1628 }
1629
docase(Testable * fn,uint32 * args)1630 void 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 auxiliary
1784 * 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 auxiliary
1795 * 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 extra
1891 * precision in result[2], then there's a tricky
1892 * judgment call about what we require in the
1893 * output. If we output the extra bits and set
1894 * errstr="?underflow" then mathtest will tolerate
1895 * the function under test rounding down to zero
1896 * _or_ up to the minimum denormal; whereas if we
1897 * suppress the extra bits and set
1898 * errstr="underflow", then mathtest will enforce
1899 * that the function really does underflow to zero.
1900 *
1901 * But where to draw the line? It seems clear to
1902 * me that numbers along the lines of
1903 * 00000000.00000000.7ff should be treated
1904 * similarly to 00000000.00000000.801, but on the
1905 * other hand, we must surely be prepared to
1906 * enforce a genuine underflow-to-zero in _some_
1907 * case where the true mathematical output is
1908 * nonzero but absurdly tiny.
1909 *
1910 * I think a reasonable place to draw the
1911 * 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 the
1914 * minimum denormal, that must mean the function
1915 * under test has managed to make an error of an
1916 * entire factor of two, and that's something we
1917 * should fix. Above that, you can misround within
1918 * the limits of your accuracy bound if you have
1919 * to.
1920 */
1921 if (result[2] < 0x40000000) {
1922 /* Total underflow (ERANGE + UFL) is required,
1923 * and we suppress the extra bits to make
1924 * mathtest enforce that the output is really
1925 * zero. */
1926 errstr = "underflow";
1927 printextra = 0;
1928 } else {
1929 /* Total underflow is not required, but if the
1930 * function rounds down to zero anyway, then
1931 * 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 a
1937 * UFL exception, warning the user of partial
1938 * underflow. The exception is if the denormal
1939 * being returned is just one of the input values,
1940 * unchanged even in principle. I bodgily handle
1941 * this by just special-casing the functions in
1942 * 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 negative
1956 * 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 on
1998 * just how close to zero the number is. See the
1999 * big comment in the double-precision case for
2000 * 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 overflow
2011 * are listed here as special cases. (See the
2012 * corresponding double case above for a fuller
2013 * 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 negative
2027 * 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
gencases(Testable * fn,int number)2140 void 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
doubletop(int x,int scale)2156 static 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
floatval(int x,int scale)2163 static 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