How to calculate the ratio of convergence for Euler's, Gauss' and Viète's approximation of $\pi$?

Question

Let $\sqrt{6\sum_{k=1}^\infty{\frac{1}{k^2}}}$ be Euler's approximation of $\pi$; $\lim_{n\rightarrow\infty}\frac{2}{g_n}$ Gauss approximation of $\pi$; and $2\cdot\frac{2}{\sqrt{2}}\cdot\frac{2}{\sqrt{2+\sqrt{2}}}\cdot\frac{2}{\sqrt{2+\sqrt{2+\sqrt{2}}}}...$ Viète's approximation of $\pi$ where:

• $\{g_n\}_{n=1}^\infty$ is the succession given by:

$1$ for $n=1$;

$1/2$ for $n=2$;

$\sqrt{g_n-2g_{n-1}}$ for $n>2$ odd;

$\frac{g_{n-2}+g_{n-1}}{2}$ for $n>2$ even;

• $\{v_n\}_{n=1}^\infty$ is the succession given by:

$a_1=\sqrt{2}$;

$a_n=\sqrt{2+a_{n-1}}$;

$V_1= 2\cdot\frac{2}{\sqrt{2}}$;

$V_n=V_{n-1}\cdot\frac{2}{a_n}$

How to calculate $\lim_{n\rightarrow\infty}\frac{\lvert\pi_{n+1}^E-\pi_n^E\rvert}{\lvert\pi_{n}^E-\pi_{n-1}^E\rvert}$, $\lim_{n\rightarrow\infty}\frac{\lvert\pi_{n+1}^G-\pi_n^G\rvert}{\lvert\pi_{n}^G-\pi_{n-1}^G\rvert}$ and $\lim_{n\rightarrow\infty}\frac{\lvert\pi_{n+1}^V-\pi_n^V\rvert}{\lvert\pi_{n}^V-\pi_{n-1}^V\rvert}$ where $\pi_n^E, \pi_n^G, \pi_n^V$ are, respectively, Euler's, Gauss' and Viète's approximations of $\pi$?

I have a MATLAB code that calculates this ratios ($1$ for Euler's, $0.5$ Gauss' and $0.25$ for Viète's) but I wonder how would it be done "by hand", and which methods do we have to solve those limits.

Answer 1

I'd suggest not "by hand" but using some simple programming language (an interpreter language will do, so you do not need to compile your routine before running it). It should offer trace options or simple output to a cmd window (console). Myself I use REXX (I am a mainframe dinosaur). I am not up to date what else is a better option nowadays.

With it you may program the algorithm and observe how by each iteration the result gets closer to the expected value, and how many digits will be certain by each loop. Vieta is know to be quite slow. But (!) usually the progress per loop is indicated what is only half of the story: how much drudgery (or effort) is executed within one single loop? If we regard the procedures to get digits of $\pi$ as kind of a root finder maybe the Ostrowski index is a more meaningful performance indicator since it also takes into acount how much is investigated in every iteration.

Some time ago I tested several procedures to get many digits of $\pi$. I attach a REXX program I used -- if it looks too complicated then just take the single sections to get an idea how to program the algorithms in question.

BTW, following program includes also two fast algorithms, Brent-Salamin and Borwein&Borwein, just for fun.

/* some Pi */

-- trace '?R'
if arg() then numeric digits arg(1)
numeric fuzz format(RxCalcLog10(digits()), , 0)
e = 10 ** -(digits() - fuzz() - 1)

d = time("R")
if digits() > 999 then signal fast

/* Pi mit Newton, konvergiert > 3 loops per digit (keine Wurzel) */
z = 1
m = 1
n = 1
p = 1
DO i = 1 UNTIL k < e
   z = z * i
   m = m + 2
   n = n * m
   k = z / n
   p = p + k
end
say '...' || right(2 * p, 33) time("R") i "Newton"

/* Pi mit Madhava-Leibniz, konvergiert ~ 2 loops per digit (eine Wurzel ausserhalb der Schleife) */
s = 1
n = 1
DO m = 3 by 2 UNTIL m = m + n
   n = n / -3
   s = s + n / m
end
say '...' || right(2 * s * Sqrt(3), 33) time("R") (m - 1) / 2 "Madhava-Leibniz"

/*
** -- Pi mit Madhava-Leibniz mod, konvergiert ~ 1 loops per digit (eine Wurzel ausserhalb der Schleife)
** -- ca. 1/3 schneller als das Original
*/
s = 1
n = 1
DO m = 3 by 4 UNTIL m = m + n
   n = n / 9
   s = s + n / (m + 2) - 3 * n / m
end
say '...' || right(2 * s * Sqrt(3), 33) time("R") (m - 3) / 4 "Madhava-Leibniz mod"

/* Pi mit Vieta, konvergiert ~ 1,7 loops per digit (Wurzel in der Schleife) */
p = 2
t = 0
DO i = 1 UNTIL t = 2
   t = Sqrt(2 + t)
   p = 2 * p / t
end
say '...' || right(p, 33) time("R") i "Vieta"

if digits() > 666 then signal fast
/* s. atan nach Acton, Pi = 6 * atan(sqrt(1/3) (Wurzel in der Schleife, ggf. davor und danach) */
a = 1 / sqrt(1+(1/3))
b = 1
f = a
do i = 1 until a = b
   a = (a + b) / 2
   b = sqrt(a * b)
end
say '...' || right(6*sqrt(1/3)*f/a, 33) time("R") i "x*Atan(y)"

/*                                        */
/*            */  fast:  /*               */
/*                                        */

/* Pi mit Brent-Salamin, konvergiert 10..600 digits per loop, jenachdem (Wurzel in der Schleife) */
a = 1
b = 1 / 2
r = 1
s = b
x = a
DO i = 1 UNTIL k < e
   b = Sqrt(x * b)
   x = a
   a = (a + b) / 2
   c = (a - b) ** 2
   r = r + r
   k = c * r
   s = s - k
end
say '...' || right(2 * a**2 / s, 33) time("R") i "Brent-Salamin"

/*  Pi mit Borwein&Borwein (eine Quadratwurzel vor, eine vierte Wurzel in der Schleife) */
e = 10 ** -((digits() - fuzz()) % 4 + 1)
y = sqrt(2)
a = 6 - 4 * y
y = y - 1
do i = 3 by 2 until y < e
--    y = vwz(1 - y ** 4)
   y = sqrt(sqrt(1 - y ** 4))    /* so schneller als 4. Wurzel ermitteln */
   y = (1 - y) / (1 + y)
   a = a * (1 + y) ** 4 - 2 ** i * y * (1 + y + y ** 2)
end
say '...' || right(1 / a, 33) time("R") (i - 1) / 2 "Borwein&Borwein 1985"

exit

/* * * * *  Subroutinen  * * * * */

Sqrt: procedure; arg z
nd = digits()
numeric digits 22    /* start with 'few' digits only */
-- say z
-- a = z / 2
-- a = 1 + 0.5 * (z - 1)      /* erste Näherung wenn z ~ 1 */
a = RxCalcSqrt(z)             /* auch eine Näherung */
f = 0
do until b = a & f = 1
   b = a
   if digits() = nd then f = 1
   a = .5 * (a + z / a)
--    say 'Delta:' format(a - b, 2, 3, ,0) || ', Digits:' digits()
   numeric digits min(nd, 8 * digits(), digits() + 2 * word(translate(format(a - b, , 0, ,0), , '-') nd, 2))
--   say digits()
end
return a

/* Vierte Wurzel ziehen, Sekantenverfahren */
vwz: procedure expose e; arg z
nd = digits()
numeric digits 22             /* start with 'few' digits only */
-- a = (z + 3) / 4            /* erste Näherung wenn z ~ 1 */
a = RxCalcSqrt(RxCalcSqrt(z))
f = 0
do until b = a & f = 1
   b = a
   if digits() = nd then f = 1
   a = 0.25 * (3 * a + z / a ** 3)   /* a - f(a)/f'(a) */
   say 'Delta:' format(a - b, 2, 3, ,0) || ', Digits:' digits()
   numeric digits min(nd, digits() || 0, digits() + 2 * word(translate(format(a - b, , 0, ,0), , '-') nd, 2))
end
return a

::requires 'rxmath' LIBRARY

The program by itself is not of much help, thus I show an output of a sample run. Input is number of digits to compute, output per procedure is i) last 33 digits of $\pi$, ii) time elapsed in seconds, iii) number of iterations, iv) type of procedure. Note that for example "Newton" takes about twice as many loops than "Vieta" but is anyway a little bit faster. That is why I suggest to use the Ostrowski index.

C:\PRGM\rexx\Pi>piite.rx 660
...134275778960917363717872146842686 6.052000 2174 Newton
...134275778960917363717872146844191 1.276000 1370 Madhava-Leibniz
...134275778960917363717872146844191 0.911000 684 Madhava-Leibniz mod
...134275778960917363717872146843491 6.340000 1091 Vieta
...134275778960917363717872146843625 7.857000 1091 x*Atan(y)
...134275778960917363717872146844093 0.081000 9 Brent-Salamin
...134275778960917363717872146844338 0.109000 4 Borwein&Borwein 1985

Answer 2

Typically, you can try to obtain an asymptotic expansion for each one of the sequences listed in your question, which tells you a lot about how fast each sequence approaches $\pi$.

1. Start by considering the difference $$\frac{\pi^2}{6} -\sum_{k=1}^n\frac{1}{k^2}=\sum_{k=n+1}^{\infty} \frac{1}{k^2}= \frac{1}{n}+o(n^{-1}),$$ where the last step can be obtained by considering the difference $$\sum_{k=n+1}^{\infty} \frac{1}{k^2} - \int_{n}^{\infty}\frac{1}{x^2} dx$$ Therefore, we can provide the following rough estimation for $\pi^E_n$ (using Taylor expansion of $\sqrt{1+x}=1+x/2+o(x)$ $$\pi^E_n= \sqrt{\pi^2 -6n^{-1}} +o(n^{-1}) = \pi +\frac{3}{\pi n}+o(n^{-1}).$$ Which means that $\frac{\pi^E_{n+1}-\pi^E_{n}}{\pi^E_n-\pi^E_{n-1}}$ converges to $1$. Roughly speaking, $\pi^E_n$ converges rather slowly to $\pi$: If you wish to approximate $\pi$ to six decimal places using $\pi^E_n$, you need $n\sim 10^6$. Note: You can obtain more terms in the asymptotic expansion of $\pi^E_n$ by using Euler-Maclaurin formulas.

2. There is a [general form][1] that you can use to derive an asymptotic expansion for $\pi^V_{n}$. I could be mistaken, but it should be something like $\pi^V_{n}=\frac{2}{\pi}(1+c4^{-n}) +o(4^{-n})$ where $c$ is some constant. So the answer to your question should be $4^{-1}$.

3. I am not sure, but this seems to be equivalent to the Gauss-Legendre algorithm, which has a quadratic convergence, so you expect the limit to be zero.