1 /* -----------------------------------------------------------------------------
3 * (c) The GHC Team, 1998-2000
5 * Miscellaneous support for floating-point primitives
7 * ---------------------------------------------------------------------------*/
9 #include "PosixSource.h"
15 * Encoding and decoding Doubles. Code based on the HBC code
20 #define SIZEOF_LIMB_T SIZEOF_UNSIGNED_INT
22 #ifdef _LONG_LONG_LIMB
23 #define SIZEOF_LIMB_T SIZEOF_UNSIGNED_LONG_LONG
25 #define SIZEOF_LIMB_T SIZEOF_UNSIGNED_LONG
29 #if SIZEOF_LIMB_T == 4
30 #define GMP_BASE 4294967296.0
31 #elif SIZEOF_LIMB_T == 8
32 #define GMP_BASE 18446744073709551616.0
34 #error Cannot cope with SIZEOF_LIMB_T -- please add definition of GMP_BASE
37 #define DNBIGIT ((SIZEOF_DOUBLE+SIZEOF_LIMB_T-1)/SIZEOF_LIMB_T)
38 #define FNBIGIT ((SIZEOF_FLOAT +SIZEOF_LIMB_T-1)/SIZEOF_LIMB_T)
40 #if IEEE_FLOATING_POINT
41 #define MY_DMINEXP ((DBL_MIN_EXP) - (DBL_MANT_DIG) - 1)
42 /* DMINEXP is defined in values.h on Linux (for example) */
43 #define DHIGHBIT 0x00100000
44 #define DMSBIT 0x80000000
46 #define MY_FMINEXP ((FLT_MIN_EXP) - (FLT_MANT_DIG) - 1)
47 #define FHIGHBIT 0x00800000
48 #define FMSBIT 0x80000000
51 #if defined(WORDS_BIGENDIAN) || defined(FLOAT_WORDS_BIGENDIAN)
59 #define __abs(a) (( (a) >= 0 ) ? (a) : (-(a)))
62 __encodeDouble (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */
65 const mp_limb_t *const arr = (const mp_limb_t *)ba;
68 /* Convert MP_INT to a double; knows a lot about internal rep! */
69 for(r = 0.0, i = __abs(size)-1; i >= 0; i--)
70 r = (r * GMP_BASE) + arr[i];
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 __2Int_encodeDouble (I_ j_high, I_ j_low, I_ e)
88 /* assuming 32 bit ints */
89 ASSERT(sizeof(int ) == 4 );
91 r = (StgDouble)((unsigned int)j_high);
92 r *= 4294967296.0; /* exp2f(32); */
93 r += (StgDouble)((unsigned int)j_low);
95 /* Now raise to the exponent */
96 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
99 /* sign is encoded in the size */
106 /* Special version for small Integers */
108 __int_encodeDouble (I_ j, I_ e)
112 r = (StgDouble)__abs(j);
114 /* Now raise to the exponent */
115 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
118 /* sign is encoded in the size */
126 __encodeFloat (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */
129 const mp_limb_t *arr = (const mp_limb_t *)ba;
132 /* Convert MP_INT to a float; knows a lot about internal rep! */
133 for(r = 0.0, i = __abs(size)-1; i >= 0; i--)
134 r = (r * GMP_BASE) + arr[i];
136 /* Now raise to the exponent */
137 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
140 /* sign is encoded in the size */
147 /* Special version for small Integers */
149 __int_encodeFloat (I_ j, I_ e)
153 r = (StgFloat)__abs(j);
155 /* Now raise to the exponent */
156 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
159 /* sign is encoded in the size */
166 /* This only supports IEEE floating point */
169 __decodeDouble (MP_INT *man, I_ *exp, StgDouble dbl)
171 /* Do some bit fiddling on IEEE */
172 unsigned int low, high; /* assuming 32 bit ints */
174 union { double d; unsigned int i[2]; } u; /* assuming 32 bit ints, 64 bit double */
176 ASSERT(sizeof(unsigned int ) == 4 );
177 ASSERT(sizeof(dbl ) == SIZEOF_DOUBLE);
178 ASSERT(sizeof(man->_mp_d[0]) == SIZEOF_LIMB_T);
179 ASSERT(DNBIGIT*SIZEOF_LIMB_T >= SIZEOF_DOUBLE);
181 u.d = dbl; /* grab chunks of the double */
185 /* we know the MP_INT* passed in has size zero, so we realloc
188 man->_mp_alloc = DNBIGIT;
190 if (low == 0 && (high & ~DMSBIT) == 0) {
194 man->_mp_size = DNBIGIT;
195 iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
199 if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
203 /* A denorm, normalize the mantissa */
204 while (! (high & DHIGHBIT)) {
214 man->_mp_d[0] = (mp_limb_t)low;
215 man->_mp_d[1] = (mp_limb_t)high;
218 man->_mp_d[0] = ((mp_limb_t)high) << 32 | (mp_limb_t)low;
220 #error Cannot cope with DNBIGIT
224 man->_mp_size = -man->_mp_size;
229 __decodeDouble_2Int (I_ *man_high, I_ *man_low, I_ *exp, StgDouble dbl)
231 /* Do some bit fiddling on IEEE */
232 unsigned int low, high; /* assuming 32 bit ints */
234 union { double d; unsigned int i[2]; } u; /* assuming 32 bit ints, 64 bit double */
236 ASSERT(sizeof(unsigned int ) == 4 );
237 ASSERT(sizeof(dbl ) == 8 );
238 ASSERT(sizeof(dbl ) == SIZEOF_DOUBLE);
240 u.d = dbl; /* grab chunks of the double */
244 if (low == 0 && (high & ~DMSBIT) == 0) {
249 iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
253 if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
257 /* A denorm, normalize the mantissa */
258 while (! (high & DHIGHBIT)) {
270 *man_high = - *man_high;
276 __decodeFloat (MP_INT *man, I_ *exp, StgFloat flt)
278 /* Do some bit fiddling on IEEE */
279 int high, sign; /* assuming 32 bit ints */
280 union { float f; int i; } u; /* assuming 32 bit float and int */
282 ASSERT(sizeof(int ) == 4 );
283 ASSERT(sizeof(flt ) == SIZEOF_FLOAT );
284 ASSERT(sizeof(man->_mp_d[0]) == SIZEOF_LIMB_T);
285 ASSERT(FNBIGIT*SIZEOF_LIMB_T >= SIZEOF_FLOAT );
287 u.f = flt; /* grab the float */
290 /* we know the MP_INT* passed in has size zero, so we realloc
293 man->_mp_alloc = FNBIGIT;
295 if ((high & ~FMSBIT) == 0) {
299 man->_mp_size = FNBIGIT;
300 *exp = ((high >> 23) & 0xff) + MY_FMINEXP;
304 if (*exp != MY_FMINEXP) /* don't add hidden bit to denorms */
308 /* A denorm, normalize the mantissa */
309 while (! (high & FHIGHBIT)) {
315 man->_mp_d[0] = (mp_limb_t)high;
317 #error Cannot cope with FNBIGIT
320 man->_mp_size = -man->_mp_size;
324 /* Convenient union types for checking the layout of IEEE 754 types -
325 based on defs in GNU libc <ieee754.h>
329 __decodeFloat_Int (I_ *man, I_ *exp, StgFloat flt)
331 /* Do some bit fiddling on IEEE */
332 int high, sign; /* assuming 32 bit ints */
333 union { float f; int i; } u; /* assuming 32 bit float and int */
335 ASSERT(sizeof(int ) == 4 );
336 ASSERT(sizeof(flt ) == 4 );
337 ASSERT(sizeof(flt ) == SIZEOF_FLOAT );
339 u.f = flt; /* grab the float */
342 if ((high & ~FMSBIT) == 0) {
346 *exp = ((high >> 23) & 0xff) + MY_FMINEXP;
350 if (*exp != MY_FMINEXP) /* don't add hidden bit to denorms */
354 /* A denorm, normalize the mantissa */
355 while (! (high & FHIGHBIT)) {
366 union stg_ieee754_flt
372 unsigned int negative:1;
373 unsigned int exponent:8;
374 unsigned int mantissa:23;
376 unsigned int mantissa:23;
377 unsigned int exponent:8;
378 unsigned int negative:1;
384 unsigned int negative:1;
385 unsigned int exponent:8;
386 unsigned int quiet_nan:1;
387 unsigned int mantissa:22;
389 unsigned int mantissa:22;
390 unsigned int quiet_nan:1;
391 unsigned int exponent:8;
392 unsigned int negative:1;
399 To recap, here's the representation of a double precision
400 IEEE floating point number:
402 sign 63 sign bit (0==positive, 1==negative)
403 exponent 62-52 exponent (biased by 1023)
404 fraction 51-0 fraction (bits to right of binary point)
407 union stg_ieee754_dbl
413 unsigned int negative:1;
414 unsigned int exponent:11;
415 unsigned int mantissa0:20;
416 unsigned int mantissa1:32;
418 #if FLOAT_WORDS_BIGENDIAN
419 unsigned int mantissa0:20;
420 unsigned int exponent:11;
421 unsigned int negative:1;
422 unsigned int mantissa1:32;
424 unsigned int mantissa1:32;
425 unsigned int mantissa0:20;
426 unsigned int exponent:11;
427 unsigned int negative:1;
431 /* This format makes it easier to see if a NaN is a signalling NaN. */
435 unsigned int negative:1;
436 unsigned int exponent:11;
437 unsigned int quiet_nan:1;
438 unsigned int mantissa0:19;
439 unsigned int mantissa1:32;
441 #if FLOAT_WORDS_BIGENDIAN
442 unsigned int mantissa0:19;
443 unsigned int quiet_nan:1;
444 unsigned int exponent:11;
445 unsigned int negative:1;
446 unsigned int mantissa1:32;
448 unsigned int mantissa1:32;
449 unsigned int mantissa0:19;
450 unsigned int quiet_nan:1;
451 unsigned int exponent:11;
452 unsigned int negative:1;
459 * Predicates for testing for extreme IEEE fp values. Used
460 * by the bytecode evaluator and the Prelude.
464 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
465 #ifdef IEEE_FLOATING_POINT
468 isDoubleNaN(StgDouble d)
470 union stg_ieee754_dbl u;
475 u.ieee.exponent == 2047 /* 2^11 - 1 */ && /* Is the exponent all ones? */
476 (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)
477 /* and the mantissa non-zero? */
482 isDoubleInfinite(StgDouble d)
484 union stg_ieee754_dbl u;
488 /* Inf iff exponent is all ones, mantissa all zeros */
490 u.ieee.exponent == 2047 /* 2^11 - 1 */ &&
491 u.ieee.mantissa0 == 0 &&
492 u.ieee.mantissa1 == 0
497 isDoubleDenormalized(StgDouble d)
499 union stg_ieee754_dbl u;
503 /* A (single/double/quad) precision floating point number
506 - mantissa is non-zero.
507 - (don't care about setting of sign bit.)
511 u.ieee.exponent == 0 &&
512 (u.ieee.mantissa0 != 0 ||
513 u.ieee.mantissa1 != 0)
519 isDoubleNegativeZero(StgDouble d)
521 union stg_ieee754_dbl u;
524 /* sign (bit 63) set (only) => negative zero */
527 u.ieee.negative == 1 &&
528 u.ieee.exponent == 0 &&
529 u.ieee.mantissa0 == 0 &&
530 u.ieee.mantissa1 == 0);
533 /* Same tests, this time for StgFloats. */
536 To recap, here's the representation of a single precision
537 IEEE floating point number:
539 sign 31 sign bit (0 == positive, 1 == negative)
540 exponent 30-23 exponent (biased by 127)
541 fraction 22-0 fraction (bits to right of binary point)
546 isFloatNaN(StgFloat f)
548 union stg_ieee754_flt u;
551 /* Floating point NaN iff exponent is all ones, mantissa is
552 non-zero (but see below.) */
554 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
555 u.ieee.mantissa != 0);
559 isFloatInfinite(StgFloat f)
561 union stg_ieee754_flt u;
564 /* A float is Inf iff exponent is max (all ones),
565 and mantissa is min(all zeros.) */
567 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
568 u.ieee.mantissa == 0);
572 isFloatDenormalized(StgFloat f)
574 union stg_ieee754_flt u;
577 /* A (single/double/quad) precision floating point number
580 - mantissa is non-zero.
581 - (don't care about setting of sign bit.)
585 u.ieee.exponent == 0 &&
586 u.ieee.mantissa != 0);
590 isFloatNegativeZero(StgFloat f)
592 union stg_ieee754_flt u;
595 /* sign (bit 31) set (only) => negative zero */
598 u.ieee.exponent == 0 &&
599 u.ieee.mantissa == 0);
602 #else /* ! IEEE_FLOATING_POINT */
604 /* Dummy definitions of predicates - they all return false */
605 StgInt isDoubleNaN(d) StgDouble d; { return 0; }
606 StgInt isDoubleInfinite(d) StgDouble d; { return 0; }
607 StgInt isDoubleDenormalized(d) StgDouble d; { return 0; }
608 StgInt isDoubleNegativeZero(d) StgDouble d; { return 0; }
609 StgInt isFloatNaN(f) StgFloat f; { return 0; }
610 StgInt isFloatInfinite(f) StgFloat f; { return 0; }
611 StgInt isFloatDenormalized(f) StgFloat f; { return 0; }
612 StgInt isFloatNegativeZero(f) StgFloat f; { return 0; }
614 #endif /* ! IEEE_FLOATING_POINT */