X-Git-Url: http://git.megacz.com/?a=blobdiff_plain;f=rts%2FStgPrimFloat.c;h=e523f328c33b14f804acc29fb5ca7088b025e6fd;hb=7d9eb2e45b4a9ff4cb053b1ec37602be88528b62;hp=95210948872edf01f00d379dafbeb635db21caff;hpb=51f782f416a8aa05153a9f2dcef34d8bcb79f1ed;p=ghc-hetmet.git diff --git a/rts/StgPrimFloat.c b/rts/StgPrimFloat.c index 9521094..e523f32 100644 --- a/rts/StgPrimFloat.c +++ b/rts/StgPrimFloat.c @@ -1,5 +1,6 @@ /* ----------------------------------------------------------------------------- * + * (c) Lennart Augustsson * (c) The GHC Team, 1998-2000 * * Miscellaneous support for floating-point primitives @@ -9,34 +10,18 @@ #include "PosixSource.h" #include "Rts.h" +#include "StgPrimFloat.h" + #include +#include + +#define IEEE_FLOATING_POINT 1 /* * Encoding and decoding Doubles. Code based on the HBC code * (lib/fltcode.c). */ -#ifdef _SHORT_LIMB -#define SIZEOF_LIMB_T SIZEOF_UNSIGNED_INT -#else -#ifdef _LONG_LONG_LIMB -#define SIZEOF_LIMB_T SIZEOF_UNSIGNED_LONG_LONG -#else -#define SIZEOF_LIMB_T SIZEOF_UNSIGNED_LONG -#endif -#endif - -#if SIZEOF_LIMB_T == 4 -#define GMP_BASE 4294967296.0 -#elif SIZEOF_LIMB_T == 8 -#define GMP_BASE 18446744073709551616.0 -#else -#error Cannot cope with SIZEOF_LIMB_T -- please add definition of GMP_BASE -#endif - -#define DNBIGIT ((SIZEOF_DOUBLE+SIZEOF_LIMB_T-1)/SIZEOF_LIMB_T) -#define FNBIGIT ((SIZEOF_FLOAT +SIZEOF_LIMB_T-1)/SIZEOF_LIMB_T) - #if IEEE_FLOATING_POINT #define MY_DMINEXP ((DBL_MIN_EXP) - (DBL_MANT_DIG) - 1) /* DMINEXP is defined in values.h on Linux (for example) */ @@ -59,25 +44,41 @@ #define __abs(a) (( (a) >= 0 ) ? (a) : (-(a))) StgDouble -__encodeDouble (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */ +__2Int_encodeDouble (I_ j_high, I_ j_low, I_ e) { - StgDouble r; - const mp_limb_t *const arr = (const mp_limb_t *)ba; - I_ i; - - /* Convert MP_INT to a double; knows a lot about internal rep! */ - for(r = 0.0, i = __abs(size)-1; i >= 0; i--) - r = (r * GMP_BASE) + arr[i]; - - /* Now raise to the exponent */ - if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */ - r = ldexp(r, e); + StgDouble r; + + /* assuming 32 bit ints */ + ASSERT(sizeof(int ) == 4 ); - /* sign is encoded in the size */ - if (size < 0) - r = -r; + r = (StgDouble)((unsigned int)j_high); + r *= 4294967296.0; /* exp2f(32); */ + r += (StgDouble)((unsigned int)j_low); + + /* Now raise to the exponent */ + if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */ + r = ldexp(r, e); + + /* sign is encoded in the size */ + if (j_high < 0) + r = -r; + + return r; +} - return r; +/* Special version for words */ +StgDouble +__word_encodeDouble (W_ j, I_ e) +{ + StgDouble r; + + r = (StgDouble)j; + + /* Now raise to the exponent */ + if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */ + r = ldexp(r, e); + + return r; } /* Special version for small Integers */ @@ -99,28 +100,6 @@ __int_encodeDouble (I_ j, I_ e) return r; } -StgFloat -__encodeFloat (I_ size, StgByteArray ba, I_ e) /* result = s * 2^e */ -{ - StgFloat r; - const mp_limb_t *arr = (const mp_limb_t *)ba; - I_ i; - - /* Convert MP_INT to a float; knows a lot about internal rep! */ - for(r = 0.0, i = __abs(size)-1; i >= 0; i--) - r = (r * GMP_BASE) + arr[i]; - - /* Now raise to the exponent */ - if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */ - r = ldexp(r, e); - - /* sign is encoded in the size */ - if (size < 0) - r = -r; - - return r; -} - /* Special version for small Integers */ StgFloat __int_encodeFloat (I_ j, I_ e) @@ -140,10 +119,25 @@ __int_encodeFloat (I_ j, I_ e) return r; } +/* Special version for small positive Integers */ +StgFloat +__word_encodeFloat (W_ j, I_ e) +{ + StgFloat r; + + r = (StgFloat)j; + + /* Now raise to the exponent */ + if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */ + r = ldexp(r, e); + + return r; +} + /* This only supports IEEE floating point */ void -__decodeDouble (MP_INT *man, I_ *exp, StgDouble dbl) +__decodeDouble_2Int (I_ *man_sign, W_ *man_high, W_ *man_low, I_ *exp, StgDouble dbl) { /* Do some bit fiddling on IEEE */ unsigned int low, high; /* assuming 32 bit ints */ @@ -151,24 +145,18 @@ __decodeDouble (MP_INT *man, I_ *exp, StgDouble dbl) union { double d; unsigned int i[2]; } u; /* assuming 32 bit ints, 64 bit double */ ASSERT(sizeof(unsigned int ) == 4 ); + ASSERT(sizeof(dbl ) == 8 ); ASSERT(sizeof(dbl ) == SIZEOF_DOUBLE); - ASSERT(sizeof(man->_mp_d[0]) == SIZEOF_LIMB_T); - ASSERT(DNBIGIT*SIZEOF_LIMB_T >= SIZEOF_DOUBLE); u.d = dbl; /* grab chunks of the double */ low = u.i[L]; high = u.i[H]; - /* we know the MP_INT* passed in has size zero, so we realloc - no matter what. - */ - man->_mp_alloc = DNBIGIT; - if (low == 0 && (high & ~DMSBIT) == 0) { - man->_mp_size = 0; + *man_low = 0; + *man_high = 0; *exp = 0L; } else { - man->_mp_size = DNBIGIT; iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP; sign = high; @@ -187,46 +175,34 @@ __decodeDouble (MP_INT *man, I_ *exp, StgDouble dbl) } } *exp = (I_) iexp; -#if DNBIGIT == 2 - man->_mp_d[0] = (mp_limb_t)low; - man->_mp_d[1] = (mp_limb_t)high; -#else -#if DNBIGIT == 1 - man->_mp_d[0] = ((mp_limb_t)high) << 32 | (mp_limb_t)low; -#else -#error Cannot cope with DNBIGIT -#endif -#endif - if (sign < 0) - man->_mp_size = -man->_mp_size; + *man_low = low; + *man_high = high; + *man_sign = (sign < 0) ? -1 : 1; } } +/* Convenient union types for checking the layout of IEEE 754 types - + based on defs in GNU libc +*/ + void -__decodeFloat (MP_INT *man, I_ *exp, StgFloat flt) +__decodeFloat_Int (I_ *man, I_ *exp, StgFloat flt) { /* Do some bit fiddling on IEEE */ int high, sign; /* assuming 32 bit ints */ union { float f; int i; } u; /* assuming 32 bit float and int */ ASSERT(sizeof(int ) == 4 ); + ASSERT(sizeof(flt ) == 4 ); ASSERT(sizeof(flt ) == SIZEOF_FLOAT ); - ASSERT(sizeof(man->_mp_d[0]) == SIZEOF_LIMB_T); - ASSERT(FNBIGIT*SIZEOF_LIMB_T >= SIZEOF_FLOAT ); u.f = flt; /* grab the float */ high = u.i; - /* we know the MP_INT* passed in has size zero, so we realloc - no matter what. - */ - man->_mp_alloc = FNBIGIT; - if ((high & ~FMSBIT) == 0) { - man->_mp_size = 0; + *man = 0; *exp = 0; } else { - man->_mp_size = FNBIGIT; *exp = ((high >> 23) & 0xff) + MY_FMINEXP; sign = high; @@ -241,266 +217,9 @@ __decodeFloat (MP_INT *man, I_ *exp, StgFloat flt) (*exp)--; } } -#if FNBIGIT == 1 - man->_mp_d[0] = (mp_limb_t)high; -#else -#error Cannot cope with FNBIGIT -#endif + *man = high; if (sign < 0) - man->_mp_size = -man->_mp_size; + *man = - *man; } } -/* Convenient union types for checking the layout of IEEE 754 types - - based on defs in GNU libc -*/ - -union stg_ieee754_flt -{ - float f; - struct { - -#if WORDS_BIGENDIAN - unsigned int negative:1; - unsigned int exponent:8; - unsigned int mantissa:23; -#else - unsigned int mantissa:23; - unsigned int exponent:8; - unsigned int negative:1; -#endif - } ieee; - struct { - -#if WORDS_BIGENDIAN - unsigned int negative:1; - unsigned int exponent:8; - unsigned int quiet_nan:1; - unsigned int mantissa:22; -#else - unsigned int mantissa:22; - unsigned int quiet_nan:1; - unsigned int exponent:8; - unsigned int negative:1; -#endif - } ieee_nan; -}; - -/* - - To recap, here's the representation of a double precision - IEEE floating point number: - - sign 63 sign bit (0==positive, 1==negative) - exponent 62-52 exponent (biased by 1023) - fraction 51-0 fraction (bits to right of binary point) -*/ - -union stg_ieee754_dbl -{ - double d; - struct { - -#if WORDS_BIGENDIAN - unsigned int negative:1; - unsigned int exponent:11; - unsigned int mantissa0:20; - unsigned int mantissa1:32; -#else -#if FLOAT_WORDS_BIGENDIAN - unsigned int mantissa0:20; - unsigned int exponent:11; - unsigned int negative:1; - unsigned int mantissa1:32; -#else - unsigned int mantissa1:32; - unsigned int mantissa0:20; - unsigned int exponent:11; - unsigned int negative:1; -#endif -#endif - } ieee; - /* This format makes it easier to see if a NaN is a signalling NaN. */ - struct { - -#if WORDS_BIGENDIAN - unsigned int negative:1; - unsigned int exponent:11; - unsigned int quiet_nan:1; - unsigned int mantissa0:19; - unsigned int mantissa1:32; -#else -#if FLOAT_WORDS_BIGENDIAN - unsigned int mantissa0:19; - unsigned int quiet_nan:1; - unsigned int exponent:11; - unsigned int negative:1; - unsigned int mantissa1:32; -#else - unsigned int mantissa1:32; - unsigned int mantissa0:19; - unsigned int quiet_nan:1; - unsigned int exponent:11; - unsigned int negative:1; -#endif -#endif - } ieee_nan; -}; - -/* - * Predicates for testing for extreme IEEE fp values. Used - * by the bytecode evaluator and the Prelude. - * - */ - -/* In case you don't suppport IEEE, you'll just get dummy defs.. */ -#ifdef IEEE_FLOATING_POINT - -StgInt -isDoubleNaN(StgDouble d) -{ - union stg_ieee754_dbl u; - - u.d = d; - - return ( - u.ieee.exponent == 2047 /* 2^11 - 1 */ && /* Is the exponent all ones? */ - (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0) - /* and the mantissa non-zero? */ - ); -} - -StgInt -isDoubleInfinite(StgDouble d) -{ - union stg_ieee754_dbl u; - - u.d = d; - - /* Inf iff exponent is all ones, mantissa all zeros */ - return ( - u.ieee.exponent == 2047 /* 2^11 - 1 */ && - u.ieee.mantissa0 == 0 && - u.ieee.mantissa1 == 0 - ); -} - -StgInt -isDoubleDenormalized(StgDouble d) -{ - union stg_ieee754_dbl u; - - u.d = d; - - /* A (single/double/quad) precision floating point number - is denormalised iff: - - exponent is zero - - mantissa is non-zero. - - (don't care about setting of sign bit.) - - */ - return ( - u.ieee.exponent == 0 && - (u.ieee.mantissa0 != 0 || - u.ieee.mantissa1 != 0) - ); - -} - -StgInt -isDoubleNegativeZero(StgDouble d) -{ - union stg_ieee754_dbl u; - - u.d = d; - /* sign (bit 63) set (only) => negative zero */ - - return ( - u.ieee.negative == 1 && - u.ieee.exponent == 0 && - u.ieee.mantissa0 == 0 && - u.ieee.mantissa1 == 0); -} - -/* Same tests, this time for StgFloats. */ - -/* - To recap, here's the representation of a single precision - IEEE floating point number: - - sign 31 sign bit (0 == positive, 1 == negative) - exponent 30-23 exponent (biased by 127) - fraction 22-0 fraction (bits to right of binary point) -*/ - - -StgInt -isFloatNaN(StgFloat f) -{ - union stg_ieee754_flt u; - u.f = f; - - /* Floating point NaN iff exponent is all ones, mantissa is - non-zero (but see below.) */ - return ( - u.ieee.exponent == 255 /* 2^8 - 1 */ && - u.ieee.mantissa != 0); -} - -StgInt -isFloatInfinite(StgFloat f) -{ - union stg_ieee754_flt u; - u.f = f; - - /* A float is Inf iff exponent is max (all ones), - and mantissa is min(all zeros.) */ - return ( - u.ieee.exponent == 255 /* 2^8 - 1 */ && - u.ieee.mantissa == 0); -} - -StgInt -isFloatDenormalized(StgFloat f) -{ - union stg_ieee754_flt u; - u.f = f; - - /* A (single/double/quad) precision floating point number - is denormalised iff: - - exponent is zero - - mantissa is non-zero. - - (don't care about setting of sign bit.) - - */ - return ( - u.ieee.exponent == 0 && - u.ieee.mantissa != 0); -} - -StgInt -isFloatNegativeZero(StgFloat f) -{ - union stg_ieee754_flt u; - u.f = f; - - /* sign (bit 31) set (only) => negative zero */ - return ( - u.ieee.negative && - u.ieee.exponent == 0 && - u.ieee.mantissa == 0); -} - -#else /* ! IEEE_FLOATING_POINT */ - -/* Dummy definitions of predicates - they all return false */ -StgInt isDoubleNaN(d) StgDouble d; { return 0; } -StgInt isDoubleInfinite(d) StgDouble d; { return 0; } -StgInt isDoubleDenormalized(d) StgDouble d; { return 0; } -StgInt isDoubleNegativeZero(d) StgDouble d; { return 0; } -StgInt isFloatNaN(f) StgFloat f; { return 0; } -StgInt isFloatInfinite(f) StgFloat f; { return 0; } -StgInt isFloatDenormalized(f) StgFloat f; { return 0; } -StgInt isFloatNegativeZero(f) StgFloat f; { return 0; } - -#endif /* ! IEEE_FLOATING_POINT */