htswanalysis/MACS/lib/gsl/gsl-1.11/ieee-utils/test.c

   1 /* ieee-utils/test.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <math.h>
  23 #include <float.h>
  24 #include <string.h>
  25 #include <gsl/gsl_ieee_utils.h>
  26 #include <gsl/gsl_test.h>
  27
  28 #if HAVE_IRIX_IEEE_INTERFACE
  29 /* don't test denormals on IRIX */
  30 #else
  31 #if HAVE_IEEE_DENORMALS
  32 #define TEST_DENORMAL 1
  33 #endif
  34 #endif
  35
  36 #ifndef FLT_MIN
  37 #define FLT_MIN 1.17549435e-38f
  38 #endif
  39
  40 #ifndef FLT_MAX
  41 #define FLT_MAX 3.40282347e+38f
  42 #endif
  43
  44 #ifndef DBL_MIN
  45 #define DBL_MIN 2.2250738585072014e-308
  46 #endif
  47
  48 #ifndef DBL_MAX
  49 #define DBL_MAX 1.7976931348623157e+308
  50 #endif
  51
  52 int
  53 main (void)
  54 {
  55   float zerof = 0.0f, minus_onef = -1.0f ;
  56   double zero = 0.0, minus_one = -1.0 ;
  57
  58   /* Check for +ZERO (float) */
  59
  60   {
  61     float f = 0.0f;
  62     const char mantissa[] = "00000000000000000000000";
  63     gsl_ieee_float_rep r;
  64     gsl_ieee_float_to_rep (&f, &r);
  65
  66     gsl_test_int (r.sign, 0, "float x = 0, sign is +");
  67     gsl_test_int (r.exponent, -127, "float x = 0, exponent is -127");
  68     gsl_test_str (r.mantissa, mantissa, "float x = 0, mantissa");
  69     gsl_test_int (r.type, GSL_IEEE_TYPE_ZERO, "float x = 0, type is ZERO");
  70   }
  71
  72   /* Check for -ZERO (float) */
  73
  74   {
  75     float f = minus_onef;
  76     const char mantissa[] = "00000000000000000000000";
  77     gsl_ieee_float_rep r;
  78
  79     while (f < 0) {
  80       f *= 0.1f;
  81     }
  82
  83     gsl_ieee_float_to_rep (&f, &r);
  84
  85     gsl_test_int (r.sign, 1, "float x = -1*0, sign is -");
  86     gsl_test_int (r.exponent, -127, "float x = -1*0, exponent is -127");
  87     gsl_test_str (r.mantissa, mantissa, "float x = -1*0, mantissa");
  88     gsl_test_int (r.type, GSL_IEEE_TYPE_ZERO, "float x = -1*0, type is ZERO");
  89   }
  90
  91   /* Check for a positive NORMAL number (e.g. 2.1) (float) */
  92
  93   {
  94     float f = 2.1f;
  95     const char mantissa[] = "00001100110011001100110";
  96
  97     gsl_ieee_float_rep r;
  98     gsl_ieee_float_to_rep (&f, &r);
  99
 100     gsl_test_int (r.sign, 0, "float x = 2.1, sign is +");
 101     gsl_test_int (r.exponent, 1, "float x = 2.1, exponent is 1");
 102     gsl_test_str (r.mantissa, mantissa, "float x = 2.1, mantissa");
 103     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "float x = 2.1, type is NORMAL");
 104   }
 105
 106
 107   /* Check for a negative NORMAL number (e.g. -1.3304...) (float) */
 108
 109   {
 110     float f = -1.3303577090924210f ;
 111     const char mantissa[] = "01010100100100100101001";
 112
 113     gsl_ieee_float_rep r;
 114     gsl_ieee_float_to_rep (&f, &r);
 115
 116     gsl_test_int (r.sign, 1, "float x = -1.3304..., sign is -");
 117     gsl_test_int (r.exponent, 0, "float x = -1.3304..., exponent is 0");
 118     gsl_test_str (r.mantissa, mantissa, "float x = -1.3304..., mantissa");
 119     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
 120                   "float x = -1.3304..., type is NORMAL");
 121   }
 122
 123   /* Check for a large positive NORMAL number (e.g. 3.37e31) (float) */
 124
 125   {
 126     float f = 3.37e31f;
 127     const char mantissa[] = "10101001010110101001001";
 128     gsl_ieee_float_rep r;
 129     gsl_ieee_float_to_rep (&f, &r);
 130
 131     gsl_test_int (r.sign, 0, "float x = 3.37e31, sign is +");
 132     gsl_test_int (r.exponent, 104, "float x = 3.37e31, exponent is 104");
 133     gsl_test_str (r.mantissa, mantissa, "float x = 3.37e31, mantissa");
 134     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "float x = 3.37e31, type is NORMAL");
 135   }
 136
 137   /* Check for a small positive NORMAL number (e.g. 3.37e-31) (float) */
 138
 139   {
 140     float f = 3.37e-31f;
 141     const char mantissa[] = "10110101011100110111011";
 142
 143     gsl_ieee_float_rep r;
 144     gsl_ieee_float_to_rep (&f, &r);
 145
 146     gsl_test_int (r.sign, 0, "float x = 3.37e-31, sign is +");
 147     gsl_test_int (r.exponent, -102, "float x = 3.37e-31, exponent is -102");
 148     gsl_test_str (r.mantissa, mantissa, "float x = 3.37e-31, mantissa");
 149     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
 150                   "float x = 3.37e-31, type is NORMAL");
 151   }
 152
 153   /* Check for FLT_MIN (smallest possible number that is not denormal) */
 154
 155   {
 156     float f = FLT_MIN;
 157     const char mantissa[] = "00000000000000000000000";
 158     gsl_ieee_float_rep r;
 159     gsl_ieee_float_to_rep (&f, &r);
 160
 161     gsl_test_int (r.sign, 0, "float x = FLT_MIN, sign is +");
 162     gsl_test_int (r.exponent, -126, "float x = FLT_MIN, exponent is -126");
 163     gsl_test_str (r.mantissa, mantissa, "float x = FLT_MIN, mantissa");
 164     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "float x = FLT_MIN, type is NORMAL");
 165   }
 166
 167   /* Check for FLT_MAX (largest possible number that is not Inf) */
 168
 169   {
 170     float f = FLT_MAX;
 171     const char mantissa[] = "11111111111111111111111";
 172
 173     gsl_ieee_float_rep r;
 174     gsl_ieee_float_to_rep (&f, &r);
 175
 176     gsl_test_int (r.sign, 0, "float x = FLT_MAX, sign is +");
 177     gsl_test_int (r.exponent, 127, "float x = FLT_MAX, exponent is 127");
 178     gsl_test_str (r.mantissa, mantissa, "float x = FLT_MAX, mantissa");
 179     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "float x = FLT_MAX, type is NORMAL");
 180   }
 181
 182
 183   /* Check for DENORMAL numbers (e.g. FLT_MIN/2^n) */
 184
 185 #ifdef TEST_DENORMAL
 186   {
 187     float f = FLT_MIN;
 188     char mantissa[] = "10000000000000000000000";
 189
 190     int i;
 191     gsl_ieee_float_rep r;
 192
 193     for (i = 0; i < 23; i++)
 194       {
 195         float x = f / (float)pow (2.0, 1 + (float) i);
 196         mantissa[i] = '1';
 197         gsl_ieee_float_to_rep (&x, &r);
 198
 199         gsl_test_int (r.sign, 0, "float x = FLT_MIN/2^%d, sign is +", i + 1);
 200         gsl_test_int (r.exponent, -127,
 201                       "float x = FLT_MIN/2^%d, exponent is -127", i + 1);
 202         gsl_test_str (r.mantissa, mantissa,
 203                       "float x = FLT_MIN/2^%d, mantissa", i + 1);
 204         gsl_test_int (r.type, GSL_IEEE_TYPE_DENORMAL,
 205                       "float x = FLT_MIN/2^%d, type is DENORMAL", i + 1);
 206         mantissa[i] = '0';
 207       }
 208   }
 209 #endif
 210
 211   /* Check for positive INFINITY (e.g. 2*FLT_MAX) */
 212
 213   {
 214     float f = FLT_MAX;
 215     const char mantissa[] = "00000000000000000000000";
 216
 217     gsl_ieee_float_rep r;
 218
 219     float x;
 220     x = 2 * f;
 221     gsl_ieee_float_to_rep (&x, &r);
 222
 223     gsl_test_int (r.sign, 0, "float x = 2*FLT_MAX, sign is +");
 224     gsl_test_int (r.exponent, 128, "float x = 2*FLT_MAX, exponent is 128");
 225     gsl_test_str (r.mantissa, mantissa, "float x = 2*FLT_MAX, mantissa");
 226     gsl_test_int (r.type, GSL_IEEE_TYPE_INF, "float x = -2*FLT_MAX, type is INF");
 227   }
 228
 229   /* Check for negative INFINITY (e.g. -2*FLT_MAX) */
 230
 231   {
 232     float f = FLT_MAX;
 233     const char mantissa[] = "00000000000000000000000";
 234
 235     gsl_ieee_float_rep r;
 236
 237     float x;
 238     x = -2 * f;
 239     gsl_ieee_float_to_rep (&x, &r);
 240
 241     gsl_test_int (r.sign, 1, "float x = -2*FLT_MAX, sign is -");
 242     gsl_test_int (r.exponent, 128, "float x = -2*FLT_MAX, exponent is 128");
 243     gsl_test_str (r.mantissa, mantissa, "float x = -2*FLT_MAX, mantissa");
 244     gsl_test_int (r.type, GSL_IEEE_TYPE_INF, "float x = -2*FLT_MAX, type is INF");
 245   }
 246
 247   /* Check for NAN (e.g. Inf - Inf) (float) */
 248
 249   {
 250     gsl_ieee_float_rep r;
 251     float x = 1.0f, y = 2.0f, z = zerof;
 252
 253     x = x / z;
 254     y = y / z;
 255     z = y - x;
 256
 257     gsl_ieee_float_to_rep (&z, &r);
 258
 259     /* We don't check the sign and we don't check the mantissa because
 260        they could be anything for a NaN */
 261
 262     gsl_test_int (r.exponent, 128, "float x = NaN, exponent is 128");
 263     gsl_test_int (r.type, GSL_IEEE_TYPE_NAN, "float x = NaN, type is NAN");
 264   }
 265
 266
 267   /* Check for +ZERO */
 268
 269   {
 270     double d = 0.0;
 271     const char mantissa[]
 272       = "0000000000000000000000000000000000000000000000000000";
 273     gsl_ieee_double_rep r;
 274     gsl_ieee_double_to_rep (&d, &r);
 275
 276     gsl_test_int (r.sign, 0, "double x = 0, sign is +");
 277     gsl_test_int (r.exponent, -1023, "double x = 0, exponent is -1023");
 278     gsl_test_str (r.mantissa, mantissa, "double x = 0, mantissa");
 279     gsl_test_int (r.type, GSL_IEEE_TYPE_ZERO, "double x = 0, type is ZERO");
 280   }
 281
 282   /* Check for -ZERO */
 283
 284   {
 285     double d =  minus_one;
 286     const char mantissa[]
 287       = "0000000000000000000000000000000000000000000000000000";
 288     gsl_ieee_double_rep r;
 289
 290     while (d < 0) {
 291       d *= 0.1;
 292     }
 293
 294     gsl_ieee_double_to_rep (&d, &r);
 295
 296     gsl_test_int (r.sign, 1, "double x = -1*0, sign is -");
 297     gsl_test_int (r.exponent, -1023, "double x = -1*0, exponent is -1023");
 298     gsl_test_str (r.mantissa, mantissa, "double x = -1*0, mantissa");
 299     gsl_test_int (r.type, GSL_IEEE_TYPE_ZERO, "double x = -1*0, type is ZERO");
 300   }
 301
 302   /* Check for a positive NORMAL number (e.g. 2.1) */
 303
 304   {
 305     double d = 2.1;
 306     const char mantissa[]
 307       = "0000110011001100110011001100110011001100110011001101";
 308     gsl_ieee_double_rep r;
 309     gsl_ieee_double_to_rep (&d, &r);
 310
 311     gsl_test_int (r.sign, 0, "double x = 2.1, sign is +");
 312     gsl_test_int (r.exponent, 1, "double x = 2.1, exponent is 1");
 313     gsl_test_str (r.mantissa, mantissa, "double x = 2.1, mantissa");
 314     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL, "double x = 2.1, type is NORMAL");
 315   }
 316
 317
 318   /* Check for a negative NORMAL number (e.g. -1.3304...) */
 319
 320   {
 321     double d = -1.3303577090924210146738460025517269968986511230468750;
 322     const char mantissa[]
 323       = "0101010010010010010100101010010010001000100011101110";
 324     gsl_ieee_double_rep r;
 325     gsl_ieee_double_to_rep (&d, &r);
 326
 327     gsl_test_int (r.sign, 1, "double x = -1.3304..., sign is -");
 328     gsl_test_int (r.exponent, 0, "double x = -1.3304..., exponent is 0");
 329     gsl_test_str (r.mantissa, mantissa, "double x = -1.3304..., mantissa");
 330     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
 331                   "double x = -1.3304..., type is NORMAL");
 332   }
 333
 334   /* Check for a large positive NORMAL number (e.g. 3.37e297) */
 335
 336   {
 337     double d = 3.37e297;
 338     const char mantissa[]
 339       = "0100100111001001100101111001100000100110011101000100";
 340     gsl_ieee_double_rep r;
 341     gsl_ieee_double_to_rep (&d, &r);
 342
 343     gsl_test_int (r.sign, 0, "double x = 3.37e297, sign is +");
 344     gsl_test_int (r.exponent, 988, "double x = 3.37e297, exponent is 998");
 345     gsl_test_str (r.mantissa, mantissa, "double x = 3.37e297, mantissa");
 346     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
 347                   "double x = 3.37e297, type is NORMAL");
 348   }
 349
 350   /* Check for a small positive NORMAL number (e.g. 3.37e-297) */
 351
 352   {
 353     double d = 3.37e-297;
 354     const char mantissa[]
 355     = "0001101000011011101011100001110010100001001100110111";
 356     gsl_ieee_double_rep r;
 357     gsl_ieee_double_to_rep (&d, &r);
 358
 359     gsl_test_int (r.sign, 0, "double x = 3.37e-297, sign is +");
 360     gsl_test_int (r.exponent, -985, "double x = 3.37e-297, exponent is -985");
 361     gsl_test_str (r.mantissa, mantissa, "double x = 3.37e-297, mantissa");
 362     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
 363                   "double x = 3.37e-297, type is NORMAL");
 364   }
 365
 366   /* Check for DBL_MIN (smallest possible number that is not denormal) */
 367
 368   {
 369     double d = DBL_MIN;
 370     const char mantissa[]
 371       = "0000000000000000000000000000000000000000000000000000";
 372     gsl_ieee_double_rep r;
 373     gsl_ieee_double_to_rep (&d, &r);
 374
 375     gsl_test_int (r.sign, 0, "double x = DBL_MIN, sign is +");
 376     gsl_test_int (r.exponent, -1022, "double x = DBL_MIN, exponent is -1022");
 377     gsl_test_str (r.mantissa, mantissa, "double x = DBL_MIN, mantissa");
 378     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
 379                   "double x = DBL_MIN, type is NORMAL");
 380   }
 381
 382   /* Check for DBL_MAX (largest possible number that is not Inf) */
 383
 384   {
 385     double d = DBL_MAX;
 386     const char mantissa[]
 387     = "1111111111111111111111111111111111111111111111111111";
 388     gsl_ieee_double_rep r;
 389     gsl_ieee_double_to_rep (&d, &r);
 390
 391     gsl_test_int (r.sign, 0, "double x = DBL_MAX, sign is +");
 392     gsl_test_int (r.exponent, 1023, "double x = DBL_MAX, exponent is 1023");
 393     gsl_test_str (r.mantissa, mantissa, "double x = DBL_MAX, mantissa");
 394     gsl_test_int (r.type, GSL_IEEE_TYPE_NORMAL,
 395                   "double x = DBL_MAX, type is NORMAL");
 396   }
 397
 398   /* Check for DENORMAL numbers (e.g. DBL_MIN/2^n) */
 399
 400 #ifdef TEST_DENORMAL
 401   {
 402     double d = DBL_MIN;
 403     char mantissa[]
 404       = "1000000000000000000000000000000000000000000000000000";
 405     int i;
 406     gsl_ieee_double_rep r;
 407
 408     for (i = 0; i < 52; i++)
 409       {
 410         double x = d / pow (2.0, 1 + (double) i);
 411         mantissa[i] = '1';
 412         gsl_ieee_double_to_rep (&x, &r);
 413
 414         gsl_test_int (r.sign, 0, "double x = DBL_MIN/2^%d, sign is +", i + 1);
 415         gsl_test_int (r.exponent, -1023,
 416                       "double x = DBL_MIN/2^%d, exponent", i + 1);
 417         gsl_test_str (r.mantissa, mantissa,
 418                       "double x = DBL_MIN/2^%d, mantissa", i + 1);
 419         gsl_test_int (r.type, GSL_IEEE_TYPE_DENORMAL,
 420                       "double x = DBL_MIN/2^%d, type is DENORMAL", i + 1);
 421         mantissa[i] = '0';
 422       }
 423   }
 424 #endif
 425
 426   /* Check for positive INFINITY (e.g. 2*DBL_MAX) */
 427
 428   {
 429     double d = DBL_MAX;
 430     const char mantissa[]
 431       = "0000000000000000000000000000000000000000000000000000";
 432     gsl_ieee_double_rep r;
 433
 434     double x;
 435     x = 2.0 * d;
 436     gsl_ieee_double_to_rep (&x, &r);
 437
 438     gsl_test_int (r.sign, 0, "double x = 2*DBL_MAX, sign is +");
 439     gsl_test_int (r.exponent, 1024, "double x = 2*DBL_MAX, exponent is 1024");
 440     gsl_test_str (r.mantissa, mantissa, "double x = 2*DBL_MAX, mantissa");
 441     gsl_test_int (r.type, GSL_IEEE_TYPE_INF, "double x = 2*DBL_MAX, type is INF");
 442   }
 443
 444   /* Check for negative INFINITY (e.g. -2*DBL_MAX) */
 445
 446   {
 447     double d = DBL_MAX;
 448     const char mantissa[]
 449       = "0000000000000000000000000000000000000000000000000000";
 450     gsl_ieee_double_rep r;
 451
 452     double x;
 453     x = -2.0 * d;
 454     gsl_ieee_double_to_rep (&x, &r);
 455
 456     gsl_test_int (r.sign, 1, "double x = -2*DBL_MAX, sign is -");
 457     gsl_test_int (r.exponent, 1024, "double x = -2*DBL_MAX, exponent is 1024");
 458     gsl_test_str (r.mantissa, mantissa, "double x = -2*DBL_MAX, mantissa");
 459     gsl_test_int (r.type, GSL_IEEE_TYPE_INF,"double x = -2*DBL_MAX, type is INF");
 460   }
 461
 462   /* Check for NAN (e.g. Inf - Inf) */
 463
 464   {
 465     gsl_ieee_double_rep r;
 466     double x = 1.0, y = 2.0, z = zero;
 467
 468     x = x / z;
 469     y = y / z;
 470     z = y - x;
 471
 472     gsl_ieee_double_to_rep (&z, &r);
 473
 474     /* We don't check the sign and we don't check the mantissa because
 475        they could be anything for a NaN */
 476
 477     gsl_test_int (r.exponent, 1024, "double x = NaN, exponent is 1024");
 478     gsl_test_int (r.type, GSL_IEEE_TYPE_NAN, "double x = NaN, type is NAN");
 479   }
 480
 481   exit (gsl_test_summary ());
 482 }