An Algolib-aided Version of Apery's Proof of the Irrationality of

Bruno Salvy

(March 4, 2003)

Apery proved in 1978 that  is irrational. We give a short version of Apery's proof that uses several tools from Algolib: gfun, Mgfun and equivalent. We only prove irrationality here and do not compute irrationality measures.

The starting point is the definition of three sequences:

 > c[n,k]:=binomial(n,k)^2*binomial(n+k,k)^2;

 > a[n]:=Sum(c[n,k],k=0..n);

 > b[n]:='a[n]'*Sum(1/k^3,k=1..n)+Sum(Sum((-1)^(m+1)*c[n,k]/2/m^3/binomial(n,m)/binomial(n+m,m),m=1..k),k=1..n);

1.

The first part of the right-hand side of the definition of  provides the desired limit for / . The proof consists in showing that the second part tends to 0, when divided by . Let

 > U:=m^3*binomial(n,m)*binomial(n+m,m);

What we are going to show is that . Thus the second sum in  is bounded by , itself bounded by , whence the desired limit.

The following simple way of proving this was suggested to us by Philippe Dumas:

For m=n, we have the desired inequality:

 > eval(U,m=n);

Otherwise, we have that , and the conclusion follows since for  between 1 and , one has .

2.  is a positive integer and  is an integer, where  is the lcm of 1,...,

The first assertion comes from binomials being integers.

The second one is obtained by showing that each summand in the double sum defining  has a denominator dividing . This is done in two steps as follows:

 > (binomial(n,k)^2*binomial(n+k,k)^2/m^3/binomial(n,m)/binomial(n+m,m))=binomial(n,k)*binomial(n+k,k)*binomial(n-m,n-k)*binomial(n+k,k-m)/m^3/binomial(k,m)^2;

Proof:

 > convert(lhs(%)/rhs(%),GAMMA);

and then it is sufficient to observe that  divides .

3. Both  and  satisfy a simple linear recurrence

This is where the Mgfun  package will be useful.

Recurrence satisfied by

 > CT:=Mgfun[creative_telescoping](c[n,k],n::shift,k::shift);

This result means that applying the recurrence operator given by the first element of this list to   , one gets the difference , where  is the second element of this list. Explicitly,

 > eval(CT[1],_f=proc(N,K) 'c[N,K]' end)=g(n,k)-g(n,k+1);

with

 > g(n,k)=subs(_f(n,k)='c[n,k]',CT[2]);

We check this formula by computing the ratio left-hand side divided by right-hand side:

 > ct:=eval(CT,_f=unapply(c[n,k],n,k)):

 > normal(expand(ct[1]/(ct[2]-subs(k=k+1,ct[2]))));

Now, summing this identity for  from 0 to +2, we get that the right-hand side telescopes, leaving out . This in turn is 0:

 > eval(ct[2],k=0);

 > normal(limit(ct[2],k=n+3));

The sequence  is 0 for  and , and thus its sum up to  is equal to .

Thus we have proved that the sequence  satisfies the following recurrence (with )

 > collect(eval(CT[1],_f=proc(n,k) A(n) end),A,factor)=0;

 > rec:=op(1,%):

We check this recurrence by computing the first few values of :

 > lista:=[seq(value(subs(n=i,a[n])),i=1..20)];

and then substituting these in the left-hand side of the recurrence:

 > seq(eval(rec,[n=i,A=proc(k) lista[k] end]),i=1..18);

Recurrence satisfied by

First, we observe that the recurrence above seems to be satisfied by the sequence :

 > seq(eval(rec,[n=i,A=proc(k) listb[k] end]),i=2..7);

In order to get a proof  that the recurrence is actually satisfied by , we use a similar process based on creative telescoping. The situation is slightly more complicated, since it is necessary to obtain a system of recurrences satisfied by the summand first.

Indefinite sum inside

We start with the inner summand, for the indefinite sum of which we need to find description making creative telescoping possible:

 > d[n,m]:=(-1)^(m+1)/2/m^3/binomial(n,m)/binomial(n+m,m);

We first get operators by definition of the indefinite sum:

 > Mgfun[dfinite_expr_to_sys](d[n,m],f(n::shift,m::shift));

 > firstoperators:=subs(k=k+1,subs(m=k,eval(%,f=proc(n,m) g(n,m)-g(n,m-1) end)));

We check that these operators cancel our sequence when evaluated at a special point (here ):

 > eval(firstoperators,g=proc(N,K) Sum(subs(n=N,d[n,m]),m=1..K) end):

 > eval(%,n=10):

 > value(subs(k=6,%));

These operators do not generate a D-finite ideal

In order to ensure the success of creative telescoping, we need to compute "sufficiently many" operators annihilating our sequence. This "sufficiently many" is formally defined in terms of the dimension of the ideal our operators generate. A D-finite ideal is an ideal such that the quotient of the algebra of operators with the ideal is finitely dimensional as a vector space. Informally, this amounts to saying that all the shifts of a sequence annihilated by the operators can be rewritten as linear combination of a finite number of them. Technically, D-finiteness can be checked by computing the Hilbert dimension of the ideal, which has to be 0.

 > sys:=subs([seq(seq(g(n+i,k+j)=Sn^i*Sk^j,j=0..2),i=0..1)],firstoperators);

 > Groebner[hilbertdim](sys,Groebner[termorder](Ore_algebra[shift_algebra]([Sn,n],[Sk,k]),tdeg(Sn,Sk)));

This would be zero in the case of a D-finite system.

More generators are obtained by creative telescoping:

 > CT2:=Mgfun[creative_telescoping](d[n,m],n::shift,m::shift);

Then, in order to sum  for m from 1 to k, we first find annihilators of the right-hand side:

 > Mgfun[dfinite_expr_to_sys](subs(_f(n,m)=d[n,m],m=k+1,CT2[2])-subs(_f(n,m)=d[n,m],m=1,CT2[2]),f(n::shift,k::shift));

and then we get annihilators for the indefinite sum by substitution in the left-hand side:

 > ann:=eval(%,f=proc(N,K) unapply(CT2[1],n,m)(N,K) end);

Again, we check that these operators are correct at :

 > eval(ann,_f=proc(N,K) Sum(subs(n=N,d[n,m]),m=1..K) end):

 > subs(n=10,%):

 > value(subs(k=6,%));

Finally, we put both sets of generators together:

 > ann2:=map(collect,ann union subs(g=_f,firstoperators),_f,factor);

Now, this is D-finite

 > sys:=subs([seq(seq(_f(n+i,k+j)=Sn^i*Sk^j,j=0..2),i=0..3)],ann2);

 > Groebner[hilbertdim](sys,Groebner[termorder](Ore_algebra[shift_algebra]([Sn,n],[Sk,k]),tdeg(Sn,Sk)));

Building up the summand

From the (D-finite) system of operators annihilating the indefinite sum of , we are now going to construct a system of operators annihilating the summand in , before applying creative telescoping. This last step will be possible because addition and product of D-finite systems remain D-finite. Thus we now add  and multiply by . This is achieved by first constructing a system annihilating the sum:

 > inhom:=(n+1)^3*(_f(n+1,k)-_f(n,k))-1;

(the system is composed of this equation, which we render homogeneous, and another one reflecting that the sum does not depend on ). Then we construct the system satisfied by the sum:

 > Mgfun[`sys+sys`]({_f(n,k+1)-_f(n,k),subs(n=n+1,inhom)-inhom},ann2);

and similarly we multiply by :

 > Mgfun[`sys*sys`](%,Mgfun[dfinite_expr_to_sys](c[n,k],_f(n::shift,k::shift)));

Creative telescoping

This is the time-consuming part of this worksheet: we now compute a recurrence satisfied by  by applying creative telescoping to the summand defined by the system we have just computed.

The following instruction triggers the display of extra information during the computation:

 > infolevel[Mgfun]:=5:

 > res:=Mgfun[creative_telescoping](LFSol(%%),n::shift,k::shift);

`Mgfun/chyzak97:   "Suitable term order guessed"`

`Mgfun/chyzak97:   "Implicit d-finite expression recognized"`

`Mgfun/chyzak97:   "Dimension is 3"`

`Mgfun/chyzak97:   "Start of Chyzak's algorithm"`

`Mgfun/chyzak97:   "Preparation of the system: .60e-1 seconds."`

`Mgfun/uncoupling:   "Uncoupling of the LOF system"`

`Mgfun/chyzak97:   "Uncoupling of the system: 1.011 seconds."`

`Mgfun/rational_sys_solve:   "Look for rational solutions of a system"`

``Mgfun/denominator_bound`[shift]   "Intermediate bound on dispersion: infinity"`

``Mgfun/denominator_bound`[shift]   "Bound on dispersion: 0"`

``Mgfun/denominator_bound`[shift]   "Computing a resultant: .11e-1 seconds."`

``Mgfun/rational_solve`[shift]   "Computing denominator bound: .20e-1 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 1"`

`Mgfun/rational_sys_solve:   "Solving equation: .580 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 2"`

`Mgfun/rational_sys_solve:   "Solving equation: .10e-1 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 3"`

`Mgfun/rational_sys_solve:   "Solving equation: 0. seconds."`

`Mgfun/chyzak97:   "No LOFE of order 0 :-("`

`Mgfun/chyzak97:   "Preparation of the system: .70e-1 seconds."`

`Mgfun/uncoupling:   "Uncoupling of the LOF system"`

`Mgfun/chyzak97:   "Uncoupling of the system: 1.160 seconds."`

`Mgfun/rational_sys_solve:   "Look for rational solutions of a system"`

``Mgfun/denominator_bound`[shift]   "Intermediate bound on dispersion: infinity"`

``Mgfun/denominator_bound`[shift]   "Bound on dispersion: 0"`

``Mgfun/denominator_bound`[shift]   "Computing a resultant: .19e-1 seconds."`

``Mgfun/rational_solve`[shift]   "Computing denominator bound: .39e-1 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 1"`

`Mgfun/rational_sys_solve:   "Solving equation: 1.210 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 2"`

`Mgfun/rational_sys_solve:   "Solving equation: 0. seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 3"`

`Mgfun/rational_sys_solve:   "Solving equation: 0. seconds."`

`Mgfun/chyzak97:   "No LOFE of order 1 :-("`

`Mgfun/chyzak97:   "Preparation of the system: .350 seconds."`

`Mgfun/uncoupling:   "Uncoupling of the LOF system"`

`Mgfun/chyzak97:   "Uncoupling of the system: 3.360 seconds."`

`Mgfun/rational_sys_solve:   "Look for rational solutions of a system"`

``Mgfun/denominator_bound`[shift]   "Intermediate bound on dispersion: infinity"`

``Mgfun/denominator_bound`[shift]   "Bound on dispersion: 0"`

``Mgfun/denominator_bound`[shift]   "Computing a resultant: .60e-1 seconds."`

``Mgfun/rational_solve`[shift]   "Computing denominator bound: .119 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 1"`

`Mgfun/rational_sys_solve:   "Solving equation: 9.260 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 2"`

`Mgfun/rational_sys_solve:   "Solving equation: .10e-1 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 3"`

`Mgfun/rational_sys_solve:   "Solving equation: 0. seconds."`

`Mgfun/chyzak97:   "No LOFE of order 2 :-("`

`Mgfun/chyzak97:   "Preparation of the system: .869 seconds."`

`Mgfun/uncoupling:   "Uncoupling of the LOF system"`

`Mgfun/chyzak97:   "Uncoupling of the system: 2.831 seconds."`

`Mgfun/rational_sys_solve:   "Look for rational solutions of a system"`

``Mgfun/denominator_bound`[shift]   "Intermediate bound on dispersion: infinity"`

``Mgfun/denominator_bound`[shift]   "Bound on dispersion: 0"`

``Mgfun/denominator_bound`[shift]   "Computing a resultant: .11e-1 seconds."`

``Mgfun/rational_solve`[shift]   "Computing denominator bound: .11e-1 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 1"`

`Mgfun/rational_sys_solve:   "Solving equation: 46.491 seconds."`

``Mgfun/denominator_bound`[shift]   "Intermediate bound on dispersion: infinity"`

``Mgfun/denominator_bound`[shift]   "Bound on dispersion: 0"`

``Mgfun/denominator_bound`[shift]   "Computing a resultant: .401 seconds."`

``Mgfun/rational_solve`[shift]   "Computing denominator bound: .511 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 2"`

`Mgfun/rational_sys_solve:   "Solving equation: 4.259 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 3"`

`Mgfun/rational_sys_solve:   "Solving equation: 0. seconds."`

`Mgfun/chyzak97:   "No LOFE of order 3 :-("`

`Mgfun/chyzak97:   "Preparation of the system: 1.980 seconds."`

`Mgfun/uncoupling:   "Uncoupling of the LOF system"`

`Mgfun/chyzak97:   "Uncoupling of the system: 6.250 seconds."`

`Mgfun/rational_sys_solve:   "Look for rational solutions of a system"`

``Mgfun/denominator_bound`[shift]   "Intermediate bound on dispersion: infinity"`

``Mgfun/denominator_bound`[shift]   "Bound on dispersion: 0"`

``Mgfun/denominator_bound`[shift]   "Computing a resultant: .120 seconds."`

``Mgfun/rational_solve`[shift]   "Computing denominator bound: .131 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 1"`

`Mgfun/rational_sys_solve:   "Solving equation: 584.509 seconds."`

``Mgfun/denominator_bound`[shift]   "Intermediate bound on dispersion: infinity"`

``Mgfun/denominator_bound`[shift]   "Bound on dispersion: 0"`

``Mgfun/denominator_bound`[shift]   "Computing a resultant: .609 seconds."`

``Mgfun/rational_solve`[shift]   "Computing denominator bound: .900 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 2"`

`Mgfun/rational_sys_solve:   "Solving equation: 5801.681 seconds."`

`Mgfun/rational_sys_solve:   "Equations solved: 3"`

`Mgfun/rational_sys_solve:   "Solving equation: 479.559 seconds."`

`Mgfun/chyzak97:   "LOFE of order 4 found!"`

satisfies the same recurrence as

As we did for , we now sum this for k from 0 to n+4 and we observe the following:

The right-hand side telescopes to 0

 > R:=res[2]:

 > eval(R,k=0);

 > F:=c[n,k]*(Sum(1/m^3,m=1..n)+Sum(d[n,m],m=1..k));

 > ff:=eval(R,_f=proc(N,K) subs(n=N,k=K,F) end):

 > expand(eval(ff,k=n+5));

Thus, we have found the following 4th order recurrence satisfied by :

 > eval(collect(res[1],_f,factor),_f=proc(n,k) A(n) end);

 > rec2:=%:

Check:

 > seq(eval(subs(n=i,rec2),A=proc(n) listb[n] end),i=1..6);

Now, we observe using gfun   that all solutions of the second order recurrence rec found above to be satisfied by  are also solutions of this recurrence:

 > expand(gfun[`rec+rec`](rec,rec2,A(n))-rec2);

This shows that rec2 is also satisfied by all linear combinations of solutions of rec and rec2.

Therefore, to prove that  satisfies rec, it is sufficient to check that its first 4 values satisfy it, which has been done above.

4. Asymptotic behaviour of :=

We start from the recurrence satisfied by ,  and  defined in the title of this section:

 > rec;

and consider the differential equation satisfied by generating functions of its solutions

 > gfun[rectodiffeq](rec,A(n),y(z));

we make this homogeneous:

 > deq:=collect(gfun[diffeqtohomdiffeq](op(remove(type,%,`=`)),y(z))/(5*_C[0]-_C[1]),diff,factor);

The solutions of this equation (and therefore the generating functions ,  and  of the sequences ,  and ) may have singularities only at 0, , and any of the roots of the leading coefficient of this equation:

 > polsing:=coeff(deq,diff(y(z),[z\$4]))/z^2;

Here are the roots of this polynomial

 > solve(polsing,z);

and the corresponding numerical values

 > evalf([%]);

First, it should be observed (by simple bounds on the binomial coefficients) that the sequences do not grow so fast as to make their generating series divergent. Therefore 0 is not a singularity of the functions we're interested in. Secondly, the sequence  is clearly increasing, which implies that  has a singularity of modulus smaller than 1, which has to be = . Since /  has a finite limit, this is also a singularity of . A crucial step of the proof is to observe that  is not  singular at that point and therefore its singularity of smallest modulus is 1/ = , whence an exponential growth of   in . This step is achieved by considering the local behaviour of the solutions of the differential equation at :

 > alias(alpha=RootOf(polsing)):DEtools[formal_sol](deq,y(z),T,z=alpha,terms=4);

The interpretation of this result is as follows: at , the differential equation has a vector space of dimension 4 of solutions, which is the direct sum of a subspace of analytic solutions of dimension 3 and a "line" of singular solutions, generated by a solution . Thus , where  is analytic in the neighborhood of  and similarly . Now, singularity analysis tells us that the coefficients  behave asymptotically like

 > equivalent(lambda[A]*sqrt(alpha-z),z,n);

and similarly for  with  in place of . Now, we deduce that lim( / )= = /  and consequently =  is analytic at , as was to be proved.

5.

The proof is obtained by writing this difference as the infinite sum

A recurrence satisfied by the numerator of the summand is readily found using gfun:

 > gfun[poltorec](A(k)*B(k+1)-A(k+1)*B(k),[subs(n=k,rec),subs(A=B,n=k,rec)],[A(k),B(k)],u(k));

 > rsolve(%,u(k));

and thus this is a sum of positive terms since  and the initial condition is positive:

 > value(subs(n=0,a[n])*subs(n=1,b[n])-subs(n=1,a[n])*subs(n=0,b[n]));

6. Conclusion

It is classical that the sequence =lcm(1,..., ) grows asymptotically like . This worksheet has thus proved that the sequence =  enjoys the following properties: it is positive and grows asymptotically like .  Now, this tends to 0:

 > evalf(exp(3)*(17-12*sqrt(2)));

 >

If  were rational, then for  large enough,  would be an integer and the sequence  would be a sequence of positive integers tending to 0, which is impossible. This concludes the proof. Further refinements lead to irrationality measures.