cbits/primFloat.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 "HsFFI.h"
  11 #include "Rts.h" // XXX wrong (for IEEE_FLOATING_POINT and WORDS_BIGENDIAN)
  12
  13 #define IEEE_FLOATING_POINT 1
  14
  15 union stg_ieee754_flt
  16 {
  17    float f;
  18    struct {
  19
  20 #if WORDS_BIGENDIAN
  21         unsigned int negative:1;
  22         unsigned int exponent:8;
  23         unsigned int mantissa:23;
  24 #else
  25         unsigned int mantissa:23;
  26         unsigned int exponent:8;
  27         unsigned int negative:1;
  28 #endif
  29    } ieee;
  30    struct {
  31
  32 #if WORDS_BIGENDIAN
  33         unsigned int negative:1;
  34         unsigned int exponent:8;
  35         unsigned int quiet_nan:1;
  36         unsigned int mantissa:22;
  37 #else
  38         unsigned int mantissa:22;
  39         unsigned int quiet_nan:1;
  40         unsigned int exponent:8;
  41         unsigned int negative:1;
  42 #endif
  43    } ieee_nan;
  44 };
  45
  46 /*
  47
  48  To recap, here's the representation of a double precision
  49  IEEE floating point number:
  50
  51  sign         63           sign bit (0==positive, 1==negative)
  52  exponent     62-52        exponent (biased by 1023)
  53  fraction     51-0         fraction (bits to right of binary point)
  54 */
  55
  56 union stg_ieee754_dbl
  57 {
  58    double d;
  59    struct {
  60
  61 #if WORDS_BIGENDIAN
  62         unsigned int negative:1;
  63         unsigned int exponent:11;
  64         unsigned int mantissa0:20;
  65         unsigned int mantissa1:32;
  66 #else
  67 #if FLOAT_WORDS_BIGENDIAN
  68         unsigned int mantissa0:20;
  69         unsigned int exponent:11;
  70         unsigned int negative:1;
  71         unsigned int mantissa1:32;
  72 #else
  73         unsigned int mantissa1:32;
  74         unsigned int mantissa0:20;
  75         unsigned int exponent:11;
  76         unsigned int negative:1;
  77 #endif
  78 #endif
  79    } ieee;
  80     /* This format makes it easier to see if a NaN is a signalling NaN.  */
  81    struct {
  82
  83 #if WORDS_BIGENDIAN
  84         unsigned int negative:1;
  85         unsigned int exponent:11;
  86         unsigned int quiet_nan:1;
  87         unsigned int mantissa0:19;
  88         unsigned int mantissa1:32;
  89 #else
  90 #if FLOAT_WORDS_BIGENDIAN
  91         unsigned int mantissa0:19;
  92         unsigned int quiet_nan:1;
  93         unsigned int exponent:11;
  94         unsigned int negative:1;
  95         unsigned int mantissa1:32;
  96 #else
  97         unsigned int mantissa1:32;
  98         unsigned int mantissa0:19;
  99         unsigned int quiet_nan:1;
 100         unsigned int exponent:11;
 101         unsigned int negative:1;
 102 #endif
 103 #endif
 104    } ieee_nan;
 105 };
 106
 107 /*
 108  * Predicates for testing for extreme IEEE fp values.
 109  */
 110
 111 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
 112 #ifdef IEEE_FLOATING_POINT
 113
 114 HsInt
 115 isDoubleNaN(HsDouble d)
 116 {
 117   union stg_ieee754_dbl u;
 118
 119   u.d = d;
 120
 121   return (
 122     u.ieee.exponent  == 2047 /* 2^11 - 1 */ &&  /* Is the exponent all ones? */
 123     (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)
 124         /* and the mantissa non-zero? */
 125     );
 126 }
 127
 128 HsInt
 129 isDoubleInfinite(HsDouble d)
 130 {
 131     union stg_ieee754_dbl u;
 132
 133     u.d = d;
 134
 135     /* Inf iff exponent is all ones, mantissa all zeros */
 136     return (
 137         u.ieee.exponent  == 2047 /* 2^11 - 1 */ &&
 138         u.ieee.mantissa0 == 0                   &&
 139         u.ieee.mantissa1 == 0
 140       );
 141 }
 142
 143 HsInt
 144 isDoubleDenormalized(HsDouble d)
 145 {
 146     union stg_ieee754_dbl u;
 147
 148     u.d = d;
 149
 150     /* A (single/double/quad) precision floating point number
 151        is denormalised iff:
 152         - exponent is zero
 153         - mantissa is non-zero.
 154         - (don't care about setting of sign bit.)
 155
 156     */
 157     return (
 158         u.ieee.exponent  == 0 &&
 159         (u.ieee.mantissa0 != 0 ||
 160          u.ieee.mantissa1 != 0)
 161       );
 162
 163 }
 164
 165 HsInt
 166 isDoubleNegativeZero(HsDouble d)
 167 {
 168     union stg_ieee754_dbl u;
 169
 170     u.d = d;
 171     /* sign (bit 63) set (only) => negative zero */
 172
 173     return (
 174         u.ieee.negative  == 1 &&
 175         u.ieee.exponent  == 0 &&
 176         u.ieee.mantissa0 == 0 &&
 177         u.ieee.mantissa1 == 0);
 178 }
 179
 180 /* Same tests, this time for HsFloats. */
 181
 182 /*
 183  To recap, here's the representation of a single precision
 184  IEEE floating point number:
 185
 186  sign         31           sign bit (0 == positive, 1 == negative)
 187  exponent     30-23        exponent (biased by 127)
 188  fraction     22-0         fraction (bits to right of binary point)
 189 */
 190
 191
 192 HsInt
 193 isFloatNaN(HsFloat f)
 194 {
 195     union stg_ieee754_flt u;
 196     u.f = f;
 197
 198    /* Floating point NaN iff exponent is all ones, mantissa is
 199       non-zero (but see below.) */
 200    return (
 201         u.ieee.exponent == 255 /* 2^8 - 1 */ &&
 202         u.ieee.mantissa != 0);
 203 }
 204
 205 HsInt
 206 isFloatInfinite(HsFloat f)
 207 {
 208     union stg_ieee754_flt u;
 209     u.f = f;
 210
 211     /* A float is Inf iff exponent is max (all ones),
 212        and mantissa is min(all zeros.) */
 213     return (
 214         u.ieee.exponent == 255 /* 2^8 - 1 */ &&
 215         u.ieee.mantissa == 0);
 216 }
 217
 218 HsInt
 219 isFloatDenormalized(HsFloat f)
 220 {
 221     union stg_ieee754_flt u;
 222     u.f = f;
 223
 224     /* A (single/double/quad) precision floating point number
 225        is denormalised iff:
 226         - exponent is zero
 227         - mantissa is non-zero.
 228         - (don't care about setting of sign bit.)
 229
 230     */
 231     return (
 232         u.ieee.exponent == 0 &&
 233         u.ieee.mantissa != 0);
 234 }
 235
 236 HsInt
 237 isFloatNegativeZero(HsFloat f)
 238 {
 239     union stg_ieee754_flt u;
 240     u.f = f;
 241
 242     /* sign (bit 31) set (only) => negative zero */
 243     return (
 244         u.ieee.negative      &&
 245         u.ieee.exponent == 0 &&
 246         u.ieee.mantissa == 0);
 247 }
 248
 249 /*
 250  There are glibc versions around with buggy rintf or rint, hence we
 251  provide our own. We always round ties to even, so we can be simpler.
 252 */
 253
 254 #define FLT_HIDDEN 0x800000
 255 #define FLT_POWER2 0x1000000
 256
 257 HsFloat
 258 rintFloat(HsFloat f)
 259 {
 260     union stg_ieee754_flt u;
 261     u.f = f;
 262     /* if real exponent > 22, it's already integral, infinite or nan */
 263     if (u.ieee.exponent > 149)  /* 22 + 127 */
 264     {
 265         return u.f;
 266     }
 267     if (u.ieee.exponent < 126)  /* (-1) + 127, abs(f) < 0.5 */
 268     {
 269         /* only used for rounding to Integral a, so don't care about -0.0 */
 270         return 0.0;
 271     }
 272     /* 0.5 <= abs(f) < 2^23 */
 273     unsigned int half, mask, mant, frac;
 274     half = 1 << (149 - u.ieee.exponent);    /* bit for 0.5 */
 275     mask = 2*half - 1;                      /* fraction bits */
 276     mant = u.ieee.mantissa | FLT_HIDDEN;    /* add hidden bit */
 277     frac = mant & mask;                     /* get fraction */
 278     mant ^= frac;                           /* truncate mantissa */
 279     if ((frac < half) || ((frac == half) && ((mant & (2*half)) == 0)))
 280     {
 281         /* this means we have to truncate */
 282         if (mant == 0)
 283         {
 284             /* f == ±0.5, return 0.0 */
 285             return 0.0;
 286         }
 287         else
 288         {
 289             /* remove hidden bit and set mantissa */
 290             u.ieee.mantissa = mant ^ FLT_HIDDEN;
 291             return u.f;
 292         }
 293     }
 294     else
 295     {
 296         /* round away from zero, increment mantissa */
 297         mant += 2*half;
 298         if (mant == FLT_POWER2)
 299         {
 300             /* next power of 2, increase exponent an set mantissa to 0 */
 301             u.ieee.mantissa = 0;
 302             u.ieee.exponent += 1;
 303             return u.f;
 304         }
 305         else
 306         {
 307             /* remove hidden bit and set mantissa */
 308             u.ieee.mantissa = mant ^ FLT_HIDDEN;
 309             return u.f;
 310         }
 311     }
 312 }
 313
 314 #define DBL_HIDDEN 0x100000
 315 #define DBL_POWER2 0x200000
 316 #define LTOP_BIT 0x80000000
 317
 318 HsDouble
 319 rintDouble(HsDouble d)
 320 {
 321     union stg_ieee754_dbl u;
 322     u.d = d;
 323     /* if real exponent > 51, it's already integral, infinite or nan */
 324     if (u.ieee.exponent > 1074) /* 51 + 1023 */
 325     {
 326         return u.d;
 327     }
 328     if (u.ieee.exponent < 1022)  /* (-1) + 1023, abs(d) < 0.5 */
 329     {
 330         /* only used for rounding to Integral a, so don't care about -0.0 */
 331         return 0.0;
 332     }
 333     unsigned int half, mask, mant, frac;
 334     if (u.ieee.exponent < 1043) /* 20 + 1023, real exponent < 20 */
 335     {
 336         /* the fractional part meets the higher part of the mantissa */
 337         half = 1 << (1042 - u.ieee.exponent);   /* bit for 0.5 */
 338         mask = 2*half - 1;                      /* fraction bits */
 339         mant = u.ieee.mantissa0 | DBL_HIDDEN;   /* add hidden bit */
 340         frac = mant & mask;                     /* get fraction */
 341         mant ^= frac;                           /* truncate mantissa */
 342         if ((frac < half) ||
 343             ((frac == half) && (u.ieee.mantissa1 == 0)  /* a tie */
 344                 && ((mant & (2*half)) == 0)))
 345         {
 346             /* truncate */
 347             if (mant == 0)
 348             {
 349                 /* d = ±0.5, return 0.0 */
 350                 return 0.0;
 351             }
 352             /* remove hidden bit and set mantissa */
 353             u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
 354             u.ieee.mantissa1 = 0;
 355             return u.d;
 356         }
 357         else    /* round away from zero */
 358         {
 359             /* zero low mantissa bits */
 360             u.ieee.mantissa1 = 0;
 361             /* increment integer part of mantissa */
 362             mant += 2*half;
 363             if (mant == DBL_POWER2)
 364             {
 365                 /* power of 2, increment exponent and zero mantissa */
 366                 u.ieee.mantissa0 = 0;
 367                 u.ieee.exponent += 1;
 368                 return u.d;
 369             }
 370             /* remove hidden bit */
 371             u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
 372             return u.d;
 373         }
 374     }
 375     else
 376     {
 377         /* 20 <= real exponent < 52, fractional part entirely in mantissa1 */
 378         half = 1 << (1074 - u.ieee.exponent);   /* bit for 0.5 */
 379         mask = 2*half - 1;                      /* fraction bits */
 380         mant = u.ieee.mantissa1;                /* no hidden bit here */
 381         frac = mant & mask;                     /* get fraction */
 382         mant ^= frac;                           /* truncate mantissa */
 383         if ((frac < half) ||
 384             ((frac == half) &&                  /* tie */
 385             (((half == LTOP_BIT) ? (u.ieee.mantissa0 & 1)  /* yuck */
 386                                 : (mant & (2*half)))
 387                                         == 0)))
 388         {
 389             /* truncate */
 390             u.ieee.mantissa1 = mant;
 391             return u.d;
 392         }
 393         else
 394         {
 395             /* round away from zero */
 396             /* increment mantissa */
 397             mant += 2*half;
 398             u.ieee.mantissa1 = mant;
 399             if (mant == 0)
 400             {
 401                 /* low part of mantissa overflowed */
 402                 /* increment high part of mantissa */
 403                 mant = u.ieee.mantissa0 + 1;
 404                 if (mant == DBL_HIDDEN)
 405                 {
 406                     /* hit power of 2 */
 407                     /* zero mantissa */
 408                     u.ieee.mantissa0 = 0;
 409                     /* and increment exponent */
 410                     u.ieee.exponent += 1;
 411                     return u.d;
 412                 }
 413                 else
 414                 {
 415                     u.ieee.mantissa0 = mant;
 416                     return u.d;
 417                 }
 418             }
 419             else
 420             {
 421                 return u.d;
 422             }
 423         }
 424     }
 425 }
 426
 427 #else /* ! IEEE_FLOATING_POINT */
 428
 429 /* Dummy definitions of predicates - they all return false */
 430 HsInt isDoubleNaN(d) HsDouble d; { return 0; }
 431 HsInt isDoubleInfinite(d) HsDouble d; { return 0; }
 432 HsInt isDoubleDenormalized(d) HsDouble d; { return 0; }
 433 HsInt isDoubleNegativeZero(d) HsDouble d; { return 0; }
 434 HsInt isFloatNaN(f) HsFloat f; { return 0; }
 435 HsInt isFloatInfinite(f) HsFloat f; { return 0; }
 436 HsInt isFloatDenormalized(f) HsFloat f; { return 0; }
 437 HsInt isFloatNegativeZero(f) HsFloat f; { return 0; }
 438
 439
 440 /* For exotic floating point formats, we can't do much */
 441 /* We suppose the format has not too many bits */
 442 /* I hope nobody tries to build GHC where this is wrong */
 443
 444 #define FLT_UPP 536870912.0
 445
 446 HsFloat
 447 rintFloat(HsFloat f)
 448 {
 449     if ((f > FLT_UPP) || (f < (-FLT_UPP)))
 450     {
 451         return f;
 452     }
 453     else
 454     {
 455         int i = (int)f;
 456         float g = i;
 457         float d = f - g;
 458         if (d > 0.5)
 459         {
 460             return g + 1.0;
 461         }
 462         if (d == 0.5)
 463         {
 464             return (i & 1) ? (g + 1.0) : g;
 465         }
 466         if (d == -0.5)
 467         {
 468             return (i & 1) ? (g - 1.0) : g;
 469         }
 470         if (d < -0.5)
 471         {
 472             return g - 1.0;
 473         }
 474         return g;
 475     }
 476 }
 477
 478 #define DBL_UPP 2305843009213693952.0
 479
 480 HsDouble
 481 rintDouble(HsDouble d)
 482 {
 483     if ((d > DBL_UPP) || (d < (-DBL_UPP)))
 484     {
 485         return d;
 486     }
 487     else
 488     {
 489         HsInt64 i = (HsInt64)d;
 490         double e = i;
 491         double r = d - e;
 492         if (r > 0.5)
 493         {
 494             return e + 1.0;
 495         }
 496         if (r == 0.5)
 497         {
 498             return (i & 1) ? (e + 1.0) : e;
 499         }
 500         if (r == -0.5)
 501         {
 502             return (i & 1) ? (e - 1.0) : e;
 503         }
 504         if (r < -0.5)
 505         {
 506             return e - 1.0;
 507         }
 508         return e;
 509     }
 510 }
 511
 512 #endif /* ! IEEE_FLOATING_POINT */