# Decimal and Exponentiation

Discussion in 'Python' started by elventear, May 19, 2006.

1. ### elventearGuest

Hi,

I am the in the need to do some numerical calculations that involve
real numbers that are larger than what the native float can handle.

I've tried to use Decimal, but I've found one main obstacle that I
don't know how to sort. I need to do exponentiation with real
exponents, but it seems that Decimal does not support non integer
exponents.

I would appreciate if anyone could recommend a solution for this
problem.

Thank you.

elventear, May 19, 2006

2. ### Tim PetersGuest

[elventear]
> I am the in the need to do some numerical calculations that involve
> real numbers that are larger than what the native float can handle.
>
> I've tried to use Decimal, but I've found one main obstacle that I
> don't know how to sort. I need to do exponentiation with real
> exponents, but it seems that Decimal does not support non integer
> exponents.
>
> I would appreciate if anyone could recommend a solution for this
> problem.

Wait <0.3 wink>. Python's Decimal module intends to be a faithful
implementation of IBM's proposed standard for decimal arithmetic:

http://www2.hursley.ibm.com/decimal/

Last December, ln, log10, exp, and exponentiation to non-integral
powers were added to the proposed standard, but nobody yet has written
implementation code for Python's module. [Python-Dev: somebody
wants to volunteer for this ]

If you're not a numeric expert, I wouldn't recommend that you try this
yourself (in particular, trying to implement x**y as exp(ln(x)*y)
using the same precision is mathematically correct but is numerically

The GNU GMP library (for which Python bindings are available) also
supports "big floats", but their power operation is also restricted to
integer powers and/or exact roots. This can be painful even to try;
e.g.,

>>> from gmpy import mpf
>>> mpf("1e10000") ** mpf("3.01")

consumed well over a minute of CPU time (on a 3.4 GHz box) before dying with

ValueError: mpq.pow fractional exponent, inexact-root

If you're working with floats outside the range of IEEE double, you
_probably_ want to be working with logarithms instead anyway; but that
depends on you app, and I don't want to know about it ;-)

Tim Peters, May 19, 2006

3. ### Dan BishopGuest

Tim Peters wrote:
....
> Wait <0.3 wink>. Python's Decimal module intends to be a faithful
> implementation of IBM's proposed standard for decimal arithmetic:
>
> http://www2.hursley.ibm.com/decimal/
>
> Last December, ln, log10, exp, and exponentiation to non-integral
> powers were added to the proposed standard, but nobody yet has written
> implementation code for Python's module. [Python-Dev: somebody
> wants to volunteer for this ]

Here's a quick-and-dirty exp function:

def exp(x):
"""
Return e raised to the power of x.
"""
if x < 0:
return 1 / exp(-x)
partial_sum = term = 1
i = 1
while True:
term *= x / i
new_sum = partial_sum + term
if new_sum == partial_sum:
return new_sum
partial_sum = new_sum
i += 1

Dan Bishop, May 20, 2006
4. ### Raymond L. BuvelGuest

elventear wrote:
> Hi,
>
> I am the in the need to do some numerical calculations that involve
> real numbers that are larger than what the native float can handle.
>
> I've tried to use Decimal, but I've found one main obstacle that I
> don't know how to sort. I need to do exponentiation with real
> exponents, but it seems that Decimal does not support non integer
> exponents.
>
> I would appreciate if anyone could recommend a solution for this
> problem.
>
> Thank you.
>

The clnum module has arbitrary precision floating point and complex
numbers with all of the standard math functions. For example, the cube
root of 2 can be computed to 40 decimal places with the following.

>>> from clnum import mpf,mpq
>>> mpf(2,40)**mpq(1,3)

mpf('1.259921049894873164767210607278228350570251464701',46)

http://calcrpnpy.sourceforge.net/clnumManual.html

Raymond L. Buvel, May 20, 2006
5. ### Raymond L. BuvelGuest

Tim Peters wrote:
<snip>

> The GNU GMP library (for which Python bindings are available) also
> supports "big floats", but their power operation is also restricted to
> integer powers and/or exact roots. This can be painful even to try;
> e.g.,
>
> >>> from gmpy import mpf
> >>> mpf("1e10000") ** mpf("3.01")

>
> consumed well over a minute of CPU time (on a 3.4 GHz box) before dying
> with
>
> ValueError: mpq.pow fractional exponent, inexact-root
>

<snip>

The clnum module handles this calculation very quickly:

>>> from clnum import mpf
>>> mpf("1e10000") ** mpf("3.01")

mpf('9.9999999999999999999999932861e30099',26)
>>> x=_
>>> x ** (1/mpf("3.01"))

mpf('9.9999999999999999999999953924e9999',26)

See http://calcrpnpy.sourceforge.net/clnumManual.html

Raymond L. Buvel, May 20, 2006
6. ### Tim PetersGuest

[Raymond L. Buvel, on
http://calcrpnpy.sourceforge.net/clnumManual.html
]
> The clnum module handles this calculation very quickly:
>
> >>> from clnum import mpf
> >>> mpf("1e10000") ** mpf("3.01")

> mpf('9.9999999999999999999999932861e30099',26)

That's probably good enough for the OP's needs -- thanks!

OTOH, it's not good enough for the decimal module:

(10**10000)**3.01 =
10**(10000*3.01) =
10**30100

exactly, and the proposed IBM standard for decimal arithmetic requires
error < 1 ULP (which implies that if the mathematical ("infinite
precision") result is exactly representable, then that's the result
you have to get). It would take some analysis to figure out how much
of clnum's error is due to using binary floats instead of decimal, and
how much due to its pow implementation.

Tim Peters, May 20, 2006
7. ### Raymond L. BuvelGuest

Tim Peters wrote:
> [Raymond L. Buvel, on
> http://calcrpnpy.sourceforge.net/clnumManual.html
> ]
>
>> The clnum module handles this calculation very quickly:
>>
>> >>> from clnum import mpf
>> >>> mpf("1e10000") ** mpf("3.01")

>> mpf('9.9999999999999999999999932861e30099',26)

>
>
> That's probably good enough for the OP's needs -- thanks!
>
> OTOH, it's not good enough for the decimal module:
>
> (10**10000)**3.01 =
> 10**(10000*3.01) =
> 10**30100
>
> exactly, and the proposed IBM standard for decimal arithmetic requires
> error < 1 ULP (which implies that if the mathematical ("infinite
> precision") result is exactly representable, then that's the result
> you have to get). It would take some analysis to figure out how much
> of clnum's error is due to using binary floats instead of decimal, and
> how much due to its pow implementation.

Indeed, it is not clear where the error is comming from especially since
you can increase the precision of the intermediate calculation and get
the exact result.

>>> mpf(mpf("1e10000",30) ** mpf("3.01",30), 20)

mpf('1.0e30100',26)

Is this the kind of thing you will need to do in the Decimal module to
meet the specification?

Raymond L. Buvel, May 20, 2006
8. ### Tim PetersGuest

[Raymond L. Buvel, on
http://calcrpnpy.sourceforge.net/clnumManual.html
]
>>> The clnum module handles this calculation very quickly:
>>>
>>> >>> from clnum import mpf
>>> >>> mpf("1e10000") ** mpf("3.01")
>>> mpf('9.9999999999999999999999932861e30099',26)

[Tim Peters]
>> That's probably good enough for the OP's needs -- thanks!
>>
>> OTOH, it's not good enough for the decimal module:
>>
>> (10**10000)**3.01 =
>> 10**(10000*3.01) =
>> 10**30100
>>
>> exactly, and the proposed IBM standard for decimal arithmetic requires
>> error < 1 ULP (which implies that if the mathematical ("infinite
>> precision") result is exactly representable, then that's the result
>> you have to get). It would take some analysis to figure out how much
>> of clnum's error is due to using binary floats instead of decimal, and
>> how much due to its pow implementation.

[Raymond]
> Indeed, it is not clear where the error is comming from especially since
> you can increase the precision of the intermediate calculation and get
> the exact result.
>
> >>> mpf(mpf("1e10000",30) ** mpf("3.01",30), 20)

> mpf('1.0e30100',26)
>
> Is this the kind of thing you will need to do in the Decimal module to
> meet the specification?

It will vary by function and the amount of effort people are willing
to put into implementations. Temporarily (inside a function's
implementation) increasing working precision is probably the easiest
way to get results provably suffering less than 1 ULP error in the
destination precision. The "provably" part is the hardest part under
any approach ;-)

Tim Peters, May 20, 2006