1 /* -----------------------------------------------------------------------------
3 * (c) Lennart Augustsson
4 * (c) The GHC Team, 1998-2000
6 * Miscellaneous support for floating-point primitives
8 * ---------------------------------------------------------------------------*/
10 #include "PosixSource.h"
16 * Encoding and decoding Doubles. Code based on the HBC code
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
26 #define MY_FMINEXP ((FLT_MIN_EXP) - (FLT_MANT_DIG) - 1)
27 #define FHIGHBIT 0x00800000
28 #define FMSBIT 0x80000000
31 #if defined(WORDS_BIGENDIAN) || defined(FLOAT_WORDS_BIGENDIAN)
39 #define __abs(a) (( (a) >= 0 ) ? (a) : (-(a)))
42 __2Int_encodeDouble (I_ j_high, I_ j_low, I_ e)
46 /* assuming 32 bit ints */
47 ASSERT(sizeof(int ) == 4 );
49 r = (StgDouble)((unsigned int)j_high);
50 r *= 4294967296.0; /* exp2f(32); */
51 r += (StgDouble)((unsigned int)j_low);
53 /* Now raise to the exponent */
54 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
57 /* sign is encoded in the size */
64 /* Special version for words */
66 __word_encodeDouble (W_ j, I_ e)
72 /* Now raise to the exponent */
73 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
79 /* Special version for small Integers */
81 __int_encodeDouble (I_ j, I_ e)
85 r = (StgDouble)__abs(j);
87 /* Now raise to the exponent */
88 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
91 /* sign is encoded in the size */
98 /* Special version for small Integers */
100 __int_encodeFloat (I_ j, I_ e)
104 r = (StgFloat)__abs(j);
106 /* Now raise to the exponent */
107 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
110 /* sign is encoded in the size */
117 /* Special version for small positive Integers */
119 __word_encodeFloat (W_ j, I_ e)
125 /* Now raise to the exponent */
126 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
132 /* This only supports IEEE floating point */
135 __decodeDouble_2Int (I_ *man_sign, W_ *man_high, W_ *man_low, I_ *exp, StgDouble dbl)
137 /* Do some bit fiddling on IEEE */
138 unsigned int low, high; /* assuming 32 bit ints */
140 union { double d; unsigned int i[2]; } u; /* assuming 32 bit ints, 64 bit double */
142 ASSERT(sizeof(unsigned int ) == 4 );
143 ASSERT(sizeof(dbl ) == 8 );
144 ASSERT(sizeof(dbl ) == SIZEOF_DOUBLE);
146 u.d = dbl; /* grab chunks of the double */
150 if (low == 0 && (high & ~DMSBIT) == 0) {
155 iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
159 if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
163 /* A denorm, normalize the mantissa */
164 while (! (high & DHIGHBIT)) {
175 *man_sign = (sign < 0) ? -1 : 1;
179 /* Convenient union types for checking the layout of IEEE 754 types -
180 based on defs in GNU libc <ieee754.h>
184 __decodeFloat_Int (I_ *man, I_ *exp, StgFloat flt)
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 */
190 ASSERT(sizeof(int ) == 4 );
191 ASSERT(sizeof(flt ) == 4 );
192 ASSERT(sizeof(flt ) == SIZEOF_FLOAT );
194 u.f = flt; /* grab the float */
197 if ((high & ~FMSBIT) == 0) {
201 *exp = ((high >> 23) & 0xff) + MY_FMINEXP;
205 if (*exp != MY_FMINEXP) /* don't add hidden bit to denorms */
209 /* A denorm, normalize the mantissa */
210 while (! (high & FHIGHBIT)) {
221 union stg_ieee754_flt
227 unsigned int negative:1;
228 unsigned int exponent:8;
229 unsigned int mantissa:23;
231 unsigned int mantissa:23;
232 unsigned int exponent:8;
233 unsigned int negative:1;
239 unsigned int negative:1;
240 unsigned int exponent:8;
241 unsigned int quiet_nan:1;
242 unsigned int mantissa:22;
244 unsigned int mantissa:22;
245 unsigned int quiet_nan:1;
246 unsigned int exponent:8;
247 unsigned int negative:1;
254 To recap, here's the representation of a double precision
255 IEEE floating point number:
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)
262 union stg_ieee754_dbl
268 unsigned int negative:1;
269 unsigned int exponent:11;
270 unsigned int mantissa0:20;
271 unsigned int mantissa1:32;
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;
279 unsigned int mantissa1:32;
280 unsigned int mantissa0:20;
281 unsigned int exponent:11;
282 unsigned int negative:1;
286 /* This format makes it easier to see if a NaN is a signalling NaN. */
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;
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;
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;
314 * Predicates for testing for extreme IEEE fp values. Used
315 * by the bytecode evaluator and the Prelude.
319 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
320 #ifdef IEEE_FLOATING_POINT
323 isDoubleNaN(StgDouble d)
325 union stg_ieee754_dbl u;
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? */
337 isDoubleInfinite(StgDouble d)
339 union stg_ieee754_dbl u;
343 /* Inf iff exponent is all ones, mantissa all zeros */
345 u.ieee.exponent == 2047 /* 2^11 - 1 */ &&
346 u.ieee.mantissa0 == 0 &&
347 u.ieee.mantissa1 == 0
352 isDoubleDenormalized(StgDouble d)
354 union stg_ieee754_dbl u;
358 /* A (single/double/quad) precision floating point number
361 - mantissa is non-zero.
362 - (don't care about setting of sign bit.)
366 u.ieee.exponent == 0 &&
367 (u.ieee.mantissa0 != 0 ||
368 u.ieee.mantissa1 != 0)
374 isDoubleNegativeZero(StgDouble d)
376 union stg_ieee754_dbl u;
379 /* sign (bit 63) set (only) => negative zero */
382 u.ieee.negative == 1 &&
383 u.ieee.exponent == 0 &&
384 u.ieee.mantissa0 == 0 &&
385 u.ieee.mantissa1 == 0);
388 /* Same tests, this time for StgFloats. */
391 To recap, here's the representation of a single precision
392 IEEE floating point number:
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)
401 isFloatNaN(StgFloat f)
403 union stg_ieee754_flt u;
406 /* Floating point NaN iff exponent is all ones, mantissa is
407 non-zero (but see below.) */
409 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
410 u.ieee.mantissa != 0);
414 isFloatInfinite(StgFloat f)
416 union stg_ieee754_flt u;
419 /* A float is Inf iff exponent is max (all ones),
420 and mantissa is min(all zeros.) */
422 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
423 u.ieee.mantissa == 0);
427 isFloatDenormalized(StgFloat f)
429 union stg_ieee754_flt u;
432 /* A (single/double/quad) precision floating point number
435 - mantissa is non-zero.
436 - (don't care about setting of sign bit.)
440 u.ieee.exponent == 0 &&
441 u.ieee.mantissa != 0);
445 isFloatNegativeZero(StgFloat f)
447 union stg_ieee754_flt u;
450 /* sign (bit 31) set (only) => negative zero */
453 u.ieee.exponent == 0 &&
454 u.ieee.mantissa == 0);
457 #else /* ! IEEE_FLOATING_POINT */
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; }
469 #endif /* ! IEEE_FLOATING_POINT */