Make conversions from long to float correctly rounded.

[mdickinson]

BPO	3166
Nosy	@mdickinson, @abalkin, @vstinner, @drj11
Files	long_as_double3.patch long_as_double4.patch

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields: ```python assignee = 'https://github.com/mdickinson' closed_at = created_at = labels = ['interpreter-core', 'type-feature'] title = 'Make conversions from long to float correctly rounded.' updated_at = user = 'https://github.com/mdickinson' ``` bugs.python.org fields: ```python activity = actor = 'mark.dickinson' assignee = 'mark.dickinson' closed = True closed_date = closer = 'mark.dickinson' components = ['Interpreter Core'] creation = creator = 'mark.dickinson' dependencies = [] files = ['12352', '13582'] hgrepos = [] issue_num = 3166 keywords = ['patch'] message_count = 24.0 messages = ['68545', '71997', '75510', '75568', '77411', '77412', '77418', '77421', '77427', '77429', '77430', '77431', '77432', '77435', '77436', '77437', '77438', '77441', '77451', '77473', '77556', '77798', '85225', '86206'] nosy_count = 5.0 nosy_names = ['mark.dickinson', 'belopolsky', 'ggenellina', 'vstinner', 'drj'] pr_nums = [] priority = 'normal' resolution = 'accepted' stage = 'resolved' status = 'closed' superseder = None type = 'enhancement' url = 'https://bugs.python.org/issue3166' versions = ['Python 3.1', 'Python 2.7'] ```

[mdickinson]

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.

[94483907-4dc9-4819-814a-424ee96e5121]

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.

[vstinner]

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

[abalkin]

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.

[vstinner]

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!?

[vstinner]

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

[vstinner]

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

[mdickinson]

Victor, what does

>> 1e16 + 2.9999

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

[vstinner]

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.

[vstinner]

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

[mdickinson]

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.

[abalkin]

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

[mdickinson]

[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.

[mdickinson]

[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.)

[abalkin]

..

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.

[mdickinson]

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.

[abalkin]

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.

[mdickinson]

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.

[mdickinson]

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.)

[mdickinson]

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.

[mdickinson]

Minor cleanup of long_as_double2.patch.

[mdickinson]

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.

[mdickinson]

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.

[mdickinson]

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