π: A 2000-Year Search Changes Direction

Victor Adamchik
 Wolfram Research, Inc. 
100 Trade Center Dr., Champaign, IL 61820
victor@wri.com
&
Stan Wagon
 Macalester College
 St. Paul, MN 55105
wagon@macalstr.edu

One of the charms of mathematics is that it is possible to make elementary discoveries about objects that have been studied for millenia. A most striking example occurred recently when David Bailey of NASA/Ames and Peter Borwein and Simon Plouffe of the Centre for Experimental and Computational Mathematics at Simon Fraser University (henceforth, BBP) discovered a remarkable, and remarkably simple new formula for π. Here is their formula:

π =(∑_(k =0)^∞) 
            (1/(16^k))((4/(8 k +1)) -(2/(8 k +4)) -(1/(8 k +5)) -
              (1/(8 k +6)))

The formula is not too hard to prove, but that misses the point: Finding it in the first place is what took some effort. BBP report that the formulaàs discovery was "a combination of inspired guessing and extensive searching using the PSLQ integer relation algorithm". The PSLQ algorithm is a method for recognizing whether a constant is a combination of other, more fundamental, constants (it has been implemented in Mathematica by Richard Crandall [Cra]). The discoverers sought such a formula because they were aware that it could be used to compute the nth digit of π (in base 16), without computing any prior digits. This goes completely against conventional wisdom, and totally eliminates high-precision requirements from a computation of, say, the billionth hexadecimal digit! The big question now is whether such a method exists for the base-10 digits of π. For a fuller treatment, including several results about other constants and open questions, see [BBP], which can be fetched from Peter Borwein's web site (http://www.cecm.sfu.ca/~pborwein).

In this paper we will discuss the BBP formula and show how Mathematica can be used to generate many other formulas of the same sort. Moreover, the Mathematica-based methods are symbolic, and so a proof of the formula is a natural consequence of the discovery.

1. A Proof of the New Formula for π

To prove the BBP formula we analyze separately the general form of the four summands. First observe that :

(∑_(k =
            0)^∞) (1/((16^k) (8 k +i))) =
        (2^(i/2)) (∑_(k =
                0)^∞)(∫_0^((1/(Sqrt[2]))))(z^(i +8 k -
                    1))\[DifferentialD]z =
          (2^(i/2)) (∫_0^(
                  (1/(Sqrt[2]))))((∑_(k =0)^∞) (z^(i +8 
                          k -1)))\[DifferentialD]z =
          (2^(i/2)) (∫_0^((1/(Sqrt[2]))))
                  ((z^(i -1))/(1 -(z^8)))\[DifferentialD]z

We can now evaluate the full formula, first introducing the new function:

g[i_]  := 
      (2^((i/2))) (∫_0^(
              (1/(Sqrt[2]))))((∑_(k =0)^∞) (z^(8 
                      k +i -1))) \[DifferentialD]z

Then the BBP formula could be proved in one line:

Simplify[4 g[1] - 2 g[4] - g[5] - g[6]]

This sort of series for π did not arise out of the blue. Of course, there is the famous Leibniz series,

π =
      4 (∑_(k =0)^∞)((((-1))^k)/(2 k +1))

but that uses powers of -1, not a useful base. There are, however, other well-known series having the same general form, that is, sums of series, each one of which has as general term a reciprocal of an integer power times a rational number over a linear function. Here are four such that are easy to derive. We call these one-term series since each involves only one linear term; the BBP π-series is a 4-term series.

log(2)  = (∑_(k =0)^∞) (1/(2^k)) 
              (1/k)                   log(2)  =
        (2/3) (∑_(k =0)^∞) (1/(9^k)) (1/(2 k +1))

log(3)  =(∑_(k =0)^∞) (1/(4^k)) 
              (1/(2 k +1))         log( (9/10))  = -
          (∑_(k =0)^∞) (1/(10^k)) \[Null] (1/k)

The first is simply the Maclaurin series for log(1 +x), at x =(1/2). The fourth, which we mention to show that base -10 formulas do exist,  is also a Maclaurin expansion;  this one 
was used by BBP to compute the 5, 000, 000, 065 th decimal digit of  log((9/10)) a world record for decimal digits. The other two come from the Gregory series  - the Maclaurin series for ((1 +x)/(1 -x)) - with x equal to (1/3) for the  (9^k) series and (1/2) for the  (4^k) case. We illustrate the log(3) case; the infinite series represented by the following output must be log(3) because ((1 +x)/(1 -x)) is 3 when x  = (1/2), and it is clearly the same series as the one for log(3) above.

Normal[Series[Log[((1 +x)/(1 -x))], {x, 0, 12}]]/.x->(1/2)

These series could be used to extract, for example, base-2 digits of log(2) or log(3), but these real numbers do not have the cachet that π does, and apparently no one, prior to BBP, thought of using such a series for digit extraction.

2. Far-Out Hex Digits

It is a fairly routine matter to use the BBP series to extract some far-out hex digits. Consider, for example, the real

s  = (∑_(k =0)^∞) (1/(16^k)) (4/(8 k +1))

The far-out digits (say five digits starting from the dth) arise as the fractional part of 
(16^d). s. Multiplying (16^d) into the sum and then splitting it into a sum from 0 to d and a sum from d+1 to infinity does the job. The second sum

s  = (∑_(k =d +1)^∞) Frac[(1/(16^(k -d))) 
            (4/(8 k +1))]

and it is easy to see that 15 terms of this is more than enough for the accuracy required. The first sum is

s  = (∑_(k =0)^d) Frac[(1/(16^(k -d))) (4/(8 k +1))]

and this is computed by reducing the numerator modulo 8k+1. Because we want the fractional part, only the residue counts, and that can be divided by 8k+1 using floating-point arithmetic. Since modular exponentation is speedy via PowerMod, high-precision is not required. Of course, in the case of π, we have four such series to consider.

The DigitExtractor code that follows produces five digits for real numbers having representations of the type we are considering: it takes as input the base b (16 in the BBP case) and the coefficient list ({4, 0,0, -2, -1, -1, 0, 0} in the BBP case) and returns the first five base-b digits, starting from the dth. The leading 3 of π is considered to be the 0th digit. The 15 is enough to guarantee that the omitted terms are negligible, at least for the numbers that arise in this article. A Do-loop is used for mainSum, as discussed in the appendix. The use of N[coeffs] guarantees that the computation is done with machine precision reals; this does limit us, however, because adding the fractional parts of, say, 10,000,000 reals each with 16 digits of precision can lead to unacceptable roundoff. The code works fine for d up to one million, buit not ten million. Serious digit hunters could use, say, 32 digits of precision, but this will really slow things down. One can use negative values of b (example in section 4). The last line of the code is to add back any leading zeros that RealDigits has stripped off. The code is quite short and not optimized for speed. The electronic version of this paper contains a twice faster version.

\(\(DigitExtractor[d_, \ b_, \ coeffs_List]\  
      :=\ Module[{
\ \ \ \ n\  =\ d\  -\ 1, \ mainSum\  =\ 0, \ s\  
            =\ Sign[b], \ base\  =\ Abs[b], 
\ \ \ \ nc\  =\ N[coeffs], 
          
\ \ \ \ \ rd, \ m\  =\ Length[coeffs]}, 
\ \ Do[mainSum\  
            =\ Mod[mainSum\  +\ s^k\ nc\ .\ Table[If[coeffs[[i]]\  ==\ 0, 
                        \ 0, 
\ \ \ \ \ \ 
                        N[PowerMod[base, \ n\  -\ k, \ (k m\  +\ i)]]\ /\ (
                            k m\  +\ i)], \ {i, \ m}], \ 1], \ 
\ \ \ \ {k, 
            \ 0, \ n}]; 
\ \ rd\  =\ RealDigits[Mod[mainSum\  
                +
\ \ \ \ \ \ \ \ \ Sum[s^k\ *\ nc\ .\ (
                        base\ ^\ (n\  -\ k)\ /\ (k m\  +\ Range[m])), 
                  \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ {k, \ n\  +\ 1, \ n\  
                      +\ 15}], \ 1], \ base]; 
\ \ Take[Join[Array[0\ &, 
              \ -rd[[2]]], \ rd[[1]]], \ 5]]
\)\)

We check by examining hexadecimal π.

DigitExtractor[0, 16, {4, 0, 0, -2, -1, -1, 0, 0}]

BaseForm[N[Pi], 16]

Good. Here are some farther-out digits, and a different approach to verification.

DigitExtractor[1000, 16, {4, 0, 0, -2, -1, -1, 0, 0}]

% == RealDigits[N[Pi, 1300], 16][[1, Range[1001, 1005]]]

And finally, we go for the millionth hex digit; this takes about six hours on a PowerMac 8100.

DigitExtractor[10^6, 16, {4, 0, 0, -2, -1, -1, 0, 0}]

These digits agree with the ones published in [BBP]. Take a moment and think about what has been done here: using no high-precision arithmetic, but only routine operations on 16-digit floating-point numbers, and getting absolutely no information about the early digits of π, we have found the millionth digit of π. In fact, Rabinowitz and Wagon [RW] had published a π-digit algorithm that uses only low-precision integer arithmetic, but one had to compute early digits to get late ones and the memory requirements went up as farther-out digits were sought. Even that was considered somewhat surprising. But if that was counterintuitive, the work of Bailey, Borwein, and Plouffe, is totally mind-blowing!

3. A General Method for Finding Formulas (and Proofs)

Off[Solve::svars, General::spell1,  $MaxPrecision::mepr];

Mathematica`s symbolic power allows us to extend the preceding ideas to get many new formulas (together with proofs!), of which the BBP formula is a special case. Here is one result, where r denotes any real (or complex) number.

π =(∑_(k =0)^∞) 
            (1/(16^k))(((4 +8 r)/(8 k +1)) -((8 r)/(8 k +
                  2)) -((4 r)/(8 k +3)) -((2 +8 r)/(8 k +
                  4)) -((1 +2 r)/(8 k +5)) -((1 +2 r)/(8 k +6)) +
              (r/(8 k +7)))

Setting r = 0 yields the BBP formula. The proof is identical to the proof in section 1.

Simplify[PowerExpand[(4 + 8 r) g[1] - 8 r g[2] - 4 r g[3] - (2 + 8 r) g[4] - (1 + 2 r) g[5] - (1 + 2r) g[6] + r g[7]]]

This unenlightening manipulation totally obscures the question of how we found this formula. Indeed, that used an interesting combination of symbolic methods. First we let (Subscript[a, i]), i = 1, . . . , 8, denote undefined constants. The first step is to look at the abstract expression obtained by using these constants on top of the 8 k + i terms in the BBP formula (i running from 1 to 8). It is simplest to use the formulation in terms of integrals (g(i)) derived in section 1.

\(aVals =Array[a\_#&, 8]; 
rawExpr =aVals.Array[g, 8]\)

This looks like a mess. But actually, it is quite beautiful! This will be clear after some simplifications. We can automate this latter simplification as it follows.

LogExpand[expr_] := PowerExpand[expr] /. {Log[n_Integer] :> Apply[#2 . Log[#1]&, Transpose[FactorInteger[n]]],
Log[v_/;v<0] -> Log[-1] + Log[-v],
ArcTan[1/Sqrt[2]] -> Pi/2 - ArcTan[Sqrt[2]]}

We now perform some simplifications on the raw expression.

simpler = Together[LogExpand[rawExpr] /. {
  Log[1 - 1/Sqrt[2]] -> -Log[2]/2 - Log[1 + Sqrt[2]],
  Log[1 + 1/Sqrt[2]] -> -Log[2]/2 + Log[1 + Sqrt[2]]}]

Only seven transcendental numbers appear, and also π(Sqrt[2]). We can find these by defining a
GetTranscentals function as follows.

GetTranscendentals[expr_] := Union[
  Cases[expr, Pi | _ArcTan | _Log, Infinity]]

The vertical bar in the Cases definition is the disjunction operator for patterns; thus the function finds all matching expressions to the specified forms at any level. Now we can collect on the transcendentals, using the third argument of Collect to map Factor, which simplifies things.

transcs = GetTranscendentals[simpler]

collected = Collect[simpler, transcs, Factor]

Note the simple form: seven transcendental numbers are multiplied by linear combinations of the ai. If we can arrange for the multipliers to be, respectively, 0, 1, 0, 0, 0, 0, and 0, then the sum will equal π. To isolate the desired equations, we gather up all the coefficients of the transcendentals, and π(Sqrt[2]) as well. CoefficientList returns a matrix; the reader might wish to examine the steps in more detail, but flattening and deleting the 0s gets us all the a-combinations that occur as coefficients.

(system = DeleteCases[Flatten[
   CoefficientList[collected,
     Append[transcs, Sqrt[2]]]], 0]) // TableForm

The first 6 entries are the logarithmic and arctangent coefficients, the next-to-last entry is the coefficient of π, and the last entry is the coefficient of π(Sqrt[2]). If we want the entire sum to equal π, then we want all but the seventh entry to be 0, with the seventh, the coefficient of  π, being 1.

Solve[system == {0, 0, 0, 0, 0, 0, 1, 0}]

This shows that the solution space has dimension 1; we let (Subscript[a, 4]) be a free parameter, -2 - 8 r.

(Subscript[a, 4]) =-2 -8 r; solution =Simplify[aVals /.First[%%]]

This yields the formula for π at the beginning of this section. If we set r to 0, out pop the BBP coefficients.

solution /. r -> 0

Another natural choice is r = -1/2; the resulting formula for π was also known to Bailey, Borwein, and Plouffe.

solution /. r -> -1/2

In fact, as pointed out to us by P. Borwein, one can easily go from the two specific formulas to the general one by examining (Subscript[formula, 1]) + r ((Subscript[formula, 2]) - (Subscript[formula, 1])). As a curiosity, we note that r =(1/8) yields a series for π with first term (22/7). Thus we can get a series expansion of (22/7) - π. By arranging that 1 be the value of the coefficient of (tan^(-1))(2) we get a series representation of (tan^(-1))(2):

aVals /. First[
  Solve[system == {0, 0, 0, 0, 0, 1, 0, 0}]]

Letting r = - (1/4) leads to several vanishing values:

% /.r ->-1/4

thus we get the following representation of (tan^(-1))(2):

(tan^(-1))(2) =(∑_(k =0)^∞) 
            (1/(16^k))((1/(8 k +1)) +(1/(2 (8 k +3))) -
              (1/(4 (8 k +5))) -(1/(8 (8 k +7))))

A quick numerical check is comforting; speedy convergence means that 10 terms give agreement to machine precision.

Sum[1/16^k {1, 1/2, -1/4, -1/8} . (1 / (8 k + {1, 3, 5, 7})),{k, 0, 15}] == ArcTan[2.]

The formulas of this section are not particularly useful vis-a-vis π. But the (tan^(-1))(2) result seems to be new. More important, perhaps further investigations using the symbolic manipulations presented here will be useful in extending the incipient theory of digit extraction.

Exercise: Use this method to derive the following two formulas.

log (3) =(∑_(k =0)^∞) 
            (1/(16^(k +1)))((16/(4 k +1)) +(4/(4 k +3)))

log (5) =(∑_(k =0)^∞) 
            (1/(16^(k +1)))((16/(4 k +1)) +(16/(4 k +2)) +(4/(4 
                    k +3)))

4. And Still More!

It turns out to be fruitful to consider generalizations. For example, we may consider alternating series, or bases other than 16. For the former we can vary e (the sign); for the latter we vary n ( 2n gives the linear coefficient and (2^n) gives the base; thus n = 4 yields 8 and 16, respectively, as in section 4). Here is a general version of the function g of section 1.

g[i_, n_, 
        e_]  := 
      (2^((i/2))) (∫_0^(
              (1/(Sqrt[2]))))((∑_(k =0)^∞) (e^k) 
                  (z^(2 n k +i -1))) \[DifferentialD]z

Throughout this and the next section we will use the following functions, which perform various simplifications. It is natural to use such specialized simplification routines in a specific project, because the built-in simplifications cannot anticipate all possibilities. Of course, it takes a bit of work, in a given context, to discover exactly what sort of simplifications will be useful; but once that is done, they may as well be gathered together in functional form for ease of use.

LogOfRadicals[expr_] := expr /. 
Module[{tmp},
  tmp = First /@ Union[
     Cases[Level[expr, -2], Log[e_] /; !RationalQ[e]]];
  tmp = Select[Flatten[Table[{tmp[[i]], tmp[[j]]},
     {i,Length[tmp] - 1}, {j, i + 1, Length[tmp]}], 1],
  RationalQ[Expand[#[[1]] #[[2]]]]& ];
  	Apply[Log[#1] ->-Log[#2] + Log[Expand[Times[##]]]&,	tmp, {1}]
]

FullExpand[e_, fun_:Identity] := 
	Collect[LogExpand[LogOfRadicals[e]], 
	{Pi, _Log, _ArcTan}, fun]

RationalQ[_Integer | _Rational] := True
RationalQ[_] := False

LogOfRadicals simplifies combinations of logarithms of radicals

LogOfRadicals[Log[3 - 1/Sqrt[3]] + Log[3 + 1/Sqrt[3]]]

while FullExpand applies several of these simplifications at once.

FullExpand[Log[3 - 1/Sqrt[3]] + Log[3 + 1/Sqrt[3]]]

n = 2, Alternating: log(5), π, (tan^(-1))(2)

The n = 1 case does not yield anything interesting, so let us start with n = 2, with alternating signs (e = -1). There are four terms to add up in this case; FullExpand performs major simplifications. The reader might wish to examine these and the later computations in unsimplified form.

Clear[a, Subscript]; n =2; \(aVals =Array[a\_#&, 2 n]; 
rawExpr =FullExpand[aVals.Array[g[#1, 
                n, -1]&, 2 n]]\)

Things are somewhat simpler than the base-16 case in that the multipliers yield a nonsingular system. Thus we can get four representations at once by extracting and then inverting the matrix of coefficients of the (Subscript[a, i]).

transcs = GetTranscendentals[rawExpr];
Coefficient[rawExpr, #] & /@ transcs

Inverse[Transpose[Table[Coefficient[Expand[%], (Subscript[a, i])], {i, 4}]]]

The columns give representations for π, (tan^(-1))(2), log(2) and log(5), respectively. Note that the π-formula involves only 3 terms. Because log(2) has the well known series we omit its representation from the following list.

π =(∑_(k =
            0)^∞)((((-1))^k)/
            (4^k)) ((2/(4 k +1)) +(2/(4 k +2)) +(1/(4 k +3)))

(tan^(-1))(2) =(∑_(k =0)^∞)((((-1)
                )^k)/(4^k)) ((1/(4 k +1)) +(1/(2 (4 k +3))))

log(5) =(∑_(k =
            0)^∞)((((-1))^k)/
            (4^k)) ((2/(4 k +1)) -(1/(4 k +3)))

The π represention can be used to extract base-4 digits, but the base-16 extractor does this anyway, so it is not a computationally important simplification. However, it should be noted that the base-4 formula for π can be split into the case of k even and k odd, in which case a base-16 formula (a six-termer that is a case of the general formula in section 3) falls right out!

In the nonalternating case with n = 2, nothing new arises. We do get a series for log(3), but it is identical to the Gregory series mentioned in section 1.

n = 3, Nonalternating: log(7)

n =3; \(aVals\  =Array[a\_#&, 2 n]; 
rawExpr =FullExpand[aVals.\((g[#1, n, 
                  1]&)\)/@Range[2 n]]\)

transcs = GetTranscendentals[rawExpr];

(system = DeleteCases[Flatten[
  CoefficientList[rawExpr,
  Join[transcs, {Sqrt[3], Sqrt[2]}]]], 0]) // TableForm

aVals /. Solve[system == {0,0,0,0,1,0,0,0,0,0,0,0,0}]

So now we have a proved two-term formula for log(7) (this formula can also be deduced by the methods of [BBP]). We rearrange a little to get the form below.

log(7) =(∑_(k =0)^∞) 
            (1/(8^(k +1)))((12/(3 k +1)) +(6/(3 k +2)))

A numerical check is comforting.

Sum[1/8^(k +1) (12/(3 k +1)  + 6/(3 k +2)),  {k,  0, 
           20}]  == Log[7.]

To illustrate the general nature of DigitExtractor , we generate some digits of 
log(7). The reader can check that the digits obtained are correct. We must reinterpret the series as involving(8^(-k)) as opposed to (8^(-k -1)).

DigitExtractor[100, 8, {12/8, 6/8, 0}]

Exercise:  Show that the alternating case with n = 3 leads to

π (Sqrt[2]) =(∑_(k =
            0)^∞)((((-1))^k)/
            (8^k)) ((4/(6 k +1)) +(1/(6 k +3)) +(1/(6 k +5)))

5. A Powerful Closed Form

While we could examine more general bases by replacing 2 with b in the integrals defined by g, a slightly different approach to the general case is more enlightening, both in finding formulas and in trying to get a theoretical handle on the big picture. In this section we indicate this alternative method, and provide a general recipe for searching for series representations for logarithms and π. As discussed, the starting point is the following type of sum:

(∑_(k =0)^∞) 
          (1/(b^(n k)))(∑_(i =1)^n) ((Subscript[a, i])/(i +nk))

where n is a positive integer and b ≥ 1 is an algebraic number. The constants (Subscript[a, i]) (from the computational point of view we are interested only in rational values) have to be found so that the  sum yields a useful combination of logarithms and π. The key point is to represent the above sum as a sum of integrals:

(∑_(k =
          0)^∞)(( (Subscript[a, i]) )/
          (b^n))(∫_0^1) 
                ((x^(i -1))/((b^n) -(x^n)))\[DifferentialD]x

This identity is easily proved by expanding the right-side using the geometric series formula. Now, instead of relying on a symbolic integrator to evaluate the right-hand side, we can do it ourselves, which will add both speed and understanding to the investigation. Using the well-known algorithm for integrating rational functions (the details are somewhat lengthy and are omitted; the idea is to use the complex roots of (b^n) -(a^n), which come in conjugate pairs), we obtain the following complicated, but nevertheless nicely closed, expression for the right-hand side of (*).

GoldenGoose[b_, n_] :=(∑_(i =1)^n) 
          (1/n)((Subscript[a, i]) 
            (b^i) (Log[b] (
                  1 +2 (∑_(k =1)^(Floor[((n -1)/2)]))Cos[((2 
                              π i 
                              k)/n)]) -Log[b -1] -(((-1))^i) If[EvenQ[
                    n], Log[b +1] -Log[b], 
                  0] -(∑_(k =1)^(Floor[((n -1)/2)]))Cos[((2 π i 
                        k)/n)] Log[(b^2) -2 b Cos[((2 k π)/n)] +1] +
              2 (∑_(k =1)^(Floor[((n -1)/2)]))Sin[((2 i k 
                          π)/n)] ArcTan[Sin[((2 k 
                              π)/n)]/(b -Cos[((2 k 
                                π)/n)])]))

The reader might wish to compare specific instances of the formula to results of symbolic integration to provide evidence of correctness. This formula is pretty rich for further investigation. For example, the arctangent is the only place that could yield π. Let see what we can learn by trying to get the argument to arctan to be 1. We would need

sin(((2 k π)/n))/(b -cos(((2 k π)/n))) =
      1

or

b =cos(((2 k π)/n)) +sin(((2 k π)/n)) =(Sqrt[2]) sin(((2 π 
                  k)/n) +(π/4))

Since we want to have our base, b, be an integer, we need to find an integer n so that

((((Sqrt[2]) sin(((2 π k)/n) +(π/4))))^n)

is a positive integer. Simple experiments show that the first two possibilities are: n = 4 and n = 8. The second case (with the base (b^n) = 16) leads to the BBP formula. For n = 4 the base is one. The unit base is a singular case of (*): the integrals are divergent. We can avoid the resulting singularity by making (Subscript[a, 1]) +(Subscript[a, 2]) +(Subscript[a, 3]) +(Subscript[a, 4]) = 0 (one can then view the right side of (*) as

∫((p (x))/(1 -
            (x^n)))\[DifferentialD]x

where p(x) is a polynomial and p(1) = 0; thus the troublesome 1 - x in the denominator is cancelled, since 1 - x divides p(x)). This yields the following expression for the main sum.

π (((Subscript[a, 1])/8) -((Subscript[a, 3])/8)) +
      2 (((3 (Subscript[a, 1]))/4) +((Subscript[a, 2])/2) +((3 (Subscript[a, 3]))/4)) log(2)

Setting it equal to π produces a one-parameter formula for π:

π =(∑_(k =
            0)^∞)(((r +4)/(4 k +1)) -((3 r)/(4 k +
                2)) +((r -4)/(4 k +3)) +(r/(4 k +4)))

If we specialize to r = 0, we simply get a rewritten version of the Leibniz series of π. So really this should be viewed as a new formula for 0:

0 =(∑_(k =
            0)^∞)((1/(4 k +1)) -(3/(4 k +2)) +(1/(4 k +3)) +
            (1/(4 k +4)))

since that is how the Leibniz series becomes a one-parameter family of formulas. Despite 
the fact that nothing new about π was learned here, it seems to us quite possible that this methodology might lead to nonexistence proofs. The outstanding such development would be a proof that π cannot be represented with this sort of series using powers of 10.

Of course, we can also use this closed form to explore some logarithms. Earlier we obtained several representations of natural logarithms of integers (2, 3, 5, and 7) in base 2. Here we consider other bases. Setting b to 3 and n to 6 yields a representation of log(13). The procedure is similar enough to earlier work that we show only partial output.

Clear[Subscript];
{b, n} = {3, 6};
rawExpr = FullExpand[GoldenGoose[b, n]];
transcs = GetTranscendentals[rawExpr];
Collect[rawExpr, transcs][[-2]]

system = DeleteCases[Flatten[CoefficientList[rawExpr, transcs]], 0];
First[%]

\(Array[a\_#&, n] /.Solve[system =={1, 0, 0, 0, 0, 0}]\)

This gives a series for log(13), which we rearrange slightly to clear fractions.

log(13) =(∑_(k =0)^∞) 
            (1/(3^(6 k +
                5)))((567/(6 k +1)) +(81/(6 k +2)) +(36/(6 k +3)) +
              (9/(6 k +4)) +(7/(6 k +5)))

Finally, recall that a most intriguing question is whether there are any interesting formulas in base 10. So we let b = 10.

rawExpr = Collect[
  GoldenGoose[10, 2], _Log, Factor]

Setting (Subscript[a, 1]) = 10 (Subscript[a, 2]) and (Subscript[a, 2]) =-(1/100) yields a series that, after simplification, reduces to the Maclaurin expansion of log((9/10)) mentioned in section 1. Setting (Subscript[a, 1 ])=-(1/5) and (Subscript[a, 2]) = 0 yields a series that is a Gregory expansion of log((9/11)) (see section 1). So, nothing really new comes from setting b = 10.

6. Conclusion

We believe that our investigations show the power of Mathematica in serious symbolic investigations. Not only have we found some new formulas, but their discovery comes automatically with a proof. It is, of course, very satisfying to discover new formulas, and perhaps readers will be inspired to carry these methods farther, perhaps by looking at some new forms. It seems promising to investigate alternating series. Our preliminary investigations yielded very quickly a new formula for π (section 4).
             
Clearly these investigations can be carried farther. Of course, there are many more open problems than solved ones! Some of them are listed in [BBP]. We'll mention here only that it would be nice if some negative results could be obtained, such as the nonexistence of a two-term series for π, or a base-10 series for π. Our goal was to develop a technique that might yield an approach to such questions. Also, we think we have clearly shown exactly where the BBP and related formulas come from, and how π arises.

Moreover, these ideas might well yield new series representations for other special numbers or functions. Consider the polylogarithmic function

(Subscript[Li, n])(z) =(∑_(k =1)^∞) ((z^k)/(k^n))

This form is close to the ones we have looked at and could be considered as a starting point for the next generalization. Itàs natural to consider the following:

(∑_(k =1)^∞) (1/(b^(n k)))
              (1/(k^m))(∑_(i =1)^(n +m))((Subscript[a, i])/(i +nk))

This form of series can be obtained from formula (*) in section 6, integrating the latter m times with respect to x. Since the repeated integration produces m free parameters we add them to the inner sum. Then our algebraic techniques yielded:

(Subscript[Li, 2])((1/2)) =(∑_(k =1)^∞) 
            (1/((2^k) (k^2)))(-(8/(k +1)) -(10/(4 k +1)) -(54/(2 
                    k +3)) +(12/(2 k +1)) +(125/(4 k +5)))

(Subscript[Li, 3])((1/2)) =(∑_(k =1)^∞) 
            (1/((2^k) (k^3)))(-(24/(2 k +1)) -(625/(2 (4 
                        k +5))) -(2401/(6 (4 k +7))) +
              (5/(4 k +1)) +(324/(2 k +3)) +(63/(4 k +
                  3)) +(32/(3 (k +1))))

From the computational perspective these are not that valuable, except to provide slightly faster convergence than the definition. But from the theoretical standpoint this gives another method of generating various classes of new series representations. Here are examples:

ζ( 3) =(∑_(k =1)^∞) 
            (1/(k^3))(-(1/(k +1)) +(16/(k +2)) +
              (3/(4 (2 k +1))) -(81/(4 (2 k +3))))

(π^2) =(∑_(k =1)^∞) 
            (1/(k^3))(-(12/(k +1)) +(384/(k +2)) +(45/(2 (2 
                        k +1))) -(1215/(2 (2 k +3))))

π =(∑_(k =1)^∞) 
            (1/(k^3))(-(238/(k +1)) +(285/(2 (2 
                        k +1))) -(667/(32 (4 
                        k +1))) -(5103/(16 (4 
                        k +3))) +(35625/(32 (4 k +5))))

We are grateful to David Bailey, Peter Borwein, and Simon Plouffe for sharing with us many of their insights and unpublished work in this area,and to Mark Sofroniou and
Michael Trott (WRI) for some helpful observations.

Appendix: Sum Subtleties

In order to compute some far-out base-16 digits using the BBP type of formula, one must pay attention to some subtleties of Sum. There are three basic ways to form a sum: (1) use Sum, (2) apply Plus, or (3) use a Do-loop. Method (2), an example of which is Apply[Plus, Table[i, {i, 1000}]], is useful for small sets, but it is clearly costly in terms of memory since the entire table must be computed and stored. So we will restrict our discussion here to a comparison of mathods (1) and (3). Here they are, with calls to MemoryInUse to compare the memory required.

CheckMemory := Print[-mOld + (mOld = MemoryInUse[])]

mOld = MemoryInUse[];
Sum[i, {i, 1000}]; CheckMemory;
s = 0; Do[s += i, {i, 1000}]; CheckMemory;

This output shows a similar memory consumption, but in fact this is not the whole truth. A more careful comparison checks the memory as the sum is being computed, but prints out results only for the first few values of the index (the omitted values are all 0).

s1 = Sum[If[i < 4, CheckMemory]; i, {i, 1000}];
CheckMemory; Print[""];

s3 = 0; Do[If[i < 4, CheckMemory]; s3 += i, {i, 1000}];
CheckMemory

Big surprise! Method (1) has a fairly large internal consumption; the negative value indicates that memory is being freed up after the sum is computed. And the Do-loop uses up hardly any memory at all. If the number of terms is increased, the first memory profile increases accordingly, but the second stays where it is, at virtually no consumption. This means that Sum cannot be used to add up, say, a million numbers. A Do-loop is essential in such a case.

In fact, Sum was designed mostly for symbolic calulations. And Sum is very powerful working with symbols rather then performing numerical calculations, as shown by the following example.

sum1  = Sum[f[k],  {k,  500,  1,  -1}];  //Timing

\(sum2\  =\ 0; 
\(Do[sum2\  +=\ f[k], \ {k, \ 500, \ 1, \ -1}]; 
      \)\  //Timing\)

The Do-loop is much slower in this case because on each iteration, function Plus sorts its arguments, while with Sum arguments are sorted only one time (since Plus is applied only at the end).

Clear[s1, s3, sum1,  sum2,  mOld];

Addendum for Electronic Version of Paper

Here is a faster version of DigitExtractor, contributed by WRI's Mark Sofroniou. It is about twice as fast as the version given in the paper, but a little (though only a little) more complicated.

DigitExtractor::usage = "DigitExtractor[d, b, coeffs] computes five base-Abs[b] digits starting from the dth of the real number given by 
Sum[1/b^k coeffs . (1/(m k + Range[m])), {k, 1, Infinity}], where m = Length[coeffs]. The variable term controls the truncation error; 15 should be adequate. All computations are done in machine precision so d should not be greater than about ten million, or else roundoff error will be excessive for the goal of 5 digits.";

\(\(DigitExtractor[d_, b_, coeffs_List] :=Module[
        {m, nb, ncoeffs, offsets, posns, pow, reald, sum, terms}, ter
            m s =15; base =Abs[b]; sign =Sign[b]; pow =d -1; sum =0; 
        n b =N[base]; m =Length[coeffs]; ncoef
            f s =N[DeleteCases[coeffs, 0]]; pos
            n s =Flatten[Position[coeffs, _?\((# !=0&)\)]]; offsets =posns; 
        
		Do[sum =Mod[sum +sign^\((pow -k)\)*
			\((ncoeffs/offsets)
                      \).\((PowerMod[base, k, offsets])\), 1]; offse
              t s +=m, {k, pow, 0, -1}]; 
        
		offsets =posns +\((pow +1)\)\ m; 
        
		Do[sum =Mod[sum +sign^\((pow +k)\)\ \((ncoeffs)\)
                      .\((nb^\((-k)\)/offsets)\), 1]; offse
              t s +=m, {k, 1, terms}]; reald =RealDigits[sum, base]; Take
          [Join[Table[0, {-reald[\([2]\)]}], reald[\([1]\)]], 5]
	]; \)\)

References

[BBP] D. Bailey, P. Borwein, and S. Plouffe, On the rapid computation of various polylogarithmic constants, Mathematics of Computation (forthcoming).

[Cra] R. E. Crandall, Topics in Advanced Scientific Computation, Springer/TELOS, New York, 1995.

[RW] S. Rabinowitz and S. Wagon, A spigot algorithm for π, American Mathemat ical Monthly 102 (1995) 195-203.