Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > Python > Decimal and Exponentiation

Reply
Thread Tools

Decimal and Exponentiation

 
 
elventear
Guest
Posts: n/a
 
      05-19-2006
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.

 
Reply With Quote
 
 
 
 
Tim Peters
Guest
Posts: n/a
 
      05-19-2006
[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
 
Reply With Quote
 
 
 
 
Dan Bishop
Guest
Posts: n/a
 
      05-20-2006
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

 
Reply With Quote
 
Raymond L. Buvel
Guest
Posts: n/a
 
      05-20-2006
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.2599210498948731647672106072782283505702514 64701',46)

For more information see

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

 
Reply With Quote
 
Raymond L. Buvel
Guest
Posts: n/a
 
      05-20-2006
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
 
Reply With Quote
 
Tim Peters
Guest
Posts: n/a
 
      05-20-2006
[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.
 
Reply With Quote
 
Raymond L. Buvel
Guest
Posts: n/a
 
      05-20-2006
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?
 
Reply With Quote
 
Tim Peters
Guest
Posts: n/a
 
      05-20-2006
[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
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
Operator precedence of exponentiation (**) and complement (~) Andrew Savige Ruby 2 03-26-2009 08:03 PM
exponentiation operator (lack of) carlos@colorado.edu C Programming 67 01-04-2006 05:27 AM
An exponentiation function for int? Steven T. Hatton C++ 14 10-16-2004 12:23 AM
RE: strange exponentiation performance Tim Peters Python 1 11-24-2003 06:35 AM
strange exponentiation performance Jeff Davis Python 0 11-23-2003 11:31 AM



Advertisments