Re: ANN: BigDecimal - decimal arithmetic on very large intergers

Discussion in 'Python' started by M.-A. Lemburg, Apr 1, 2005.

  1. wrote:
    > BigDecimal is a Python class that supports decimal arithmetic on very large integers. BigDecimal was inspired by the posting of BigDec to c.l.py by Tim Peters. BigDecimal implements all the commonly used integer methods. (It doesn't implement any of the binary/shifting operations.)
    >
    > It has been optimized for performance. It uses a 4x4 Toom-Cook algorithm for multiplication and a new, very fast, division algorithm. If GMPY is available, it will be automatically used.
    >
    > Performance examples, computing the decimal represendation of the 42nd Mersenne prime:
    > 2**25964951 - 1
    >
    > Tim Peter's posting to c.l.py: 13 minutes 41 seconds
    > BigDecimal: 59 seconds
    > BigDecimal w/gmpy: 10 seconds


    You might want to look into mxNumber (part of the egenix-mx-experimental
    package):

    http://www.egenix.com/files/python/mxNumber.html

    There's a C type called Integer in that package that basically
    wraps the GMP integer type.

    --
    Marc-Andre Lemburg
    eGenix.com

    Professional Python Services directly from the Source (#1, Apr 01 2005)
    >>> Python/Zope Consulting and Support ... http://www.egenix.com/
    >>> mxODBC.Zope.Database.Adapter ... http://zope.egenix.com/
    >>> mxODBC, mxDateTime, mxTextTools ... http://python.egenix.com/

    ________________________________________________________________________

    ::: Try mxODBC.Zope.DA for Windows,Linux,Solaris,FreeBSD for free ! ::::
    M.-A. Lemburg, Apr 1, 2005
    #1
    1. Advertising

  2. M.-A. Lemburg

    Guest

    M.-A. Lemburg wrote:
    > wrote:
    > > BigDecimal is a Python class that supports decimal arithmetic on

    very large integers. BigDecimal was inspired by the posting of BigDec
    to c.l.py by Tim Peters. BigDecimal implements all the commonly used
    integer methods. (It doesn't implement any of the binary/shifting
    operations.)
    > >
    > > It has been optimized for performance. It uses a 4x4 Toom-Cook

    algorithm for multiplication and a new, very fast, division algorithm.
    If GMPY is available, it will be automatically used.
    > >
    > > Performance examples, computing the decimal represendation of the

    42nd Mersenne prime:
    > > 2**25964951 - 1
    > >
    > > Tim Peter's posting to c.l.py: 13 minutes 41 seconds
    > > BigDecimal: 59 seconds
    > > BigDecimal w/gmpy: 10 seconds

    >
    > You might want to look into mxNumber (part of the

    egenix-mx-experimental
    > package):
    >
    > http://www.egenix.com/files/python/mxNumber.html
    >
    > There's a C type called Integer in that package that basically
    > wraps the GMP integer type.
    >
    > --
    > Marc-Andre Lemburg
    > eGenix.com
    >
    > Professional Python Services directly from the Source (#1, Apr 01

    2005)
    > >>> Python/Zope Consulting and Support ...

    http://www.egenix.com/
    > >>> mxODBC.Zope.Database.Adapter ...

    http://zope.egenix.com/
    > >>> mxODBC, mxDateTime, mxTextTools ...

    http://python.egenix.com/
    >

    ________________________________________________________________________
    >
    > ::: Try mxODBC.Zope.DA for Windows,Linux,Solaris,FreeBSD for free !

    ::::

    I have used mxNumber in the past. This library uses "super" digits with
    hundreds / thousands of decimal digits and then builds builds
    multiplication and division on those. The original motivation was that
    the conversion time to / from decimal string format is O(n) instead of
    O(n^2).

    Using just native Python long support, the 4-way Toom-Cook
    multiplication algorithm is faster than the built-in multiplication
    when the numbers are several hundred thousand digits long. On my
    machine, it is roughly twice as fast when multiplying one million digit
    numbers. (The 4-way Toom-Cook is O(n^~1.4) versus O(n^~1.585) for the
    Karatsuba multiplication in Python.)

    The division algortihm in BigDecimal is effectively O(n^~1.4) also.
    Using just native Python long support, the division algorithm is faster
    than the built-in division algorithm when the numbers are several tens
    of thousands digits long.

    Interestingly, BigDecimal can do division faster than GMP 3.1.x with
    numbers approximately 10 million digits in length. BigDecimal is faster
    than GMP 4.1.4 with numbers of approximately 1 million digits in
    length. (GMP 4 is faster for small, ~10,000 digits, than GMP 3, but
    grows more quickly.)

    casevh
    , Apr 1, 2005
    #2
    1. Advertising

  3. M.-A. Lemburg

    Kay Schluehr Guest

    Some notes about float approximations in mxNumber

    Hi Mark,

    I was a bit surprised to find the very slow Farey approximation by
    means of the <FareyRational> class in the mxNumber package. If the goal
    was to reconstruct a rational from a float it is not a good choice and
    should be replaced by a continued fractions approximation. Some time
    ago I implemented it by myself so it can be published here:


    def cfrac(z,bound=10**-5,n=10):
    '''
    creates a continued fraction from some number z.
    @ bound - terminate cf if the rest is lower than the provided bound
    @ n - terminate cf after n steps
    '''
    l = []
    while 1:
    a = int(z)
    l.append(a)
    y = z-a
    if y<=bound or len(l)==n:
    return l
    z = 1./y


    def fold(cf):
    '''
    create u/v = cfrac(a0,a1,...,an) using following rules:

    1. / u1 u0 \ / a1 1 \ / a0 1 \
    | | = | | * | |
    \ v1 v0 / \ 1 0 / \ 1 0 /

    2. The recursion rules
    v(n+1) = v(n)*a(n+1)+v(n-1)
    u(n+1) = u(n)*a(n+1)+u(n-1)
    '''
    if len(cf)<2:
    return Rational(0,1)
    un = cf[0]*cf[1]+1
    vn = cf[1]
    un_1 = cf[0]
    vn_1 = 1
    for a in cf[2:]:
    b = un
    un = un*a+un_1
    un_1 = b
    b = vn
    vn = vn*a+vn_1
    vn_1 = b
    return Rational(un,vn)


    >>> fold(cfrac(1./3))

    1/3
    >>> fold(cfract(1525/42.))

    1525/42

    import math

    >>> fold(cfract(math.sqrt(2)))

    3363/2378


    It is possible to provide this functionality in an even more efficient
    manner because it is usefull to bound only the approximation error of
    the rational, not the size of the continued fraction.

    def float2ratio(z, bound=10**-5):
    '''
    convert a float into a Rational.
    The approximation is bound by the error-limit 'bound'.
    '''
    a = int(z)
    y = z-a
    z = 1./y
    b = int(z)
    un = a*b+1
    vn = b
    un_1 = a
    vn_1 = 1
    a = b
    while 1:
    y = z-a
    if y<bound:
    return Rational(un,vn),k
    z = 1./y
    a = int(z)
    x = un
    un = un*a+un_1
    un_1 = x
    x = vn
    vn = vn*a+vn_1
    vn_1 = x
    xn = float(un)/vn
    yn = float(un_1)/vn_1
    if abs(xn-yn)<=bound:
    return Rational(un,vn)

    >>> math.sqrt(2)

    1.4142135623730951

    >>> float2ratio(math.sqrt(2))

    1393/985
    >>> 1393./985

    1.4142131979695431
    ^

    >>> float2ratio(math.sqrt(2),10**-10)

    275807/195025
    >>> 275807./195025

    1.4142135623637995
    ^

    >>> math.pi

    3.1415926535897931

    >>> float2ratio(math.pi,10**-14)

    245850922/78256779
    >>> 245850922./78256779

    3.1415926535897931

    Note that the algorithm needed only 13 iteration steps to approximate
    pi
    in this accuracy.

    Regards,
    Kay
    Kay Schluehr, Apr 1, 2005
    #3
  4. Re: Some notes about float approximations in mxNumber

    Kay Schluehr wrote:
    > Hi Marc,
    >
    > I was a bit surprised to find the very slow Farey approximation by
    > means of the <FareyRational> class in the mxNumber package. If the goal
    > was to reconstruct a rational from a float it is not a good choice and
    > should be replaced by a continued fractions approximation.


    The idea was to be able to create a Rational() from a float
    using a given upper bound on the denominator.

    The standard Rational() constructor already provides a way
    to construct a rational out of a Python or mxNumber float,
    but always uses maximum precision - which may not always be
    what the user wants (e.g. to work around rounding errors).

    Thanks for posting your faster version. I think it would be
    a good candidate for a Python Cookbook Entry and I'll
    see whether I can add something like this to one of the next
    mxNumber releases.

    > Some time
    > ago I implemented it by myself so it can be published here:
    >
    >
    > def cfrac(z,bound=10**-5,n=10):
    > '''
    > creates a continued fraction from some number z.
    > @ bound - terminate cf if the rest is lower than the provided bound
    > @ n - terminate cf after n steps
    > '''
    > l = []
    > while 1:
    > a = int(z)
    > l.append(a)
    > y = z-a
    > if y<=bound or len(l)==n:
    > return l
    > z = 1./y
    >
    >
    > def fold(cf):
    > '''
    > create u/v = cfrac(a0,a1,...,an) using following rules:
    >
    > 1. / u1 u0 \ / a1 1 \ / a0 1 \
    > | | = | | * | |
    > \ v1 v0 / \ 1 0 / \ 1 0 /
    >
    > 2. The recursion rules
    > v(n+1) = v(n)*a(n+1)+v(n-1)
    > u(n+1) = u(n)*a(n+1)+u(n-1)
    > '''
    > if len(cf)<2:
    > return Rational(0,1)
    > un = cf[0]*cf[1]+1
    > vn = cf[1]
    > un_1 = cf[0]
    > vn_1 = 1
    > for a in cf[2:]:
    > b = un
    > un = un*a+un_1
    > un_1 = b
    > b = vn
    > vn = vn*a+vn_1
    > vn_1 = b
    > return Rational(un,vn)
    >
    >
    >
    >>>>fold(cfrac(1./3))

    >
    > 1/3
    >
    >>>>fold(cfract(1525/42.))

    >
    > 1525/42
    >
    > import math
    >
    >
    >>>>fold(cfract(math.sqrt(2)))

    >
    > 3363/2378
    >
    >
    > It is possible to provide this functionality in an even more efficient
    > manner because it is usefull to bound only the approximation error of
    > the rational, not the size of the continued fraction.
    >
    > def float2ratio(z, bound=10**-5):
    > '''
    > convert a float into a Rational.
    > The approximation is bound by the error-limit 'bound'.
    > '''
    > a = int(z)
    > y = z-a
    > z = 1./y
    > b = int(z)
    > un = a*b+1
    > vn = b
    > un_1 = a
    > vn_1 = 1
    > a = b
    > while 1:
    > y = z-a
    > if y<bound:
    > return Rational(un,vn),k
    > z = 1./y
    > a = int(z)
    > x = un
    > un = un*a+un_1
    > un_1 = x
    > x = vn
    > vn = vn*a+vn_1
    > vn_1 = x
    > xn = float(un)/vn
    > yn = float(un_1)/vn_1
    > if abs(xn-yn)<=bound:
    > return Rational(un,vn)
    >
    >
    >>>>math.sqrt(2)

    >
    > 1.4142135623730951
    >
    >
    >>>>float2ratio(math.sqrt(2))

    >
    > 1393/985
    >
    >>>>1393./985

    >
    > 1.4142131979695431
    > ^
    >
    >
    >>>>float2ratio(math.sqrt(2),10**-10)

    >
    > 275807/195025
    >
    >>>>275807./195025

    >
    > 1.4142135623637995
    > ^
    >
    >
    >>>>math.pi

    >
    > 3.1415926535897931
    >
    >
    >>>>float2ratio(math.pi,10**-14)

    >
    > 245850922/78256779
    >
    >>>>245850922./78256779

    >
    > 3.1415926535897931
    >
    > Note that the algorithm needed only 13 iteration steps to approximate
    > pi
    > in this accuracy.
    >
    > Regards,
    > Kay
    >


    --
    Marc-Andre Lemburg
    eGenix.com

    Professional Python Services directly from the Source (#1, Apr 04 2005)
    >>> Python/Zope Consulting and Support ... http://www.egenix.com/
    >>> mxODBC.Zope.Database.Adapter ... http://zope.egenix.com/
    >>> mxODBC, mxDateTime, mxTextTools ... http://python.egenix.com/

    ________________________________________________________________________

    ::: Try mxODBC.Zope.DA for Windows,Linux,Solaris,FreeBSD for free ! ::::
    M.-A. Lemburg, Apr 4, 2005
    #4
    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. Replies:
    0
    Views:
    292
  2. MLyngsie
    Replies:
    2
    Views:
    862
    MLyngsie
    Oct 26, 2006
  3. Replies:
    3
    Views:
    1,220
  4. Stanimir Stamenkov
    Replies:
    4
    Views:
    2,578
    Stanimir Stamenkov
    Jul 18, 2008
  5. Ramon
    Replies:
    12
    Views:
    1,341
Loading...

Share This Page