1 /* -----------------------------------------------------------------------------
2 * $Id: StgPrimFloat.c,v 1.6 2000/11/07 13:30:41 simonmar Exp $
4 * (c) The GHC Team, 1998-2000
6 * Miscellaneous support for floating-point primitives
8 * ---------------------------------------------------------------------------*/
13 * Encoding and decoding Doubles. Code based on the HBC code
17 #define GMP_BASE 4294967296.0
18 #define DNBIGIT 2 /* mantissa of a double will fit in two longs */
19 #define FNBIGIT 1 /* for float, one long */
21 #if IEEE_FLOATING_POINT
22 #define MY_DMINEXP ((DBL_MIN_EXP) - (DBL_MANT_DIG) - 1)
23 /* DMINEXP is defined in values.h on Linux (for example) */
24 #define DHIGHBIT 0x00100000
25 #define DMSBIT 0x80000000
27 #define MY_FMINEXP ((FLT_MIN_EXP) - (FLT_MANT_DIG) - 1)
28 #define FHIGHBIT 0x00800000
29 #define FMSBIT 0x80000000
32 #ifdef WORDS_BIGENDIAN
40 #define __abs(a) (( (a) >= 0 ) ? (a) : (-(a)))
43 __encodeDouble (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */
49 /* Convert MP_INT to a double; knows a lot about internal rep! */
50 for(r = 0.0, i = __abs(size)-1; i >= 0; i--)
51 r = (r * GMP_BASE) + arr[i];
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 small Integers */
66 __int_encodeDouble (I_ j, I_ e)
70 r = (StgDouble)__abs(j);
72 /* Now raise to the exponent */
73 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
76 /* sign is encoded in the size */
84 __encodeFloat (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */
90 /* Convert MP_INT to a float; knows a lot about internal rep! */
91 for(r = 0.0, i = __abs(size)-1; i >= 0; i--)
92 r = (r * GMP_BASE) + arr[i];
94 /* Now raise to the exponent */
95 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
98 /* sign is encoded in the size */
105 /* Special version for small Integers */
107 __int_encodeFloat (I_ j, I_ e)
111 r = (StgFloat)__abs(j);
113 /* Now raise to the exponent */
114 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
117 /* sign is encoded in the size */
124 /* This only supports IEEE floating point */
127 __decodeDouble (MP_INT *man, I_ *exp, StgDouble dbl)
129 /* Do some bit fiddling on IEEE */
130 nat low, high; /* assuming 32 bit ints */
132 union { double d; int i[2]; } u; /* assuming 32 bit ints, 64 bit double */
134 u.d = dbl; /* grab chunks of the double */
138 /* we know the MP_INT* passed in has size zero, so we realloc
141 man->_mp_alloc = DNBIGIT;
143 if (low == 0 && (high & ~DMSBIT) == 0) {
147 man->_mp_size = DNBIGIT;
148 iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
152 if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
156 /* A denorm, normalize the mantissa */
157 while (! (high & DHIGHBIT)) {
168 man->_mp_d[1] = high;
171 man->_mp_d[0] = ((unsigned long)high) << 32 | (unsigned long)low;
173 error : error : error : Cannae cope with DNBIGIT
177 man->_mp_size = -man->_mp_size;
182 __decodeFloat (MP_INT *man, I_ *exp, StgFloat flt)
184 /* Do some bit fiddling on IEEE */
185 int high, sign; /* assuming 32 bit ints */
186 union { float f; int i; } u; /* assuming 32 bit float and int */
188 u.f = flt; /* grab the float */
191 /* we know the MP_INT* passed in has size zero, so we realloc
194 man->_mp_alloc = FNBIGIT;
196 if ((high & ~FMSBIT) == 0) {
200 man->_mp_size = FNBIGIT;
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)) {
216 man->_mp_d[0] = high;
218 error : error : error : Cannae cope with FNBIGIT
221 man->_mp_size = -man->_mp_size;
225 /* Convenient union types for checking the layout of IEEE 754 types -
226 based on defs in GNU libc <ieee754.h>
229 union stg_ieee754_flt
235 unsigned int negative:1;
236 unsigned int exponent:8;
237 unsigned int mantissa:23;
239 unsigned int mantissa:23;
240 unsigned int exponent:8;
241 unsigned int negative:1;
247 unsigned int negative:1;
248 unsigned int exponent:8;
249 unsigned int quiet_nan:1;
250 unsigned int mantissa:22;
252 unsigned int mantissa:22;
253 unsigned int quiet_nan:1;
254 unsigned int exponent:8;
255 unsigned int negative:1;
262 To recap, here's the representation of a double precision
263 IEEE floating point number:
265 sign 63 sign bit (0==positive, 1==negative)
266 exponent 62-52 exponent (biased by 1023)
267 fraction 51-0 fraction (bits to right of binary point)
270 union stg_ieee754_dbl
276 unsigned int negative:1;
277 unsigned int exponent:11;
278 unsigned int mantissa0:20;
279 unsigned int mantissa1:32;
281 unsigned int mantissa1:32;
282 unsigned int mantissa0:20;
283 unsigned int exponent:11;
284 unsigned int negative:1;
287 /* This format makes it easier to see if a NaN is a signalling NaN. */
291 unsigned int negative:1;
292 unsigned int exponent:11;
293 unsigned int quiet_nan:1;
294 unsigned int mantissa0:19;
295 unsigned int mantissa1:32;
297 unsigned int mantissa1:32;
298 unsigned int mantissa0:19;
299 unsigned int quiet_nan:1;
300 unsigned int exponent:11;
301 unsigned int negative:1;
307 * Predicates for testing for extreme IEEE fp values. Used
308 * by the bytecode evaluator and the Prelude.
312 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
313 #ifdef IEEE_FLOATING_POINT
316 isDoubleNaN(StgDouble d)
318 union stg_ieee754_dbl u;
323 u.ieee.exponent == 2047 /* 2^11 - 1 */ && /* Is the exponent all ones? */
324 (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)
325 /* and the mantissa non-zero? */
330 isDoubleInfinite(StgDouble d)
332 union stg_ieee754_dbl u;
336 /* Inf iff exponent is all ones, mantissa all zeros */
338 u.ieee.exponent == 2047 /* 2^11 - 1 */ &&
339 u.ieee.mantissa0 == 0 &&
340 u.ieee.mantissa1 == 0
345 isDoubleDenormalized(StgDouble d)
347 union stg_ieee754_dbl u;
351 /* A (single/double/quad) precision floating point number
354 - mantissa is non-zero.
355 - (don't care about setting of sign bit.)
359 u.ieee.exponent == 0 &&
360 (u.ieee.mantissa0 != 0 ||
361 u.ieee.mantissa1 != 0)
367 isDoubleNegativeZero(StgDouble d)
369 union stg_ieee754_dbl u;
372 /* sign (bit 63) set (only) => negative zero */
375 u.ieee.negative == 1 &&
376 u.ieee.exponent == 0 &&
377 u.ieee.mantissa0 == 0 &&
378 u.ieee.mantissa1 == 0);
381 /* Same tests, this time for StgFloats. */
384 To recap, here's the representation of a single precision
385 IEEE floating point number:
387 sign 31 sign bit (0 == positive, 1 == negative)
388 exponent 30-23 exponent (biased by 127)
389 fraction 22-0 fraction (bits to right of binary point)
394 isFloatNaN(StgFloat f)
396 union stg_ieee754_flt u;
399 /* Floating point NaN iff exponent is all ones, mantissa is
400 non-zero (but see below.) */
402 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
403 u.ieee.mantissa != 0);
407 isFloatInfinite(StgFloat f)
409 union stg_ieee754_flt u;
412 /* A float is Inf iff exponent is max (all ones),
413 and mantissa is min(all zeros.) */
415 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
416 u.ieee.mantissa == 0);
420 isFloatDenormalized(StgFloat f)
422 union stg_ieee754_flt u;
425 /* A (single/double/quad) precision floating point number
428 - mantissa is non-zero.
429 - (don't care about setting of sign bit.)
433 u.ieee.exponent == 0 &&
434 u.ieee.mantissa != 0);
438 isFloatNegativeZero(StgFloat f)
440 union stg_ieee754_flt u;
443 /* sign (bit 31) set (only) => negative zero */
446 u.ieee.exponent == 0 &&
447 u.ieee.mantissa == 0);
450 #else /* ! IEEE_FLOATING_POINT */
452 /* Dummy definitions of predicates - they all return false */
453 StgInt isDoubleNaN(d) StgDouble d; { return 0; }
454 StgInt isDoubleInfinite(d) StgDouble d; { return 0; }
455 StgInt isDoubleDenormalized(d) StgDouble d; { return 0; }
456 StgInt isDoubleNegativeZero(d) StgDouble d; { return 0; }
457 StgInt isFloatNaN(f) StgFloat f; { return 0; }
458 StgInt isFloatInfinite(f) StgFloat f; { return 0; }
459 StgInt isFloatDenormalized(f) StgFloat f; { return 0; }
460 StgInt isFloatNegativeZero(f) StgFloat f; { return 0; }
462 #endif /* ! IEEE_FLOATING_POINT */