mpmath puzzle

[Fredrik Johansson]

Dick Moores wrote:

> For 1234 ** 10.9, why the wrong result from mpmath.power()?
>
> ========================================
> #!/usr/bin/env python
> #coding=utf-8
> from mpmath import *
>
> mpf.dps = 32
>
> x = mpf(1234)
> y = mpf(10.9)
>
> print power(x,y)
> print "4.9583278648155041477415234438717e+33" # from Windows calculator
>
> """
> output:
> 4.9583278648155166864966558721921e+33
> 4.9583278648155041477415234438717e+33
> """
> ========================================
>  (Code is also at <http://python.pastebin.com/m72a277b8>)
>
>  Thanks,
>
>  Dick Moores

Hi,

When you create y, you first create an inexact Python float whose
value is actually 10.900000000000000355... and then pass it to mpf,
which just copies that value along with the error. To create an
accurate representation of a fractional number (integers are safe),
pass a string instead:

>>> from mpmath import *
>>>
>>> mpf.dps = 32
>>> x = mpf(1234)
>>> y = mpf('10.9')
>>> print power(x,y)
4.9583278648155041477415234438719e+33

The last printed digit is wrong, which is to be expected in binary
arithmetic since 10.9 still cannot be represented exactly (with
mpf.dps = 32, the error is about 5e-033); computing a power amplifies
that error. You can temporarily increase the precision by a couple of
digits during calculations and then reduce it before printing the
final value if you want all digits to be right.

I should update the mpmath documentation to be much clearer about all this.

Fredrik