Here's a nice little Mathematica routine for evaluating Tito's continued fraction with precision prec
:
prec = 10^4;
y = N[4, prec];
c = y; d = 0; k = 1;
u = 1; v = y;
While[True,
c = 1 + u/c; d = 1/(1 + u d);
h = c*d; y *= h;
v += 96 k^2 + 8;
c = v + u/c; d = 1/(v + u d);
h = c*d; y *= h;
If[Abs[h - 1] <= 10^-prec, Break[]];
u += 3 k (k + 1) + 1;
k++];
6/y
where I used the Lentz-Thompson-Barnett method for the evaluation.
For prec = 10^4
, the thing evaluates in 120 seconds (via AbsoluteTiming[]
), giving a result that agrees with $\zeta(3)$ to 10,000 digits.
One can consider the even part of Tito's CF, which converges at twice the rate of the original:
$$\cfrac{6}{5-\cfrac{u_1}{v_1-\cfrac{u_2}{v_2-\cfrac{u_3}{v_3-\cdots}}}}$$
where
$$\begin{align*}
u_k&=k^6\\
v_k&=(17k^2+17k+5)(2k+1)
\end{align*}$$
Here's Mathematica code corresponding to this CF:
prec = 10^4;
y = N[5, prec];
c = y; d = 0; k = 1;
While[True,
u = k^6;
v = (2 k + 1) ((17 k + 17) k + 5);
c = v - u/c; d = 1/(v - u d);
h = c*d; y *= h;
If[Abs[h - 1] <= 10^-prec, Break[]];
k++];
6/y
For prec = 10^4
, the thing evaluates in 70 seconds (via AbsoluteTiming[]
). There may be further ways to accelerate the convergence of the CF, but I have yet to look into them.
Added, quite a bit later:
As it turns out, the even part I derived is precisely Apéry's CF for $\zeta(3)$ (thanks Américo!). Conversely put, Tito's CF is an extension of Apéry's CF. Here's how to derive Apéry's CF from Tito's CF (while proving convergence along the way).
We start from an equivalence transformation of Tito's CF. A general equivalence transformation of a CF
$$b_0+\cfrac{a_1}{b_1+\cfrac{a_2}{b_2+\cfrac{a_3}{b_3+\cdots}}}$$
with some sequence $\mu_k, k>0$ looks like this:
$$b_0+\cfrac{\mu_1 a_1}{\mu_1 b_1+\cfrac{\mu_1 \mu_2 a_2}{\mu_2 b_2+\cfrac{\mu_2 \mu_3 a_3}{\mu_3 b_3+\cdots}}}$$
Now, given a CF
$$b_0+\cfrac{a_1}{b_1+\cfrac{a_2}{b_2+\cdots}}$$
one can transform this into a CF of the form
$$b_0+\cfrac{w_1}{1+\cfrac{w_2}{1+\cdots}}$$
where $w_1=\dfrac{a_1}{b_1}$ and $w_k=\dfrac{a_k}{b_k b_{k-1}}$ for $k > 1$, where we used $\mu_k=\dfrac1{b_k}$. Applying this transformation to Tito's CF yields the CF
$$\cfrac{\frac32}{1+\cfrac{w_2}{1+\cfrac{w_3}{1+\cdots}}}$$
where $w_{2k}=\dfrac{k^3}{4(2k-1)^3}$ and $w_{2k+1}=\dfrac{k^3}{4(2k+1)^3}$. (You can easily demonstrate that this transformed CF and Tito's CF have identical convergents.)
At this point, we find that since the $w_k \leq\dfrac14$, we have convergence of the CF by Worpitzky's theorem.
Now, we move on to extracting the even part of this transformed CF. Recall that if a CF has the sequence of convergents
$$u_0=b_0,u_1=b_0+\cfrac{a_1}{b_1},u_2=b_0+\cfrac{a_1}{b_1+\cfrac{a_2}{b_2}},\dots$$
then the even part is the CF whose convergents are $u_0,u_2,u_4,\dots$ (Analogously, there is the odd part with the sequence of convergents $u_1,u_3,u_5,\dots$)
Now, given a CF of the form
$$b_0+\cfrac{w_1}{1+\cfrac{w_2}{1+\cdots}}$$
its even part is the CF
$$b_0+\cfrac{w_1}{1+w_2-\cfrac{w_2 w_3}{1+w_3+w_4-\cfrac{w_4 w_5}{1+w_5+w_6-\cdots}}}$$
Thus, the even part of the previously transformed CF is given by
$$\cfrac{\frac32}{\frac54-\cfrac{\beta_1}{\delta_1-\cfrac{\beta_2}{\delta_2-\cdots}}}$$
where
$$\begin{align*}
\beta_k&=\frac{k^3}{4(2k-1)^3}\frac{k^3}{4(2k+1)^3}=\frac{k^6}{16(2k-1)^3(2k+1)^3}\\
\delta_k&=1+\frac{k^3}{4(2k+1)^3}+\frac{(k+1)^3}{4(2k+1)^3}=\frac{17k^2+17k+5}{4(2k+1)^2}
\end{align*}$$
We're almost there! We only need to perform another equivalence transformation, which I'll split into two steps to ease understanding. First, the easy one with $\mu_k=4$, which yields the CF
$$\cfrac{6}{5-\cfrac{16\beta_1}{4\delta_1-\cfrac{16\beta_2}{4\delta_2-\cdots}}}$$
The last step is to cancel out the odd integer denominators of the $\beta_k$ and $\delta_k$; to do this, we take $\mu_k=(2k+1)^3$; this finally yields the CF
$$\cfrac{6}{5-\cfrac{u_1}{v_1-\cfrac{u_2}{v_2-\cfrac{u_3}{v_3-\cdots}}}}$$
where
$$\begin{align*}
u_k&=k^6\\
v_k&=(17k^2+17k+5)(2k+1)
\end{align*}$$
and this is Apéry's CF.
For completeness, I present a formula for the odd part of Tito's CF, after some post-processing with a few equivalence transformations:
$$\zeta(3)=\frac32-\cfrac{81}{\lambda_1-\cfrac{\eta_1}{\lambda_2-\cfrac{\eta_2}{\lambda_3-\ddots}}}$$
where
$$\begin{align*}
\eta_k&=4\times(4k^4+8k^3+k^2-3k)^3=4\times10^3,\,4\times126^3,\dots\\
\lambda_k&=8\times(68k^6-45k^4+12k^2-1)=8\times34,\,8\times3679,\dots
\end{align*}$$
The formula is somewhat more complicated, and converges at the same rate as the even part.
UPDATE 2
Sanchez proposed the example of $\frac 1{1-x}=1+x+x^2+\cdots$ with $\,x\in\mathbb{R}\,$ that becomes $\frac 1{1-z}=1+z+z^2+\cdots$ with $\,z\in\mathbb{C}\,$. Note that this is not defined for $z=1$.
Practically this means that the computations over reals have to be replaced by the corresponding computations over complex values.
Let's start on the real axis near $x_0$ with the more general expansion :
$$f(x)=a_0+a_1(x-x_0)+a_2(x-x_0)^2+\cdots$$
To simplify a little analytic continuation is about getting the same expansion $$f(z)=a_0+a_1(z-z_0)+a_2(z-z_0)^2+\cdots$$ for complex numbers $z$ around a complex point $z_0$ and similar expansions $f(z)=b_0+b_1(z-z_1)+b_2(z-z_1)^2+\cdots\;$ around other points $z_1,\;z_2$ and so to cover as much as possible the complex plane.
Let's start for example with $z_0=x_0=0$ for the $\frac 1{1-z}$ function (with a unit radius convergence disk) followed by two other disks centered at the complex values $z_1$, $z_2$ (avoiding the pole at $z=1$) then we get this picture :
One advantage of complex analysis is that you'll see directly that for example $\displaystyle\frac 1{1+z^2}$ has a problem at $z=i$ and $z=-i$ (because the denominator vanishes !). Since $|i-0|=|-i-0|=1$ you may deduce that the convergence radius at $z_0=0$ of $\frac 1{1+z^2}=1-z^2+z^4-\cdots\,$ is $1$ (since the unit circle centered at $0$ will 'touch' the imaginary $i$ as well as $-i$ as illustrated above).
The complex points we need to avoid are named 'non removable singularities'. Near one of these points the smooth Taylor expansion will be replaced by a Laurent series :
$$f(z)=\cdots+\frac{a_{-2}}{(z-z_0)^2}+\frac{a_{-1}}{z-z_0}+a_0+a_1(z-z_0)+a_2(z-z_0)^2+\cdots$$
(with an infinite number of negative exponents in the case of an 'essential singularity')
Computing derivatives and integrals is done nearly as in the real case but with the second dimension we gained the possibility of integrating around circular contours (coming back to the starting point using a different path).
Let's compute $\int_C f(z)\;dz$ by evaluating its constituents : we have as in the real case $\int (z-z_0)^n\,dz=\frac {(z-z_0)^{n+1}}{n+1}$ except in the case $n=-1$ where we get the complex logarithm $\log(z-z_0)$. For $n\not =-1$ a counterclockwise contour integral will subtract two times the same value and be $0$ while for $n=-1$ the imaginary part will have gained and additional $2\pi\,i$ (see the picture from the last link).
This implies that the integral of the whole Laurent series over the contour will simply be $\ 2\pi\;i\;a_{-1}\ $ with $a_{-1}$ the famous 'residue' of $f$ at $z_0$.
Let's come to the zeta function that you probably saw defined as :
$$\zeta(z)=\sum_{n=1}^\infty \frac 1{n^z}$$
(or with $z$ written as $x$ real but at this point you should understand that it doesn't matter : the expression will be the same !)
The problem with this series is convergence : for real values of $z$ it will converge only for values larger than $1$ (from the 'integral test'). For complex values the problem is the same : it won't converge for $\Re(z)<1$ so that this series can't be used directly to find any zero of $\zeta$.
Fortunately a slight modification will allow to get convergence for $\Re(z)>0$ :
\begin{align}
\zeta(z)&=\frac 1{1^z}+\frac 1{2^z}+\frac 1{3^z}+\frac 1{4^z}+\frac 1{5^z}+\cdots\\
&=\frac 1{1^z}-\frac 1{2^z}+\frac 1{3^z}-\frac 1{4^z}+\frac 1{5^z}+\cdots+\frac 2{2^z}+\frac 2{4^z}+\cdots\\
&=\sum_{n=1}^\infty \frac {(-1)^{n-1}}{n^z}+\frac 2{2^z}\left(1+\frac 1{2^z}+\frac 1{3^z}\cdots\right)\\
\zeta(z)&=\sum_{n=1}^\infty \frac {(-1)^{n-1}}{n^z}+2^{1-z}\zeta(z)\\
\end{align}
Putting both $\zeta(z)$ at the left gives $\ \displaystyle\zeta(z)\left(1-2^{1-z}\right)=\sum_{n=1}^\infty \frac {(-1)^{n-1}}{n^z}\;$ and :
$$\zeta(z)=\frac 1{1-2^{1-z}}\sum_{n=1}^\infty \frac {(-1)^{n-1}}{n^z}$$
There is still a singularity at $z=1$ but we have now an alternate series that will converge on the real interval $(0,1)$ and in fact in the half-plane $\Re(z)>0$ (excluding only $z=1$) (proof). The (grossly schematized) analytic continuation process allows to go from the half-plane $\Re(z)>1$ to the half-plane $\Re(z)>0$ (minus $z=1$) and defines a unique function. We may for example use a convergence disk of radius $1$ centered at $z_3=2$ and then turn around $1$ using the $z_2$ from the picture and continue the process (remaining in the region $\Re(z)>0$).
Another formula for $\zeta$ that works well is obtained using Euler Maclaurin's expansion :
$$\zeta(z) \sim \sum_{k=1}^N \frac 1{k^z} \color{#bb0000}{+\frac 1{(z-1)\;N^{z-1}}}\\\color{#006600}{-\frac 1{2\;N^z}+\frac z{12\;N^{z+1}}-\frac{z(z+1)(z+2)}{720\;N^{z+3}}+\frac{z(z+1)(z+2)(z+3)(z+4)}{30240\;N^{z+5}}}$$
The neat idea here (from Euler like most of the initial work about the real zeta function) is to :
- compute the first terms of the $\zeta$ series (more of them will mean higher precision and in fact we will need $2\,\pi\,N >\ |\Im(z)|$ so that $N=10$ is enough for imaginary parts under $60$)
- replace the remaining terms by the corresponding integral $\;\color{#bb0000}{\int_{N+1}^\infty \frac {dk}{k^z}}$
- add corrective terms to take into account the previous replacement of the sum by an integral. The numerical coefficients appearing here are the famous Bernoulli numbers divided by the corresponding factorial $\frac {B_n}{n!}$ with $B_{2n+1}=0$ for $n>0$.
It is interesting to observe that this formula remains valid for complex values and works too for negative real parts of $z$ (depending on the number of Bernoulli terms at the end).
Let's show too the Laurent series of $\zeta$ at $z=1$ with the simple pole at $1$ :
$$\zeta(z)=\frac 1{z-1}+\sum_{n=0}^\infty \frac{(-1)^n}{n!}\gamma_n\;(z-1)^n$$
with $\gamma_n$ the Stieltjes constants and $\gamma_0=0.5772156649\cdots$ the famous Euler constant.
For faster evaluation see the different links from this thread or this discussion.
Best Answer
Article $[1]$ by Alfred van der Poorten is a good starting point, the best I know of.
Convergence speed
It is known that the given series converges very slowly. Apéry's rational approximation $a_n/b_n$ below improves the speed of convergence to $\zeta(3)$, while controlling the increase of the size of $b_n$.
Apéry stated the following equality:
$$\zeta (3)\overset{\mathrm{def}}{=}\sum_{n=1}^{\infty }\frac{1}{n^{3}}=\frac{5}{2}\sum_{n=1}^{\infty}\frac{(-1)^{n-1}}{n^{3}\binom{2n}{n}}.\tag{1}$$
In section 3 van der Poorten shows that
$$\sum_{n=1}^{N}\frac{1}{n^{3}}+\sum_{k=1}^{N}\frac{(-1)^{k-1}}{2k^{3}\binom{N}{k}\binom{N+k}{k}}=\frac{5}{2}\sum_{k=1}^{N }\frac{(-1)^{n-1}}{k^{3}\binom{2k}{k}},\tag{2}$$
from which $(1)$ follows by letting $N$ tend to infinity (see this answer of mine). The formula $(2)$ is not very intuitive though.
The second series in $(1)$ is a fast convergent series, faster by far than the defining series for the Apéry's constant $\zeta(3)$. There are even faster convergent series, obtained by techniques of convergence acceleration.
Proof of the irrationality
As far as the irrationality concerns, Apéry constructed two sequences $(a_n),(b_n)$ $^1$ whose ratio $a_n/b_n\to\zeta(3)$ and
This is enough to prove the irrationality of $\zeta (3)$ by contradiction. $[2]$.
Intuition
Given that in van der Poorten's words
it is natural to ask: Where did these ideas come from? My tentative explanation is based on the following fact. A few years later after his proof $[3]$, Roger Apéry derived the rational approximation $a_n/b_n$ to $\zeta(3)$ by transforming the defining series for $\zeta(3)$ into a continued fraction and applying iterated transformations to it, which improved the speed of convergence, and obtained the recurrence relations satisfied by $a_{n},b_{n}$ $[4]$.
--
$^1$ The sequences are $$\begin{equation*} a_{n}=\sum_{k=0}^{n}\binom{n}{k}^{2}\binom{n+k}{k}^{2}c_{n,k}, \qquad b_{n}=\sum_{k=0}^{n}\binom{n}{k}^{2}\binom{n+k}{k}^{2},\end{equation*}$$
where
$$\begin{equation*} c_{n,k}=\sum_{m=1}^{n}\frac{1}{m^{3}}+\sum_{m=1}^{k}\frac{\left( -1\right) ^{m-1}}{2m^{3}\binom{n}{m}\binom{n+m}{m}}\quad k\leq n. \end{equation*}$$
References.
$[1]$ Poorten, Alf., A Proof that Euler Missed…, Apéry’s proof of the irrationality of $\zeta(3)$. An informal report, Math. Intelligencer 1, nº 4, 1978/79, pp. 195-203.
$[2]$ Fischler, Stéfane, Irrationalité de valeurs de zêta (d’ après Apéry, Rivoal, …), Séminaire Bourbaki 2002-2003, exposé nº 910 (nov. 2002), Astérisque 294 (2004), 27-62
$[3]$ Apéry, Roger (1979), Irrationalité de $\zeta2$ et $\zeta3$, Astérisque 61: 11–13
$[4]$ Apéry, Roger (1981) Interpolations de Fractions Continues et Irrationalité de certaines Constantes, Bull. section des sciences du C.T.H.S., n.º3, p.37-53