Fix rounding ZZ -> Python float

[jdemeyer]

Consider the following conversions:
A. Sage Integer -> Python float (or Sage RDF)
B. Python int/long -> Python float
C. Sage Integer -> Sage RR
D. C unsigned long -> C double

Conversion A rounds to zero, while the others round to nearest (by default). We should make A consistent with the rest:

sage: int(float(ZZ(10^17-1)))
99999999999999984
sage: int(float(int(10^17-1)))
100000000000000000
sage: int(RR(10^17-1))
100000000000000000
sage: cython("""print <unsigned long>(<double> (%sUL))""" % (10^17-1))
100000000000000000

Currently, Integer -> float uses mpz_get_d() which rounds to zero. We should fix this, similar to #14416.

CC: @zimmermann6

Component: basic arithmetic

Author: Jeroen Demeyer

Branch/Commit: 04138f3

Reviewer: Christoph Lauter, Marc Mezzarobba

Issue created by migration from https://trac.sagemath.org/ticket/16385

[jdemeyer]

Branch: u/jdemeyer/ticket/16385

[jdemeyer]

Commit: 04138f3

[jdemeyer]

New commits:

04138f3

Correct rounding in ZZ -> float conversion

[mezzarobba]

Reviewer: Christoph Lauter, Marc Mezzarobba

[mezzarobba]

comment:6

Looks good to me. But Christoph Lauter sitting next to me noticed that expressions like 1.0/0.0 raise the Division by zero FP exception, while IEEE-754 conversions are supposed to raise Overflow when the result is too large to be represented in the target format.

Another minor issue is that your mpz_get_d_nearest actually assumes that the current hardware rounding mode is round-to-nearest, because the final call to ldexp can overflow:

sage: b = 2^1024-1
sage: float(b)
inf
sage: from ctypes import cdll                                                                 
sage: from ctypes.util import find_library                                                    
sage: libm = cdll.LoadLibrary(find_library('m'))                                              
sage: FE_TOWARDZERO = int(0xc00)                                                              
sage: libm.fesetround(FE_TOWARDZERO)
0
sage: float(b)
1.7976931348623157e+308

But I don't know if we really want to go into this kind of business in sage, and your patch clearly improves the previous implementation. So I'll give it positive review and let you decide if you want to support FP flags and non-default rounding modes.

Christoph also suggests to simplify the code a bit by remplacing the part that rounds from 64 to 63 bits by something like

d = <double> ((q64 << 1 ) | !remainder_is_zero)

(Of course, this version also assumes that the FPU rounds to nearest.)

[mezzarobba]

comment:7

Another minor point in case you decide to further improve the implementation: why test it with ZZ((2^53 - 3/4) * 2^971) rather than ZZ((2^53 - 1/2) * 2^971) - 1?

[vbraun]

Changed branch from u/jdemeyer/ticket/16385 to 04138f3