ghc/rts/StgPrimFloat.c

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