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