/* Test mpn_get_d. Copyright 2002-2004 Free Software Foundation, Inc. This file is part of the GNU MP Library test suite. The GNU MP Library test suite is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. The GNU MP Library test suite is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ #include "config.h" #include #include #include #include #include "gmp-impl.h" #include "tests.h" #ifndef _GMP_IEEE_FLOATS #define _GMP_IEEE_FLOATS 0 #endif /* Exercise various 2^n values, with various exponents and positive and negative. */ void check_onebit (void) { static const int bit_table[] = { 0, 1, 2, 3, GMP_NUMB_BITS - 2, GMP_NUMB_BITS - 1, GMP_NUMB_BITS, GMP_NUMB_BITS + 1, GMP_NUMB_BITS + 2, 2 * GMP_NUMB_BITS - 2, 2 * GMP_NUMB_BITS - 1, 2 * GMP_NUMB_BITS, 2 * GMP_NUMB_BITS + 1, 2 * GMP_NUMB_BITS + 2, 3 * GMP_NUMB_BITS - 2, 3 * GMP_NUMB_BITS - 1, 3 * GMP_NUMB_BITS, 3 * GMP_NUMB_BITS + 1, 3 * GMP_NUMB_BITS + 2, 4 * GMP_NUMB_BITS - 2, 4 * GMP_NUMB_BITS - 1, 4 * GMP_NUMB_BITS, 4 * GMP_NUMB_BITS + 1, 4 * GMP_NUMB_BITS + 2, 5 * GMP_NUMB_BITS - 2, 5 * GMP_NUMB_BITS - 1, 5 * GMP_NUMB_BITS, 5 * GMP_NUMB_BITS + 1, 5 * GMP_NUMB_BITS + 2, 6 * GMP_NUMB_BITS - 2, 6 * GMP_NUMB_BITS - 1, 6 * GMP_NUMB_BITS, 6 * GMP_NUMB_BITS + 1, 6 * GMP_NUMB_BITS + 2, }; static const int exp_table[] = { 0, -100, -10, -1, 1, 10, 100, }; /* FIXME: It'd be better to base this on the float format. */ #if defined (__vax) || defined (__vax__) int limit = 127; /* vax fp numbers have limited range */ #else int limit = 511; #endif int bit_i, exp_i, i; double got, want; mp_size_t nsize, sign; long bit, exp, want_bit; mp_limb_t np[20]; for (bit_i = 0; bit_i < numberof (bit_table); bit_i++) { bit = bit_table[bit_i]; nsize = BITS_TO_LIMBS (bit+1); refmpn_zero (np, nsize); np[bit/GMP_NUMB_BITS] = CNST_LIMB(1) << (bit % GMP_NUMB_BITS); for (exp_i = 0; exp_i < numberof (exp_table); exp_i++) { exp = exp_table[exp_i]; want_bit = bit + exp; if (want_bit >= limit || want_bit <= -limit) continue; want = 1.0; for (i = 0; i < want_bit; i++) want *= 2.0; for (i = 0; i > want_bit; i--) want *= 0.5; for (sign = 0; sign >= -1; sign--, want = -want) { got = mpn_get_d (np, nsize, sign, exp); if (got != want) { printf ("mpn_get_d wrong on 2^n\n"); printf (" bit %ld\n", bit); printf (" exp %ld\n", exp); printf (" want_bit %ld\n", want_bit); printf (" sign %ld\n", (long) sign); mpn_trace (" n ", np, nsize); printf (" nsize %ld\n", (long) nsize); d_trace (" want ", want); d_trace (" got ", got); abort(); } } } } } /* Exercise values 2^n+1, while such a value fits the mantissa of a double. */ void check_twobit (void) { int i, mant_bits; double got, want; mp_size_t nsize, sign; mp_ptr np; mant_bits = tests_dbl_mant_bits (); if (mant_bits == 0) return; np = refmpn_malloc_limbs (BITS_TO_LIMBS (mant_bits)); want = 3.0; for (i = 1; i < mant_bits; i++) { nsize = BITS_TO_LIMBS (i+1); refmpn_zero (np, nsize); np[i/GMP_NUMB_BITS] = CNST_LIMB(1) << (i % GMP_NUMB_BITS); np[0] |= 1; for (sign = 0; sign >= -1; sign--) { got = mpn_get_d (np, nsize, sign, 0); if (got != want) { printf ("mpn_get_d wrong on 2^%d + 1\n", i); printf (" sign %ld\n", (long) sign); mpn_trace (" n ", np, nsize); printf (" nsize %ld\n", (long) nsize); d_trace (" want ", want); d_trace (" got ", got); abort(); } want = -want; } want = 2.0 * want - 1.0; } free (np); } /* Expect large negative exponents to underflow to 0.0. Some systems might have hardware traps for such an underflow (though usually it's not the default), so watch out for SIGFPE. */ void check_underflow (void) { static const long exp_table[] = { -999999L, LONG_MIN, }; static const mp_limb_t np[1] = { 1 }; static long exp; mp_size_t nsize, sign; double got; int exp_i; nsize = numberof (np); if (tests_setjmp_sigfpe() == 0) { for (exp_i = 0; exp_i < numberof (exp_table); exp_i++) { exp = exp_table[exp_i]; for (sign = 0; sign >= -1; sign--) { got = mpn_get_d (np, nsize, sign, exp); if (got != 0.0) { printf ("mpn_get_d wrong, didn't get 0.0 on underflow\n"); printf (" nsize %ld\n", (long) nsize); printf (" exp %ld\n", exp); printf (" sign %ld\n", (long) sign); d_trace (" got ", got); abort (); } } } } else { printf ("Warning, underflow to zero tests skipped due to SIGFPE (exp=%ld)\n", exp); } tests_sigfpe_done (); } /* Expect large values to result in +/-inf, on IEEE systems. */ void check_inf (void) { static const long exp_table[] = { 999999L, LONG_MAX, }; static const mp_limb_t np[4] = { 1, 1, 1, 1 }; long exp; mp_size_t nsize, sign, got_sign; double got; int exp_i; if (! _GMP_IEEE_FLOATS) return; for (nsize = 1; nsize <= numberof (np); nsize++) { for (exp_i = 0; exp_i < numberof (exp_table); exp_i++) { exp = exp_table[exp_i]; for (sign = 0; sign >= -1; sign--) { got = mpn_get_d (np, nsize, sign, exp); got_sign = (got >= 0 ? 0 : -1); if (! tests_isinf (got)) { printf ("mpn_get_d wrong, didn't get infinity\n"); bad: printf (" nsize %ld\n", (long) nsize); printf (" exp %ld\n", exp); printf (" sign %ld\n", (long) sign); d_trace (" got ", got); printf (" got sign %ld\n", (long) got_sign); abort (); } if (got_sign != sign) { printf ("mpn_get_d wrong sign on infinity\n"); goto bad; } } } } } /* Check values 2^n approaching and into IEEE denorm range. Some systems might not support denorms, or might have traps setup, so watch out for SIGFPE. */ void check_ieee_denorm (void) { static long exp; mp_limb_t n = 1; long i; mp_size_t sign; double want, got; if (! _GMP_IEEE_FLOATS) return; if (tests_setjmp_sigfpe() == 0) { exp = -1020; want = 1.0; for (i = 0; i > exp; i--) want *= 0.5; for ( ; exp > -1500 && want != 0.0; exp--) { for (sign = 0; sign >= -1; sign--) { got = mpn_get_d (&n, (mp_size_t) 1, sign, exp); if (got != want) { printf ("mpn_get_d wrong on denorm\n"); printf (" n=1\n"); printf (" exp %ld\n", exp); printf (" sign %ld\n", (long) sign); d_trace (" got ", got); d_trace (" want ", want); abort (); } want = -want; } want *= 0.5; FORCE_DOUBLE (want); } } else { printf ("Warning, IEEE denorm tests skipped due to SIGFPE (exp=%ld)\n", exp); } tests_sigfpe_done (); } /* Check values 2^n approaching exponent overflow. Some systems might trap on overflow, so watch out for SIGFPE. */ void check_ieee_overflow (void) { static long exp; mp_limb_t n = 1; long i; mp_size_t sign; double want, got; if (! _GMP_IEEE_FLOATS) return; if (tests_setjmp_sigfpe() == 0) { exp = 1010; want = 1.0; for (i = 0; i < exp; i++) want *= 2.0; for ( ; exp < 1050; exp++) { for (sign = 0; sign >= -1; sign--) { got = mpn_get_d (&n, (mp_size_t) 1, sign, exp); if (got != want) { printf ("mpn_get_d wrong on overflow\n"); printf (" n=1\n"); printf (" exp %ld\n", exp); printf (" sign %ld\n", (long) sign); d_trace (" got ", got); d_trace (" want ", want); abort (); } want = -want; } want *= 2.0; FORCE_DOUBLE (want); } } else { printf ("Warning, IEEE overflow tests skipped due to SIGFPE (exp=%ld)\n", exp); } tests_sigfpe_done (); } /* ARM gcc 2.95.4 was seen generating bad code for ulong->double conversions, resulting in for instance 0x81c25113 incorrectly converted. This test exercises that value, to see mpn_get_d has avoided the problem. */ void check_0x81c25113 (void) { #if GMP_NUMB_BITS >= 32 double want = 2176995603.0; double got; mp_limb_t np[4]; mp_size_t nsize; long exp; if (tests_dbl_mant_bits() < 32) return; for (nsize = 1; nsize <= numberof (np); nsize++) { refmpn_zero (np, nsize-1); np[nsize-1] = CNST_LIMB(0x81c25113); exp = - (nsize-1) * GMP_NUMB_BITS; got = mpn_get_d (np, nsize, (mp_size_t) 0, exp); if (got != want) { printf ("mpn_get_d wrong on 2176995603 (0x81c25113)\n"); printf (" nsize %ld\n", (long) nsize); printf (" exp %ld\n", exp); d_trace (" got ", got); d_trace (" want ", want); abort (); } } #endif } void check_rand (void) { gmp_randstate_ptr rands = RANDS; int rep, i; unsigned long mant_bits; long exp, exp_min, exp_max; double got, want, d; mp_size_t nalloc, nsize, sign; mp_limb_t nhigh_mask; mp_ptr np; mant_bits = tests_dbl_mant_bits (); if (mant_bits == 0) return; /* Allow for vax D format with exponent 127 to -128 only. FIXME: Do something to probe for a valid exponent range. */ exp_min = -100 - mant_bits; exp_max = 100 - mant_bits; /* space for mant_bits */ nalloc = BITS_TO_LIMBS (mant_bits); np = refmpn_malloc_limbs (nalloc); nhigh_mask = MP_LIMB_T_MAX >> (GMP_NAIL_BITS + nalloc * GMP_NUMB_BITS - mant_bits); for (rep = 0; rep < 200; rep++) { /* random exp_min to exp_max, inclusive */ exp = exp_min + (long) gmp_urandomm_ui (rands, exp_max - exp_min + 1); /* mant_bits worth of random at np */ if (rep & 1) mpn_random (np, nalloc); else mpn_random2 (np, nalloc); nsize = nalloc; np[nsize-1] &= nhigh_mask; MPN_NORMALIZE (np, nsize); if (nsize == 0) continue; sign = (mp_size_t) gmp_urandomb_ui (rands, 1L) - 1; /* want = {np,nsize}, converting one bit at a time */ want = 0.0; for (i = 0, d = 1.0; i < mant_bits; i++, d *= 2.0) if (np[i/GMP_NUMB_BITS] & (CNST_LIMB(1) << (i%GMP_NUMB_BITS))) want += d; if (sign < 0) want = -want; /* want = want * 2^exp */ for (i = 0; i < exp; i++) want *= 2.0; for (i = 0; i > exp; i--) want *= 0.5; got = mpn_get_d (np, nsize, sign, exp); if (got != want) { printf ("mpn_get_d wrong on random data\n"); printf (" sign %ld\n", (long) sign); mpn_trace (" n ", np, nsize); printf (" nsize %ld\n", (long) nsize); printf (" exp %ld\n", exp); d_trace (" want ", want); d_trace (" got ", got); abort(); } } free (np); } int main (void) { tests_start (); mp_trace_base = -16; check_onebit (); check_twobit (); check_inf (); check_underflow (); check_ieee_denorm (); check_ieee_overflow (); check_0x81c25113 (); #if ! (defined (__vax) || defined (__vax__)) check_rand (); #endif tests_end (); exit (0); }