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

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

**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 11These 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.

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

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 16Good. 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!

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.

**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]]) / 8Only 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

**transcs = GetTranscendentals[simpler]**

{Pi, ArcTan[2], ArcTan[Sqrt[2]], Log[2], Log[3], Log[5], Log[1 + Sqrt[2]]}

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.

**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 2The 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} 2In 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 4Letting 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.

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

*LogOfRadical*s 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

e a Cos[r] + (1 - -) Log[1 - s] + (1 + -) Tan[x] 2 3It is useful to map

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

(2 - e) Log[1 - s] (3 + a) Tan[x] Cos[r] + ------------------ + -------------- 2 3

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

Log[2] - Log[3] + Log[13]

**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] ------------------------------- 4Things 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

1 1 1 {{2, 1, 1, 2}, {2, 0, 0, 0}, {1, -, -(-), -1}, {0, 0, -(-), 0}} 2 2 2The 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.

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

**(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

3 {{0, 3, 0, -, 0, 0}} 2So 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.]**

TrueTo illustrate the general nature of

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

{3, 0, 0, 5, 7}

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 4Setting 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

a[1] 3 a[2] 9 a[3] 27 a[4] 81 a[5] 243 a[6] ---- + ------ - ------ + ------- + ------- - -------- 4 4 2 4 4 2

7 1 4 1 7 {{-, -, --, --, ---, 0}} 3 3 27 27 243This 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.

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.

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