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