Re: A Formula for Pi

Mensanator wrote:
On Jun 20, 2:35 pm, r...@xxxxxxxxxxxxxx (Rob Johnson) wrote:
In article <4f348778-0c41-45e4-bb51-9fcc3dca9...@xxxxxxxxxxxxxxxxxxxxxxxxxxxx>,
Mensanator <mensana...@xxxxxxx> wrote:
On Jun 20, 10:57 am, Maury Barbato <mauriziobarb...@xxxxxxxx> wrote:
Jose Carlos Santos wrote:
On 20-06-2008 7:16, Maury Barbato wrote:
I found in the book "The Penguin Dictionary of
Curious and Interesting Numbers" by Wells the
following formula involving pi
....
oo
--- k+1 1
> (-1) --------------
--- 2k(2k+1)(2k+2)
k=1
....
pi - 3
= ------
4

Being an alternating series with monotonically decreasing terms, the
error is less than 1/(8n^3) after n terms.

Some of us don't know how to do this. But given the series, I can
write a program to apply it.

But I would prefer a series that converges in ~300 terms to the
one that converges in ~10**34 terms.
....

The magic sumalt algorithm (based on Cohen, Villegas & Zagier) will do
here, it is implemented in PARI and here is one in Maple (taken from
Arndt Joerg's book 'Algorithms for programmers').

It uses rational arithmetics only (and since Maple uses GMP for that
may fit your needs).

For the task it needs only 132 steps for 101 decimal and 1307 for 1001.

sumalt:= proc(A,n)
# computes S(n) = Sum( A(k), k=0..n) with error smaller then
# abs( S - S(n) ) <= 2*abs(S)/(3+sqrt(8))^n, S = S(infinity)
local b,c,s,k;

b:=2^(2*n-1);
c:=b;
s:=0;

for k from n-1 to 0 by -1 do
s:= s + c* A(k); # A(k) = a(k)*(-1)^k;
b:= b * ((2*k+1)*(k+1)) / (2*(n+k)*(n-k));
c:= c + b;
end do;

return s/c;
end proc;

For summing from 0 we need A := k -> (-1)^k / ((k+1)*(2*k+3)*(k+2))
as function for the coefficients by re-indexing.

We need n = (ln(2)-ln(epsilon)+ln(abs(S)))/ln(3+2*2^(1/2)) for an
error epsilon, so theN:=132 will do for epsilon = 10^(-theDigits),
theDigits:=101, S:= Pi - 3 the fractional part of Pi.

Then

R0:=sumalt(A, theN);
R1:=evalf[theDigits](R0);
R2:=evalf[theDigits+1](Pi-3);

shows exactness up to the last decimal place (one need an additional
Digit in Maple, since Pi has 1 place before the decimal point).

Execution is very fast (based on GMP, even if there are layer between
Maple and that library).
.

Relevant Pages

  • Re: Overflow in MMA in Maple
    ... I take it that Maple ... How does it do the modular reductions so quickly?? ... If it's using GMP, then it does whatever GMP uses, obviously, which ... They no longer do my traditional winks tournament lunch - liver and bacon. ...
    (sci.math)
  • Re: Overflow in MMA in Maple = need larger universe?
    ... Since you have GMP, why do you care? ... cannot be stored, digit by digit, explicitly, ... If you used the corresponding program in Maple or Mma you ...
    (sci.math.symbolic)
  • Re: Overflow in MMA in Maple
    ... > I was wondering why MMA and Maple overflow on this example: ... > output conversion took 8 ms ... > (This is computed by demos/pexpr.c from the GMP distribution) ... I thought that from Maple 8 on GMP was included. ...
    (sci.math.symbolic)
  • Re: Overflow in MMA in Maple
    ... > I was wondering why MMA and Maple overflow on this example: ... As for Maple, you are just using the wrong symbol. ... computation by modular exponentiation. ... Maple does indeed use GMP, and it was used for the above computation. ...
    (sci.math.symbolic)
Loading