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
21 #define SIZEOF_LIMB_T SIZEOF_UNSIGNED_INT
23 #ifdef _LONG_LONG_LIMB
24 #define SIZEOF_LIMB_T SIZEOF_UNSIGNED_LONG_LONG
26 #define SIZEOF_LIMB_T SIZEOF_UNSIGNED_LONG
30 #if SIZEOF_LIMB_T == 4
31 #define GMP_BASE 4294967296.0
32 #elif SIZEOF_LIMB_T == 8
33 #define GMP_BASE 18446744073709551616.0
35 #error Cannot cope with SIZEOF_LIMB_T -- please add definition of GMP_BASE
38 #define DNBIGIT ((SIZEOF_DOUBLE+SIZEOF_LIMB_T-1)/SIZEOF_LIMB_T)
39 #define FNBIGIT ((SIZEOF_FLOAT +SIZEOF_LIMB_T-1)/SIZEOF_LIMB_T)
41 #if IEEE_FLOATING_POINT
42 #define MY_DMINEXP ((DBL_MIN_EXP) - (DBL_MANT_DIG) - 1)
43 /* DMINEXP is defined in values.h on Linux (for example) */
44 #define DHIGHBIT 0x00100000
45 #define DMSBIT 0x80000000
47 #define MY_FMINEXP ((FLT_MIN_EXP) - (FLT_MANT_DIG) - 1)
48 #define FHIGHBIT 0x00800000
49 #define FMSBIT 0x80000000
52 #if defined(WORDS_BIGENDIAN) || defined(FLOAT_WORDS_BIGENDIAN)
60 #define __abs(a) (( (a) >= 0 ) ? (a) : (-(a)))
63 __encodeDouble (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */
66 const mp_limb_t *const arr = (const mp_limb_t *)ba;
69 /* Convert MP_INT to a double; knows a lot about internal rep! */
70 for(r = 0.0, i = __abs(size)-1; i >= 0; i--)
71 r = (r * GMP_BASE) + arr[i];
73 /* Now raise to the exponent */
74 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
77 /* sign is encoded in the size */
85 __2Int_encodeDouble (I_ j_high, I_ j_low, I_ e)
89 /* assuming 32 bit ints */
90 ASSERT(sizeof(int ) == 4 );
92 r = (StgDouble)((unsigned int)j_high);
93 r *= 4294967296.0; /* exp2f(32); */
94 r += (StgDouble)((unsigned int)j_low);
96 /* Now raise to the exponent */
97 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
100 /* sign is encoded in the size */
107 /* Special version for words */
109 __word_encodeDouble (W_ j, I_ e)
115 /* Now raise to the exponent */
116 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
122 /* Special version for small Integers */
124 __int_encodeDouble (I_ j, I_ e)
128 r = (StgDouble)__abs(j);
130 /* Now raise to the exponent */
131 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
134 /* sign is encoded in the size */
142 __encodeFloat (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */
145 const mp_limb_t *arr = (const mp_limb_t *)ba;
148 /* Convert MP_INT to a float; knows a lot about internal rep! */
149 for(r = 0.0, i = __abs(size)-1; i >= 0; i--)
150 r = (r * GMP_BASE) + arr[i];
152 /* Now raise to the exponent */
153 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
156 /* sign is encoded in the size */
163 /* Special version for small Integers */
165 __int_encodeFloat (I_ j, I_ e)
169 r = (StgFloat)__abs(j);
171 /* Now raise to the exponent */
172 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
175 /* sign is encoded in the size */
182 /* Special version for small positive Integers */
184 __word_encodeFloat (W_ j, I_ e)
190 /* Now raise to the exponent */
191 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
197 /* This only supports IEEE floating point */
200 __decodeDouble (MP_INT *man, I_ *exp, StgDouble dbl)
202 /* Do some bit fiddling on IEEE */
203 unsigned int low, high; /* assuming 32 bit ints */
205 union { double d; unsigned int i[2]; } u; /* assuming 32 bit ints, 64 bit double */
207 ASSERT(sizeof(unsigned int ) == 4 );
208 ASSERT(sizeof(dbl ) == SIZEOF_DOUBLE);
209 ASSERT(sizeof(man->_mp_d[0]) == SIZEOF_LIMB_T);
210 ASSERT(DNBIGIT*SIZEOF_LIMB_T >= SIZEOF_DOUBLE);
212 u.d = dbl; /* grab chunks of the double */
216 /* we know the MP_INT* passed in has size zero, so we realloc
219 man->_mp_alloc = DNBIGIT;
221 if (low == 0 && (high & ~DMSBIT) == 0) {
225 man->_mp_size = DNBIGIT;
226 iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
230 if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
234 /* A denorm, normalize the mantissa */
235 while (! (high & DHIGHBIT)) {
245 man->_mp_d[0] = (mp_limb_t)low;
246 man->_mp_d[1] = (mp_limb_t)high;
249 man->_mp_d[0] = ((mp_limb_t)high) << 32 | (mp_limb_t)low;
251 #error Cannot cope with DNBIGIT
255 man->_mp_size = -man->_mp_size;
260 __decodeDouble_2Int (I_ *man_sign, W_ *man_high, W_ *man_low, I_ *exp, StgDouble dbl)
262 /* Do some bit fiddling on IEEE */
263 unsigned int low, high; /* assuming 32 bit ints */
265 union { double d; unsigned int i[2]; } u; /* assuming 32 bit ints, 64 bit double */
267 ASSERT(sizeof(unsigned int ) == 4 );
268 ASSERT(sizeof(dbl ) == 8 );
269 ASSERT(sizeof(dbl ) == SIZEOF_DOUBLE);
271 u.d = dbl; /* grab chunks of the double */
275 if (low == 0 && (high & ~DMSBIT) == 0) {
280 iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
284 if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
288 /* A denorm, normalize the mantissa */
289 while (! (high & DHIGHBIT)) {
300 *man_sign = (sign < 0) ? -1 : 1;
304 /* Convenient union types for checking the layout of IEEE 754 types -
305 based on defs in GNU libc <ieee754.h>
309 __decodeFloat_Int (I_ *man, I_ *exp, StgFloat flt)
311 /* Do some bit fiddling on IEEE */
312 int high, sign; /* assuming 32 bit ints */
313 union { float f; int i; } u; /* assuming 32 bit float and int */
315 ASSERT(sizeof(int ) == 4 );
316 ASSERT(sizeof(flt ) == 4 );
317 ASSERT(sizeof(flt ) == SIZEOF_FLOAT );
319 u.f = flt; /* grab the float */
322 if ((high & ~FMSBIT) == 0) {
326 *exp = ((high >> 23) & 0xff) + MY_FMINEXP;
330 if (*exp != MY_FMINEXP) /* don't add hidden bit to denorms */
334 /* A denorm, normalize the mantissa */
335 while (! (high & FHIGHBIT)) {
346 union stg_ieee754_flt
352 unsigned int negative:1;
353 unsigned int exponent:8;
354 unsigned int mantissa:23;
356 unsigned int mantissa:23;
357 unsigned int exponent:8;
358 unsigned int negative:1;
364 unsigned int negative:1;
365 unsigned int exponent:8;
366 unsigned int quiet_nan:1;
367 unsigned int mantissa:22;
369 unsigned int mantissa:22;
370 unsigned int quiet_nan:1;
371 unsigned int exponent:8;
372 unsigned int negative:1;
379 To recap, here's the representation of a double precision
380 IEEE floating point number:
382 sign 63 sign bit (0==positive, 1==negative)
383 exponent 62-52 exponent (biased by 1023)
384 fraction 51-0 fraction (bits to right of binary point)
387 union stg_ieee754_dbl
393 unsigned int negative:1;
394 unsigned int exponent:11;
395 unsigned int mantissa0:20;
396 unsigned int mantissa1:32;
398 #if FLOAT_WORDS_BIGENDIAN
399 unsigned int mantissa0:20;
400 unsigned int exponent:11;
401 unsigned int negative:1;
402 unsigned int mantissa1:32;
404 unsigned int mantissa1:32;
405 unsigned int mantissa0:20;
406 unsigned int exponent:11;
407 unsigned int negative:1;
411 /* This format makes it easier to see if a NaN is a signalling NaN. */
415 unsigned int negative:1;
416 unsigned int exponent:11;
417 unsigned int quiet_nan:1;
418 unsigned int mantissa0:19;
419 unsigned int mantissa1:32;
421 #if FLOAT_WORDS_BIGENDIAN
422 unsigned int mantissa0:19;
423 unsigned int quiet_nan:1;
424 unsigned int exponent:11;
425 unsigned int negative:1;
426 unsigned int mantissa1:32;
428 unsigned int mantissa1:32;
429 unsigned int mantissa0:19;
430 unsigned int quiet_nan:1;
431 unsigned int exponent:11;
432 unsigned int negative:1;
439 * Predicates for testing for extreme IEEE fp values. Used
440 * by the bytecode evaluator and the Prelude.
444 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
445 #ifdef IEEE_FLOATING_POINT
448 isDoubleNaN(StgDouble d)
450 union stg_ieee754_dbl u;
455 u.ieee.exponent == 2047 /* 2^11 - 1 */ && /* Is the exponent all ones? */
456 (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)
457 /* and the mantissa non-zero? */
462 isDoubleInfinite(StgDouble d)
464 union stg_ieee754_dbl u;
468 /* Inf iff exponent is all ones, mantissa all zeros */
470 u.ieee.exponent == 2047 /* 2^11 - 1 */ &&
471 u.ieee.mantissa0 == 0 &&
472 u.ieee.mantissa1 == 0
477 isDoubleDenormalized(StgDouble d)
479 union stg_ieee754_dbl u;
483 /* A (single/double/quad) precision floating point number
486 - mantissa is non-zero.
487 - (don't care about setting of sign bit.)
491 u.ieee.exponent == 0 &&
492 (u.ieee.mantissa0 != 0 ||
493 u.ieee.mantissa1 != 0)
499 isDoubleNegativeZero(StgDouble d)
501 union stg_ieee754_dbl u;
504 /* sign (bit 63) set (only) => negative zero */
507 u.ieee.negative == 1 &&
508 u.ieee.exponent == 0 &&
509 u.ieee.mantissa0 == 0 &&
510 u.ieee.mantissa1 == 0);
513 /* Same tests, this time for StgFloats. */
516 To recap, here's the representation of a single precision
517 IEEE floating point number:
519 sign 31 sign bit (0 == positive, 1 == negative)
520 exponent 30-23 exponent (biased by 127)
521 fraction 22-0 fraction (bits to right of binary point)
526 isFloatNaN(StgFloat f)
528 union stg_ieee754_flt u;
531 /* Floating point NaN iff exponent is all ones, mantissa is
532 non-zero (but see below.) */
534 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
535 u.ieee.mantissa != 0);
539 isFloatInfinite(StgFloat f)
541 union stg_ieee754_flt u;
544 /* A float is Inf iff exponent is max (all ones),
545 and mantissa is min(all zeros.) */
547 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
548 u.ieee.mantissa == 0);
552 isFloatDenormalized(StgFloat f)
554 union stg_ieee754_flt u;
557 /* A (single/double/quad) precision floating point number
560 - mantissa is non-zero.
561 - (don't care about setting of sign bit.)
565 u.ieee.exponent == 0 &&
566 u.ieee.mantissa != 0);
570 isFloatNegativeZero(StgFloat f)
572 union stg_ieee754_flt u;
575 /* sign (bit 31) set (only) => negative zero */
578 u.ieee.exponent == 0 &&
579 u.ieee.mantissa == 0);
582 #else /* ! IEEE_FLOATING_POINT */
584 /* Dummy definitions of predicates - they all return false */
585 StgInt isDoubleNaN(d) StgDouble d; { return 0; }
586 StgInt isDoubleInfinite(d) StgDouble d; { return 0; }
587 StgInt isDoubleDenormalized(d) StgDouble d; { return 0; }
588 StgInt isDoubleNegativeZero(d) StgDouble d; { return 0; }
589 StgInt isFloatNaN(f) StgFloat f; { return 0; }
590 StgInt isFloatInfinite(f) StgFloat f; { return 0; }
591 StgInt isFloatDenormalized(f) StgFloat f; { return 0; }
592 StgInt isFloatNegativeZero(f) StgFloat f; { return 0; }
594 #endif /* ! IEEE_FLOATING_POINT */