In article <kv0dvx....@spenarnc.xs4all.nl>,
Albert van der Horst <alb...@spenarnc.xs4all.nl> wrote:
>In article <a266b155-02df-4294-8b91-41be4a1e4...@a32g2000yqm.googlegroups.com>,
>Mark Dickinson <dicki...@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
>>Taylor series for arcsin(x), applied to x =3D 1/2 to get asin(1/2) =3D pi/
>>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
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<-----------------------
>
>--
>--
>Albert van der Horst, UTRECHT,THE NETHERLANDS
>Economic growth -- being exponential -- ultimately falters.
>alb...@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst
>
--
--
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
alb...@spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst
--
http://mail.python.org/mailman/listinfo/python-list
