Prev: Python tricks with applescript in OS-X
Next: Which graph library is best suited for large graphs?
From: Albert van der Horst on 21 Dec 2009 15:40
In article <kv0dvx.7k5(a)spenarnc.xs4all.nl>,
Albert van der Horst <albert(a)spenarnc.xs4all.nl> wrote:
>In article <a266b155-02df-4294-8b91-41be4a1e4f81(a)a32g2000yqm.googlegroups.com>,
>Mark Dickinson <dickinsm(a)gmail.com> wrote:
>>On Dec 11, 10:30=A0am, Mark Dickinson <dicki...(a)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.
>albert(a)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.
albert(a)spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst

From: Steven D'Aprano on 21 Dec 2009 16:05
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?


--
Steven
From: Albert van der Horst on 21 Dec 2009 23:04
In article <00b967e1$0$15623$c3e8da3(a)news.astraweb.com>,
Steven D'Aprano <steve(a)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.

>--
>Steven

Groetjes Albert

--
--
Albert van der Horst, UTRECHT,THE NETHERLANDS
Economic growth -- being exponential -- ultimately falters.
albert(a)spe&ar&c.xs4all.nl &=n http://home.hccnet.nl/a.w.m.van.der.horst

From: Gabriel Genellina on 21 Dec 2009 23:02
En Mon, 21 Dec 2009 17:40:40 -0300, Albert van der Horst
<albert(a)spenarnc.xs4all.nl> escribi�:
> In article <kv0dvx.7k5(a)spenarnc.xs4all.nl>,
> Albert van der Horst <albert(a)spenarnc.xs4all.nl> wrote:
>> In article
>> <a266b155-02df-4294-8b91-41be4a1e4f81(a)a32g2000yqm.googlegroups.com>,
>> Mark Dickinson <dickinsm(a)gmail.com> wrote:
>>>
>>> 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.

> 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

From: Steven D'Aprano on 22 Dec 2009 01:03
On Tue, 22 Dec 2009 04:04:56 +0000, Albert van der Horst wrote:

> In article <00b967e1$0$15623$c3e8da3(a)news.astraweb.com>, Steven D'Aprano
> <steve(a)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
First  |  Prev  | 
Pages: 1 2
Prev: Python tricks with applescript in OS-X
Next: Which graph library is best suited for large graphs?