which pi formula is given in the decimal module documentation?
Anh Hai Trinh
docs.python.org/library/decimal.html#recipes>?
def pi():
"""Compute Pi to the current precision.
>>> print pi()
3.141592653589793238462643383
"""
getcontext().prec += 2 # extra digits for intermediate steps
three = Decimal(3) # substitute "three=3.0" for regular
floats
lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
while s != lasts:
lasts = s
n, na = n+na, na+8
d, da = d+da, da+32
t = (t * n) / d
s += t
getcontext().prec -= 2
return +s # unary plus applies the new precision
It looks like an infinite series with term `t`, where`n` = (2k-1)^2
and `d` = d = 4k(4k+2) for k = 1... Does it have a name?
Mark Dickinson
Interesting. So the general term here is
3 * (2k choose k) / (16**k * (2*k+1)), k >= 0.
Python 3.2a0 (py3k:76334:76335, Nov 18 2009, 11:34:44)
[GCC 4.0.1 (Apple Inc. build 5490)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from math import factorial as f
>>> 3*sum(f(2*k)/f(k)/f(k)/(2*k+1)/16**k for k in range(50))
3.141592653589794
I've no idea what its name is or where it comes from, though. I
expect Raymond Hettinger would know.
--
Mark
Mark Dickinson
> > It looks like an infinite series with term `t`, where`n` = (2k-1)^2
> > and `d` = d = 4k(4k+2) for k = 1... Does it have a name?
>
> Interesting. So the general term here is
> 3 * (2k choose k) / (16**k * (2*k+1)), k >= 0.
>
> expect Raymond Hettinger would know.
After a cup of coffee, it's much clearer: this just comes from the
Taylor series for arcsin(x), applied to x = 1/2 to get asin(1/2) = pi/
6.
--
Mark
Mark Dickinson
> docs.python.org/library/decimal.html#recipes>?
>
> def pi():
> """Compute Pi to the current precision.
>
> >>> print pi()
> 3.141592653589793238462643383
>
> """
> getcontext().prec += 2 # extra digits for intermediate steps
> three = Decimal(3) # substitute "three=3.0" for regular
> floats
> lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
> while s != lasts:
> lasts = s
> n, na = n+na, na+8
> d, da = d+da, da+32
> t = (t * n) / d
> s += t
> getcontext().prec -= 2
> return +s # unary plus applies the new precision
And just for fun, here's a one-liner (well, okay, two lines including
the import) that uses Decimal to print the first 28 digits of pi:
Python 2.6.4 (r264:75706, Nov 16 2009, 15:42:08)
[GCC 4.0.1 (Apple Inc. build 5490)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from decimal import Decimal as D
>>> print reduce(lambda x,k:2+k/2*x/k,range(999,1,-2),D())
3.141592653589793238462643383
--
Mark
Albert van der Horst
Mark Dickinson <dick...@gmail.com> wrote:
>On Dec 11, 10:30=A0am, Mark Dickinson <dicki...@gmail.com> wrote:
>> > It looks like an infinite series with term `t`, where`n` =3D (2k-1)^2
>> > and `d` =3D d =3D 4k(4k+2) for k =3D 1... Does it have a name?
>>
>> Interesting. =A0So the general term here is
>> 3 * (2k choose k) / (16**k * (2*k+1)), =A0k >=3D 0.
>>
>> I've no idea what its name is or where it comes from, though. =A0I
>> expect Raymond Hettinger would know.
>
>After a cup of coffee, it's much clearer: this just comes from the
>6.
Curious. It seems better to calculate the zero of sin(pi/6)-1/2.
Not that we can forego the need of a Taylor series, but sin
converges much faster than arcsin.
The derivative is known analytically, and we have a 5th order process
before we know it.
It would be a big win for large precisions.
(Especially if we remember a previous value of pi to start up.)
The trick with temporarily increasing precision could be superfluous.
(I implemented this once in FORTRAN and was much disappointed that
double precision wasn't enough to show off the 5th order convergence. )
>--
>Mark
Groetjes Albert
--
--
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
albert@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst
Albert van der Horst
What the heck. I tried it all out. It turns out that sin(pi/6)=.5
is not favourable because the derivative contains sqrt(3)
So I tried to find the zero of cos(x) near pi/2 which is pi/2.
The derivative of the cos is -sin. So a better approximation than
x is x+cos(x). The second derivative is cos which is zero.
The third derivative is again -cos(x).
So a still better approximation is x+cos(x)+cos(x)^3/6.
This can be iterated. [The remaining error is 3/40.cos(x)^5.
I found that experimentally and a proof is left to the reader.]
Below you see cos which just calculates cosine with a
Taylor series.
Then there is pi2() that uses it to calculate pi, for the
normal fp precision.
pi3() shows the algorithm in its glory and should work for
any floating point package.
pi4() does the same, if precision is given.
And last but not least pi5() that uses the Decimal package
to advantage. It precalculates a starting point in 1/5 of
the precision. Then it does one more iteration in the full precision.
For 1000 digits it is about 5 times faster than pi(), for
a moderate increase in complexity.
# ----------- 8<----------------8<-----------------------
# $Id: pi.py,v 1.3 2009/12/21 19:01:15 albert Exp albert $
# Copyright (2008): Albert van der Horst {by GNU Public License}
# <http:// docs.python.org/library/decimal.html#recipes>
from decimal import getcontext,Decimal
def pi():
"""Compute Pi to the current precision.
>>> print pi()
3.141592653589793238462643383
"""
getcontext().prec += 2 # extra digits for intermediate steps
three = Decimal(3) # substitute "three=3.0" for regular floats
lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
while s != lasts:
print s
lasts = s
n, na = n+na, na+8
d, da = d+da, da+32
t = (t * n) / d
s += t
getcontext().prec -= 2
return +s # unary plus applies the new precision
def cos(halfpi):
"""Compute cos of halfpi
"""
x = halfpi**2
t = 1
lasts = 1
s = 0 # First term is actually 1
n=1
while s != lasts:
print s
lasts = s
t = -t*x / (n*(n+1))
n += 2
s += t
# Add 1 now, this saves iterations that don't contribute to precision.
return 1+s
def pi2():
' Calculate pi by a correction based on derivatives '
x=1.57
q=cos(1.57)
# Deviation 3/40.q^5
return 2*(x+q+q**3/6)
def pi3():
' Calculate pi by a 5th order process '
x=1.5
xold =1.
while x != xold:
xold = x
q = cos(x)
x += q*(1+q*q/6)
return 2*x
def pi4():
' Calculate pi by a 5th order process, with favorable stop criterion'
precision = 10e-20
rp = precision ** .2 # Required precision with room to spare.
print rp
x=1.5
q=1
while q>rp:
q = cos(x)
x += q*(1+q*q/6)
return 2*x
def pi5():
' Calculate pi by a 5th order process, adjusting precision'
oldprec = getcontext().prec
getcontext().prec = oldprec/4+1
rp = Decimal(10)**(-getcontext().prec+1)
x=Decimal("1.5")
q = x
while q >rp:
print x,q,rp
q = cos(x)
x += q*(1+q*q/6)
# One more iteration with full precision
getcontext().prec = oldprec + 2
q = cos(x)
x += q*(1+q*q/6)
getcontext().prec = oldprec
return 2*x
# ----------- 8<----------------8<-----------------------
Steven D'Aprano
On Mon, 21 Dec 2009 20:40:40 +0000, Albert van der Horst wrote:
> def pi4():
> ' Calculate pi by a 5th order process, with favorable stop
> criterion'
> precision = 10e-20
Why do you say 10e-20 instead of 1e-19?
--
Steven
Albert van der Horst
No thought went into that.
Note that the error jumps from 1e-5 to 1e-25 between iterations,
so 1e-20 or 1e-19 hardly makes a difference.
>--
>Steven
Groetjes Albert
Gabriel Genellina
<alb...@spenarnc.xs4all.nl> escribi�:
> In article <kv0dv...@spenarnc.xs4all.nl>,
> Albert van der Horst <alb...@spenarnc.xs4all.nl> wrote:
>> In article
>> <a266b155-02df-4294...@a32g2000yqm.googlegroups.com>,
>> Mark Dickinson <dick...@gmail.com> wrote:
>>>
>>> After a cup of coffee, it's much clearer: this just comes from the
>>> pi/6.
> Below you see cos which just calculates cosine with a
> Taylor series.
> Then there is pi2() that uses it to calculate pi, for the
> normal fp precision.
> pi3() shows the algorithm in its glory and should work for
> any floating point package.
> pi4() does the same, if precision is given.
> And last but not least pi5() that uses the Decimal package
> to advantage. It precalculates a starting point in 1/5 of
> the precision. Then it does one more iteration in the full precision.
> For 1000 digits it is about 5 times faster than pi(), for
> a moderate increase in complexity.
You may try Demo/scripts/pi.py in the source distribution; it uses long
integers to compute a continued fraction approximation to pi. I wonder how
does it compare to those other algorithms.
--
Gabriel Genellina
Steven D'Aprano
> In article <00b967e1$0$15623$c3e...@news.astraweb.com>, Steven D'Aprano
> <st...@REMOVE-THIS-cybersource.com.au> wrote:
>>Nice work! But I have a question...
>>
>>On Mon, 21 Dec 2009 20:40:40 +0000, Albert van der Horst wrote:
>>
>>> def pi4():
>>> ' Calculate pi by a 5th order process, with favorable stop
>>> criterion'
>>> precision = 10e-20
>>
>>
>>Why do you say 10e-20 instead of 1e-19?
>
> No thought went into that.
> Note that the error jumps from 1e-5 to 1e-25 between iterations, so
> 1e-20 or 1e-19 hardly makes a difference.
Ah, then it's a typo -- you've written TEN e-20 instead of ONE e-20 in
your code.
--
Steven