t-get_d.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:13k
源码类别:

数学计算

开发平台:

Unix_Linux

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