X-Git-Url: https://oss.titaniummirror.com/gitweb?a=blobdiff_plain;f=gmp%2Ftests%2Fmpn%2Ft-get_d.c;fp=gmp%2Ftests%2Fmpn%2Ft-get_d.c;h=c3999e2a7c5eb09664f84cb645890e1a8d1cce48;hb=6fed43773c9b0ce596dca5686f37ac3fc0fa11c0;hp=0000000000000000000000000000000000000000;hpb=27b11d56b743098deb193d510b337ba22dc52e5c;p=msp430-gcc.git diff --git a/gmp/tests/mpn/t-get_d.c b/gmp/tests/mpn/t-get_d.c new file mode 100644 index 00000000..c3999e2a --- /dev/null +++ b/gmp/tests/mpn/t-get_d.c @@ -0,0 +1,502 @@ +/* Test mpn_get_d. + +Copyright 2002, 2003, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser 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 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 Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ + +/* Note that we don't use for LONG_MIN, but instead our own + definition in gmp-impl.h. In gcc 2.95.4 (debian 3.0) under + -mcpu=ultrasparc, limits.h sees __sparc_v9__ defined and assumes that + means long is 64-bit long, but it's only 32-bits, causing fatal compile + errors. */ + +#include "config.h" + +#include +#include +#include +#include + +#include "gmp.h" +#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. */ +#ifdef __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 (); + check_rand (); + + tests_end (); + exit (0); +}