StudyWebAward

Pi: A 2000-Year Search Changes Direction

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


1. A New Formula for Pi

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 Pi. Here is their formula:

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 Pi (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 Pi. 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.

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.

To prove the BBP formula we analyze separately the general form of the four summands. First observe that f(i) = g(i) = h(i), where f, g, and h are as defined below.

f[i_] := Sum[1/(16^k (8 k + i)), {k, 0, Infinity}]

g[i_] := 2^(i/2) Integrate[Sum[z^(8 k + i - 1),{k, 0, Infinity}], {z, 0, 1/Sqrt[2]}]

h[i_] := 2^(i/2) Sum[Integrate[z^(8 k + i - 1), {z, 0, 1/Sqrt[2]}],{k, 0, Infinity}]

The equality of h and g is just a switch of the integral and sum; and f = h by a simple integration, as follows.

2^(i/2) * Integrate[z^(8 k + i - 1), z] /. z -> 1/Sqrt[2] // Simplify

        1
--------------
 4 k
2    (i + 8 k)
We can now evaluate the full formula, first loading SymbolicSum in order to get a closed form for the infinite sum in g. PowerExpand expands out some logarithms, which then disappear.

Needs["Algebra`SymbolicSum`"]
Simplify[PowerExpand[4 g[1] - 2 g[4] - g[5] - g[6]]]
Pi

This sort of series for Pi did not arise out of the blue. Of course, there is the famous Leibniz series, Pi = Sum[(-1)^k 4/(2k+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 Pi-series is a 4-term series.

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,065th 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 -> HoldForm[1/2]

         1 3      1 5      1 7      1 9      1 11
      2 (-)    2 (-)    2 (-)    2 (-)    2 (-)
  1      2        2        2        2        2
2 - + ------ + ------ + ------ + ------ + -------
  2     3        5        7        9        11
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 Pi does, and apparently no one, prior to BBP, thought of using such a series for digit extraction.

2. Sum Subtleties

Before discussing how the BBP formula can produce some far-out base-16 digits, we need to 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;
628
648

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
4716
0
0
-3584

744
0
0
372

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
{0.05 Second, Null}

sum2 = 0;
Do[sum2 += f[k], {k, 500, 1, -1}]; //Timing
{4.81667 Second, Null}

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).

3. 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 =Sum[1/16^k 4/(8k+1),{k,0,Infinity}]. 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 is Sum[ Frac[1/16^(k-d) 4/(8k+1)],{k,d+1,Infinity}], and it is easy to see that 15 terms of this is more than enough for the accuracy required. The first sum is Sum[ Frac[1/16^(k-d) 4/(8k+1)],{k,0,d}]], 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 Pi, 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 Pi 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 section 2. 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 5). The last line of the code is to add back any leading zeros that RealDigits has stripped off. The code isquite short and not optimized for speed. The electronic version of this paper contains a much faster version (due to Mark Sofroniou of WRI).

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[s^k nc . Table[If[coeffs[[i]] == 0, 0, 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 Pi

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

{3, 2, 4, 3, 15}

BaseForm[N[Pi], 16]

3.243f
      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}]

{3, 4, 9, 15, 1}

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

True

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}]

{2, 6, 12, 6, 5}

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 Pi, we have found the millionth digit of Pi. In fact, Rabinowitz and Wagon [RW] had published a 9-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!

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

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.

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]]]

Pi

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 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]

                             3        5                                     1        3
                  a[4] (-Log[-] + Log[-])                -2 ArcTan[2] - Log[-] + Log[-]
            15               4        4             Pi                      2        2
-2 a[8] Log[--] + ----------------------- + 2 a[2] (-- + ------------------------------) + 
            16               2                      8                  8
 
                                  1        3
                2 ArcTan[2] - Log[-] + Log[-]
          -Pi                     2        2
  8 a[6] (--- + -----------------------------) + 
           8                  8
 
                                           1                   1                5
  (a[3] (2 Sqrt[2] ArcTan[2] - 4 ArcTan[-------] + Sqrt[2] Log[-] - Sqrt[2] Log[-] - 
                                        Sqrt[2]                2                2
 
                    1                    1
       2 Log[1 - -------] + 2 Log[1 + -------])) / (4 Sqrt[2]) + 
                 Sqrt[2]              Sqrt[2]
 
                                            1                   1                5
  (a[5] (-2 Sqrt[2] ArcTan[2] + 4 ArcTan[-------] + Sqrt[2] Log[-] - Sqrt[2] Log[-] - 
                                         Sqrt[2]                2                2
 
                    1                    1
       2 Log[1 - -------] + 2 Log[1 + -------])) / (2 Sqrt[2]) + 
                 Sqrt[2]              Sqrt[2]
 
                                            1                   1                5
  (a[7] (-2 Sqrt[2] ArcTan[2] - 4 ArcTan[-------] - Sqrt[2] Log[-] + Sqrt[2] Log[-] - 
                                         Sqrt[2]                2                2
 
                    1                    1
       2 Log[1 - -------] + 2 Log[1 + -------])) / Sqrt[2] + 
                 Sqrt[2]              Sqrt[2]
 
                                           1                   1                5
  (a[1] (2 Sqrt[2] ArcTan[2] + 4 ArcTan[-------] - Sqrt[2] Log[-] + Sqrt[2] Log[-] - 
                                        Sqrt[2]                2                2
 
                    1                    1
       2 Log[1 - -------] + 2 Log[1 + -------])) / (8 Sqrt[2])
                 Sqrt[2]              Sqrt[2]
This looks like a mess. But actually, it is quite beautiful! This will be clear after some simplifications. PowerExpand breaks down terms of the form Log[5/3], but it is also useful to break down, say, Log[4] to 2 Log[2]. We can automate this lattter simplification, including PowerExpand as well, as follows. We also define a function to pull out certain types of transcendental numbers from an expression.

LogExpand[expr_] := PowerExpand[expr] /. Log[n_Integer] :> Apply[#2 . Log[#1]&, Transpose[FactorInteger[n]]]

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

We now perform some simplifications on the raw expression.

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

(Sqrt[2] Pi a[1] + 2 Pi a[2] - 2 Sqrt[2] Pi a[3] + 4 Sqrt[2] Pi a[5] - 
 
    8 Pi a[6] - 8 Sqrt[2] Pi a[7] + 2 a[1] ArcTan[2] - 4 a[2] ArcTan[2] + 
 
    4 a[3] ArcTan[2] - 8 a[5] ArcTan[2] + 16 a[6] ArcTan[2] - 
 
    16 a[7] ArcTan[2] - 2 Sqrt[2] a[1] ArcTan[Sqrt[2]] + 
 
    4 Sqrt[2] a[3] ArcTan[Sqrt[2]] - 8 Sqrt[2] a[5] ArcTan[Sqrt[2]] + 
 
    16 Sqrt[2] a[7] ArcTan[Sqrt[2]] + 64 a[8] Log[2] + 2 a[2] Log[3] - 
 
    4 a[4] Log[3] + 8 a[6] Log[3] - 16 a[8] Log[3] + a[1] Log[5] - 
 
    2 a[3] Log[5] + 4 a[4] Log[5] - 4 a[5] Log[5] + 8 a[7] Log[5] - 
 
    16 a[8] Log[5] + 2 Sqrt[2] a[1] Log[1 + Sqrt[2]] + 
 
    4 Sqrt[2] a[3] Log[1 + Sqrt[2]] + 8 Sqrt[2] a[5] Log[1 + Sqrt[2]] + 
 
    16 Sqrt[2] a[7] Log[1 + Sqrt[2]]) / 8
Only seven transcendental numbers appear, and also Sqrt[2]. We can find these semi-automatically by pattern-matching, and use them to collect, which simplifies things a lot (mapping Factor cleans it up a little). The vertical bar in Cases is the disjunction operator for patterns; thus this selects expressions at any level satisfying any of the three patterns.

transcs = GetTranscendentals[simpler]

{Pi, ArcTan[2], ArcTan[Sqrt[2]], Log[2], Log[3], Log[5], Log[1 + Sqrt[2]]}
collected = Map[Factor, Collect[simpler, transcs]]
Pi (Sqrt[2] a[1] + 2 a[2] - 2 Sqrt[2] a[3] + 4 Sqrt[2] a[5] - 8 a[6] - 8 Sqrt[2] a[7])
-------------------------------------------------------------------------------------- + 
                                          8
 
  (a[1] - 2 a[2] + 2 a[3] - 4 a[5] + 8 a[6] - 8 a[7]) ArcTan[2]
  ------------------------------------------------------------- + 
                                4
 
  (-a[1] + 2 a[3] - 4 a[5] + 8 a[7]) ArcTan[Sqrt[2]]
  -------------------------------------------------- + 8 a[8] Log[2] + 
                      2 Sqrt[2]
 
  (a[2] - 2 a[4] + 4 a[6] - 8 a[8]) Log[3]
  ---------------------------------------- + 
                     4
 
  (a[1] - 2 a[3] + 4 a[4] - 4 a[5] + 8 a[7] - 16 a[8]) Log[5]
  ----------------------------------------------------------- + 
                               8
 
  (a[1] + 2 a[3] + 4 a[5] + 8 a[7]) Log[1 + Sqrt[2]]
  --------------------------------------------------
                      2 Sqrt[2]
Note the simple form: seven transcendental numbers are multiplied by linear combinations of the a[i]. If we can arrange for the multipliers to be, respectively, 1, 0, 0, 0, 0, 0, and 0, then the sum will equal Pi. 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

a[1] + 2 a[3] + 4 a[5] + 8 a[7]
-------------------------------
               4

a[1] - 2 a[3] + 4 a[4] - 4 a[5] + 8 a[7] - 16 a[8]
--------------------------------------------------
                        8

a[2] - 2 a[4] + 4 a[6] - 8 a[8]
-------------------------------
               4

8 a[8]

-a[1] + 2 a[3] - 4 a[5] + 8 a[7]
--------------------------------
               4

a[1] - 2 a[2] + 2 a[3] - 4 a[5] + 8 a[6] - 8 a[7]
-------------------------------------------------
                        4

a[2]
---- - a[6]
 4

a[1]   a[3]   a[5]
---- - ---- + ---- - a[7]
 8      4      2
The first 6 entries are the logarithmic and arctangent coefficients, the next-to-last entry is the coefficient of Pi, and the last entry is the coefficient of Pi Sqrt[2]. If we want the entire sum to equal Pi, then we want all but the seventh entry to be 0, with the seventh, the coefficient of Pi, being 1.

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

{{a[4] -> -2 - 8 a[7], a[2] -> -8 a[7], a[6] -> -1 - 2 a[7], a[8] -> 0, 
 
   a[1] -> 4 + 8 a[7], a[3] -> -4 a[7], a[5] -> -1 - 2 a[7]}}
This shows that the solution space has dimension 1; we let a[7] be a free parameter, r.

a[7] = r;
solution = aVals /. First[%%]

{4 + 8 r, -8 r, -4 r, -2 - 8 r, -1 - 2 r, -1 - 2 r, r, 0}
This yields the formula for Pi at the beginning of this section. If we set r to 0, out pop the BBP coefficients.

solution /. r -> 0

{4, 0, 0, -2, -1, -1, 0, 0}
Another natural choice is r = Ð1/2; the resulting formula for Pi was also known to Bailey, Borwein, and Plouffe.

solution /. r -> -1/2

                     1
{0, 4, 2, 2, 0, 0, -(-), 0}
                     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 formula1 + r (formula2 - formula1). As a curiosity, we note that r = 1/8 yields a series for Pi with first term 22/7. Thus we can get a series expansion of 22/7- Pi.

By arranging that 1 be the value of the coefficient of arctan 2 we get a series representation of arctan 2.

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

                                      1           1
{2 + 8 r, -1 - 8 r, -4 r, -1 - 8 r, -(-) - 2 r, -(-) - 2 r, r, 0}
                                      2           4
Letting r = -1/8 leads to several vanishing values; thus we get the following representation of arctan 2.

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.]
True

The formulas of this section are not particularly useful vis-a-vis Pi. But the arctan 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.

5. 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) Integrate[Sum[e^k z^(2 n k + i - 1),{k, 0, Infinity}], {z, 0, 1/Sqrt[2]}]

The integral that follows tells us that g(i, n, e) = Sum[e^k/(2^(kn) (2 n k + i)),{k,0,Infinity}]; so we may use g to further our explorations.

2^(i/2) * Integrate[e^k z^(2 n k + i - 1), z] /. z -> 1/Sqrt[2] // Simplify

        k
       e
----------------
 k n
2    (i + 2 k n)
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.

CollectPatterns[expr_, p_List, fun_:Identity] :=
Function[l, Fold[Map[fun, Collect[##]]&, expr,
Union[Cases[l, #]]& /@ p]][Level[expr, -1]]

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] := Map[Factor, CollectPatterns[
LogExpand[LogOfRadicals[e]], {Pi, _Log, _ArcTan}, fun]]

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

LogOfRadicals simplifies combinations of logarithms of radicals, while CollectPatterns collects the patterns.

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

    26
Log[--]
    3
CollectPatterns[Cos[r] + Log[1 - s] - e/2 Log[1 - s] + Tan[x] + a/3 Tan[x], {_Log, _Tan}]
              e                    a
Cos[r] + (1 - -) Log[1 - s] + (1 + -) Tan[x]
              2                    3
It is useful to map Factor onto the summands of an expressions, and CollectPatterns allows a third argument to map such a function.

CollectPatterns[%, {_Log, _Tan}, Factor]

         (2 - e) Log[1 - s]   (3 + a) Tan[x]
Cos[r] + ------------------ + --------------
                 2                  3
FullExpand applies several of these simplifications at once.

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

Log[2] - Log[3] + Log[13]
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.

n = 2;
aVals = Array[a, 2 n];
collected = FullExpand[aVals . Array[g[#, n, -1]&, 2n]]

Pi a[2]   (a[1] - 2 a[2] + 2 a[3]) ArcTan[2]
------- + ---------------------------------- - 2 a[4] Log[2] + 
   2                      2
 
  (a[1] - 2 a[3] + 4 a[4]) Log[5]
  -------------------------------
                 4
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 a[i]

Coefficient[collected, #] & /@ {Pi, ArcTan[2], Log[2], Log[5]}

 a[2]  a[1] - 2 a[2] + 2 a[3]           a[1] - 2 a[3] + 4 a[4]
{----, ----------------------, -2 a[4], ----------------------}
  2              2                                4
Inverse[Transpose[Table[Coefficient[Expand[%], a[i]], {i, 4}]]]
                                 1    1                 1
{{2, 1, 1, 2}, {2, 0, 0, 0}, {1, -, -(-), -1}, {0, 0, -(-), 0}}
                                 2    2                 2
The columns give representations for Pi, arctan 2, log 2, and log 5, respectively. Note that the Pi-formula involves only 3 terms. Because log 2 has the well known form Sum[1/(k 2^k), {k,1,Infinity}], we omit its new representation from the following list.

The Pi 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 Pi 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 4) 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;
aVals = Map[a, Range[2 n]];
rawExpr = FullExpand[aVals . Map[g[#, n, 1] &, Range[2 n]]]
      a[2]       2 a[4]        a[1]      a[2]     2 a[4]           2
Pi (--------- - ---------) + (------- + ------- - ------- - 2 Sqrt[-] a[5]) 
    3 Sqrt[3]   3 Sqrt[3]     Sqrt[6]   Sqrt[3]   Sqrt[3]          3
 
          -1 + Sqrt[2]      a[1]      a[2]     2 a[4]           2               1 + Sqrt[2]
   ArcTan[------------] + (------- - ------- + ------- - 2 Sqrt[-] a[5]) ArcTan[-----------] + 
            Sqrt[3]        Sqrt[6]   Sqrt[3]   Sqrt[3]          3                 Sqrt[3]
 
     a[1]      Sqrt[2] a[3]   2 Sqrt[2] a[5]
  (--------- - ------------ + -------------- + 4 a[6]) Log[2] + 
   3 Sqrt[2]        3               3
 
     -a[1]     a[2]   Sqrt[2] a[3]   a[4]   Sqrt[2] a[5]   4 a[6]
  (--------- + ---- - ------------ + ---- - ------------ - ------) Log[7] + 
   6 Sqrt[2]    6          3          3          3           3
 
   Sqrt[2] a[1]   4 Sqrt[2] a[5]             1
  (------------ + --------------) Log[1 + -------] + 
        3               3                 Sqrt[2]
 
     a[1]      2 Sqrt[2] a[5]                     2 Sqrt[2] a[3] Log[4 + Sqrt[2]]
  (--------- + --------------) Log[3 + Sqrt[2]] + -------------------------------
   3 Sqrt[2]         3                                           3
transcs = GetTranscendentals[rawExpr];

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

2 a[3]
------
  3

a[1] + 4 a[5]
-------------
      6

a[1] + 4 a[5]
-------------
      3

a[2]   a[4]   4 a[6]
---- + ---- - ------
 6      3       3

-a[1]   a[3]   a[5]
----- - ---- - ----
 12      3      3

4 a[6]

a[1]   a[3]   2 a[5]
---- - ---- + ------
 6      3       3

-a[2]   2 a[4]
----- + ------
  3       3

a[1]   2 a[5]
---- - ------
 6       3

a[2]   2 a[4]
---- - ------
 3       3

a[1]   2 a[5]
---- - ------
 6       3

a[2] - 2 a[4]
-------------
      9
aVals /. Solve[system == {0,0,0,1,0,0,0,0,0,0,0,0}]
           3
{{0, 3, 0, -, 0, 0}}
           2
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.

A numerical check is comforting.

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

True
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 1/8^kk as opposed to 1/8^(k+1).

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

{3, 0, 0, 5, 7}
Exercise: Show that the alternating case with n = 3 leads to

6. 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 Pi. As discussed, the starting point is the following type of sum:

where n is a positive integer and b >= 1 is an algebraic number. The constants 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 Pi. The key point is to represent the above sum as a sum of integrals:

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 - x^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_] :=
Sum[a[i] b^i/n * (
Log[b] (1 + 2 Sum[Cos[2 Pi i k/n], {k, 1, Floor[(n-1)/2]}]) -
Log[b - 1] -
(-1)^i If[ EvenQ[n], Log[b+1] - Log[b], 0] -
Sum[Cos[2 Pi i k/n] Log[b^2 - 2 b Cos[2 k Pi/n] + 1], {k, 1, Floor[(n-1)/2]}] +
2 Sum[ Sin[2 i k Pi/n] ArcTan[Sin[2 k Pi/n]/(b - Cos[2 k Pi/n])] , {k, 1, Floor[(n-1)/2]}]),
{i, 1, 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 Pi. Let see what we can learn by trying to get the argument to arctan to be 1. We would need sin(2 k Pi/n)/(b - cos(2 k Pi/n)) = 1, or b = sin(2 k Pi/n) + cos(2 k Pi/n) = Sqrt[2] sin(1/4 + 2 k Pi/n). Since we want to have our base, bn, be an integer, we need to find an integer n so that 2^(n/2) sin(1/4 + 2 k Pi/n)^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 a1 + a2 + a3 + a4 = 0 (one can then view the right side of (*) as Integrate[p(x)/(1 - x^n), x], where p 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.

    a[1]   a[3]     3 a[1]   a[2]   3 a[3]
Pi (---- - ----) + (------ + ---- + ------) Log[2]
     8      8         4       2       4
Setting it equal to Pi produces a one-parameter formula for Pi:

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

since that is how the Leibniz series becomes one-parameter family of formulas. Despite the fact that nothing new about Pi 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 Pi 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.

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

 a[1]   3 a[2]   9 a[3]   27 a[4]   81 a[5]   243 a[6]
(---- + ------ - ------ + ------- + ------- - --------) Log[13]
  4       4        2         4         4         2
system = DeleteCases[Flatten[CoefficientList[rawExpr, transcs]], 0];
First[%]
a[1]   3 a[2]   9 a[3]   27 a[4]   81 a[5]   243 a[6]
---- + ------ - ------ + ------- + ------- - --------
 4       4        2         4         4         2
Array[a, n] /. Solve[system == {1, 0, 0, 0, 0, 0}]
  
  7  1  4   1    7
{{-, -, --, --, ---, 0}}
  3  3  27  27  243
This gives a series for log 13, which we rearrange slightly to clear fractions.

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

rawExpr = CollectPatterns[ GoldenGoose[10, 2] /. Log[9] :> Log[9/10] + Log[10], {_Log}]

 
                        9
(-5 a[1] - 50 a[2]) Log[--] + (-5 a[1] + 50 a[2]) Log[10] + 
                        10
 
  (5 a[1] - 50 a[2]) Log[11]
Setting a[1] = -1/10 and a[2] = -1/100 yields a series that, after simplification, reduces to the Maclaurin expansion of log 9/10 mentioned in section 1. Setting a[1] = -1/5 and a2 = 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.

7. 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 Pi (section 5).

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 Pi, or a base-10 series for Pi. 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 Pi arises.

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

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:

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:

and

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:

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.

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 Pi, American Mathemat ical Monthly 102 (1995) 195Ð203.