If n is a Python long, then one might expect float(n) to return the 
closest float to n.  Currently it doesn't do this.  For example (with 
Python 2.6, on OS X 10.5.2/Intel):

>>> n = 295147905179352891391L

The closest float to n is equal to n+1.  But float(n) returns the 
further of the two floats bracketing n, equal to n-65535:

>>> float(n)
2.9514790517935283e+20
>>> long(float(n))
295147905179352825856L
>>> n - long(float(n))
65535L

It's fairly straightforward to fix PyLong_AsDouble to return the closest 
double to a given long integer n (using the round-half-to-even rule in 
the case of a tie).  The attached patch does this.

Having a correctly rounded float(n) can be useful for testing other 
floating-point routines that are supposed to be correctly rounded.

I agree, longs should be correctly rounded when coerced to floats.

There is an ugly (but amusing) workaround while people wait for this 
patch:  Go via a string:

int(float(repr(295147905179352891391)[:-1]))

Though I assume this relies on the platform's strtod working correctly.  
Which it does for me.

You may use "if (nbits == (size_t)-1 && PyErr_Occurred())" to check 
_PyLong_NumBits() error (overflow). Well, "if (numbits > DBL_MAX_EXP)" 
should already catch overflow, but I prefer explicit test to check the 
error case.

Anyway, interresting patch! Python3 vanilla:
>>> n = 295147905179352891391; int(float(n)) - n
-65535

Python3 + your patch:
>>> int(float(n)) - n
1

Mark,

I noticed that you replaced a call to _PyLong_AsScaledDouble with your 
round to nearest algorithm.  I wonder if _PyLong_AsScaledDouble itself 
would benefit from your change.  Currently it is used in PyLong_AsDouble 
and long_true_divide.  I would think that long_true_divide would be more 
accurate if longs were rounded to the nearest float.

I also wonder whether round to nearest float can be implemented without 
floating point arithmetics.  I would think round towards zero should be 
a simple matter of extracting an appropriate number of bits from the 
long and round to nearest would at most require a long addition.

I believe _PyLong_AsScaledDouble is written the way it is to support 
non-IEEE floating formats, but I am not sure that your algorithm would 
always return the nearest float on an arbitrary non-IEEE platform.

Maybe it would be worthwhile to provide a simple IEEE specific code with   
well specified semantics for both PyLong_AsDouble and long_true_divide, 
but fall back to the current code on non-IEEE platforms.

float(295147905179352891391L) gives different result on Python 2.5 and 
Python 2.6:
- 2.9514790517935289e+20   # Python 2.5.1
- 2.9514790517935283e+20   # 2.7a0

whereas the code is the same!?

Python 2.5.1 (r251:54863, Jul 31 2008, 23:17:40)
>>> reduce(lambda x,y: x*32768.0 + y, [256, 0, 0, 1, 32767])
2.9514790517935283e+20
>>> float(295147905179352891391L)
2.9514790517935289e+20

Python 2.7a0 (trunk:67679M, Dec  9 2008, 14:29:12)
>>> reduce(lambda x,y: x*32768.0 + y, [256, 0, 0, 1, 32767])
2.9514790517935283e+20
>>> float(295147905179352891391L)
2.9514790517935283e+20

Python 3.1a0 (py3k:67652M, Dec  9 2008, 13:08:19)
>>> float(295147905179352891391)
2.9514790517935283e+20
>>> digits=[256, 0, 0, 1, 32767]; x=0
>>> for d in digits:
...  x*=32768.0
...  x+= d
...
>>> x
2.9514790517935283e+20

All results are the same, except float(295147905179352891391L) in 
Python 2.5!? Python 2.5 rounds correctly:

Python 2.5.1 (r251:54863, Jul 31 2008, 23:17:40)
>>> x=295147905179352891391L
>>> long(float(long(x))) - x
1L

Ok, I understand why different versions of the same code gives 
different results: compiler flags! Python 2.5.1 is my Ubuntu version 
(should be compiled with -O3) whereas Python 2.7 and 3.1a0 are 
compiled by me with -00.

Results with Python 2.5.1:
 - with -O0, float(295147905179352891391L) gives 
2.9514790517935283e+20
 - with -O1, float(295147905179352891391L) gives 
2.9514790517935289e+20

I'm unable to isolate the exact compiler flag which changes the 
result. I tried all options listed in the gcc doc for the -O1 option:
http://gcc.gnu.org/onlinedocs/gcc-4.1.2/gcc/Optimize-Options.html#Optimize-Options

Victor, what does

>>> 1e16 + 2.9999

give on your Ubuntu 2.5 machine?
(Humor me. :) )

About -O0 vs -O1, I think that I understood (by reading the 
assembler).

pseudocode of the -O0 version:
  while (....)
  {
     load x from the stack
     x = x * ... + ...
     write x to the stack
  }

pseudocode of the -O1 version:
  while (....)
  {
     x = x * ... + ...
  }

Intel uses 80 bits float in internals, but load/store uses 64 bits 
float. Load/store looses least significant bits.

Hey, floating point numbers are funny :-)

---

Python 2.5.1 (r251:54863, Jul 31 2008, 23:17:40)
>>> 1e16 + 2.999
10000000000000002.0
>>> 1e16 + 2.9999
10000000000000004.0

Same result with python 2.7/3.1.

An interresting document: "Request for Comments: Rounding in PHP"
http://wiki.php.net/rfc/rounding

> Intel uses 80 bits float in internals, but load/store uses 64 bits 
> float. Load/store looses least significant bits.

Exactly.  If your Intel machine is Pentium 4 or newer, you can get
around this by using the SSE2 extensions, which work with 64-bit doubles
throughout.  I don't remember offhand precisely which gcc flags you need
for this.

> Hey, floating point numbers are funny :-)

Yup.

On Tue, Dec 9, 2008 at 11:02 AM, Mark Dickinson <report@bugs.python.org> wrote:
...
> If your Intel machine is Pentium 4 or newer, you can get
> around this by using the SSE2 extensions, which work with 64-bit doubles
> throughout.  I don't remember offhand precisely which gcc flags you need
> for this.

The flags you may be looking for are -msse2 -mfpmath=sse

[Alexander]
> The flags you may be looking for are -msse2 -mfpmath=sse

Thanks, Alexander!

[Alexander again, from an earlier post...]
> I noticed that you replaced a call to _PyLong_AsScaledDouble with your 
> round to nearest algorithm.  I wonder if _PyLong_AsScaledDouble itself 
> would benefit from your change.  Currently it is used in 
PyLong_AsDouble 
> and long_true_divide.  I would think that long_true_divide would be 
more 
> accurate if longs were rounded to the nearest float.

You read my mind!  I've got another issue open about making long 
division round correctly, somewhere.  And indeed I'd like to make 
_PyLong_AsScaledDouble do correct rounding.  (I'd also like to make it 
return the exponent in bits, rather than digits, so that mathmodule.c 
doesn't have to know about the long int representation, but that's 
another story...)

> I believe _PyLong_AsScaledDouble is written the way it is to support 
> non-IEEE floating formats, but I am not sure that your algorithm would 
> always return the nearest float on an arbitrary non-IEEE platform.

Well, I had other possible formats in mind when I wrote the code, and I 
hope it's good whenever FLT_RADIX is 2.  If you can think of explicit 
cases where it's not going to work, please let me know.  My belief that 
it's safe rests on two facts: (1) it entirely avoids IEEE 754 oddities 
like negative zero, denormals and NaNs, and (2) all the floating-point 
operations actually performed should have exact results, so rounding 
doesn't come into play anywhere.

When FLT_RADIX is some other power of 2 (FLT_RADIX=16 is the only 
example I know of that exists in the wild) the code probably doesn't 
produce correctly rounded results in all cases, but it shouldn't do 
anything catastrophic either---I think the error still should't be more 
than 1ulp in this case.  When FLT_RADIX is not a power of 2 then so much 
else is going to be broken anyway that it's not worth worrying about.

[Alexander]
> I also wonder whether round to nearest float can be implemented 
without
> floating point arithmetics.  I would think round towards zero should 
be
> a simple matter of extracting an appropriate number of bits from the
> long and round to nearest would at most require a long addition.

The idea's attractive.  The problem is finding an integer type that's 
guaranteed to have enough bits to store the mantissa for the float 
(probably plus one or two bits more for comfort);  for IEEE 754 this 
means a 64-bit integer type, and there aren't any of those in C89.

(One could use two 32-bit integer variables, but that's going to get 
messy.)

..
> The idea's attractive.  The problem is finding an integer type that's
> guaranteed to have enough bits to store the mantissa for the float
> (probably plus one or two bits more for comfort);  for IEEE 754 this
> means a 64-bit integer type, and there aren't any of those in C89.

But Python already has an arbitrary precision integer type, why not
use it?  Performance may suffer, but optimization can be considered
later possibly first on the most popular platforms.

As you say, performance would suffer.

What would using Python's integer type solve, that isn't already solved by  
the patch?

I know the code isn't terribly readable;  I'll add some comments 
explaining clearly what's going on.

On Tue, Dec 9, 2008 at 12:39 PM, Mark Dickinson <report@bugs.python.org> wrote:
..
> What would using Python's integer type solve, that isn't already solved by
> the patch?
>

Speaking for myself, it would alleviate the irrational fear of
anything involving the FPU. :-)

Seriously, it is not obvious that your algorithm is correct and does
not depend on x being stored in an extended precision register.
Floating point experts are in short supply, so it may take a while for
your patch to be thoroughly reviewed by one of them.  If you could
produce a formal proof of correctness of your algorithm, that would be
more helpful than code comments.

On the other hand, an implementation that uses only integer
arithmetics plus bit shifts is likely to be much easier to understand.
 I don't think the performance will suffer much.  We can even start
with a pure python prototype using struct.pack or ctypes to produce
the double result.

By the way, the algorithm here is essentially the same as the algorithm that I 
implemented for the float.fromhex method, except that the float.fromhex method is more 
complicated in that it may have to deal with signed zeros or subnormals.

So any mathematical defects that you find in this patch probably indicate a defect in 
float.fromhex too.

In fact, the code *does* do integer arithmetic, except that one of the integers happens 
to be stored as a double.  If you look at the code you'll find that at every stage, the 
floating-point variable "x" has an exact nonnegative integer value between 0 and 
2**DBL_MANT_DIG.  All such values are exactly representable as a double, and all the 
arithmetic operations involved are exact.  This holds right up to the ldexp call at the 
end.  So the arithmetic with x is essentially integer arithmetic.

I accept the code needs extra documentation;  I was planning to put the equivalent 
Python code into the comments to make things clearer.

> floating-point variable "x" has an exact nonnegative integer value 
> between 0 and 2**DBL_MANT_DIG. 

Hmm.  On closer inspection that's not quite true.  After the line

x = x * PyLong_BASE + (dig & (PyLong_BASE - pmask));

x has a value of the form n * pmask, where pmask is a power of 2 and
n is in the range [0, 2**DBL_MANT_DIG).  It's still exactly represented, 
provided that FLT_RADIX is 2.  (It's the multiplications by powers of 2 
that get hairy when FLT_RADIX is 16, since they *can* lose information.)

Thanks for your comments, Alexander.

Here's a rewritten version of the patch that's better commented and 
somewhat less convoluted;  I think it should be easier to verify the 
correctness of this one.

Minor cleanup of long_as_double2.patch.

Updated patch;  cleanup of comments and slight refactoring of code.

Int->float conversions are even a speck faster than the current code, for 
small inputs.  (On my machine, on a Friday night, during a full moon.  
Your results may differ.  :))

Also, retarget this for 2.7 and 3.1.

Updated patch;  applies cleanly to current trunk.  No significant changes.

Note that there's now a new reason to apply this patch: it ensures that 
the result of a long->float conversion is independent of whether we're 
using 30-bit digits or 15-bit digits for longs.

One problem:  if long->float conversions are correctly rounded, then
int->float conversions should be correctly rounded as well.  (And ideally, 
we should have float(int(n)) == float(long(n)) for any integer n.)

This problem only affects 64-bit machines:  on 32-bit machines, all 
integers are exactly representable as floats, and the C99 standard 
specifies that in that case the conversion should be exact.

(Slightly updated version of) patch applied in r71772 (trunk),
r71773 (py3k).