rts/StgPrimFloat.c

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