rounding difference between Python 3 and SymPy

[Chris Smith]

Python 3 implements "round to even on tie" logic so `round(12.5)` -> 12,  not 13. I have updated, in #16608, the round function but there is a difference in how ties are detected. I shift the desired position to the ones position and then check for a tie so 12.345 is shifted to 1234.5 and rounded to 1234 then is divided by 100 to give 12.34. Python doesn't do this. I suspect it adds 0.05 and then detects that12.395 > 12395/1000 and rounds up to 12.35

>>> Rational(*.345.as_integer_ratio())-Rational(345,1000)  # .345 < 345/1000
-3/112589990684262400
>>> Rational(*.395.as_integer_ratio())-Rational(395,1000)  # .395 > 395/1000
1/56294995342131200


>>> round(12.345,2)  # 12.345 rounds up b/c a tie is not detected
12.35

Does anyone have objections to the proposed rounding?

/c

[Aaron Meurer]

Doesn't Python do rounding based on the binary representation of the float?

I'm a little confused what "round to even" means in that case.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/a84085c6-aa90-437c-b063-a87f909beac4%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Aaron Meurer]

Here is the Python implementation
https://github.com/python/cpython/blob/a10d426bab66a4e1f20d5e1b9aee3dbb435cf309/Objects/floatobject.c#L917

It's possible you have found a bug in Python's round. I'm unclear how
it is supposed to work.

Aaron Meurer

[Aaron Meurer]

The round() function used in that function comes from the C math
standard library. Here is an implementation that is used if it isn't
available. https://github.com/python/cpython/blob/a10d426bab66a4e1f20d5e1b9aee3dbb435cf309/Python/pymath.c#L72

Aaron Meurer

[Chris Smith]

That floats are stored in binary is an implementation detail which need not prevent base-10 rounding to still work. The 2nd argument to round is intended to tell at which base-10 digit the rounding is to take place. By shifting that digit to the ones position and then doing the rounding one can easily detect whether what follows is above or below a half.

>>> .345.as_integer_ratio()
(1553741871442821, 4503599627370496)

>>> (_*100).as_integer_ratio()
(69, 2)

"round to even" means "if what follows the digit to which you are rounding is 5 then round so your digit is even". As is shown in the integer ratios above, what follows the 4 is exactly 1/2 so we can round to the even 4 (as in .34) instead of the odd 5 (as in .35) just like round(0.75,1) rounds to 0.8 while round(0.25) rounds to 0.2.

On Wednesday, April 10, 2019 at 7:29:18 PM UTC-5, Aaron Meurer wrote:

Doesn't Python do rounding based on the binary representation of the float?

I'm a little confused what "round to even" means in that case.

Aaron Meurer

On Wed, Apr 10, 2019 at 6:16 PM Chris Smith <smi...@gmail.com> wrote:
>
> Python 3 implements "round to even on tie" logic so `round(12.5)` -> 12,  not 13. I have updated, in #16608, the round function but there is a difference in how ties are detected. I shift the desired position to the ones position and then check for a tie so 12.345 is shifted to 1234.5 and rounded to 1234 then is divided by 100 to give 12.34. Python doesn't do this. I suspect it adds 0.05 and then detects that12.395 > 12395/1000 and rounds up to 12.35
>
>
> >>> Rational(*.345.as_integer_ratio())-Rational(345,1000)  # .345 < 345/1000
> -3/112589990684262400
> >>> Rational(*.395.as_integer_ratio())-Rational(395,1000)  # .395 > 395/1000
> 1/56294995342131200
>
>
> >>> round(12.345,2)  # 12.345 rounds up b/c a tie is not detected
> 12.35
>
>
>
> Does anyone have objections to the proposed rounding?
>
> /c
>
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to sy...@googlegroups.com.

[Oscar Benjamin]

The fact that the numbers are stored in binary is significant:

In [16]: nums = [eval('1.%d5' % n) for n in range(10)]

In [17]: nums
Out[17]: [1.05, 1.15, 1.25, 1.35, 1.45, 1.55, 1.65, 1.75, 1.85, 1.95]

In [18]: [round(n, 1) for n in nums]
Out[18]: [1.1, 1.1, 1.2, 1.4, 1.4, 1.6, 1.6, 1.8, 1.9, 1.9]

Proper decimal rounding might look like:

In [20]: from decimal import Decimal

In [21]: nums = [Decimal('1.%d5' % n) for n in range(10)]

In [22]: nums
Out[22]:
[Decimal('1.05'),
Decimal('1.15'),
Decimal('1.25'),
Decimal('1.35'),
Decimal('1.45'),
Decimal('1.55'),
Decimal('1.65'),
Decimal('1.75'),
Decimal('1.85'),
Decimal('1.95')

In [23]: [round(n, 1) for n in nums]
Out[23]:
[Decimal('1.1'),
Decimal('1.2'),
Decimal('1.3'),
Decimal('1.4'),
Decimal('1.5'),
Decimal('1.6'),
Decimal('1.7'),
Decimal('1.8'),
Decimal('1.9'),
Decimal('2.0')]

Or with half-even rounding:

In [24]: from decimal import getcontext

In [25]: getcontext().rounding = 'ROUND_HALF_EVEN'

In [26]: [round(n, 1) for n in nums]
Out[26]:
[Decimal('1.0'),
Decimal('1.2'),
Decimal('1.2'),
Decimal('1.4'),
Decimal('1.4'),
Decimal('1.6'),
Decimal('1.6'),
Decimal('1.8'),
Decimal('1.8'),
Decimal('2.0')]

The binary floats don't work out correct because some are above and
some are below the number suggested by the original float literal.

> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/7873bfe6-6f91-4d39-af6f-9f4039714fa8%40googlegroups.com.

[Chris Smith]

I am aware of how the numbers are stored but was overly optimistic that the shift could resolve this in all cases. It can't (and thanks for the correction). 

But my suggested alternative makes a significant difference in how often the problem arises:

>>> bad=[]
>>> for i in range(1,1000):
...  n = str(i)+'5'
...  if int(str(round(int(n)/10**len(n),len(n)-1))[-1])%2!=0:bad.append(i)  # e.g. round(0.1235, 3)
...
>>> len(bad)
546


>>> bad=[]
>>> for i in range(1,1000):
...  n = str(i)+'5'
...  if round(int(n)/10**len(n)*10**(len(n)-1))%2!=0:bad.append(i)  # e.g. round(0.1235*1000)
...
>>> len(bad)
8
>>> bad  # e.g. 0.545*100 != 54.5
[54, 57, 501, 503, 505, 507, 509, 511]

So the question is whether we want to do better and keep the SymPy algorithm.

[Oscar Benjamin]

I think that Python's float.__round__ is correct. AIUI it rounds
correctly based on the true value represented by the float:

In [4]: round(1.05, 1)
Out[4]: 1.1

In [5]: import decimal

In [6]: decimal.Decimal(1.1)
Out[6]: Decimal('1.100000000000000088817841970012523233890533447265625')

That's surprising just because the true exact value isn't what we
thought it should be. There is a rounding error but it happens in the
literal 1.1 that doesn't give us the float we think it should. It's
extra confusing because str(1.1) gives '1.1' making it seem like the
exact number we requested. The behaviour of round is mathematically
correct here though given the input received: it is computing the true
exact result correctly rounded (where "correctly rounded" refers to
the decimal to binary rounding at the end).

For decimals that can be exactly represented in binary the results are
as we would intuitively expect:

In [7]: round(1.25, 1)
Out[7]: 1.2

In [8]: round(1.75, 1)
Out[8]: 1.8

(Of course 1.2 and 1.8 are not exactly representable in binary
floating point but we have the floats that result from correctly
rounding the true numbers 1.2 and 1.8.)

--
Oscar

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/b50c0504-ec96-4583-84b1-d5a3902277c2%40googlegroups.com.

[Chris Smith]

The question is whether we can, with limitations of binary representation, give a result that is consistent with what we would expect if using base-10 notation. The advantage that SymPy and Decimal have is that they know, by merit of values given at instantiation, what the last significant digit is. So if we round to that digit and use only those digits when making the decision when rounding to fewer places we can get it right. e.g. 0.345 is stored as something like '0.344999999999999' If we know we are working with a precision of 3 then we can multiply by 1000, round to int and now use that number (345) to round to the second digit and use tie-breaking logic to give the result as 340 instead of 350 and return a final result of 0.340.

 >>> def r(n,p,prec):
 ... if p > prec:
 ... return n
 ... i = int(round(n*10**prec))
 ... m = 10**(prec-p)
 ... i, r = divmod(i, m)
 ... if i%2 and 2*r==m:
 ...     i -= 1
 ... return i/float(10**p)
 ...
 >>> bad=[]
 ... for i in range(1,1000):

 ...   n = str(i)+'5'

 ... if int(str(r(eval('.'+n),len(n)-1,len(n)))[:len(n)+1][-1])%2!=0:bad.append(i)
 ...
 >>> len(bad)
 0
 >>> r(.545, 2, 3)
 0.540000000000000
 >>> round(.545,2)
 0.55

So I guess the answer to the original question is that it's ok for Decimal or SymPy to give a different result than Python when dealing with non-native number objects.

I will double check that SymPy's round is doing what it should and not generating any "incorrect/bad" results.

/c

> To unsubscribe from this group and stop receiving emails from it, send an email to sy...@googlegroups.com.