1 /* Test mpn_get_d.
2
3 Copyright 2002-2004 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library test suite.
6
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15 Public License for more details.
16
17 You should have received a copy of the GNU General Public License along with
18 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
19
20 #include "config.h"
21
22 #include <setjmp.h>
23 #include <signal.h>
24 #include <stdio.h>
25 #include <stdlib.h>
26
27 #include "gmp-impl.h"
28 #include "tests.h"
29
30
31 #ifndef _GMP_IEEE_FLOATS
32 #define _GMP_IEEE_FLOATS 0
33 #endif
34
35
36 /* Exercise various 2^n values, with various exponents and positive and
37 negative. */
38 void
39 check_onebit (void)
40 {
41 static const int bit_table[] = {
42 0, 1, 2, 3,
43 GMP_NUMB_BITS - 2, GMP_NUMB_BITS - 1,
44 GMP_NUMB_BITS,
45 GMP_NUMB_BITS + 1, GMP_NUMB_BITS + 2,
46 2 * GMP_NUMB_BITS - 2, 2 * GMP_NUMB_BITS - 1,
47 2 * GMP_NUMB_BITS,
48 2 * GMP_NUMB_BITS + 1, 2 * GMP_NUMB_BITS + 2,
49 3 * GMP_NUMB_BITS - 2, 3 * GMP_NUMB_BITS - 1,
50 3 * GMP_NUMB_BITS,
51 3 * GMP_NUMB_BITS + 1, 3 * GMP_NUMB_BITS + 2,
52 4 * GMP_NUMB_BITS - 2, 4 * GMP_NUMB_BITS - 1,
53 4 * GMP_NUMB_BITS,
54 4 * GMP_NUMB_BITS + 1, 4 * GMP_NUMB_BITS + 2,
55 5 * GMP_NUMB_BITS - 2, 5 * GMP_NUMB_BITS - 1,
56 5 * GMP_NUMB_BITS,
57 5 * GMP_NUMB_BITS + 1, 5 * GMP_NUMB_BITS + 2,
58 6 * GMP_NUMB_BITS - 2, 6 * GMP_NUMB_BITS - 1,
59 6 * GMP_NUMB_BITS,
60 6 * GMP_NUMB_BITS + 1, 6 * GMP_NUMB_BITS + 2,
61 };
62 static const int exp_table[] = {
63 0, -100, -10, -1, 1, 10, 100,
64 };
65
66 /* FIXME: It'd be better to base this on the float format. */
67 #if defined (__vax) || defined (__vax__)
68 int limit = 127; /* vax fp numbers have limited range */
69 #else
70 int limit = 511;
71 #endif
72
73 int bit_i, exp_i, i;
74 double got, want;
75 mp_size_t nsize, sign;
76 long bit, exp, want_bit;
77 mp_limb_t np[20];
78
79 for (bit_i = 0; bit_i < numberof (bit_table); bit_i++)
80 {
81 bit = bit_table[bit_i];
82
83 nsize = BITS_TO_LIMBS (bit+1);
84 refmpn_zero (np, nsize);
85 np[bit/GMP_NUMB_BITS] = CNST_LIMB(1) << (bit % GMP_NUMB_BITS);
86
87 for (exp_i = 0; exp_i < numberof (exp_table); exp_i++)
88 {
89 exp = exp_table[exp_i];
90
91 want_bit = bit + exp;
92 if (want_bit >= limit || want_bit <= -limit)
93 continue;
94
95 want = 1.0;
96 for (i = 0; i < want_bit; i++)
97 want *= 2.0;
98 for (i = 0; i > want_bit; i--)
99 want *= 0.5;
100
101 for (sign = 0; sign >= -1; sign--, want = -want)
102 {
103 got = mpn_get_d (np, nsize, sign, exp);
104 if (got != want)
105 {
106 printf ("mpn_get_d wrong on 2^n\n");
107 printf (" bit %ld\n", bit);
108 printf (" exp %ld\n", exp);
109 printf (" want_bit %ld\n", want_bit);
110 printf (" sign %ld\n", (long) sign);
111 mpn_trace (" n ", np, nsize);
112 printf (" nsize %ld\n", (long) nsize);
113 d_trace (" want ", want);
114 d_trace (" got ", got);
115 abort();
116 }
117 }
118 }
119 }
120 }
121
122
123 /* Exercise values 2^n+1, while such a value fits the mantissa of a double. */
124 void
125 check_twobit (void)
126 {
127 int i, mant_bits;
128 double got, want;
129 mp_size_t nsize, sign;
130 mp_ptr np;
131
132 mant_bits = tests_dbl_mant_bits ();
133 if (mant_bits == 0)
134 return;
135
136 np = refmpn_malloc_limbs (BITS_TO_LIMBS (mant_bits));
137 want = 3.0;
138 for (i = 1; i < mant_bits; i++)
139 {
140 nsize = BITS_TO_LIMBS (i+1);
141 refmpn_zero (np, nsize);
142 np[i/GMP_NUMB_BITS] = CNST_LIMB(1) << (i % GMP_NUMB_BITS);
143 np[0] |= 1;
144
145 for (sign = 0; sign >= -1; sign--)
146 {
147 got = mpn_get_d (np, nsize, sign, 0);
148 if (got != want)
149 {
150 printf ("mpn_get_d wrong on 2^%d + 1\n", i);
151 printf (" sign %ld\n", (long) sign);
152 mpn_trace (" n ", np, nsize);
153 printf (" nsize %ld\n", (long) nsize);
154 d_trace (" want ", want);
155 d_trace (" got ", got);
156 abort();
157 }
158 want = -want;
159 }
160
161 want = 2.0 * want - 1.0;
162 }
163
164 free (np);
165 }
166
167
168 /* Expect large negative exponents to underflow to 0.0.
169 Some systems might have hardware traps for such an underflow (though
170 usually it's not the default), so watch out for SIGFPE. */
171 void
172 check_underflow (void)
173 {
174 static const long exp_table[] = {
175 -999999L, LONG_MIN,
176 };
177 static const mp_limb_t np[1] = { 1 };
178
179 static long exp;
180 mp_size_t nsize, sign;
181 double got;
182 int exp_i;
183
184 nsize = numberof (np);
185
186 if (tests_setjmp_sigfpe() == 0)
187 {
188 for (exp_i = 0; exp_i < numberof (exp_table); exp_i++)
189 {
190 exp = exp_table[exp_i];
191
192 for (sign = 0; sign >= -1; sign--)
193 {
194 got = mpn_get_d (np, nsize, sign, exp);
195 if (got != 0.0)
196 {
197 printf ("mpn_get_d wrong, didn't get 0.0 on underflow\n");
198 printf (" nsize %ld\n", (long) nsize);
199 printf (" exp %ld\n", exp);
200 printf (" sign %ld\n", (long) sign);
201 d_trace (" got ", got);
202 abort ();
203 }
204 }
205 }
206 }
207 else
208 {
209 printf ("Warning, underflow to zero tests skipped due to SIGFPE (exp=%ld)\n", exp);
210 }
211 tests_sigfpe_done ();
212 }
213
214
215 /* Expect large values to result in +/-inf, on IEEE systems. */
216 void
217 check_inf (void)
218 {
219 static const long exp_table[] = {
220 999999L, LONG_MAX,
221 };
222 static const mp_limb_t np[4] = { 1, 1, 1, 1 };
223 long exp;
224 mp_size_t nsize, sign, got_sign;
225 double got;
226 int exp_i;
227
228 if (! _GMP_IEEE_FLOATS)
229 return;
230
231 for (nsize = 1; nsize <= numberof (np); nsize++)
232 {
233 for (exp_i = 0; exp_i < numberof (exp_table); exp_i++)
234 {
235 exp = exp_table[exp_i];
236
237 for (sign = 0; sign >= -1; sign--)
238 {
239 got = mpn_get_d (np, nsize, sign, exp);
240 got_sign = (got >= 0 ? 0 : -1);
241 if (! tests_isinf (got))
242 {
243 printf ("mpn_get_d wrong, didn't get infinity\n");
244 bad:
245 printf (" nsize %ld\n", (long) nsize);
246 printf (" exp %ld\n", exp);
247 printf (" sign %ld\n", (long) sign);
248 d_trace (" got ", got);
249 printf (" got sign %ld\n", (long) got_sign);
250 abort ();
251 }
252 if (got_sign != sign)
253 {
254 printf ("mpn_get_d wrong sign on infinity\n");
255 goto bad;
256 }
257 }
258 }
259 }
260 }
261
262 /* Check values 2^n approaching and into IEEE denorm range.
263 Some systems might not support denorms, or might have traps setup, so
264 watch out for SIGFPE. */
265 void
266 check_ieee_denorm (void)
267 {
268 static long exp;
269 mp_limb_t n = 1;
270 long i;
271 mp_size_t sign;
272 double want, got;
273
274 if (! _GMP_IEEE_FLOATS)
275 return;
276
277 if (tests_setjmp_sigfpe() == 0)
278 {
279 exp = -1020;
280 want = 1.0;
281 for (i = 0; i > exp; i--)
282 want *= 0.5;
283
284 for ( ; exp > -1500 && want != 0.0; exp--)
285 {
286 for (sign = 0; sign >= -1; sign--)
287 {
288 got = mpn_get_d (&n, (mp_size_t) 1, sign, exp);
289 if (got != want)
290 {
291 printf ("mpn_get_d wrong on denorm\n");
292 printf (" n=1\n");
293 printf (" exp %ld\n", exp);
294 printf (" sign %ld\n", (long) sign);
295 d_trace (" got ", got);
296 d_trace (" want ", want);
297 abort ();
298 }
299 want = -want;
300 }
301 want *= 0.5;
302 FORCE_DOUBLE (want);
303 }
304 }
305 else
306 {
307 printf ("Warning, IEEE denorm tests skipped due to SIGFPE (exp=%ld)\n", exp);
308 }
309 tests_sigfpe_done ();
310 }
311
312
313 /* Check values 2^n approaching exponent overflow.
314 Some systems might trap on overflow, so watch out for SIGFPE. */
315 void
316 check_ieee_overflow (void)
317 {
318 static long exp;
319 mp_limb_t n = 1;
320 long i;
321 mp_size_t sign;
322 double want, got;
323
324 if (! _GMP_IEEE_FLOATS)
325 return;
326
327 if (tests_setjmp_sigfpe() == 0)
328 {
329 exp = 1010;
330 want = 1.0;
331 for (i = 0; i < exp; i++)
332 want *= 2.0;
333
334 for ( ; exp < 1050; exp++)
335 {
336 for (sign = 0; sign >= -1; sign--)
337 {
338 got = mpn_get_d (&n, (mp_size_t) 1, sign, exp);
339 if (got != want)
340 {
341 printf ("mpn_get_d wrong on overflow\n");
342 printf (" n=1\n");
343 printf (" exp %ld\n", exp);
344 printf (" sign %ld\n", (long) sign);
345 d_trace (" got ", got);
346 d_trace (" want ", want);
347 abort ();
348 }
349 want = -want;
350 }
351 want *= 2.0;
352 FORCE_DOUBLE (want);
353 }
354 }
355 else
356 {
357 printf ("Warning, IEEE overflow tests skipped due to SIGFPE (exp=%ld)\n", exp);
358 }
359 tests_sigfpe_done ();
360 }
361
362
363 /* ARM gcc 2.95.4 was seen generating bad code for ulong->double
364 conversions, resulting in for instance 0x81c25113 incorrectly converted.
365 This test exercises that value, to see mpn_get_d has avoided the
366 problem. */
367 void
368 check_0x81c25113 (void)
369 {
370 #if GMP_NUMB_BITS >= 32
371 double want = 2176995603.0;
372 double got;
373 mp_limb_t np[4];
374 mp_size_t nsize;
375 long exp;
376
377 if (tests_dbl_mant_bits() < 32)
378 return;
379
380 for (nsize = 1; nsize <= numberof (np); nsize++)
381 {
382 refmpn_zero (np, nsize-1);
383 np[nsize-1] = CNST_LIMB(0x81c25113);
384 exp = - (nsize-1) * GMP_NUMB_BITS;
385 got = mpn_get_d (np, nsize, (mp_size_t) 0, exp);
386 if (got != want)
387 {
388 printf ("mpn_get_d wrong on 2176995603 (0x81c25113)\n");
389 printf (" nsize %ld\n", (long) nsize);
390 printf (" exp %ld\n", exp);
391 d_trace (" got ", got);
392 d_trace (" want ", want);
393 abort ();
394 }
395 }
396 #endif
397 }
398
399
400 void
401 check_rand (void)
402 {
403 gmp_randstate_ptr rands = RANDS;
404 int rep, i;
405 unsigned long mant_bits;
406 long exp, exp_min, exp_max;
407 double got, want, d;
408 mp_size_t nalloc, nsize, sign;
409 mp_limb_t nhigh_mask;
410 mp_ptr np;
411
412 mant_bits = tests_dbl_mant_bits ();
413 if (mant_bits == 0)
414 return;
415
416 /* Allow for vax D format with exponent 127 to -128 only.
417 FIXME: Do something to probe for a valid exponent range. */
418 exp_min = -100 - mant_bits;
419 exp_max = 100 - mant_bits;
420
421 /* space for mant_bits */
422 nalloc = BITS_TO_LIMBS (mant_bits);
423 np = refmpn_malloc_limbs (nalloc);
424 nhigh_mask = MP_LIMB_T_MAX
425 >> (GMP_NAIL_BITS + nalloc * GMP_NUMB_BITS - mant_bits);
426
427 for (rep = 0; rep < 200; rep++)
428 {
429 /* random exp_min to exp_max, inclusive */
430 exp = exp_min + (long) gmp_urandomm_ui (rands, exp_max - exp_min + 1);
431
432 /* mant_bits worth of random at np */
433 if (rep & 1)
434 mpn_random (np, nalloc);
435 else
436 mpn_random2 (np, nalloc);
437 nsize = nalloc;
438 np[nsize-1] &= nhigh_mask;
439 MPN_NORMALIZE (np, nsize);
440 if (nsize == 0)
441 continue;
442
443 sign = (mp_size_t) gmp_urandomb_ui (rands, 1L) - 1;
444
445 /* want = {np,nsize}, converting one bit at a time */
446 want = 0.0;
447 for (i = 0, d = 1.0; i < mant_bits; i++, d *= 2.0)
448 if (np[i/GMP_NUMB_BITS] & (CNST_LIMB(1) << (i%GMP_NUMB_BITS)))
449 want += d;
450 if (sign < 0)
451 want = -want;
452
453 /* want = want * 2^exp */
454 for (i = 0; i < exp; i++)
455 want *= 2.0;
456 for (i = 0; i > exp; i--)
457 want *= 0.5;
458
459 got = mpn_get_d (np, nsize, sign, exp);
460
461 if (got != want)
462 {
463 printf ("mpn_get_d wrong on random data\n");
464 printf (" sign %ld\n", (long) sign);
465 mpn_trace (" n ", np, nsize);
466 printf (" nsize %ld\n", (long) nsize);
467 printf (" exp %ld\n", exp);
468 d_trace (" want ", want);
469 d_trace (" got ", got);
470 abort();
471 }
472 }
473
474 free (np);
475 }
476
477
478 int
479 main (void)
480 {
481 tests_start ();
482 mp_trace_base = -16;
483
484 check_onebit ();
485 check_twobit ();
486 check_inf ();
487 check_underflow ();
488 check_ieee_denorm ();
489 check_ieee_overflow ();
490 check_0x81c25113 ();
491 #if ! (defined (__vax) || defined (__vax__))
492 check_rand ();
493 #endif
494
495 tests_end ();
496 exit (0);
497 }