Math errors in python

[Bengt Richter]

On Sun, 19 Sep 2004 14:41:53 +0200, Peter Otten <__peter__ at web.de> wrote:

>Paul Rubin wrote:
>
>> Peter Otten <__peter__ at web.de> writes:
>>> Starting with Python 2.4 there will be the 'decimal' module supporting
>>> "arithmetic the way you know it":
>>> 
>>> >>> from decimal import *
>>> >>> Decimal("12.10") + Decimal("8.30")
>> 
>> I haven't tried 2.4 yet.  After
>
>The auther is currently working on an installer, but just dropping it into
>2.3's site-packages should work, too.
>
>>     a = Decimal("1") / Decimal("3")
>>     b = a * Decimal("3")
>>     print b
>> 
>> What happens?  Is that arithmetic as the way I know it?
>
>Decimal as opposed to rational:
>
>>>> from decimal import *
>>>> Decimal(1)/Decimal(3)
>Decimal("0.3333333333333333333333333333")
>>>> 3*_
>Decimal("0.9999999999999999999999999999")
>
>Many people can cope with the inaccuracy induced by base 10 representations
>and are taken by surprise by base 2 errors.
>But you are right I left too much room for interpretation.
>
I hacked a little rational + decimal exponent representation based toy a while
back. The original post had a bug, which someone pointed out and I posted a
followup fix for, but the revised version was not posted. But I can if someone
is interested.

 >>> from ut.exactdec import ED
 >>> ED(1)/ED(3)
 ED('1 / 3')
 >>> 3*_
 ED('1')

If you give it a float, it wants to know how many decimals you mean:
 >>> ED(1./3)
 Traceback (most recent call last):
   File "<stdin>", line 1, in ?
   File "c:\pywk\ut\exactdec.py", line 93, in __init__
     raise ValueError(
 ValueError: Specify decimals for least significant digit of 10**(-decimals)
 (decimals may also be specified as 'all' to capture all bits of float)

 >>> ED(1./3, 'all')
 ED('0.333333333333333314829616256247390992939472198486328125')

If you give it a string literal, it takes it as accurate, but you can round it
to create a new accurate number:
 >>> ED('1/3', 54)
 ED('0.333333333333333333333333333333333333333333333333333333')
 >>> ED('1/3', 60)
 ED('0.333333333333333333333333333333333333333333333333333333333333')

That's an accurate number that has all zeroes to the right of those 60 3's
 >>> ED('1/3', 60)*3
 ED('0.999999999999999999999999999999999999999999999999999999999999')

If you don't round, you get a fully accurate result"
 >>> ED('1/3')*3
 ED('1')

It's interesting to look at pi:
 >>> import math
 >>> math.pi
 3.1415926535897931
 >>> ED(math.pi, 'all')
 ED('3.141592653589793115997963468544185161590576171875')
 >>> ED(3.1415926535897931, 'all')
 ED('3.141592653589793115997963468544185161590576171875')

Same actual exact decimal value gets created from 
 >>>
 >>> repr(math.pi)
 '3.1415926535897931'

meaning they both have the same floating point hardware representation,
but the short version decimal literal is sufficient to set all the bits right
even though it doesn't represent the fully exact value in decimal.
Economy courtesy of the Timbot I think ;-)

I don't know what the rules in Decimal are for stage-wise rounding vs keeping
accuracy, but I imagine you could get the same kind of surprises that are
available in binary from floating point, e.g., 

 >>> from ut.exactdec import ED

Floating point:
 >>> acc = 1.0
 >>> for i in xrange(100): acc += 1e-300
 ...
 >>> acc
 1.0

That really is exactly 1.0
 >>> ED(acc,'all')

Now the calculation accurately:
 ED('1')
 >>> ecc = ED(1)
 >>> for i in xrange(100): ecc += ED('1e-300')
 ...
 >>> ecc
 ED('1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
 0000000000000001')
 >>> ecc-1
 ED('1.0e-298')
 
If you add a small Decimal delta repeatedly, will it get rounded away like the floating point
version, or will accuracy get promoted, or what? Sorry, I haven't read the docs yet ;-/

Regards,
Bengt Richter