Prev: Python tricks with applescript in OSX
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 <a266b15502df42948b9141be4a1e4f81(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 (2k1)^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 = 10e20 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 = 10e20 Why do you say 10e20 instead of 1e19?  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)REMOVETHIScybersource.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 = 10e20 > > >Why do you say 10e20 instead of 1e19? No thought went into that. Note that the error jumps from 1e5 to 1e25 between iterations, so 1e20 or 1e19 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 >> <a266b15502df42948b9141be4a1e4f81(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)REMOVETHIScybersource.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 = 10e20 >> >> >>Why do you say 10e20 instead of 1e19? > > No thought went into that. > Note that the error jumps from 1e5 to 1e25 between iterations, so > 1e20 or 1e19 hardly makes a difference. Ah, then it's a typo  you've written TEN e20 instead of ONE e20 in your code.  Steven
First

Prev

Pages: 1 2 Prev: Python tricks with applescript in OSX Next: Which graph library is best suited for large graphs? 