Velocity Reviews > Re: ANN: BigDecimal - decimal arithmetic on very large intergers

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

M.-A. Lemburg
Guest
Posts: n/a

 04-01-2005
http://www.velocityreviews.com/forums/(E-Mail Removed) 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, mxDateTime, mxTextTools ... http://python.egenix.com/

__________________________________________________ ______________________

casevh@comcast.net
Guest
Posts: n/a

 04-01-2005
M.-A. Lemburg wrote:
> (E-Mail Removed) 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/

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

http://python.egenix.com/
>

__________________________________________________ ______________________
>

::::

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

Kay Schluehr
Guest
Posts: n/a

 04-01-2005
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

M.-A. Lemburg
Guest
Posts: n/a

 04-04-2005
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, mxDateTime, mxTextTools ... http://python.egenix.com/

__________________________________________________ ______________________