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 } 498