From: Anh Hai Trinh on 11 Dec 2009 03:16 I'm just curious which formula for pi is given here: ? 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? From: Mark Dickinson on 11 Dec 2009 05:30 On Dec 11, 8:16 am, Anh Hai Trinh wrote:> I'm just curious which formula for pi is given here: 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? 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 From: Mark Dickinson on 11 Dec 2009 05:43 On Dec 11, 10:30 am, Mark Dickinson wrote:> > 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. > > I've no idea what its name is or where it comes from, though.  I > 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 From: Mark Dickinson on 11 Dec 2009 07:49 On Dec 11, 8:16 am, Anh Hai Trinh wrote:> I'm just curious which formula for pi is given here: 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 From: Albert van der Horst on 21 Dec 2009 10:33 In article , Mark Dickinson wrote:>On Dec 11, 10:30=A0am, Mark Dickinson 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 -- -- 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  |  Next  |  Last