+/*
+ There are glibc versions around with buggy rintf or rint, hence we
+ provide our own. We always round ties to even, so we can be simpler.
+*/
+
+#define FLT_HIDDEN 0x800000
+#define FLT_POWER2 0x1000000
+
+HsFloat
+rintFloat(HsFloat f)
+{
+ union stg_ieee754_flt u;
+ u.f = f;
+ /* if real exponent > 22, it's already integral, infinite or nan */
+ if (u.ieee.exponent > 149) /* 22 + 127 */
+ {
+ return u.f;
+ }
+ if (u.ieee.exponent < 126) /* (-1) + 127, abs(f) < 0.5 */
+ {
+ /* only used for rounding to Integral a, so don't care about -0.0 */
+ return 0.0;
+ }
+ /* 0.5 <= abs(f) < 2^23 */
+ unsigned int half, mask, mant, frac;
+ half = 1 << (149 - u.ieee.exponent); /* bit for 0.5 */
+ mask = 2*half - 1; /* fraction bits */
+ mant = u.ieee.mantissa | FLT_HIDDEN; /* add hidden bit */
+ frac = mant & mask; /* get fraction */
+ mant ^= frac; /* truncate mantissa */
+ if ((frac < half) || ((frac == half) && ((mant & (2*half)) == 0)))
+ {
+ /* this means we have to truncate */
+ if (mant == 0)
+ {
+ /* f == ±0.5, return 0.0 */
+ return 0.0;
+ }
+ else
+ {
+ /* remove hidden bit and set mantissa */
+ u.ieee.mantissa = mant ^ FLT_HIDDEN;
+ return u.f;
+ }
+ }
+ else
+ {
+ /* round away from zero, increment mantissa */
+ mant += 2*half;
+ if (mant == FLT_POWER2)
+ {
+ /* next power of 2, increase exponent an set mantissa to 0 */
+ u.ieee.mantissa = 0;
+ u.ieee.exponent += 1;
+ return u.f;
+ }
+ else
+ {
+ /* remove hidden bit and set mantissa */
+ u.ieee.mantissa = mant ^ FLT_HIDDEN;
+ return u.f;
+ }
+ }
+}
+
+#define DBL_HIDDEN 0x100000
+#define DBL_POWER2 0x200000
+#define LTOP_BIT 0x80000000
+
+HsDouble
+rintDouble(HsDouble d)
+{
+ union stg_ieee754_dbl u;
+ u.d = d;
+ /* if real exponent > 51, it's already integral, infinite or nan */
+ if (u.ieee.exponent > 1074) /* 51 + 1023 */
+ {
+ return u.d;
+ }
+ if (u.ieee.exponent < 1022) /* (-1) + 1023, abs(d) < 0.5 */
+ {
+ /* only used for rounding to Integral a, so don't care about -0.0 */
+ return 0.0;
+ }
+ unsigned int half, mask, mant, frac;
+ if (u.ieee.exponent < 1043) /* 20 + 1023, real exponent < 20 */
+ {
+ /* the fractional part meets the higher part of the mantissa */
+ half = 1 << (1042 - u.ieee.exponent); /* bit for 0.5 */
+ mask = 2*half - 1; /* fraction bits */
+ mant = u.ieee.mantissa0 | DBL_HIDDEN; /* add hidden bit */
+ frac = mant & mask; /* get fraction */
+ mant ^= frac; /* truncate mantissa */
+ if ((frac < half) ||
+ ((frac == half) && (u.ieee.mantissa1 == 0) /* a tie */
+ && ((mant & (2*half)) == 0)))
+ {
+ /* truncate */
+ if (mant == 0)
+ {
+ /* d = ±0.5, return 0.0 */
+ return 0.0;
+ }
+ /* remove hidden bit and set mantissa */
+ u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
+ u.ieee.mantissa1 = 0;
+ return u.d;
+ }
+ else /* round away from zero */
+ {
+ /* zero low mantissa bits */
+ u.ieee.mantissa1 = 0;
+ /* increment integer part of mantissa */
+ mant += 2*half;
+ if (mant == DBL_POWER2)
+ {
+ /* power of 2, increment exponent and zero mantissa */
+ u.ieee.mantissa0 = 0;
+ u.ieee.exponent += 1;
+ return u.d;
+ }
+ /* remove hidden bit */
+ u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
+ return u.d;
+ }
+ }
+ else
+ {
+ /* 20 <= real exponent < 52, fractional part entirely in mantissa1 */
+ half = 1 << (1074 - u.ieee.exponent); /* bit for 0.5 */
+ mask = 2*half - 1; /* fraction bits */
+ mant = u.ieee.mantissa1; /* no hidden bit here */
+ frac = mant & mask; /* get fraction */
+ mant ^= frac; /* truncate mantissa */
+ if ((frac < half) ||
+ ((frac == half) && /* tie */
+ (((half == LTOP_BIT) ? (u.ieee.mantissa0 & 1) /* yuck */
+ : (mant & (2*half)))
+ == 0)))
+ {
+ /* truncate */
+ u.ieee.mantissa1 = mant;
+ return u.d;
+ }
+ else
+ {
+ /* round away from zero */
+ /* increment mantissa */
+ mant += 2*half;
+ u.ieee.mantissa1 = mant;
+ if (mant == 0)
+ {
+ /* low part of mantissa overflowed */
+ /* increment high part of mantissa */
+ mant = u.ieee.mantissa0 + 1;
+ if (mant == DBL_HIDDEN)
+ {
+ /* hit power of 2 */
+ /* zero mantissa */
+ u.ieee.mantissa0 = 0;
+ /* and increment exponent */
+ u.ieee.exponent += 1;
+ return u.d;
+ }
+ else
+ {
+ u.ieee.mantissa0 = mant;
+ return u.d;
+ }
+ }
+ else
+ {
+ return u.d;
+ }
+ }
+ }
+}
+