Decimal and Exponentiation

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

  1. elventear

    elventear Guest

    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
    #1
    1. Advertising

  2. elventear

    Tim Peters Guest

    [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 ;-)
     
    Tim Peters, May 19, 2006
    #2
    1. Advertising

  3. elventear

    Dan Bishop Guest

    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
    #3
  4. 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)

    For more information see

    http://calcrpnpy.sourceforge.net/clnumManual.html
     
    Raymond L. Buvel, May 20, 2006
    #4
  5. 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
    #5
  6. elventear

    Tim Peters Guest

    [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
    #6
  7. 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
    #7
  8. elventear

    Tim Peters Guest

    [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
    #8
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Steven T. Hatton

    An exponentiation function for int?

    Steven T. Hatton, Oct 13, 2004, in forum: C++
    Replies:
    14
    Views:
    741
    JXStern
    Oct 16, 2004
  2. Jeff Davis

    strange exponentiation performance

    Jeff Davis, Nov 23, 2003, in forum: Python
    Replies:
    0
    Views:
    340
    Jeff Davis
    Nov 23, 2003
  3. Tim Peters
    Replies:
    1
    Views:
    432
    Jeff Davis
    Nov 24, 2003
  4. exponentiation operator (lack of)

    , Dec 22, 2005, in forum: C Programming
    Replies:
    67
    Views:
    1,508
    Dave Thompson
    Jan 4, 2006
  5. Andrew Savige
    Replies:
    2
    Views:
    241
    Andrew Savige
    Mar 26, 2009
Loading...

Share This Page