tests/mpn/t-get_d.c

   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.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library.  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 #ifdef __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   check_rand ();
 499
 500   tests_end ();
 501   exit (0);
 502 }