ghc/rts/StgPrimFloat.c

   1 /* -----------------------------------------------------------------------------
   2  * $Id: StgPrimFloat.c,v 1.4 1999/02/18 12:26:12 simonm Exp $
   3  *
   4  * (c) The GHC Team, 1998-1999
   5  *
   6  * Miscellaneous support for floating-point primitives
   7  *
   8  * ---------------------------------------------------------------------------*/
   9
  10 #include "Rts.h"
  11
  12 /*
  13  * Encoding and decoding Doubles.  Code based on the HBC code
  14  * (lib/fltcode.c).
  15  */
  16
  17 #define GMP_BASE 4294967296.0
  18 #if FLOATS_AS_DOUBLES /* defined in StgTypes.h */
  19 #define DNBIGIT 1   /* mantissa of a double will fit in one long */
  20 #else
  21 #define DNBIGIT  2  /* mantissa of a double will fit in two longs */
  22 #endif
  23 #define FNBIGIT  1  /* for float, one long */
  24
  25 #if IEEE_FLOATING_POINT
  26 #define MY_DMINEXP  ((DBL_MIN_EXP) - (DBL_MANT_DIG) - 1)
  27 /* DMINEXP is defined in values.h on Linux (for example) */
  28 #define DHIGHBIT 0x00100000
  29 #define DMSBIT   0x80000000
  30
  31 #define MY_FMINEXP  ((FLT_MIN_EXP) - (FLT_MANT_DIG) - 1)
  32 #define FHIGHBIT 0x00800000
  33 #define FMSBIT   0x80000000
  34 #endif
  35
  36 #ifdef WORDS_BIGENDIAN
  37 #define L 1
  38 #define H 0
  39 #else
  40 #define L 0
  41 #define H 1
  42 #endif
  43
  44 #define __abs(a)                (( (a) >= 0 ) ? (a) : (-(a)))
  45
  46 StgDouble
  47 __encodeDouble (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */
  48 {
  49     StgDouble r;
  50     W_ *arr = (W_ *)ba;
  51     I_ i;
  52
  53     /* Convert MP_INT to a double; knows a lot about internal rep! */
  54     i = __abs(size)-1;
  55     if (i < 0) {
  56       r = 0.0;
  57     } else {
  58       for (r = arr[i], i--; i >= 0; i--)
  59         r = r * GMP_BASE + arr[i];
  60     }
  61
  62     /* Now raise to the exponent */
  63     if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
  64         r = ldexp(r, e);
  65
  66     /* sign is encoded in the size */
  67     if (size < 0)
  68         r = -r;
  69
  70     return r;
  71 }
  72
  73 /* Special version for small Integers */
  74 StgDouble
  75 __int_encodeDouble (I_ j, I_ e)
  76 {
  77   StgDouble r;
  78
  79   r = (StgDouble)__abs(j);
  80
  81   /* Now raise to the exponent */
  82   if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
  83     r = ldexp(r, e);
  84
  85   /* sign is encoded in the size */
  86   if (j < 0)
  87     r = -r;
  88
  89   return r;
  90 }
  91
  92 #if ! FLOATS_AS_DOUBLES
  93 StgFloat
  94 __encodeFloat (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */
  95 {
  96     StgFloat r;
  97     W_ *arr = (W_ *)ba;
  98     I_ i;
  99
 100     /* Convert MP_INT to a float; knows a lot about internal rep! */
 101     for(r = 0.0, i = __abs(size); i >= 0; i--)
 102         r = (r * GMP_BASE) + arr[i];
 103
 104     /* Now raise to the exponent */
 105     if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
 106         r = ldexp(r, e);
 107
 108     /* sign is encoded in the size */
 109     if (size < 0)
 110         r = -r;
 111
 112     return r;
 113 }
 114
 115 /* Special version for small Integers */
 116 StgFloat
 117 __int_encodeFloat (I_ j, I_ e)
 118 {
 119   StgFloat r;
 120
 121   r = (StgFloat)__abs(j);
 122
 123   /* Now raise to the exponent */
 124   if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
 125     r = ldexp(r, e);
 126
 127   /* sign is encoded in the size */
 128   if (j < 0)
 129     r = -r;
 130
 131   return r;
 132 }
 133
 134 #endif  /* FLOATS_AS_DOUBLES */
 135
 136 /* This only supports IEEE floating point */
 137
 138 void
 139 __decodeDouble (MP_INT *man, I_ *exp, StgDouble dbl)
 140 {
 141     /* Do some bit fiddling on IEEE */
 142     nat low, high;              /* assuming 32 bit ints */
 143     int sign, iexp;
 144     union { double d; int i[2]; } u;    /* assuming 32 bit ints, 64 bit double */
 145
 146     u.d = dbl;      /* grab chunks of the double */
 147     low = u.i[L];
 148     high = u.i[H];
 149
 150     /* we know the MP_INT* passed in has size zero, so we realloc
 151         no matter what.
 152     */
 153     man->_mp_alloc = DNBIGIT;
 154
 155     if (low == 0 && (high & ~DMSBIT) == 0) {
 156         man->_mp_size = 0;
 157         *exp = 0L;
 158     } else {
 159         man->_mp_size = DNBIGIT;
 160         iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
 161         sign = high;
 162
 163         high &= DHIGHBIT-1;
 164         if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
 165             high |= DHIGHBIT;
 166         else {
 167             iexp++;
 168             /* A denorm, normalize the mantissa */
 169             while (! (high & DHIGHBIT)) {
 170                 high <<= 1;
 171                 if (low & DMSBIT)
 172                     high++;
 173                 low <<= 1;
 174                 iexp--;
 175             }
 176         }
 177         *exp = (I_) iexp;
 178 #if DNBIGIT == 2
 179         man->_mp_d[0] = low;
 180         man->_mp_d[1] = high;
 181 #else
 182 #if DNBIGIT == 1
 183         man->_mp_d[0] = ((unsigned long)high) << 32 | (unsigned long)low;
 184 #else
 185         error : error : error : Cannae cope with DNBIGIT
 186 #endif
 187 #endif
 188         if (sign < 0)
 189             man->_mp_size = -man->_mp_size;
 190     }
 191 }
 192
 193 #if ! FLOATS_AS_DOUBLES
 194 void
 195 __decodeFloat (MP_INT *man, I_ *exp, StgFloat flt)
 196 {
 197     /* Do some bit fiddling on IEEE */
 198     int high, sign;                 /* assuming 32 bit ints */
 199     union { float f; int i; } u;    /* assuming 32 bit float and int */
 200
 201     u.f = flt;      /* grab the float */
 202     high = u.i;
 203
 204     /* we know the MP_INT* passed in has size zero, so we realloc
 205         no matter what.
 206     */
 207     man->_mp_alloc = FNBIGIT;
 208
 209     if ((high & ~FMSBIT) == 0) {
 210         man->_mp_size = 0;
 211         *exp = 0;
 212     } else {
 213         man->_mp_size = FNBIGIT;
 214         *exp = ((high >> 23) & 0xff) + MY_FMINEXP;
 215         sign = high;
 216
 217         high &= FHIGHBIT-1;
 218         if (*exp != MY_FMINEXP) /* don't add hidden bit to denorms */
 219             high |= FHIGHBIT;
 220         else {
 221             (*exp)++;
 222             /* A denorm, normalize the mantissa */
 223             while (! (high & FHIGHBIT)) {
 224                 high <<= 1;
 225                 (*exp)--;
 226             }
 227         }
 228 #if FNBIGIT == 1
 229         man->_mp_d[0] = high;
 230 #else
 231         error : error : error : Cannae cope with FNBIGIT
 232 #endif
 233         if (sign < 0)
 234             man->_mp_size = -man->_mp_size;
 235     }
 236 }
 237 #endif  /* FLOATS_AS_DOUBLES */
 238
 239 /* Convenient union types for checking the layout of IEEE 754 types -
 240    based on defs in GNU libc <ieee754.h>
 241 */
 242
 243 union stg_ieee754_flt
 244 {
 245    float f;
 246    struct {
 247
 248 #if WORDS_BIGENDIAN
 249         unsigned int negative:1;
 250         unsigned int exponent:8;
 251         unsigned int mantissa:23;
 252 #else
 253         unsigned int mantissa:23;
 254         unsigned int exponent:8;
 255         unsigned int negative:1;
 256 #endif
 257    } ieee;
 258    struct {
 259
 260 #if WORDS_BIGENDIAN
 261         unsigned int negative:1;
 262         unsigned int exponent:8;
 263         unsigned int quiet_nan:1;
 264         unsigned int mantissa:22;
 265 #else
 266         unsigned int mantissa:22;
 267         unsigned int quiet_nan:1;
 268         unsigned int exponent:8;
 269         unsigned int negative:1;
 270 #endif
 271    } ieee_nan;
 272 };
 273
 274 /*
 275
 276  To recap, here's the representation of a double precision
 277  IEEE floating point number:
 278
 279  sign         63           sign bit (0==positive, 1==negative)
 280  exponent     62-52        exponent (biased by 1023)
 281  fraction     51-0         fraction (bits to right of binary point)
 282 */
 283
 284 union stg_ieee754_dbl
 285 {
 286    double d;
 287    struct {
 288
 289 #if WORDS_BIGENDIAN
 290         unsigned int negative:1;
 291         unsigned int exponent:11;
 292         unsigned int mantissa0:20;
 293         unsigned int mantissa1:32;
 294 #else
 295         unsigned int mantissa1:32;
 296         unsigned int mantissa0:20;
 297         unsigned int exponent:11;
 298         unsigned int negative:1;
 299 #endif
 300    } ieee;
 301     /* This format makes it easier to see if a NaN is a signalling NaN.  */
 302    struct {
 303
 304 #if WORDS_BIGENDIAN
 305         unsigned int negative:1;
 306         unsigned int exponent:11;
 307         unsigned int quiet_nan:1;
 308         unsigned int mantissa0:19;
 309         unsigned int mantissa1:32;
 310 #else
 311         unsigned int mantissa1:32;
 312         unsigned int mantissa0:19;
 313         unsigned int quiet_nan:1;
 314         unsigned int exponent:11;
 315         unsigned int negative:1;
 316 #endif
 317    } ieee_nan;
 318 };
 319
 320 /*
 321  * Predicates for testing for extreme IEEE fp values. Used
 322  * by the bytecode evaluator and the Prelude.
 323  *
 324  */
 325
 326 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
 327 #ifdef IEEE_FLOATING_POINT
 328
 329 StgInt
 330 isDoubleNaN(d)
 331 StgDouble d;
 332 {
 333   union stg_ieee754_dbl u;
 334
 335   u.d = d;
 336
 337   return (
 338     u.ieee.exponent  == 2047 /* 2^11 - 1 */ &&  /* Is the exponent all ones? */
 339     (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)
 340         /* and the mantissa non-zero? */
 341     );
 342 }
 343
 344 StgInt
 345 isDoubleInfinite(d)
 346 StgDouble d;
 347 {
 348     union stg_ieee754_dbl u;
 349
 350     u.d = d;
 351
 352     /* Inf iff exponent is all ones, mantissa all zeros */
 353     return (
 354         u.ieee.exponent  == 2047 /* 2^11 - 1 */ &&
 355         u.ieee.mantissa0 == 0                   &&
 356         u.ieee.mantissa1 == 0
 357       );
 358 }
 359
 360 StgInt
 361 isDoubleDenormalized(d)
 362 StgDouble d;
 363 {
 364     union stg_ieee754_dbl u;
 365
 366     u.d = d;
 367
 368     /* A (single/double/quad) precision floating point number
 369        is denormalised iff:
 370         - exponent is zero
 371         - mantissa is non-zero.
 372         - (don't care about setting of sign bit.)
 373
 374     */
 375     return (
 376         u.ieee.exponent  == 0 &&
 377         (u.ieee.mantissa0 != 0 ||
 378          u.ieee.mantissa1 != 0)
 379       );
 380
 381 }
 382
 383 StgInt
 384 isDoubleNegativeZero(d)
 385 StgDouble d;
 386 {
 387     union stg_ieee754_dbl u;
 388
 389     u.d = d;
 390     /* sign (bit 63) set (only) => negative zero */
 391
 392     return (
 393         u.ieee.negative  == 1 &&
 394         u.ieee.exponent  == 0 &&
 395         u.ieee.mantissa0 == 0 &&
 396         u.ieee.mantissa1 == 0);
 397 }
 398
 399 /* Same tests, this time for StgFloats. */
 400
 401 /*
 402  To recap, here's the representation of a single precision
 403  IEEE floating point number:
 404
 405  sign         31           sign bit (0 == positive, 1 == negative)
 406  exponent     30-23        exponent (biased by 127)
 407  fraction     22-0         fraction (bits to right of binary point)
 408 */
 409
 410
 411 StgInt
 412 isFloatNaN(f)
 413 StgFloat f;
 414 {
 415 # ifdef FLOATS_AS_DOUBLES
 416     return (isDoubleNaN(f));
 417 # else
 418     union stg_ieee754_flt u;
 419     u.f = f;
 420
 421    /* Floating point NaN iff exponent is all ones, mantissa is
 422       non-zero (but see below.) */
 423    return (
 424         u.ieee.exponent == 255 /* 2^8 - 1 */ &&
 425         u.ieee.mantissa != 0);
 426
 427 # endif /* !FLOATS_AS_DOUBLES */
 428 }
 429
 430 StgInt
 431 isFloatInfinite(f)
 432 StgFloat f;
 433 {
 434 # ifdef FLOATS_AS_DOUBLES
 435     return (isDoubleInfinite(f));
 436 # else
 437     union stg_ieee754_flt u;
 438     u.f = f;
 439
 440     /* A float is Inf iff exponent is max (all ones),
 441        and mantissa is min(all zeros.) */
 442     return (
 443         u.ieee.exponent == 255 /* 2^8 - 1 */ &&
 444         u.ieee.mantissa == 0);
 445 # endif /* !FLOATS_AS_DOUBLES */
 446 }
 447
 448 StgInt
 449 isFloatDenormalized(f)
 450 StgFloat f;
 451 {
 452 # ifdef FLOATS_AS_DOUBLES
 453     return (isDoubleDenormalized(f));
 454 # else
 455     union stg_ieee754_flt u;
 456     u.f = f;
 457
 458     /* A (single/double/quad) precision floating point number
 459        is denormalised iff:
 460         - exponent is zero
 461         - mantissa is non-zero.
 462         - (don't care about setting of sign bit.)
 463
 464     */
 465     return (
 466         u.ieee.exponent == 0 &&
 467         u.ieee.mantissa != 0);
 468 #endif /* !FLOATS_AS_DOUBLES */
 469 }
 470
 471 StgInt
 472 isFloatNegativeZero(f)
 473 StgFloat f;
 474 {
 475 #ifdef FLOATS_AS_DOUBLES
 476     return (isDoubleNegativeZero(f));
 477 # else
 478     union stg_ieee754_flt u;
 479     u.f = f;
 480
 481     /* sign (bit 31) set (only) => negative zero */
 482     return (
 483         u.ieee.negative      &&
 484         u.ieee.exponent == 0 &&
 485         u.ieee.mantissa == 0);
 486 # endif /* !FLOATS_AS_DOUBLES */
 487 }
 488
 489 #else /* ! IEEE_FLOATING_POINT */
 490
 491 /* Dummy definitions of predicates - they all return false */
 492 StgInt isDoubleNaN(d) StgDouble d; { return 0; }
 493 StgInt isDoubleInfinite(d) StgDouble d; { return 0; }
 494 StgInt isDoubleDenormalized(d) StgDouble d; { return 0; }
 495 StgInt isDoubleNegativeZero(d) StgDouble d; { return 0; }
 496 StgInt isFloatNaN(f) StgFloat f; { return 0; }
 497 StgInt isFloatInfinite(f) StgFloat f; { return 0; }
 498 StgInt isFloatDenormalized(f) StgFloat f; { return 0; }
 499 StgInt isFloatNegativeZero(f) StgFloat f; { return 0; }
 500
 501 #endif /* ! IEEE_FLOATING_POINT */