SOME NEW FORMULAS FOR π

Gert Almkvist, Christian Krattenthaler^{†}^{∗}^{§} and Joakim Petersson

Matematiska Institutionen, Lunds Universitet, SE-221 00 Lund, Sweden.

     E-mail: Gert.Almkvist@math.lu.se, Joakim.Petersson@math.lu.se

†Institut für Mathematik der Universität Wien, Strudlhofgasse 4, A-1090 Wien, Austria.

                     E-mail: KRATT@Ap.Univie.Ac.At
              WWW: http://www.mat.univie.ac.at/~kratt

Abstract. We show how to find series expansions for π of the form π = ∑_{n=0}^{∞} S(n)/\binom{mn}{pn}a^{n}, where S(n) is some polynomial in n (depending on m, p, a). We prove that there exist such expansions for m = 8k, p = 4k, a = (−4)^{k}, for any k, and give explicit examples for such expansions for small values of m, p and a.

1. Introduction

Using the formula (due to Bill Gosper [3])

π = ∑_{n=0}^{∞} \frac{50n − 6}{\binom{3n}{n}2^{n}},

Fabrice Bellard [1, file pi1.c] found an algorithm for computing the n-th decimal of π without computing the earlier ones. Thus he improved an earlier algorithm due to Simon Plouffe [6].

This formula can be proved in the following way. We have

\frac{1}{\binom{3n}{n}} = (3n + 1) ∫_{0}^{1} x^{2n}(1 − x)^{n}dx.

∗Research partially supported by the Austrian Science Foundation FWF, grant P13190-MAT, and by EC’s IHRP Programme, grant HPRN-CT-2001-00272.

§Current Address: Institut Girard Desargues, Université Claude Bernard Lyon-I, 21, avenue Claude Bernard, F-69622 Villeurbanne Cedex, France.

SOME NEW FORMULAS FOR π 2

Hence the right hand side of the formula will be

∫_{0}^{1} ∑_{n=0}^{∞}(50n − 6)(3n + 1) (\frac{x^{2}(1 − x)}{2})^{n} dx.

But

∑_{n=0}^{∞}(50n − 6)(3n + 1)y^{n} = \frac{2(56y^{2} + 97y − 3)}{(1 − y)^{3}},

so we get

RHS = 8 ∫_{0}^{1} \frac{28x^{6} − 56x^{5} + 28x^{4} − 97x^{3} + 97x^{2} − 6}{(x^{3} − x^{2} + 2)^{3}}dx

= [\frac{4x(x − 1)(x^{3} − 28x^{2} + 9x + 8)}{(x^{3} − x^{2} + 2)^{2}} + 4 arctan(x − 1)]_{0}^{1} = π.

We then asked if there are other such formulas of type

π = ∑_{n=0}^{∞} \frac{S(n)}{\binom{mn}{pn}a^{n}} ,

where S(n) is a polynomial in n. Using the same trick as before we have

\frac{1}{\binom{mn}{pn}} = (mn + 1) ∫_{0}^{1} x^{pn}(1 − x)^{(m−p)n} dx,

and hence

RHS = ∫_{0}^{1} ∑_{n=0}^{∞}(mn + 1)S(n) (\frac{x^{p}(1 − x)^{m−p}}{a} )^{n} dx.  (1.1)

If S(n) has degree d then

∑_{n=0}^{∞}(mn + 1)S(n)y^{n} = \frac{T (y)}{(1 − y)^{d+2}} ,  (1.2)

where T has degree d + 1. It follows with

y = \frac{x^{p}(1 − x)^{m−p}}{a}

that

RHS = ∫_{0}^{1} \frac{P (x)}{(x^{p}(1 − x)^{m−p} − a)^{d+2}} dx,  (1.3)

where P (x) is a polynomial in x of degree m(d+1). We want this integral to be equal to π. A good way to get π is to have arctan(x) or arctan(x − 1) after integration. This means that

x^{p}(1 − x)^{m−p} − a

must have the factor x^{2}+1 or (x−1)^{2}+1, that is have a zero at i or 1+i. This restricts m and p and gives the value of a. After experimenting with the LLL-algorithm we found formulas for π in the following cases (there could be many more):

m p  a deg(S)
3 1 2  1 Gosper [3]
7 2 2  5 Bellard [1]^{1}
8 4 −4  4
10 4 4  8
12 4 −4  8
16 8 16  8
24 12 −64  12
32 16 256 16
40 20 −4^{5}  20
48 24 4^{6}  24
56 28 −4^{7}  28
64 32 4^{8}  32
72 36 −4^{9}  36
80 40 4^{10}  40

E.g., to find the formula in the last case we computed

s(k) = ∑_{n=0}^{∞} \frac{n^{k}}{\binom{80n}{40n}4^{10n}}

for k = 0, 1, 2, …, 40 with 6000 digits. Then the LLL-algorithm found (in about three days) the linear combination between π and s(0), s(1), …, s(40). A good check is that the coefficient of π is a product of small primes (actually all less than 80).

1 This formula was used in Bellard’s world record setting computation of the 1000 billionth binary digit of π, being based on the algorithm in [2].

2. Proving the formulas

For m ⩽ 16 the integrals can be computed by brute force using Maple. The higher cases are all symmetric, i.e., m = 2p. Using the symmetry

x ⟷ 1 − x

one can assume that arctan(x) and arctan(x − 1) each contribute π/2. So we make the wild assumption that

∫ \frac{P (x)}{Q(x)^{d+2}} dx = \frac{R(x)}{Q(x)^{d+1}} + 2 arctan(x) + 2 arctan(x − 1),

where R(x) is a polynomial with R(0) = R(1) = 0. Differentiation with respect to x yields

\frac{P}{Q^{d+2}} = \frac{R′}{Q^{d+1}} − (d + 1)\frac{Q′R}{Q^{d+2}} + 2 (\frac{1}{x^{2} + 1} + \frac{1}{x^{2} − 2x + 2}) ,

or

QR′ − (d + 1)Q′R = P − 2Q^{d+2} (\frac{1}{x^{2} + 1} + \frac{1}{x^{2} − 2x + 2}) .

This is a differential equation of first order in R with polynomial coefficients. Taking the case \binom{32n}{16n}, we have deg(P) = 544 and d = 16, so the denominator Q^{d+2} has degree 576. We solve the differential equation using Maple’s

                       dsolve(*,R(0)=0,series)

setting the order to 600. Then we find that R is a polynomial of degree 543. In practice one replaces the coefficients of \frac{1}{x^{2}+1} and \frac{1}{x^{2}−2x+2} by 2 times the denominator in the formula found for π. Then R will have (huge) integer coefficients.

Remark 1. Computing the integral in the case \binom{8n}{4n}, Maple got the terms

arctan (\frac{2x^{3} − 3x^{2} + 7x − 3}{5}) + arctan (\frac{x}{2} − \frac{1}{4}) ,

which can be shown to be equal to

arctan(x) + arctan(x − 1).

In the nonsymmetric case \binom{12n}{4n} one gets arctan(x − 1) and arctan(x^{3} − 2x^{2} + x − 1), but the latter term makes no contribution since it has the same values at 0 and 1.

Remark 2. When we set this up then, in principle, we could let Maple do the evaluation of T (y) (and thus of P (x)) through (1.2), i.e., let Maple’s formal summation routines compute the sum

∑_{n=0}^{∞}(mn + 1)S(n)y^{n},

or, to make the point clearly visible, compute a sum of the form

∑_{n=0}^{∞} S̃(n)y^{n},

where S̃(n) is a certain polynomial in n. However, it turns out that this is very time-consuming. It is a much better idea to expand S̃(n) as S̃(n) = ∑_{j=0}^{d+1} a(j)\binom{n+j}{d+1}, and then use that ∑_{n=0}^{∞} \binom{n+j}{d+1}y^{n} = y^{d+1−j}/(1 − y)^{d+2} for j = 0, 1, …, d + 1. See also the approach that we follow in Section 4 when we prove our main result.

3. The case \binom{8kn}{4kn}, k = 1, 2, 3, …

Solving the differential equation for R(x) one finds that R(x) has the factor x(1−x)(2x−1). If one divides out the factor 2x−1, the remaining factor is invariant under the substitution x ⟷ 1 − x. Hence we make the substitution

R(x) = (2x − 1)Ř(z),

where

z = x(1 − x).

Then

\frac{dz}{dx} = 1 − 2x and (1 − 2x)^{2} = 1 − 4z.

Now

Q = z^{4k} − (−4)^{k},

so

Q′ = 4kz^{4k−1}(1 − 2x)

and

R′ = 2Ř(z) + (2x − 1)\frac{dŘ}{dz}(1 − 2x) = 2Ř − (1 − 4z)\frac{dŘ}{dz}.

It follows that the left hand side of the differential equation is

QR′ − (4k + 1)Q′R = Q (2Ř − (1 − 4z)\frac{dŘ}{dz}) + (4k + 1)4kz^{4k−1}(1 − 2x)^{2}Ř

= (2Q + 4k(4k + 1)(1 − 4z)z^{4k−1})Ř − (1 − 4z)Q\frac{dŘ}{dz}.

On the right hand side we have

\frac{1}{x^{2} + 1} + \frac{1}{x^{2} − 2x + 2} = \frac{3 − 2z}{z^{2} − 2z + 2},

and Q is for all k divisible by

z^{4} + 4 = (z^{2} − 2z + 2)(z^{2} + 2z + 2).

Hence P and Q and

Q (\frac{1}{x^{2} + 1} + \frac{1}{x^{2} − 2x + 2})

are all polynomials in z. Finally we get the following differential equation for Ř:

−(1−4z)Q\frac{dŘ}{dz} +(2Q+4k(4k+1)(1−4z)z^{4k−1})Ř−P +2(3−2z)\frac{Q^{4k+2}}{z^{2} − 2z + 2} = 0. (3.1)

If P is known then one solves this equation for Ř just as before (but the degrees are cut in half).

But there is a possibility to find the formula for π and prove it in one stroke. Let

N(k) = 4k(4k + 1),

and assume that

Ř(z) = ∑_{j=1}^{N(k)−1} a(j)z^{j}  (3.2)

and

π = ∑_{n=0}^{∞} \frac{S_{k}(n)}{\binom{8kn}{4kn}(−4)^{kn}} ,  (3.3)

with

S_{k}(n) = ∑_{j=0}^{4k} a(N(k) + j)n^{j}.  (3.4)

Substituting this into the differential equation we will get a system of linear equations for the a(j)’s of size (N(k) + 4k) × (N(k) + 4k). If this system is nonsingular, we can be sure to be able to solve it and thus find and prove an expansion for π of the form (3.3).

We went (again) to the computer and generated the system of equations for small values of k. Let A(k) be the matrix of the system. Let further r(k) be the least common denominator of the coefficients of S_{k}(n). We obtained the following tables:

\begin{matrix}k & det(A(k)) \\ 1 & 2^{91}3^{8}5^{7}7^{2} \\ 2 & −2^{523}3^{52}5^{17}7^{14}11^{4}13^{3} \\ 3 & 2^{1367}3^{177}5^{41}7^{25}11^{20}13^{19}17^{5}19^{4}23^{2} \\ 4 & −2^{3231}3^{167}5^{83}7^{53}11^{28}13^{27}17^{25}19^{8}23^{6}29^{3}31^{2} \\ 5 & 2^{5399}3^{290}5^{345}7^{93}11^{41}13^{37}17^{33}19^{32}23^{10}29^{7}31^{6}37^{3} \end{matrix}

\begin{matrix}k & r(k) \\ 1 & 3^{2}5^{2}7^{2} \\ 2 & 3^{6}5^{3}7^{2}11^{2}13^{2} \\ 3 & 2^{5}3^{3}5^{2}7^{2}11^{2}13^{2}17^{2}19^{2}23^{2} \\ 4 & 2^{3}3^{10}5^{6}7^{3}11 · 13^{2}17^{2}19^{2}23^{2}29^{2}31^{2} \\ 5 & 2^{11}3^{13}5^{3}7^{4}11^{3}13^{3}17^{2}19^{2}23^{2}29^{2}31^{2}37^{2} \\ 6 & 2^{10}3^{7}5^{8}7^{3}11^{3}13^{3}17^{2}19^{2}23^{2}29^{2}31^{2}37^{2}41^{2}43^{2}47^{2} \\ 7 & 2^{17}3^{16}5^{9}7^{3}11^{4}13^{3}17^{3}19^{2}23^{2}29^{2}31^{2}37^{2}41^{2}43^{2}47^{2}53^{2} \\ 8 & 2^{10}3^{19}5^{8}7^{8}11^{4}13^{3}17^{3}19^{3}23^{2}29^{2}31^{2}37^{2}41^{2}43^{2}47^{2}53^{2}59^{2}61^{2} \\ 9 & 2^{23}3^{4}5^{10}7^{8}11^{4}13^{4}17^{3}19^{3}23^{3}29^{2}31^{2}37^{2}41^{2}43^{2}47^{2}53^{2}59^{2}61^{2}67^{2}71^{2} \\ 10 & 2^{20}3^{20}5^{4}7^{9}11^{5}13^{4}17^{3}19^{3}23^{2}29^{2}31^{2}37^{2}41^{2}43^{2}47^{2}53^{2}59^{2}61^{2}67^{2}71^{2}73^{2}79^{2} \end{matrix}

From the table it seems “obvious” that the determinant of A(k) will never vanish as it grows rather quickly in absolute value with k. We also see that the same primes occur in the factorizations of r(k) and det(A(k)). Even more striking, it seems that the largest prime factors occurring grow only slowly (namely approximately linearly) when k increases. The last fact strongly indicates that there may even be a closed form formula for det(A(k)). As it turns out, this is indeed the case. In the next section we will explicitly compute the determinant of an equivalent system of equations, from which it follows that (see the remark after Theorem 6)

det(A(k)) = (−1)^{k−1}2^{32k^{3}+24k^{2}+2k−1}k^{8k^{2}+2k}((4k + 1)!)^{4k} \frac{(8k)!}{(4k)!} ∏_{j=1}^{4k} \frac{(2j)!}{j!}.  (3.5)

Hence we have the following theorem, which is the main theorem of our paper.

Theorem 1. For all k ≥ 1 there is a formula

π = ∑_{n=0}^{∞} \frac{S_{k}(n)}{\binom{8kn}{4kn}(−4)^{kn}} ,

where S_{k}(n) is a polynomial in n of degree 4k with rational coefficients. The polynomial S_{k}(n) can be found by solving the system of linear equations generated by (3.1) and the Ansatz (3.2) and (3.4).

The denominators of all a(j)’s divide r(k) which is much smaller than det(A(k)). This means that there must be some miracle occurring at the end when solving the system. E.g., when k = 5 then det(A(5)) has about 2400 digits but r(5) only 40. Unfortunately, we are not able to offer an explanation for that.

In practice we are only interested in the coefficients of S_{k} so we try to eliminate a(1), a(2), …, a(N(k) − 1) first. This can be done by first avoiding all equations containing a(N(k)), …, a(N(k)+4k), i.e., the equations coming from the coefficients of t^{4kv}, v = 0, 1, …, 4k + 1. This looks very nice theoretically, in particular as the system for a(1), …, a(N(k)) is triangular, but in practice the computer breaks down since the rational numbers occuring become very large.

We close this section by listing a few explicit examples.

Example 2. We have

π = \frac{1}{r} ∑_{n=0}^{∞} \frac{S(n)}{\binom{8n}{4n}(−4)^{n}} ,

where

r = 3^{2}5^{2}7^{2}

and

S(n) = −89286 + 3875948n − 34970134n^{2} + 110202472n^{3} − 115193600n^{4}.

Example 3. We have

π = \frac{1}{r} ∑_{n=0}^{∞} \frac{S(n)}{\binom{16n}{8n}16^{n}} ,

where

r = 3^{6}5^{3}7^{2}11^{2}13^{2}

and

S(n) = −869897157255 − 3524219363487888n + 112466777263118189n^{2} − 1242789726208374386n^{3} + 6693196178751930680n^{4} − 19768094496651298112n^{5} +32808347163463348736n^{6}−28892659596072587264n^{7}+10530503748472012800n^{8}.

SOME NEW FORMULAS FOR π 9

Example 4. We have

π = \frac{1}{r} ∑_{n=0}^{∞} \frac{S(n)}{\binom{32n}{16n}256^{n}} ,

where

r = 2^{3}3^{10}5^{6}7^{3}11 · 13^{2}17^{2}19^{2}23^{2}29^{2}31^{2}

and

S(n) = − 2062111884756347479085709280875

+ 1505491740302839023753569717261882091900n − 112401149404087658213839386716211975291975n^{2} + 3257881651942682891818557726225840674110002n^{3} − 51677309510890630500607898599463036267961280n^{4} + 517337977987354819322786909541179043148522720n^{5} − 3526396494329560718758086392841258152390245120n^{6} + 171145766235995166227501216110074805943799363584n^{7} − 60739416613228219940886539658145904402068029440n^{8} + 159935882563435860391195903248596461569183580160n^{9} − 313951952615028230229958218839819183812205608960n^{10} + 457341091673257198565533286493831205566468325376n^{11} − 486846784774707448105420279985074159657397780480n^{12} + 367314505118245777241612044490633887668208926720n^{13} − 185647326591648164598342857319777582801297080320n^{14} + 56224688035707015687999128994324690418467340288n^{15} − 7687255778816557786073977795149360408612044800n^{16}.

Example 5. (nonsymmetric). We have

π = \frac{1}{r} ∑_{n=0}^{∞} \frac{S(n)}{\binom{10n}{4n}4^{n}} ,

where

r = 3^{2}5 · 7^{2}11 · 13 · 17 · 19

and

S(n) = −4843934523072 − 1008341177146848n + 23756198610834352n^{2} − 242873913552020704n^{3} + 1195813551184400032n^{4} − 3272960363556054592n^{5} + 4909379167837011328n^{6} − 3816399750842818816n^{7} + 1190182007407360000n^{8}.

4. Proof of the Theorem

We want to prove that, by making the Ansatz (3.2) and (3.4) and substituting this into the differential equation (3.1) (the polynomial P being given by S_{k} through (1.1)– (1.3), Q being given by (x(1 − x))^{4k} − (−4)^{k}), the resulting system of linear equations will always have a solution. In fact, we aim at finding an explicit formula for the determinant of the corresponding matrix of coefficients that allows us to conclude that it can never vanish.

It turns out that for that purpose it is more convenient to set up the system of linear equations in a different, but equivalent way. This equivalent system will have coefficient matrix M (see (4.4)). The evaluation of its determinant will be accomplished through Eqs. (4.6), (4.7), (4.8), (4.10), and Theorem 6.

To be precise, we encode the polynomial S_{k}(x) (and, thus, T (y) and P (x)) differently. We claim that T (y) has an expansion of the form

T (y) = ∑_{j=0}^{4k+1} c_{j}y^{j},  (4.1)

subject to the single constraint

∏_{i=1}^{4k+1}(4ik − 1)c_{0} + ∑_{j=1}^{4k+1}(−1)^{j} ( ∏_{i=1}^{4k+1−j}(4ik − 1)) (∏_{i=1}^{j−1}(4ik + 1)) c_{j} = 0.  (4.2)

(As usual, empty poducts have to be interpreted as 1.) This is seen as follows. The polynomial S_{k}(x) can be written in the form

S_{k}(n) = ∑_{j=0}^{4k} \binom{n + j}{4k}s_{j},

for some coefficients s_{j}. Hence, we have

\frac{T (y)}{(1 − y)^{4k+2}} = ∑_{n=0}^{∞}(4kn + 1)S_{k}(n)y^{n}

= ∑_{n=0}^{∞}(4kn + 1) ∑_{j=0}^{4k} \binom{n + j}{4k}s_{j}y^{n}

= ∑_{j=0}^{4k} s_{j} ∑_{n=0}^{∞}(4k(n + j + 1) − 4k(j + 1) + 1)\binom{n + j}{4k}y^{n}

= ∑_{j=0}^{4k} s_{j} ∑_{n=0}^{∞} (4k(4k + 1)\binom{n + j + 1}{4k + 1}y^{n} − (4k(j + 1) − 1)\binom{n + j}{4k}y^{n})

= ∑_{j=0}^{4k} s_{j} (4k(4k + 1)\frac{y^{4k−j}}{(1 − y)^{4k+2}} − (4k(j + 1) − 1)\frac{y^{4k−j}}{(1 − y)^{4k+1}} )

= \frac{1}{(1 − y)^{4k+2}} ∑_{j=0}^{4k+1} y^{j} ((4jk + 1)s_{4k−j} + (4k(4k − j + 2) − 1)s_{4k−j+1}) .

In the last line, s_{−1} and s_{4k+1} have to be read as 0. It is now a trivial exercise to substitute the coefficients of y^{j} in the sum in the last line into the left-hand side of (4.2) and verify the truth of (4.2).

The above implies that

P (z) = ∑_{j=0}^{4k+1} c_{j}z^{4jk}/(−4)^{jk},

where the coefficients c_{j} obey (4.2).

Now we are ready to set up the system of linear equations. We make again the Ansatz (3.2), but we replace (3.4) by

P (z) = ∑_{j=0}^{4k+1} a(N(k) + j)z^{4jk}/(−4)^{jk},  (4.3)

where the a(N(k)+j), j = 0, 1, …, 4k +1, are subject to (4.2) (i.e., the relation (4.2) holds when c_{j} is replaced by a(N(k) + j)). Clearly, we have to add (4.2) to the set of equations that result from the differential equation (3.1).

The coefficient matrix of the system looks a follows:

M = (\begin{matrix}x & y \\ U & V \end{matrix} ) ,  (4.4)

where x is a line vector of N(k) − 1 zeroes, y = (y_{0}, y_{1}, …, y_{4k+1}) is the vector of coefficients of (4.2), i.e., y_{0} = ∏_{i=1}^{4k+1}(4ik − 1), and

y_{ℓ} = (−1)^{ℓ} ( ∏_{i=1}^{4k+1−ℓ}(4ik − 1)) (∏_{i=1}^{ℓ−1}(4ik + 1)) ,

ℓ = 1, 2, …, 4k + 1, where U is an (N(k) + 4k) × (N(k) − 1) matrix and V is an (N(k) + 4k) × (4k + 2) matrix, both of which we define below.

We consider the top-most line of M (which is formed out of x and y) as row 0 of M. We label the rows of U and V by i running from 1 to N(k) + 4k. Furthermore, we label the columns of M by j running from 1 to N(k) + 4k + 1.

Following this labelling scheme, the matrix U has nonzero entries only in the four diagonals i = j, i = j + 1, i = j + 4k, i = j + 4k + 1. We denote the entries in column j on these four diagonals in order f_{0}(j), f_{1}(j), g_{0}(j), and g_{1}(j), where

f_{0}(j) = j(−4)^{k},
f_{1}(j) = −(4j + 2)(−4)^{k},
g_{0}(j) = (N(k) − j),
g_{1}(j) = −(4N(k) − 4j − 2).

To be precise, the (j, j)-entry is f_{0}(j), the (j +1, j)-entry is f_{1}(j), the (j +4k, j)-entry is g_{0}(j), the (j + 4k + 1, j)-entry is g_{1}(j), j = 1, 2, …, N(k) − 1.

On the other hand, the matrix V is composed out of columns, labelled N(k), N(k) + 1, …, N(k) + 4k + 1, each of which containing just one nonzero entry. To be precise, the nonzero entry of column N(k) + j is located in the (4jk + 1)-st row (according to our labelling scheme), and it is equal to (−4)^{−jk}, j = 0, 1, …, 4k + 1.

We will now compute the determinant of M and show that it does not vanish for any k.

We perform some row operations on M, with the effect that the entries of y get eliminated. This is achieved by subtracting ∏_{i=1}^{4k+1}(4ik − 1) times row 1 from row 0, and

(−1)^{ℓ}(−4)^{ℓk} ( ∏_{i=1}^{4k+1−ℓ}(4ik − 1)) (∏_{i=1}^{ℓ−1}(4ik + 1))

times row 4ℓk+1 from row 0, ℓ = 1, 2, …, 4k+1. Doing this, we must expect changes in row 0 in columns 1, 4k, 4k + 1, 8k, 8k + 1, …,N(k) − 4k = 16k^{2}. However, at this point a miracle occurs: the new entries in row 0 in columns 4ℓk+1, ℓ = 0, 1, …, 4k−1, are still 0. On the other hand, the values of the new entries in row 0 in columns 4ℓk, ℓ = 1, 2, …, 4k are

(−1)^{ℓ−1}(−4)^{(ℓ+1)k}8k(4k + 1) ( ∏_{i=1}^{4k−ℓ}(4ik − 1)) (∏_{i=1}^{ℓ−1}(4ik + 1)) .  (4.5)

After these manipulations we obtain a matrix of the form

M′ = (\begin{matrix}x′ & y′ \\ U & V \end{matrix}) ,

where x′ and y′ have the same dimensions as before x and y, respectively, where the only nonzero entries of x′ are in columns labelled by numbers which are divisible by 4k, with the entry in column 4ℓk given by (4.5), and where y′ consists only of zeroes. We have

det M = det M′.  (4.6)

The next step consists in expanding the determinant of M′ with respect to columns N(k), N(k) + 1, …, N(k) + 4k + 1 (i.e., the last 4k + 2 columns). Since each of these columns contains just one nonzero entry (which is a power of (−4)^{−k}), we have

det M′ = ±(−4)^{−k(2k+1)(4k+1)} det M′′,  (4.7)

where M′′ is the matrix arising from M′ by deleting the last 4k + 2 columns and the rows 4ℓk + 1, ℓ = 0, 1, …, 4k + 1. More precisely, the matrix M′′ has the following form:

(\begin{matrix}x_{1} & x_{2} & x_{3} & . . . . . . . . & x_{4k+1} \\ F_{1} & 0 & 0 & . . . . . . . . & 0 \\ G_{1} & F_{2} & 0 & . . . . . . . . & 0 \\ 0 & G_{2} & F_{3} & . . . . . . . . & 0 \\ 0 & 0 & G_{3} & ... & ⋮ \\ & & ... & ... & 0 \\ 0 & G_{4k} & F_{4k+1} \\ … & 0 & G_{4k+1} \end{matrix}) ,

where x_{ℓ}, ℓ = 1, 2, …, 4k, is a line vector with 4k entries, all of them being zero except for the last, which is equal to (4.5), where x_{4k+1} is a line vector of 4k−1 zeroes, where F_{ℓ} and G_{ℓ}, ℓ = 1, 2, …, 4k, are (4k − 1) × (4k) matrices with nonzero entries only in the two main diagonals, and where F_{4k+1} and G_{4k+1} are (4k − 1) × (4k − 1) matrices, G_{4k+1} being upper triangular. To be precise, for ℓ = 1, 2, …, 4k we have

F_{ℓ} = ( \begin{matrix}f_{1}(4(ℓ−1)k+1) & f_{0}(4(ℓ−1)k+2) & 0 & ... \\ 0 & f_{1}(4(ℓ−1)k+2) & f_{0}(4(ℓ−1)k+3) & 0 & ... \\ & ... & ... \\ & & ... & ... \\ & & 0 & f_{1}(4ℓk−2) & f_{0}(4ℓk−1) & 0 \\ & & & 0 & f_{1}(4ℓk−1) & f_{0}(4ℓk) \end{matrix} )

and

G_{ℓ} = ( \begin{matrix}g_{1}(4(ℓ−1)k+1) & g_{0}(4(ℓ−1)k+2) & 0 & ... \\ 0 & g_{1}(4(ℓ−1)k+2) & g_{0}(4(ℓ−1)k+3) & 0 & ... \\ & ... & ... \\ & & ... & ... \\ & & 0 & g_{1}(4ℓk−2) & g_{0}(4ℓk−1) & 0 \\ & & & 0 & g_{1}(4ℓk−1) & g_{0}(4ℓk) \end{matrix} ) ,

and we have

G_{4k+1} = ( \begin{matrix}g_{1}(16k^{2}+1) & g_{0}(16k^{2}+2) & 0 & ... \\ 0 & g_{1}(16k^{2}+2) & g_{0}(16k^{2}+3) & ... \\ & ... & ... \\ & & ... & ... \\ & & & g_{1}(16k^{2}+4k−2) & g_{0}(16k^{2}+4k−1) \\ & & & 0 & g_{1}(16k^{2}+4k−1) \end{matrix} ) .

The precise form of F_{4k+1} is without relevance for us. We do a Laplace expansion with respect to the last 4k − 1 rows. Because of the triangular form of G_{4k+1} we obtain

det M′′ = ( ∏_{i=16k^{2}+1}^{16k^{2}+4k−1} g_{1}(i)) det M′′′,  (4.8)

where

M′′′ = (\begin{matrix}x_{1} & x_{2} & x_{3} & … & x_{4k} \\ F_{1} & 0 & 0 & … \\ G_{1} & F_{2} & 0 & … \\ 0 & G_{2} & F_{3} & … \\ 0 & 0 & G_{3} & ... & ⋮ \\ & ... & ... & ... & 0 \\ & & & G_{4k−1} & F_{4k} \\ & & & 0 & G_{4k} \end{matrix}) .

Instead of M′′′ we consider a more general matrix. Define the functions

f_{0}(t, j) = ((N(k) + j)Y_{t} − X_{2,t})(−4)^{k},
f_{1}(t, j) = −((4N(k) + 4j + 2)Y_{t} − 4X_{1,t})(−4)^{k},
g_{0}(t, j) = (X_{2,t} − jY_{t}),
g_{1}(t, j) = −(4X_{1,t} − (4j + 2)Y_{t}).

It should be noted that these functions specialize to f_{0}(j), f_{1}(j), g_{0}(j), g_{1}(j), respec-

SOME NEW FORMULAS FOR π 15

tively, if X_{1,t} = X_{2,t} = N(k) and Y_{t} = 1. Now we define the matrix M^{X} by

M^{X} = (\begin{matrix}x_{1} & x_{2} & x_{3} & … & x_{4k} \\ F_{1}^{X} & 0 & 0 & … \\ G_{1}^{X} & F_{2}^{X} & 0 & … \\ 0 & G_{2}^{X} & F_{3}^{X} & … \\ 0 & 0 & G_{3}^{X} & ... & ⋮ \\ & ... & ... & ... & 0 \\ & & & G_{4k−1}^{X} & F_{4k}^{X} \\ & & & 0 & G_{4k}^{X} \end{matrix}) ,  (4.9)

where

F_{ℓ}^{X} = ( \begin{matrix}f_{1}(1,4(ℓ−1)k+1) & f_{0}(1,4(ℓ−1)k+2) & 0 & ... \\ 0 & f_{1}(2,4(ℓ−1)k+2) & f_{0}(2,4(ℓ−1)k+3) & 0 & ... \\ & ... & ... \\ & & ... & ... \\ & & f_{1}(4k−2,4ℓk−2) & f_{0}(4k−2,4ℓk−1) & 0 \\ & & 0 & f_{1}(4k−1,4ℓk−1) & f_{0}(4k−1,4ℓk) \end{matrix} )

and

G_{ℓ}^{X} = ( \begin{matrix}g_{1}(1,4(ℓ−1)k+1) & g_{0}(1,4(ℓ−1)k+2) & 0 & ... \\ 0 & g_{1}(2,4(ℓ−1)k+2) & g_{0}(2,4(ℓ−1)k+3) & 0 & ... \\ & ... & ... \\ & & ... & ... \\ & & g_{1}(4k−2,4ℓk−2) & g_{0}(4k−2,4ℓk−1) & 0 \\ & & 0 & g_{1}(4k−1,4ℓk−1) & g_{0}(4k−1,4ℓk) \end{matrix} ) ,

Clearly, we have

M′′′ = M^{X} |_{X_{1,t}=X_{2,t}=N(k), Y_{t}=1}.  (4.10)

The evaluation of det M^{X} is given in the theorem below. From the result it is obvious that det M^{X} |_{X_{1,t}=X_{2,t}=N(k), Y_{t}=1} is nonzero, and, thus, also det M. Theorem 6. We have

det M^{X} = (−1)^{k−1}4^{2k(4k^{2}+7k+2)}k^{2k(4k+1)} ∏_{i=1}^{4k}(i + 1)_{4k−i+1}

× ∏_{a=1}^{4k−1} (2X_{1,a} − (32k^{2} + 2a − 1)Y_{a})

× ∏_{1≤a≤b≤4k−1}(2X_{2,b}Y_{a} − 2X_{1,a}Y_{b} − (2b − 2a + 1)Y_{a}Y_{b}), (4.11)

where (α)_{k} is the standard notation for shifted factorials, (α)_{k} := α(α + 1) ⋯ (α + k − 1), k ≥ 1, and (α)_{0} := 1.

Remark 3. Once having found this theorem, it is not difficult to prove (3.5), by working out how the coefficients of P (z) resulting from the Ansatz (4.1)–(4.3) are related to the coefficients of P (z) resulting from the Ansatz (3.4). Since this is not essential for the proof of Theorem 1, we leave the details to the reader.

Proof of Theorem 6. We follow the “identification of factors” method as described in Section 2.4 in [5].

First we show that (2X_{1,a} − (32k^{2} + 2a − 1)Y_{a}) divides det M^{X}, a = 1, 2, …, 4k− 1. What has to be proved is that det M^{X} vanishes for X_{1,a} = (32k^{2} + 2a − 1)Y_{a}/2. This can be done by showing that for this choice of X_{1,a} there is a nontrivial linear combination of the rows of M^{X}. Indeed, if X_{1,a} = (32k^{2} + 2a − 1)Y_{a}/2 we have

\frac{2(X_{2,4k−1} − (N(k) − 1)Y_{4k−1})}{(−4)^{k(4k+1)+1}(16k^{2} + 1) ∏_{t=1}^{4k−1}(4tk + 1)} · (row 0 of M^{X})

+ ∑_{r=0}^{4k} ∑_{s=0}^{4k−a−1} (\frac{(−1)^{r(k−1)}}{4^{rk}} ∏_{t=0}^{r−1} \frac{4k − 1 + 4tk}{16k^{2} + 1 − 4tk}

· 2^{s} ∏_{t=4k−s}^{4k−1} \frac{2X_{1,t} − (32k^{2} + 2t − 1)Y_{t}}{X_{2,t−1} − (16k^{2} + t − 1)Y_{t−1}} )

· (row (16k^{2} − (4k − 1)r − s − 1) of M^{X}) = 0,

as is easy to verify.

Next we claim that (2X_{2,b}Y_{a}−2X_{1,a}Y_{b}−(2b−2a+1)Y_{a}Y_{b}) divides det M^{X}, 1 ≤ a ≤ b ≤ 4k − 1. Let us first impose the additional restriction that a < b. Using the above reasoning, the claim then follows from the fact that if X_{2,b} = \frac{Y_{b}}{Y_{a}} X_{1,a}+(2b−2a+1)Y_{b}/2 we have

∑_{r=0}^{4k} ∑_{s=4k−b−1}^{4k−a−1} \frac{1}{(−4)^{rk}} (∏_{t=1}^{r} \frac{2X_{1,a} − (32k^{2} − 8kt + 2a + 1)Y_{a}}{2X_{1,a} − (64k^{2} + 8k − 8kt + 2a + 1)Y_{a}} )

· 4^{b−4k+s+1} ( ∏_{t=4k−s}^{b} \frac{2X_{1,t}Y_{a} − 2X_{1,a}Y_{t} − (2t − 2a)Y_{a}Y_{t}}{2X_{2,t−1}Y_{a} − 2X_{1,a}Y_{t−1} − (2t − 2a − 1)Y_{a}Y_{t−1}} )

· (row (16k^{2} − (4k − 1)r − s − 1) of M^{X}) = 0,

as is again easy to verify. On the other hand, if a = b, then the same argument shows that (2X_{2,a} − 2X_{1,a} − Y_{a}) divides det M^{X}. It remains to be checked that also Y_{a} divides det M^{X}. Indeed, if Y_{a} = 0 then we have

∑_{r=0}^{4k} \frac{1}{(−4)^{rk}} · (row (16k^{2} − (4k − 1)r − 4k + a) of M^{X}) = 0,

whence Y_{a} divides det M^{X} for a = 1, 2, …, 4k − 1.

These arguments show that the product on the right-hand side of (4.11) divides det M^{X} as a polynomial in the X_{1,a}’s, X_{2,a}’s, and Y_{a}’s.

Clearly, the degree in the X_{1,a}’s, X_{2,a}’s, and Y_{a}’s of det M^{X} is at most 16k^{2} − 1. But the degree of the right-hand side of (4.11) is exactly 16k^{2} − 1. Therefore we have proved that

det M^{X} = C_{1} ∏_{i=1}^{4k}(i + 1)_{4k−i+1} ∏_{a=1}^{4k−1} (2X_{1,a} − (32k^{2} + 2a − 1)Y_{a})

× ∏_{1≤a≤b≤4k−1}(2X_{2,b}Y_{a} − 2X_{1,a}Y_{b} − (2b − 2a + 1)Y_{a}Y_{b}), (4.12)

where C_{1} is a constant independent of the X_{1,a}’s, X_{2,a}’s, and Y_{a}’s.

In order to determine C_{1}, we compare coefficients of

X_{1,1}^{4k}X_{1,2}^{4k−1} ⋯ X_{1,4k−1}^{2}Y_{1}^{1}Y_{2}^{2} ⋯ Y_{4k−1}^{4k−1}  (4.13)

on both sides of (4.12). We claim that the coefficient of this monomial in det M^{X} is equal to det M^{C}, where M^{C} is defined exactly in the same way as M^{X} (see (4.9)), except that the definitions of the functions f_{0}, f_{1}, g_{0}, g_{1} are replaced by

f_{0}(t, j) = (N(k) + j)(−4)^{k},  (4.14a)
f_{1}(t, j) = 4(−4)^{k},  (4.14b)
g_{0}(t, j) = −j,  (4.14c)
g_{1}(t, j) = −4.  (4.14d)

This is seen as follows. The monomial (4.13) does not contain any X_{2,a}. Therefore, for finding its coefficient in det M^{X}, we may set X_{2,a} = 0 in M^{X} for all a.

In which way may the monomial (4.13) appear in det M^{X} (with all X_{2,a} equal to 0)? A typical term in the expansion of det M^{X} is the product of 16k^{2} entries of M^{X}, each from a different row and column. The monomial (4.13) contains X_{1,1}^{4k}. The variable X_{1,1} is only found in columns 4ℓk + 1, ℓ = 0, 1, …, 4k − 1 (and rows labelled by numbers ≡ 1 mod 4k − 1, according to our labelling scheme). Therefore in a product of entries (each from a different row and column) which produces a term containing X_{1,1}^{4k} all the entries from columns 4ℓk + 1 must be ones containing X_{1,1}. This explains the above definitions (4.14b) and (4.14d) of f_{1}(1, 4ℓk + 1) and g_{1}(1, 4ℓk + 1), ℓ = 0, 1, …, 4k − 1, respectively. Moreover, we must generate the Y_{1} in (4.13) from an entry in a column 4ℓk + 2, for some ℓ. (The variable Y_{1} is also found in entries in columns 4ℓk + 1, but these columns are already taken by our choice of entries which contain the X_{1,1}’s.) This explains the definitions (4.14a) and (4.14c) of f_{0}(1, 4ℓk + 2) and g_{0}(1, 4ℓk + 2), ℓ = 0, 1, …, 4k − 1, respectively. Next we ask how we can find (in the remaining columns and rows) entries which contain X_{1,2}^{4k−1}. Arguing in an analogous manner, the variable X_{1,2} only appears in columns 4ℓk + 2, ℓ = 0, 1, 2 …, 4k − 1. One of these columns is already taken by the entry from which we picked Y_{1}. Therefore in all the remaining ones we must choose entries containing X_{1,2}. This explains the definitions (4.14b) and (4.14d) of f_{1}(2, 4ℓk + 2) and g_{1}(2, 4ℓk + 2), ℓ = 0, 1, …, 4k − 1, respectively. Next we consider the term Y_{2}^{2} in (4.13). It must come from two entries in columns 4ℓk + 3, for two different ℓ’s. This explains the definitions (4.14a) and (4.14c) of f_{0}(2, 4ℓk + 3) and g_{0}(2, 4ℓk + 3), ℓ = 0, 1, …, 4k − 1, respectively. Etc.

The evaluation of det M^{C} follows from Lemma 7 below with X_{a} = 1 and Z_{a} = N(k) for a = 1, 2, …, 4k − 1. □

We consider now a more general determinant than det M^{C}, the latter having been defined through the functions in (4.14). Replace these functions by

f_{0}(t, j) = (Z_{t} + j)(−4)^{k},
f_{1}(t, j) = 4(−4)^{k}X_{t},
g_{0}(t, j) = −j,
g_{1}(t, j) = −4X_{t}.

Let us denote the matrix defined by these functions in the same way as before M^{C} by M^{Z}. Clearly, M^{Z} specializes to M^{C} if all X_{t} are set equal to 1 and all Z_{t} to N(k).

The determinant of M^{Z} evaluates as follows.

Lemma 7. We have

det M^{Z} = (−1)^{k−1}2^{16k^{3}+20k^{2}+14k−1}k^{4k}(4k + 1)! ∏_{a=1}^{4k−1} (X_{a}^{4k+1−a} ∏_{b=0}^{a−1}(Z_{a} − 4bk)). (4.15)

Proof. We proceed in a similar way as in the proof of Theorem 6. In the first step we show that the product on the right-hand side of (4.15) divides det M^{Z} as a polynomial in the X_{a}’s and Z_{a}’s. Then, in the second step, we compare the degrees of the product and det M^{Z}. Since the degree of det M^{Z} turns out to be at most the degree of the product, it then follows that det M^{Z} is equal to the product times some constant which is independent of the X_{t}’s and Z_{t}’s. Finally, in the third step, this constant is found by computing the leading coefficient of det M^{Z}.

Step 1. The product ∏_{a=1}^{4k−1} (X_{a}^{4k+1−a} ∏_{b=0}^{a−1}(Z_{a} − 4bk)) divides det M^{Z}. We start by applying several row and column operations to det M^{Z}, with the final goal of reducing the size of the determinant. First, for i = 16k^{2} − 1, 16k^{2} − 2, …, 4k, in this order, we add (−4)^{k} times row i to row i − 4k + 1. (It should be recalled that, according to our labelling scheme, we number the rows of M^{Z} from 0 to 16k^{2} − 1.) Thus, we obtain the determinant of the following matrix:

(\begin{matrix}x_{1} & x_{2} & x_{3} & . . . . . . . . . . & x_{4k} \\ F′_{1} & F′_{2} & F′_{3} & . . . . . . . . . . & F′_{4k} \\ G_{1} & F′_{1} & F′_{2} & . . . . . . . . . . & F′_{4k−1} \\ 0 & G_{2} & F′_{1} & . . . . . . . . . . & F′_{4k−2} \\ 0 & 0 & G_{3} & ... & ⋮ \\ ... & ... & ... & F′_{2} \\ 0 & G_{4k−1} & F′_{1} \\ … & 0 & G_{4k} \end{matrix}) ,

where the x_{ℓ}’s, ℓ = 1, 2, …, 4k, are defined as earlier, and where the (F′_{ℓ})’s and G_{ℓ}’s, ℓ = 1, 2, …, 4k, are the (4k − 1) × (4k) matrices

F′_{ℓ} = ( \begin{matrix}0 & (−4)^{ℓk}Z_{1} & 0 & ... \\ 0 & 0 & (−4)^{ℓk}Z_{2} & 0 & ... \\ & ... & ... \\ & & ... & ... \\ & & & 0 & (−4)^{ℓk}Z_{4k−2} & 0 \\ & & & 0 & 0 & (−4)^{ℓk}Z_{4k−1} \end{matrix} )

and

G_{ℓ} = ( \begin{matrix}−4X_{1} & −4(ℓ−1)k−2 & 0 & ... \\ 0 & −4X_{2} & −4(ℓ−1)k−3 & 0 & ... \\ & ... & ... \\ & & ... & ... \\ & & 0 & −4X_{4k−2} & −4ℓk−1 & 0 \\ & & & 0 & −4X_{4k−1} & −4ℓk \end{matrix} ) .

Next we “make” the submatrices G_{ℓ}, ℓ = 1, 2, …, 4k, to diagonal matrices, by subtracting (j + 1)/4X_{j} times column j from column j + 1, j = 1, 2, …, 4k − 1, (j + 1)/4X_{j−4k} times column j from column j + 1, j = 4k + 1, 4k + 2, …, 8k − 1, …, and (j + 1)/4X_{j−16k^{2}+4k} times column j from column j + 1, j = 16k^{2} − 4k + 1, 16k^{2} − 4k + 2, …, 16k^{2} − 1. After these operations we obtain the determinant of the matrix

(\begin{matrix}x_{1} & x_{2} & x_{3} & . . . . . . . & x_{4k} \\ F_{1,1} & F_{2,2} & F_{3,3} & . . . . . . . & F_{4k,4k} \\ G & F_{1,2} & F_{2,3} & . . . . . . . & F_{4k−1,4k} \\ 0 & G & F_{1,3} & . . . . . . . & F_{4k−2,4k} \\ 0 & 0 & G & ... & ⋮ \\ ... & ... & ... & F_{2,4k} \\ 0 & G & F_{1,4k} \\ … & 0 & G \end{matrix}) ,

where G and the F_{α,β}’s, 1 ≤ α ≤ β ≤ 4k, are the (4k − 1) × (4k) matrices

G = ( \begin{matrix}−4X_{1} & 0 & 0 & ... \\ 0 & −4X_{2} & 0 & 0 & ... \\ & ... & ... & ... \\ & & ... & ... & ... \\ & & & 0 & −4X_{4k−2} & 0 & 0 \\ & & & & 0 & −4X_{4k−1} & 0 \end{matrix} )

and

F_{α,β} = ( \begin{matrix}0 & f_{1,2}^{(α,β)} & f_{1,3}^{(α,β)} & ... & & f_{1,4k}^{(α,β)} \\ 0 & 0 & f_{2,3}^{(α,β)} & f_{2,4}^{(α,β)} & ... & f_{2,4k}^{(α,β)} \\ & ... & ... & & & ⋮ \\ & & ... & ... & & ⋮ \\ & & & 0 & f_{4k−2,4k−1}^{(α,β)} & f_{4k−2,4k}^{(α,β)} \\ & & & 0 & 0 & f_{4k−1,4k}^{(α,β)} \end{matrix} ) ,

where

f_{rs}^{(α,β)} = \frac{(−4)^{αk−s+r+1}Z_{r} (4(β − 1)k + r + 2)_{s−r−1}}{X_{r+1}X_{r+2} ⋯ X_{s−1}} .

Now we eliminate the last columns in F_{αβ} for 1 ≤ α < β ≤ 4k. We start by eliminating the last column of F_{1,4k}. We do this by adding f_{r,4k}^{(1,4k)}/4X_{r} times column 16k^{2} − 8k + r to column 16k^{2}, r = 1, 2, …, 4k − 1. This makes all the entries in the last column which are in rows 16k^{2} − 8k + 2, …, 16k^{2} − 4k − 1, 16k^{2} − 4k zero, whereas the entries in the last column in rows 16k^{2} − 12k + 3, …, 16k^{2} − 8k + 1 are modified. Next we eliminate these entries in a similar fashion, by using the columns 16k^{2} − 12k + r, r = 1, 2, …, 4k − 1, etc. In the end all the entries in the last column in rows 4k, 4k + 1, …, 16k^{2} − 1 will be zero, whereas the entries in the last column in rows 1, 2, …, 4k − 1 will have been (significantly) modified. An analogous procedure is applied to eliminate the entries in the last columns of F_{α,4k−1}, α = 1, 2, …, 4k − 2. Just to mention the first step: We add f_{r,4k}^{(1,4k−1)}/4X_{r} times column 16k^{2} − 12k + r to column 16k^{2} −4k, r = 1, 2, …, 4k −1. This makes all the entries in column 16k^{2} −4k which are in rows 16k^{2} −12k +3, …, 16k^{2} −8k, 16k^{2} −8k +1 zero, whereas the entries in the last column in rows 16k^{2} − 16k + 4, …, 16k^{2} − 12k + 2 are modified. Etc.

The advantage after having done all this is that now all the entries in columns 4ℓk, ℓ = 1, 2, …, 4k, are zero except for entries in rows 0, 1, …, 4k − 1. This fact, and the fact that the submatrices G are diagonal matrices (of rectangular form) with last column consisting entirely of zeroes, makes it possible to reduce the determinant of the (new) matrix significantly. For i = 16k^{2} − 1, 16k^{2} − 2, …, 4k we may expand the determinant with respect to row i, in this order. If the details are worked out, then we see that our original determinant det M^{Z} is equal to

(−4)^{16k^{2}−4k} ∏_{a=1}^{4k−1} X_{a}^{4k} det (\begin{matrix}u_{1} & & & u_{2} & & & … & & & u_{4k} \\ e_{1,1} & & & e_{1,2} & & & … & & & e_{1,4k} \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ e_{4k−1,1} & & & e_{4k−1,2} & & & … & & & e_{4k−1,4k} \end{matrix}) ,

where u_{ℓ} is given by (4.5), and where e_{ij} is given by

∑_{r=1}^{j}_{∑_{4k=t_{0}>t_{1}>···>t_{r}=i} ; j+1=n_{0}>n_{1}>n_{2}>···>n_{r}=1}(−1)^{r−1}(−4)^{jk−4k+i+1}Z_{t_{1}}Z_{t_{2}} ⋯ Z_{t_{r}}

× \frac{∏_{ν=0}^{r−1}(4k(n_{ν} − 1) − 4k + t_{ν+1} + 2)_{t_{ν}−t_{ν+1}−1}}{X_{i+1}X_{i+2} ⋯ X_{4k−1}} .

Clearly, we may extract (−4)^{jk} from column j, j = 1, 2, …, 4k, and (−4)^{−4k+i+1}/X_{i+1}X_{i+2} ⋯ X_{4k−1} from row i, i = 1, 2, …, 4k − 1 (still using our nonstandard labelling scheme where the rows are numbered 0, 1, …, 4k − 1), so that we obtain the expression

(−4)^{16k^{2}−4k+2k^{2}(4k+1)−(2k−1)(4k−1)} ∏_{a=1}^{4k−1} X_{a}^{4k+1−a}

× det (\begin{matrix}ũ_{1} & & & ũ_{2} & & & … & & & ũ_{4k} \\ ẽ_{1,1} & & & ẽ_{1,2} & & & … & & & ẽ_{1,4k} \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ ẽ_{4k−1,1} & & & ẽ_{4k−1,2} & & & … & & & ẽ_{4k−1,4k} \end{matrix}) , (4.16)

where ũ_{ℓ} is given by

(−1)^{ℓ−1}(−4)^{k}8k(4k + 1) ( ∏_{i=1}^{4k−ℓ}(4ik − 1)) (∏_{i=1}^{ℓ−1}(4ik + 1)) ,

and where ẽ_{ij} is given by

∑_{r=1}^{j}_{∑_{4k=t_{0}>t_{1}>···>t_{r}=i} ; j+1=ñ_{0}>ñ_{1}>ñ_{2}>···>ñ_{r}=1}(−1)^{r−1}Z_{t_{1}}Z_{t_{2}} ⋯ Z_{t_{r}}

× ∏_{ν=0}^{r−1}(4k(ñ_{ν} − 1) − 4k + t_{ν+1} + 2)_{t_{ν}−t_{ν+1}−1}. (4.17)

From (4.16) it is abundantly clear that ∏_{a=1}^{4k−1} X_{a}^{4k+1−a} divides det M^{Z} (which, after all, is equal to (4.16)). It remains to show that also ∏_{a=1}^{4k−1} ∏_{b=0}^{a−1}(Z_{a} − 4bk) divides det M^{Z}.

Let a and b, 1 ≤ a ≤ 4k−1, 0 ≤ b ≤ a−1, be given. We want to show that Z_{a}−4bk divides det M^{Z}. We will show the equivalent fact that the rows of det M^{Z} |_{Z_{a}=4bk} are linearly dependent. The crucial observation, from which this claim follows easily, is that the entries ẽ_{ij}, j = 1, 2, …, 4k, in row i of the determinant in (4.16) are given by a polynomial in j, p_{i}(j) say, of degree (4k − i − 1) and with leading coefficient (i.e., coefficient of j^{4k−i−1}) equal to

∑_{r=1}^{j} ∑_{4k=t_{0}>t_{1}>···>t_{r}=i}(−1)^{r−1}Z_{t_{1}}Z_{t_{2}} ⋯ Z_{t_{r}} \frac{(4k)^{4k−i−r}}{∏_{ν=0}^{r−1}(t_{ν} − i)}.  (4.18)

This is seen as follows. The summand in (4.17) is a polynomial in j, ñ_{1}, ñ_{2}, …, ñ_{r−1} of multidegree (4k − t_{1} − 1, t_{1} − t_{2} − 1, t_{2} − t_{3} − 1, …, t_{r−1} − t_{r} − 1) (i.e., the degree in j is 4k − t_{1} − 1, the degree in ñ_{1} is t_{1} − t_{2} − 1, etc.). Because of the fact that (for fixed u and varying v) ∑_{γ=u}^{v} γ^{e} is a polynomial in v of degree e + 1 with leading coefficient (i.e., coefficient of v^{e+1}) equal to 1/(e+1), successive summation over ñ_{r−1}, ñ_{r−2}, …, ñ_{1} yields the claimed facts.

Under the specialization Z_{a} = 4bk, it is seen by “inspection” that (4.18) with i replaced by a − b vanishes, because the summand corresponding to t_{0} > ⋯ > t_{ω−1} > t_{ω} = a > t_{ω+1} > ⋯ > t_{r} cancels with the summand corresponding to t_{0} > ⋯ > t_{ω−1} > t_{ω+1} > ⋯ > t_{r}. Hence, the polynomial p_{a−b}(j) has degree (at most) 4k − a + b − 2 (instead of 4k − a + b − 1). Consequently, if Z_{a} = 4bk then the entries in rows a − b, a − b + 1, …, 4k − 1 are given by polynomials in j (to wit: j denoting the column index of the entries) of respective degrees 4k − a + b − 2 (!), 4k − a + b − 2, 4k − a + b − 3, …, 2, 1. These are 4k − a + b polynomials, all of degree at most 4k − a + b − 2. It follows that there must be a nontrivial linear combination of these polynomials that vanishes. Hence, the rows a − b, a − b + 1, …, 4k − 1 are linearly dependent, which, in turn, implies that the determinant in (4.16) (and, thus, also det M^{Z}) vanishes for Z_{a} = 4bk, 1 ≤ a ≤ 4k − 1, 0 ≤ b ≤ a − 1.

Step 2. Comparison of degrees. Clearly, the degree of det M^{Z} as a polynomial in the X_{t}’s and Z_{t}’s is at most 16k^{2} − 1, whereas the degree of the product on the right-hand side of (4.15) is exactly 16k^{2} − 1. Hence, we have

det M^{Z} = C_{2} ∏_{a=1}^{4k−1} (X_{a}^{4k+1−a} ∏_{b=0}^{a−1}(Z_{a} − 4bk)),

where C_{2} is a constant independent of the X_{t}’s and Z_{t}’s.

Step 3. Computation of the leading coefficient. In order to determine C_{2}, we determine the coefficient of ∏_{a=1}^{4k−1} X_{a}^{4k+1−a}Z_{a}^{a} in the expansion of det M^{Z}. By arguments similar to those at the end of the proof of Theorem 6, it is seen that this coefficient is given by the determinant of the following matrix, which we denote by M^{L}. It is defined exactly in the same way as M^{X} (see (4.9)), except that the definitions of the functions f_{0}, f_{1}, g_{0}, g_{1} are replaced by

f_{0}(t, j) = (−4)^{k},
f_{1}(t, j) = 4(−4)^{k},
g_{0}(t, j) = 0,
g_{1}(t, j) = −4.

By expanding this determinant with respect to row 0, we obtain

det M^{L} = ∑_{ℓ=1}^{4k}(−1)^{ℓ−1}u_{ℓ} det M_{ℓ}^{L},  (4.19)

where u_{ℓ} is, as earlier, given by (4.5), and M_{ℓ}^{L} is the matrix arising from M^{L} by deleting row 0 and column 4ℓk.

Let ℓ, 1 ≤ ℓ ≤ 4k, be fixed. We will next compute det M_{ℓ}^{L}. When we built M_{ℓ}^{L} from M^{L}, we deleted in particular row 0. Therefore we will now switch to the usual labelling scheme for rows and columns of a matrix, i.e., we will subsequently not only label the columns by 1, 2, … but also the rows.

If ℓ < 4k, then we expand det M_{ℓ}^{L} with respect to the last 4k−1 rows and then with respect to the last column. Since these rows and this column contain only one nonzero entry, we obtain some multiple of the determinant of a (16k^{2}−4k−1)×(16k^{2}−4k−1) matrix. If ℓ < 4k − 1, then we continue by expanding the (now reduced) determinant with respect to the last 4k − 2 rows and then with respect to the last 2 columns. We continue in the same manner until we have reduced det M_{ℓ}^{L} to the determinant of a (4ℓk − 1) × (4ℓk − 1) matrix, more precisely, until we arrive at

det M_{ℓ}^{L} = 4^{\binom{4k}{2}−\binom{ℓ}{2}}(−4)^{k\binom{4k−ℓ+1}{2}} det M′_{ℓ},

SOME NEW FORMULAS FOR π 24

where M′_{ℓ} is the matrix

(\begin{matrix}F & 0 & 0 & . . . . . . & 0 \\ G & F & 0 & . . . . . . & 0 \\ 0 & G & F & . . . . . . & 0 \\ ⋮ & ... & ... & ⋮ \\ 0 & … & 0 & G & F & 0 \\ 0 & . . . . . . . & 0 & G & F′ \\ 0 & . . . . . . . . . . . . & 0 & U \end{matrix}) , (ℓ occurrences of F),

with F and G the (4k − 1) × (4k) matrices

F = (\begin{matrix}4(−4)^{k} & (−4)^{k} & 0 & … \\ 0 & 4(−4)^{k} & (−4)^{k} & 0 & … \\ & ... & ... & ... \\ & & ... & ... & ... \\ & & & 0 & 4(−4)^{k} & (−4)^{k} & 0 \\ & & & & 0 & 4(−4)^{k} & (−4)^{k} \end{matrix})

and

G = (\begin{matrix}−4 & 0 & 0 & … \\ 0 & −4 & 0 & 0 & … \\ & ... & ... & ... \\ & & ... & ... & ... \\ & & & 0 & −4 & 0 & 0 \\ & & & & 0 & −4 & 0 \end{matrix}) ,

F′ the (4k − 1) × (4k − 1) matrix which arises from F by deleting its last column, and U the (ℓ − 1) × (4k − 1) matrix

U = (\begin{matrix}−4 & 0 & 0 & . . . . . . . . . . . & 0 \\ 0 & −4 & 0 & . . . . . . . . . . . & 0 \\ ⋮ & ... & ... & ⋮ \\ 0 & . . . . . . . & −4 & 0 & … & 0 \end{matrix}) .

We continue by expanding M′_{ℓ} with respect to the last ℓ − 1 rows. Thus we obtain

det M_{ℓ}^{L} = 4^{\binom{4k}{2}−\binom{ℓ}{2}+ℓ−1}(−4)^{k\binom{4k−ℓ+1}{2}} det \overline{M}_{ℓ}  (4.20)

for det M_{ℓ}^{L}, where \overline{M}_{ℓ} is the matrix

(\begin{matrix}F & 0 & 0 & . . . . . . & 0 \\ G & F & 0 & . . . . . . & 0 \\ 0 & G & F & . . . . . . & 0 \\ ⋮ & ... & ... & ⋮ \\ 0 & … & 0 & G & F & 0 \\ 0 & . . . . . . . & 0 & G & V \end{matrix} ) , (ℓ − 1 occurrences of F),

with V the (4k − 1) × (4k − ℓ) matrix

V = (\begin{matrix}0 & & 0 & & & . . . . . . . . . . . . . . . . . . . & & 0 \\ . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . & . \\ 0 & & 0 & & & . . . . . . . . . . . . . . . . . . . & & 0 \\ (−4)^{k} & & 0 & & & … \\ 4(−4)^{k} & (−4)^{k} & 0 & … \\ 0 & 4(−4)^{k} & (−4)^{k} & 0 & … \\ ... & ... & ... \\ ... & ... & ... \\ 0 & 4(−4)^{k} & (−4)^{k} \\ & & & & & 0 & & 4(−4)^{k} \end{matrix}) .  (4.21)

Next we prepare for a reduction from the top of \overline{M}_{ℓ}. We subtract 4 times column j from column j − 1, j = 4(ℓ − 1)k, 4(ℓ − 1)k − 1, …, 4k(ℓ − 2) + 2, j = 4k(ℓ − 2), 4k(ℓ − 2) − 1, …, 4k(ℓ − 3) + 2, …, j = 4k, 4k − 1, …, 2, in this order. Thus det M′_{ℓ} is converted to det M′′_{ℓ}, where M′′_{ℓ} is the matrix

(\begin{matrix}F′′ & 0 & 0 & . . . . . . . . . . . . & 0 \\ G′′ & F′′ & 0 & . . . . . . . . . . . . & 0 \\ 0 & G′′ & F′′ & . . . . . . . . . . . . & 0 \\ 0 & 0 & G′′ & ... & ⋮ \\ ⋮ & ... & ... & ... & 0 \\ 0 & . . . . . . . & 0 & G′′ & F′′ & 0 \\ 0 & . . . . . . . . . . . . & 0 & G′′ & V \end{matrix} ) ,

with F′′ and G′′ the (4k − 1) × (4k) matrices

F′′ = (\begin{matrix}0 & (−4)^{k} & 0 & … \\ 0 & 0 & (−4)^{k} & 0 & … \\ & ... & ... & ... \\ & & ... & ... & ... \\ & & & 0 & 0 & (−4)^{k} & 0 \\ & & & & 0 & 0 & (−4)^{k} \end{matrix})

and

G′′ = (\begin{matrix}−4 & 0 & 0 & … \\ (−4)^{2} & −4 & 0 & 0 & … \\ ⋮ & ... & ... & ... \\ ⋮ & ... & ... & ... \\ (−4)^{4k−2} & . . . . . . . & (−4)^{2} & −4 & 0 & 0 \\ (−4)^{4k−1} & . . . . . . . & (−4)^{3} & (−4)^{2} & −4 & 0 \end{matrix}) .

Our next goal is to “push” the (nonzero) entries in columns 4tk+1, t = 0, 1, …, ℓ− 2, down to rows (4k − 1)(ℓ − 1) + 1, …, (4k − 1)ℓ − 1, (4k − 1)ℓ. (This is similar to what we did in Step 1 when we “pushed” all the nonzero entries in columns 4tk, t = 1, 2, …, 4k up to rows 0, 1, …, 4k − 1.) In order to achieve this for the 1-st column, we add

∑_{r=1}^{ℓ−2} ∑_{s=2}^{4k}(−1)^{r}(−4)^{−rk+s−1}\binom{s − 1}{r − 1} · (column (4kr + s))

to column 1. Similarly, in order to achieve this for the (4k + 1)-st column, we add

∑_{r=2}^{ℓ−2} ∑_{s=2}^{4k}(−1)^{r−1}(−4)^{−(r−1)k+s−1}\binom{s − 1}{r − 2} · (column (4kr + s))

to column 4k + 1. Etc. As a result, the determinant det M′′_{ℓ} is converted to the determinant of the matrix

(\begin{matrix}F′′ & 0 & 0 & . . . . . . . . . . . . . . & 0 \\ G′ & F′′ & 0 & . . . . . . . . . . . . . . & 0 \\ 0 & G′ & F′′ & . . . . . . . . . . . . . . & 0 \\ 0 & 0 & G′ & ... & ⋮ \\ ⋮ & ... & ... & ... & 0 \\ 0 & . . . . . . . & 0 & G′ & F′′ & 0 \\ H_{1} & H_{2} & . . . . . . . & H_{ℓ−2} & G′′ & V \end{matrix} ) ,  (4.22)

where F′′ and G′′ are as before, G′ is the (4k − 1) × (4k) matrix

G′ = (\begin{matrix}0 & 0 & 0 & … \\ 0 & −4 & 0 & 0 & … \\ ⋮ & ... & ... & ... \\ ⋮ & & ... & ... & ... \\ 0 & (−4)^{4k−3} & … & (−4)^{2} & −4 & 0 & 0 \\ 0 & (−4)^{4k−2} & … & (−4)^{3} & (−4)^{2} & −4 & 0 \end{matrix}) ,

and H_{t}, t = 1, 2, …, ℓ − 2, is a (4k − 1) × (4k) matrix with all entries equal to 0, except for the entries in the first column. To be precise, the entry in the first column and row s of H_{t} is given by

(−1)^{ℓ−1−t}(−4)^{−(ℓ−1−t)k+s}\binom{s − 1}{ℓ − 1 − t}.  (4.23)

It should be noted that the entries in the first column of G′′ (G′′ appearing at the bottom of the matrix (4.22), as do the matrices H_{t}) is given by (4.23) with t = ℓ − 1.

Now everything is prepared for the reduction. We expand the determinant of (4.22) with respect to rows 1, 2, …, (4k − 1)(ℓ − 1). This reduces the determinant of (4.22) to

(−1)^{\binom{ℓ}{2}}(−4)^{(ℓ−1)(4k−1)k} det M̃_{ℓ},

where M̃_{ℓ} is a (4k − 1) × (4k − 1) matrix of the form

(N V ) ,

with the (s, t)-entry of N being given by (4.23), s = 1, 2, …, 4k−1, t = 1, 2, …, ℓ−1, and V the (4k − 1) × (4k − ℓ) matrix from above. If we substitute all this in (4.20), we obtain that det M_{ℓ}^{L} is equal to

(−1)^{\binom{ℓ}{2}}4^{\binom{4k}{2}−\binom{ℓ}{2}+ℓ−1}(−4)^{k\binom{4k−ℓ+1}{2}+(ℓ−1)(4k−1)k} det M̃_{ℓ}.  (4.24)

The submatrix V of M̃_{ℓ} is almost diagonal. Subtraction of 4 times row s from row s + 1 in M̃_{ℓ}, s = ℓ − 1, ℓ, …, 4k − 2, will transform it into a completely diagonal matrix (namely into the matrix on the right-hand side of (4.21) with all entries 4(−4)^{k} replaced by 0). As a side effect, this will turn the (4k − 1, 1)-entry of M̃_{ℓ} into

(−1)^{ℓ−2}(−4)^{−(ℓ−2)k+4k−1}\binom{4k − 1}{ℓ − 1}.

As is easily seen, the determinant of the in this way modified matrix, M_{ℓ}^{∗} say, is

(−1)^{4k−ℓ+\binom{ℓ−1}{2}}(−4)^{(4k−ℓ)k} · ((4k − 1, 1)-entry of M_{ℓ}^{∗})

× ∏_{s=1}^{ℓ−2} ((s, ℓ − s)-entry of M_{ℓ}^{∗}) ,

or, explicitly,

(−1)^{ℓ}(−4)^{(4k−ℓ)k−(ℓ−2)k+4k−1−\binom{ℓ−2}{2}k+\binom{ℓ−1}{2}}\binom{4k − 1}{ℓ − 1}.

Substitution of the above in (4.24) yields that the determinant det M_{ℓ}^{L} is equal to

(−1)^{ℓk}4^{8k^{3}+10k^{2}−ℓk+2k−1}\binom{4k − 1}{ℓ − 1}.

Now we substitute this in (4.19). We obtain that det M^{L} is equal to

∑_{ℓ=1}^{4k}(−1)^{k}4^{8k^{3}+10k^{2}+3k}\binom{4k − 1}{ℓ − 1}2k(4k + 1) ( ∏_{i=1}^{4k−ℓ}(4ik − 1)) (∏_{i=1}^{ℓ−1}(4ik + 1))

= (−1)^{k}2^{16k^{3}+20k^{2}+14k−1}k^{4k}(4k + 1)(4k − 1)! ∑_{ℓ=1}^{4k} \binom{\frac{1}{4k} − 1}{4k − ℓ}\binom{−\frac{1}{4k} − 1}{ℓ − 1}.

The sum is readily evaluated by means of the Chu–Vandermonde summation (see e.g. [4, Sec. 5.1, (5.27)]), so that we obtain

(−1)^{k}2^{16k^{3}+20k^{2}+14k−1}k^{4k}(4k + 1)(4k − 1)! \binom{−2}{4k − 1}

= (−1)^{k−1}2^{16k^{3}+20k^{2}+14k−1}k^{4k}(4k + 1)! .

Since the coefficient of ∏_{a=1}^{4k−1} X_{a}^{4k+1−a}Z_{a}^{a} in the expression on the right-hand side of (4.15) is exactly the same, we have completed the proof of the lemma. □

References

[1] Fabrice Bellard, π page, http://fabrice.bellard.free.fr/pi/.

[2] D. Bailey, P. Borwein, S. Plouffe, On the rapid computation of various polylogarithmic constants, Math. Comp. 66 (1997), 903–913.

[3] R. Wm. Gosper, unpublished research announcement, 1974.

[4] R. L. Graham, D. E. Knuth and O. Patashnik, Concrete Mathematics, Addison-Wesley, Reading, Massachusetts, 1989.

[5] C. Krattenthaler, Advanced determinant calculus, Séminaire Lotharingien Combin. 42 (“The Andrews Festschrift”) (1999), Article B42q, 67 pp.

[6] S. Plouffe, On the computation of the n’th decimal digit of various transcendental num-

   bers, manuscript, 1996, available at http://www.lacim.uqam.ca/~plouffe/.