t-get_d.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:13k
- /* 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 <limits.h> 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 <setjmp.h>
- #include <signal.h>
- #include <stdio.h>
- #include <stdlib.h>
- #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^nn");
- printf (" bit %ldn", bit);
- printf (" exp %ldn", exp);
- printf (" want_bit %ldn", want_bit);
- printf (" sign %ldn", (long) sign);
- mpn_trace (" n ", np, nsize);
- printf (" nsize %ldn", (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 + 1n", i);
- printf (" sign %ldn", (long) sign);
- mpn_trace (" n ", np, nsize);
- printf (" nsize %ldn", (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 underflown");
- printf (" nsize %ldn", (long) nsize);
- printf (" exp %ldn", exp);
- printf (" sign %ldn", (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 infinityn");
- bad:
- printf (" nsize %ldn", (long) nsize);
- printf (" exp %ldn", exp);
- printf (" sign %ldn", (long) sign);
- d_trace (" got ", got);
- printf (" got sign %ldn", (long) got_sign);
- abort ();
- }
- if (got_sign != sign)
- {
- printf ("mpn_get_d wrong sign on infinityn");
- 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 denormn");
- printf (" n=1n");
- printf (" exp %ldn", exp);
- printf (" sign %ldn", (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 overflown");
- printf (" n=1n");
- printf (" exp %ldn", exp);
- printf (" sign %ldn", (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 %ldn", (long) nsize);
- printf (" exp %ldn", 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 datan");
- printf (" sign %ldn", (long) sign);
- mpn_trace (" n ", np, nsize);
- printf (" nsize %ldn", (long) nsize);
- printf (" exp %ldn", 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);
- }