Prev: Python tricks with applescript in OS-X
Next: Which graph library is best suited for large graphs?
From: Anh Hai Trinh on 11 Dec 2009 03:16
I'm just curious which formula for pi is given here: <http://
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?
From: Mark Dickinson on 11 Dec 2009 05:30
On Dec 11, 8:16 am, Anh Hai Trinh <anh.hai.tr...(a)gmail.com> wrote:
> I'm just curious which formula for pi is given here: <http://
> 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 <dicki...(a)gmail.com> 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 <anh.hai.tr...(a)gmail.com> wrote:
> I'm just curious which formula for pi is given here: <http://
> 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 <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

--
--
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
Pages: 1 2
Prev: Python tricks with applescript in OS-X
Next: Which graph library is best suited for large graphs?