Decimal and Exponentiation

[elventear]

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.

[Tim Peters]

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

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

[Dan Bishop]

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

[Raymond L. Buvel]

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.25992104 989487316476721 060727822835057 0251464701',46)

For more information see

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

[Raymond L. Buvel]

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.99999999 999999999999999 32861e30099',26 ) x=_
x ** (1/mpf("3.01"))

mpf('9.99999999 999999999999999 53924e9999',26)

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

[Tim Peters]

[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.99999999 999999999999999 32861e30099',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.

[Raymond L. Buvel]

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.99999999 999999999999999 32861e30099',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("1e1000 0",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?

[Tim Peters]

[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.99999999 999999999999999 32861e30099',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("1e1000 0",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 ;-)