1 /* -----------------------------------------------------------------------------
3 * (c) Lennart Augustsson
4 * (c) The GHC Team, 1998-2000
6 * Miscellaneous support for floating-point primitives
8 * ---------------------------------------------------------------------------*/
11 #include "Rts.h" // XXX wrong (for IEEE_FLOATING_POINT and WORDS_BIGENDIAN)
13 #define IEEE_FLOATING_POINT 1
21 unsigned int negative:1;
22 unsigned int exponent:8;
23 unsigned int mantissa:23;
25 unsigned int mantissa:23;
26 unsigned int exponent:8;
27 unsigned int negative:1;
33 unsigned int negative:1;
34 unsigned int exponent:8;
35 unsigned int quiet_nan:1;
36 unsigned int mantissa:22;
38 unsigned int mantissa:22;
39 unsigned int quiet_nan:1;
40 unsigned int exponent:8;
41 unsigned int negative:1;
48 To recap, here's the representation of a double precision
49 IEEE floating point number:
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)
62 unsigned int negative:1;
63 unsigned int exponent:11;
64 unsigned int mantissa0:20;
65 unsigned int mantissa1:32;
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;
73 unsigned int mantissa1:32;
74 unsigned int mantissa0:20;
75 unsigned int exponent:11;
76 unsigned int negative:1;
80 /* This format makes it easier to see if a NaN is a signalling NaN. */
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;
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;
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;
108 * Predicates for testing for extreme IEEE fp values.
111 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
112 #ifdef IEEE_FLOATING_POINT
115 isDoubleNaN(HsDouble d)
117 union stg_ieee754_dbl u;
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? */
129 isDoubleInfinite(HsDouble d)
131 union stg_ieee754_dbl u;
135 /* Inf iff exponent is all ones, mantissa all zeros */
137 u.ieee.exponent == 2047 /* 2^11 - 1 */ &&
138 u.ieee.mantissa0 == 0 &&
139 u.ieee.mantissa1 == 0
144 isDoubleDenormalized(HsDouble d)
146 union stg_ieee754_dbl u;
150 /* A (single/double/quad) precision floating point number
153 - mantissa is non-zero.
154 - (don't care about setting of sign bit.)
158 u.ieee.exponent == 0 &&
159 (u.ieee.mantissa0 != 0 ||
160 u.ieee.mantissa1 != 0)
166 isDoubleNegativeZero(HsDouble d)
168 union stg_ieee754_dbl u;
171 /* sign (bit 63) set (only) => negative zero */
174 u.ieee.negative == 1 &&
175 u.ieee.exponent == 0 &&
176 u.ieee.mantissa0 == 0 &&
177 u.ieee.mantissa1 == 0);
180 /* Same tests, this time for HsFloats. */
183 To recap, here's the representation of a single precision
184 IEEE floating point number:
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)
193 isFloatNaN(HsFloat f)
195 union stg_ieee754_flt u;
198 /* Floating point NaN iff exponent is all ones, mantissa is
199 non-zero (but see below.) */
201 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
202 u.ieee.mantissa != 0);
206 isFloatInfinite(HsFloat f)
208 union stg_ieee754_flt u;
211 /* A float is Inf iff exponent is max (all ones),
212 and mantissa is min(all zeros.) */
214 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
215 u.ieee.mantissa == 0);
219 isFloatDenormalized(HsFloat f)
221 union stg_ieee754_flt u;
224 /* A (single/double/quad) precision floating point number
227 - mantissa is non-zero.
228 - (don't care about setting of sign bit.)
232 u.ieee.exponent == 0 &&
233 u.ieee.mantissa != 0);
237 isFloatNegativeZero(HsFloat f)
239 union stg_ieee754_flt u;
242 /* sign (bit 31) set (only) => negative zero */
245 u.ieee.exponent == 0 &&
246 u.ieee.mantissa == 0);
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.
254 #define FLT_HIDDEN 0x800000
255 #define FLT_POWER2 0x1000000
260 union stg_ieee754_flt u;
262 /* if real exponent > 22, it's already integral, infinite or nan */
263 if (u.ieee.exponent > 149) /* 22 + 127 */
267 if (u.ieee.exponent < 126) /* (-1) + 127, abs(f) < 0.5 */
269 /* only used for rounding to Integral a, so don't care about -0.0 */
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)))
281 /* this means we have to truncate */
284 /* f == ±0.5, return 0.0 */
289 /* remove hidden bit and set mantissa */
290 u.ieee.mantissa = mant ^ FLT_HIDDEN;
296 /* round away from zero, increment mantissa */
298 if (mant == FLT_POWER2)
300 /* next power of 2, increase exponent an set mantissa to 0 */
302 u.ieee.exponent += 1;
307 /* remove hidden bit and set mantissa */
308 u.ieee.mantissa = mant ^ FLT_HIDDEN;
314 #define DBL_HIDDEN 0x100000
315 #define DBL_POWER2 0x200000
316 #define LTOP_BIT 0x80000000
319 rintDouble(HsDouble d)
321 union stg_ieee754_dbl u;
323 /* if real exponent > 51, it's already integral, infinite or nan */
324 if (u.ieee.exponent > 1074) /* 51 + 1023 */
328 if (u.ieee.exponent < 1022) /* (-1) + 1023, abs(d) < 0.5 */
330 /* only used for rounding to Integral a, so don't care about -0.0 */
333 unsigned int half, mask, mant, frac;
334 if (u.ieee.exponent < 1043) /* 20 + 1023, real exponent < 20 */
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 */
343 ((frac == half) && (u.ieee.mantissa1 == 0) /* a tie */
344 && ((mant & (2*half)) == 0)))
349 /* d = ±0.5, return 0.0 */
352 /* remove hidden bit and set mantissa */
353 u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
354 u.ieee.mantissa1 = 0;
357 else /* round away from zero */
359 /* zero low mantissa bits */
360 u.ieee.mantissa1 = 0;
361 /* increment integer part of mantissa */
363 if (mant == DBL_POWER2)
365 /* power of 2, increment exponent and zero mantissa */
366 u.ieee.mantissa0 = 0;
367 u.ieee.exponent += 1;
370 /* remove hidden bit */
371 u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
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 */
384 ((frac == half) && /* tie */
385 (((half == LTOP_BIT) ? (u.ieee.mantissa0 & 1) /* yuck */
390 u.ieee.mantissa1 = mant;
395 /* round away from zero */
396 /* increment mantissa */
398 u.ieee.mantissa1 = mant;
401 /* low part of mantissa overflowed */
402 /* increment high part of mantissa */
403 mant = u.ieee.mantissa0 + 1;
404 if (mant == DBL_HIDDEN)
408 u.ieee.mantissa0 = 0;
409 /* and increment exponent */
410 u.ieee.exponent += 1;
415 u.ieee.mantissa0 = mant;
427 #else /* ! IEEE_FLOATING_POINT */
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; }
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 */
444 #define FLT_UPP 536870912.0
449 if ((f > FLT_UPP) || (f < (-FLT_UPP)))
464 return (i & 1) ? (g + 1.0) : g;
468 return (i & 1) ? (g - 1.0) : g;
478 #define DBL_UPP 2305843009213693952.0
481 rintDouble(HsDouble d)
483 if ((d > DBL_UPP) || (d < (-DBL_UPP)))
489 HsInt64 i = (HsInt64)d;
498 return (i & 1) ? (e + 1.0) : e;
502 return (i & 1) ? (e - 1.0) : e;
512 #endif /* ! IEEE_FLOATING_POINT */