Implement round-to-even compliant rounding

Author: remyoudompheng

src/flint/types/fmpz.pyx

@@ -3,7 +3,7 @@ from flint.utils.typecheck cimport typecheck
 from flint.utils.conversion cimport chars_from_str
 from flint.utils.conversion cimport str_from_chars, _str_trunc
 cimport libc.stdlib
-import math
+from libc.math cimport nextafter
 
 from flint.flintlib.types.flint cimport FMPZ_REF, FMPZ_TMP, FMPZ_UNKNOWN, COEFF_IS_MPZ
 from flint.flintlib.functions.flint cimport flint_free
@@ -106,15 +106,37 @@ cdef class fmpz(flint_scalar):
         return fmpz_get_intlong(self.val)
 
     def __float__(self):
-        if fmpz_bits(self.val) <= 1023:
+        cdef slong bits = fmpz_bits(self.val)
+        cdef long tz
+        cdef double d
+        cdef fmpz_struct xabs[1]
+        cdef int roundup = 0
+        if bits <= 1023:
             # Known to be representable by a IEEE-754 double
-            return fmpz_get_d(self.val)
-
-        cdef slong exp
-        # fmpz_get_d_2exp is always accurate
-        # math.ldexp handles overflow checks
-        cdef double d = fmpz_get_d_2exp(&exp, self.val)
-        return math.ldexp(d, exp)
+            d = fmpz_get_d(self.val)
+            # fmpz_get_d always rounds towards zero
+            # Sometimes we need to round away from zero:
+            # - if the 54-th most significant bit is 1
+            # - if further bits are not zero
+            # - or if all further bits are zero but the 53-th bit is 1 (round-to-even convention)
+            tz = fmpz_val2(self.val)
+            if fmpz_sgn(self.val) == -1:
+                fmpz_init(xabs)
+                fmpz_abs(xabs, self.val)
+                if bits >= 54 and fmpz_tstbit(xabs, bits - 54) == 1 and \
+                        (tz < bits - 54 or (tz == bits - 54 and fmpz_tstbit(xabs, bits - 53) == 1)):
+                    roundup = 1
+                fmpz_clear(xabs)
+            else:
+                if bits >= 54 and fmpz_tstbit(self.val, bits - 54) == 1 and \
+                        (tz < bits - 54 or (tz == bits - 54 and fmpz_tstbit(self.val, bits - 53) == 1)):
+                    roundup = 1
+            if roundup:
+                # increase the mantissa of d by 1
+                d = nextafter(d, 2 * d)
+            return d
+
+        return float(fmpz_get_intlong(self.val))
 
     def __floor__(self):
         return self