ON THE RAPID COMPUTATION OF
VARIOUS POLYLOGARITHMIC CONSTANTSDavid Bailey, Peter Borwein1 and Simon Plouffe
Abstract.
We give algorithms for the computation of the d-th digit of certain transcendental
numbers in various bases. These algorithms can be easily implemented (multiple
precision arithmetic is not needed), require virtually no memory, and feature runtimes that scale nearly linearly with the order of the digit desired. They makeit feasible to compute, for example, the billionth binary digit of log (2) or π on amodest work station in a few hours run time.
We demonstrate this technique by computing the ten billionth hexadecimal digit ofπ, the billionth hexadecimal digits of π 2 , log(2) and log2 (2), and the ten billionth
decimal digit of log(9/10).
These calculations rest on the observation that very special types of identities exist for
certain numbers like π, π 2 , log(2) and log2 (2). These are essentially polylogarithmic
ladders in an integer base. A number of these identities that we derive in this work
appear to be new, for example the critical identity for π:
π=

∞

1  4
2
1
1 
−
−
−
.
i 8i + 1
16
8i
+
4
8i
+
5
8i
+6
i=0

Ref: Mathematics of Computation, vol. 66, no. 218 (April 1997), pg. 903–913

1 Research supported in part by NSERC of Canada.
1991 Mathematics Subject Classiﬁcation. 11A05 11Y16 68Q25.
Key words and phrases. Computation, digits, log, polylogarithms, SC, π, algorithm.

Typeset by AMS-TEX

1

1. Introduction.
It is widely believed that computing just the d-th digit of a number like π is reallyno easier than computing all of the ﬁrst d digits. From a bit complexity point ofview this may well be true, although it is probably very hard to prove. What we willshow is that it is possible to compute just the d-th digit of many transcendentals in(essentially) linear time and logarithmic space. So while this is not of fundamentallylower complexity than the best known algorithms (for say π or log(2)), this makes
such calculations feasible on modest workstations without needing to implement
arbitrary precision arithmetic.
We illustrate this by computing the ten billionth hexadecimal digit of π, the billionth
hexadecimal digits of π 2 , log(2) and log2 (2), and the ten billionth decimal digit oflog(9/10). Details are given in Section 4. A previous result in this same spirit
is the Rabinowitz-Wagon “spigot” algorithm for π. In that scheme, however, thecomputation of the digit at position n depends on all digits preceding position n.
We are interested in computing in polynomially logarithmic space and polynomial
time. This class is usually denoted SC (space = logO(1) (d) and time = dO(1) whered is the place of the “digit” to be computed). Actually we are most interested in
the space we will denote by SC∗ of polynomially logarithmic space and (almost)linear time (here we want the time = O(d logO(1) (d))). There is always a possible
ambiguity when computing a digit string base b in distinguishing a sequence ofdigits a(b − 1)(b − 1)(b − 1) from (a + 1)000. In this particular case we consider
either representation as an acceptable computation. In practice this problem does
not arise.
It is not known whether division is possible in SC, similarly it is not known whether
base change is possible in SC. The situation is even worse in SC∗ , where it is noteven known whether multiplication is possible. If two numbers are in SC∗ (in thesame base) then their product computes in time = O(d2 logO(1) (d)) and is in SC
but not obviously in SC∗ . The d2 factor here is present because the logarithmic
space requirement precludes the usage of advanced multiplication techniques, suchas those based on FFTs.
We will not dwell on complexity issues except to point out that diﬀerent algorithms
are needed for diﬀerent bases (at least given our current ignorance about basechange) and very little closure exists on the class of numbers with d-th digit computable in SC. Various of the complexity related issues are discussed in [6,8,9,11,14].
As we will show in Section 3, the class of numbers we can compute in SC∗ in baseb includes all numbers of the form
(1.1)

∞

p(k)
,
ck
b q(k)

k=1

where p and q are polynomials with integer coeﬃcients and c is a positive integer.
Since addition is possible in SC∗ , integer linear combinations of such numbers arealso feasible (provided the base is ﬁxed).
The algorithm for the binary digits of π, which also shows that π is in SC∗ in base2, rests on the following remarkable identity:
2

Theorem 1. The following identity holds:
(1.2)

π=

∞

2
1
1 
1  4
−
−
−
.
i
16 8i + 1 8i + 4 8i + 5 8i + 6
i=0This can also be written as:
(1.3)

π=

∞


pi

 8i  i
i=1 16[pi ] = [4, 0, 0, −2, −1, −1, 0, 0],

where the overbar notation indicates that the sequence is periodic.
Proof. This identity is equivalent to:√ 1/√2 √4 2 − 8x3 − 4 2x4 − 8x5
(1.4)
π=
dx,
1 − x80√which on substituting y := 2x becomes
 1
π=

0

16 y − 16dy .
y4 − 2 y3 + 4 y − 4

The equivalence of (1.2) and (1.4) is straightforward. It follows from the identity 1/√20xk−1dx =1 − x8 1/√2 ∞0

xk−1+8i dx

i=0

∞
1 
1
=√ k
i (8i + k)
16
2 i=0

That the integral (1.4) evaluates to π is an exercise in partial fractions most easily
done in Maple or Mathematica. This proof entirely conceals the route to discovery. We found the identity (1.2) bya combination of inspired guessing and extensive searching using the PSLQ integer
relation algorithm [3,12].
Shortly after the authors originally announced the result (1.2), several colleagues,
including Helaman Ferguson, Tom Hales, Victor Adamchik, Stan Wagon, Donald
Knuth and Robert Harley, pointed out to us other formulas for π of this type. One
intriguing example is
π=

∞

2
1
1/2
1/2
1/4 
1  2
+
+
−
−
−
,
i 8i + 1
16
4i
+
2
4i
+
3
4i
+
5
4i
+
6
4i
+7
i=0

which can be written more compactly asπ=∞(−1)i 
i=0

4i

2
2
1 
+
+
.
4i + 1 4i + 2 4i + 3
3

In [2], this and some related identities are derived using Mathematica.
As it turns out, these other formulas for π can all be written as formula (1.2) plus
a rational multiple of the identity0=∞
1  −8
8
4
8
2
2
1 
+
+
+
+
+
−
.
i
16 8i + 1 8i + 2 8i + 3 8i + 4 8i + 5 8i + 6 8i + 7
i=0

The proof of this identity is similar to that of Theorem 1.
The identities of the next section and Section 5 show that, in base 2, π 2 , log2 (2)and various other constants, including {log(2), log(3), . . . , log(22)} are in SC∗ . (We
don’t know however if log(23) is even in SC.)We will describe the algorithm in the Section 3. Complexity issues are discussed
in [3,5,6,7,8,9,14,19,21] and algorithmic issues in [5,6,7,8,14]. The requisite special
function theory may be found in [1,5,15,16,17,20].
2. Identities.
As usual, we deﬁne the m-th polylogarithm Lm by(2.1)

Lm (z) :=

∞

zi
i=1

im

,

|z| < 1 .

The most basic identity is− log(1 − 2−n ) = L1 (1/2n )(2.2)

which shows that log(1 − 2−n ) is in SC∗ base 2 for integer n. (See also section 5.)
Much less obvious are the identities(2.3)

π 2 = 36L2 (1/2) − 36L2 (1/4) − 12L2 (1/8) + 6L2 (1/64)

and
(2.4)

log2 (2) = 4L2 (1/2) − 6L2 (1/4) − 2L2 (1/8) + L2 (1/64) .

These can be written as
π 2 = 36

(2.5)

∞

ai
i=1

(2.6).

log2 (2) = 2

2i i2

∞

bi
,
i i2
2
i=1

,

[ai ] = [1, −3, −2, −3, 1, 0][bi ] = [2, −10, −7, −10, 2, −1] .

Here the overline notation indicates that the sequences repeat. Thus we see that
π 2 and log2 (2) are in SC∗ in base 2. These two formulas can alternately be written
π2 =

∞


16
9 1 
24
8
6
1
−
−
−
+
i
2
2
2
2
2
8 i=0 64 (6i + 1)
(6i + 2)(6i + 3)(6i + 4)(6i + 5)4

log2 (2) =∞1  1  −16
16
40
14
10
1
.
+
−
−
−
+
i
2
2
2
2
2
2
8 i=0 64 (6i) (6i + 1) (6i + 2) (6i + 3) (6i + 4) (6i + 5)

Identities (2.3)-(2-6) are examples of polylogarithmic ladders in the base 1/2 inthe sense of [16]. As with (1.2) we found them by searching for identities of thistype using an integer relation algorithm. We have not found them directly in print.
However (2.5) follows from equation (4.70) of [15] with α = π/3, β = π/2 and γ =π/3. Identity (2.6) now follows from the well known identity
12L2 (1/2) = π 2 − 6 log2 (2) .

(2.7)

A distinct but similar formula that we have found for π 2 is
π2 =

∞


1 
16
16
8
16
4
4
2
,
−
−
−
−
−
+
i
2
2
2
2
2
2
2
16 (8i + 1) (8i + 2) (8i + 3) (8i + 4) (8i + 5) (8i + 6) (8i + 7)i=0which can be derived from the methods of section 1.
There are several ladder identities involving L3 :(2.8)

35/2ζ(3) − π 2 log(2) = 36L3 (1/2) − 18L3 (1/4) − 4L3 (1/8) + L3 (1/64) ,(2.9)

2 log3 (2) − 7ζ(3) = −24L3 (1/2) + 18L3 (1/4) + 4L3 (1/8) − L3 (1/64) ,(2.10) 10 log3 (2)−2π 2 log(2) = −48L3 (1/2)+54L3 (1/4)+12L3 (1/8)−3L3 (1/64) .
The favored algorithms for π of the last centuries involved some variant of Machin’s1706 formula:11π= 4 arctan − arctan
.
4
5
239

(2.11)

There are many related formula [15,16,17,20] but to be useful to us all the argumentsof the arctans have to be a power of a common base, and we have not discovered
any such formula for π . One can however writeπ
11= 2 arctan √ + arctan √228(2.12)
This can be written as(2.13)√

2π = 4f (1/2) + f (1/8)

where

f (x) :=∞
(−1)i xi
i=1

and allows for the calculation of√

2i + 12π in SC∗ .

Another two identities involving Catalan’s constant G, π and log(2) are:(2.14)

G−∞
π log(2) 
ci
=
,
i+1

8
2  i2i=1 2
5[ci ] = [1, 1, 1, 0, −1, −1, −1, 0]

and
∞
di
5 2 log2 (2) 
π −
=
,
i+1

96
8
2  i2i=1 2(2.15)[di ] = [1, 0, −1, −1, −1, 0, 1, 1]

These may be found in [17 p. 105, p. 151]. Thus 8G − π log(2) is also in SC∗ inbase 2, but it is open and interesting as to whether G is itself in SC∗ in base 2.
A family of base 2 ladder identities exist:
Lm (1/64)
Lm (1/8) 2 Lm (1/4) 4 Lm (1/2) 5 (− log(2))
−
−
−
+
6m−1
3m−1
2m−1
9
9 m!

(2.16)m−2+

π 2 (− log(2))
54 (m − 2)!

−m−4

π 4 (− log(2))
486 (m − 4)!

−403 ζ(5) (− log(2))
1296 (m − 5)!

m

m−5

= 0.

The above identity holds for 1 ≤ m ≤ 5; when the arguments to factorials arenegative they are taken to be inﬁnite so the corresponding terms disappear. See[16, p. 45].
As in the case of formula (1.2) for π, colleagues of the authors have subsequently
pointed out several other formulas of this type for various constants. Three examples reported by Knuth, which are based on formulas in [13, p. 17, 18, 22, 47, 139],
are
∞

√
√
1  1
1/2
1/4
1/8 
+
+
+
2 ln(1 + 2) =i 8i + 1
16
8i
+
3
8i
+
5
8i
+7
i=0
√

∞

√
1  1
1/2
1/4
1/8 
−
+
−
2 arctan(1/ 2) =i
16 8i + 1 8i + 3 8i + 5 8i + 7
i=0

arctan(1/3) =

∞

1
1/2
1/4 
1  1
−
−
−
i
16 8i + 1 8i + 2 8i + 4 8i + 5
i=0Thus these constants are also in class SC∗ . Some other examples can be found in[18].
3. The Algorithm.
Our algorithm to compute individual base-b digits of certain constants is based on
the binary scheme for exponentiation, wherein one evaluates xn rapidly by successive squaring and multiplication. This reduces the number of multiplications toless than 2 log2 (n). According to Knuth [14], where details are given, this trickgoes back at least to 200 B.C. In our application, we need to perform exponentiation modulo a positive integer c, but the overall scheme is the same — one merely
performs all operations modulo c. An eﬃcient formulation of this algorithm is asfollows.
To compute r = bn mod c, ﬁrst set t to be the largest power of two ≤ n, and set
r = 1. Then
A: if n ≥ t then r ← br mod c;n ← n − t;
6

endif

t ← t/2
if t ≥ 1 then r ← r2 mod c;

go to A;

endif

Here and in what follows, “mod” is used in the binary operator sense, namely as thebinary function deﬁned by x mod y := x − [x/y]y. Note that the above algorithm
is entirely performed with positive integers that do not exceed c2 in size. Thus it
can be correctly performed, without round-oﬀ error, provided a numeric precisionof at least 1 + 2 log2 c bits is used.
Consider now a constant deﬁned by a series of the form
∞

1
S=
,
bck p(k)k=0where b and c are positive integers and p(k) is a polynomial with integer coeﬃcients.
First observe that the digits in the base b expansion of S beginning at position n+1
can be obtained from the fractional part of bn S. Thus we can write∞ n−ckbn(3.4)
b S mod 1 =
mod 1
p(k)
k=0

n/c

=

 b

k=0

n−ck

mod p(k)mod 1 +p(k)∞
k=n/c+1

bn−ck
mod 1
p(k)

For each term of the ﬁrst summation, the binary exponentiation scheme is used
to evaluate the numerator. Then ﬂoating-point arithmetic is used to perform thedivision and add the result to the sum mod 1. The second summation, where
the exponent of b is negative, may be evaluated as written using ﬂoating-point
arithmetic. It is only necessary to compute a few terms of this second summation,
just enough to insure that the remaining terms sum to less than the “epsilon” ofthe ﬂoating-point arithmetic being used. The ﬁnal result, a fraction between 0 and
1, is then converted to the desired base b.
Since ﬂoating-point arithmetic is used here in divisions and in addition modulo 1,
the result is of course subject to round-oﬀ error. If the ﬂoating-point arithmetic
system being used has the property that the result of each individual ﬂoating-point
operation is in error by at most one bit (as in systems implementing the IEEE
arithmetic standard), then no more than log2 (2n) bits of the ﬁnal result will becorrupted. This is actually a generous estimate, since it does not assume any
cancelation of errors, which would yield a lower estimate. In any event, it is clearthat ordinary IEEE 64-bit arithmetic is suﬃcient to obtain a numerically signiﬁcantresult for even a large computation, and “quad precision” (i.e. 128-bit) arithmetic,
if available, can insure that the ﬁnal result is accurate to several digits beyond
the one desired. One can check the signiﬁcance of a computed result beginningat position n by also performing a computation at position n + 1 or n − 1 andcomparing the trailing digits produced.
The most basic interesting constant whose digits can be computed using this schemeis∞1log(2) =k2k
k=1
7

in base 2. Using this scheme to compute hexademical digits of π from identity (1.2)
is only marginally more complicated, since one can rewrite formula (1.2) using foursums of the required form. Details are given in the next section. In both cases, inorder to compute the n-th binary digit (or a ﬁxed number of binary digits at then-th place) we must sum O(n) terms of the series. Each term requires O(log(n))
arithmetic operations and the required precision is O(log(n)) digits. This gives
a total bit complexity of O(n log(n)M (log(n))) where M (j) is the complexity ofmultiplying j bit integers. So even with ordinary multiplication the bit complexity
is O(n log3 (n)).
This algorithm is, by a factor of log(log(log(n))), asymptotically slower than thefastest known algorithms for generating the n-th digit by generating all of the ﬁrst
n digits of log(2) or π [7]. The asymptotically fastest algorithms for all the ﬁrst
n digits known requires a Strassen-Schönhage multiplication [19]; the algorithms
actually employed use an FFT based multiplication and are marginally slower than
our algorithm, from a complexity point of view, for computing just the n-th digit. Of
course this complexity analysis is totally misleading: the strength of our algorithm
rests mostly on its easy implementation in standard precision without requiring
FFT methods to accelerate the computation.
It is clear that the above methods can easily be extended to evaluate digits ofcontstants deﬁned by a formula of the form
∞

p(k)
,
S=
ck
b q(k)
k=0

where p and q are polynomials with integer coeﬃcients and c is a positive integer.
Similarly if p and q are slowly growing analytic functions of various types themethod extends.
4. Computations.
We report here computations of π, log(2), log2 (2), π 2 and log(9/10), based on
the formulas (1.1), (2.2), (2.5), (2.6) and the identity log(9/10) = −L1 (1/10),
respectively.
Each of our computations employed quad precision ﬂoating-point arithmetic fordivision and sum mod 1 operations. Quad precision is supported from Fortran on
the IBM RS6000/590 and the SGI Power Challenge (R8000), which were employedby the authors in these computations. We were able to avoid the usage of explicit
quad precision in the exponentiation scheme by exploiting a hardware feature common to these two systems, namely the 106-bit internal registers in the multiply-add
operation. This saved considerable time, because quad precision operations are
signiﬁcantly more expensive than 64-bit operations.
Computation of π 2 and log2 (2) presented a special challenge, because one must
perform the exponentiation algorithm modulo k 2 instead of k. When n is larger than
only 213 , some terms of the series (2.5) and (2.6) must be computed with a modulus
k 2 that is greater than 226 . Squares that appear in the exponentiation algorithm willthen exceed 252 , which is the nearly the maximum precision of IEEE 64-bit ﬂoatingpoint numbers. When n is larger than 226 , then squares in the exponentiation
algorithm will exceed 2104 , which is nearly the limit of quad precision.
8

This diﬃculty can be remedied using a method which has been employed for example in searches for Wieferich primes [10]. Represent the running value r in the
exponentiation algorithm by the ordered pair (r1 , r2 ), where r = r1 + kr2 , andwhere r1 and r2 are positive integers less than k. Then one can write
r2 = (r1 + kr2 )2 = r12 + 2r1 r2 k + r22 k 2
When this is reduced mod k 2 , the last term disappears. The remaining expression isof the required ordered pair form, provided that r12 is ﬁrst reduced mod k, the carryfrom this reduction is added to 2r1 r2 , and this sum is also reduced mod k. Notethat this scheme can be implemented with integers of size not exceeding 2k 2 . Since
the computation of r2 mod k 2 is the key operation of the binary exponentiation
algorithm, this means that ordinary IEEE 64-bit ﬂoating-point arithmetic can beused to compute the n-th hexadecimal digit of π 2 or log2 (2) for n up to about
224 . For larger n, we still used this basic scheme, but we employed the multiplyadd “trick” mentioned above to avoid the need for explicit quad precision in thissection of code.
Our results are given below. The ﬁrst entry, for example, gives the 106 -th through106 +13-th hexadecimal digits of π after the “decimal” point. In all cases we did thecalculations twice — the second calculation was similar to the ﬁrst, except shiftedback one position. Since this changes all the arithmetic performed, it is a highly
rigorous validity check. Thus we believe that all the digits shown below are correct.

9

Constant:

Base:

Position:

Digits from Position:

π

16

106
107
108
109
1010

26C65E52CB4593
17AF5863EFED8D
ECB840E21926EC
85895585A0428B
921C73C6838FB2

log(2)

16

106
107
108
109

418489A9406EC9
815F479E2B9102
E648F40940E13E
B1EEF1252297EC

π2

16

106
107
108
109

685554E1228505
9862837AD8AABF
4861AAF8F861BE
437A2BA4A13591

log2 (2)

16

106
107
108
109

2EC7EDB82B2DF733374B47882B32
3F55150F1AB3DC
8BA7C885CEFCE8

log(9/10)

10

106
107
108
109
1010

80174212190900
21093001236414
01309302330968
44066397959215
82528693381274

These computations were done at NASA Ames Research Center, using workstation
cycles that otherwise would have been idle.
5. Logs in base 2.
It is easy to compute, in base 2, the d-th binary digit oflog(1 − 2−n ) = L1 (1/2n ) .

(5.1)

So it is easy to compute log(m) for any integer m that can be written as
(5.2)

m :=

(2a1 − 1)(2a2 − 1) · · · (2ah − 1).
(2b1 − 1)(2b2 − 1) · · · (2bj − 1)

In particular the n-th cyclotomic polynomial evaluated at 2 is so computable. A
check shows that all primes less than 19 are of this form. The beginning of this listis:{2, 3, 5, 7, 11, 13, 17, 31, 43, 57, 73, 127, 151, 205, 257} .
10

Sinceand since 7,
from√

218 − 1 = 7 · 9 · 19 · 73,
9 and 73 are all on the above list we can compute log(19) in SC∗
log(19) = log(218 − 1) − log(7) − log(9) − log(73).

Note that 211 − 1 = 23 · 89 so either both log(23) and log(89) are in SC∗ or neitheris.
We would like to thank Carl Pomerance for showing that an identity of type (5.2)
does not exist for 23. This is a consequence of the fact that each cyclotomic polynomial evaluated at two has a new distinct prime factor. We would also like tothank Robert Harley for pointing out that 29 and 37 are in SC ∗ in base 2 via
consideration of the Aurefeuillian factors 22n−1 + 2n + 1 and 22n−1 − 2n + 1.
6. Relation Bounds.
One of the ﬁrst questions that arises in the wake of the above study is whether there
exists a scheme of this type to compute decimal digits of π. At present we knowof no identity like (1.2) in base 10. The chances that there is such an identity aredimmed by some numerical results that we have obtained using the PSLQ integer
relation algorithm [3, 12]. These computations establish (with the usual provisosof computer “proofs”) that there are no identities (except for the case n = 16) of
the form


∞
a2
a3
am+1
a1
1  1
+
+ ··· +
π=
+
,
a0
a0
nk mk + 1 mk + 2
mk + mk=0

where n ranges from 2 to 128, where m ranges from 1 to min(n, 32), and where
the Euclidean norm of the integer vector (a0 , a1 , · · · , am+1 ) is 1012 or less. These
results of course do not have any bearing on the possibility that there is a formula
not of this form which permits computation of π in some non-binary base.
In fact, J. P. Buhler has reported a proof that any identity for π of the above form
√ K
must have n = 2K or n = 2 . This also does not exclude more complicatedformulae for the computation of π base 10.
7. Questions.
As mentioned in the previous section, we cannot at present compute decimal digitsof π by our methods because we know of no identity like (1.2) in base 10. But it
seems unlikely that it is fundamentally impossible to do so. This raises the following
obvious problem:
1] Find an algorithm for the n-th decimal digit of π in SC∗ . It is not even clearthat π is in SC in base 10 but it ought to be possible to show this.
2] Show that π is in SC in all bases.
√
3] Are e and 2 in SC (SC∗ ) in any base?
Similarly the treatment of log is incomplete:11

4] Is log(2) in SC∗ in base 10?
5] Is log(23) in SC∗ in base 2?
8. Acknowledgments.
The authors wish to acknowledge the following for their helpful comments: V.
Adamchik, J. Borwein, J. Buhler, R. Crandall, H. Ferguson, T. Hales, R. Harley,
D. Knuth, C. Pomerance and S. Wagon.

12References
1. M. Abramowitz & I.A. Stegun, Handbook of Mathematical Functions, Dover, New York, NY,1965.
2. V. Adamchik & S. Wagon, Pi: A 2000-year search changes direction (preprint).
3. A. V. Aho, J.E. Hopcroft, & J. D. Ullman, The Design and Analysis of Computer Algorithms,
Addison-Wesley, Reading, MA, 1974.
4. D. H. Bailey, J. Borwein and R. Girgensohn, Experimental evaluation of Euler sums, Experimental Mathematics 3 (1994), 17–30.
5. J. Borwein, & P Borwein, Pi and the AGM – A Study in Analytic Number Theory andComputational Complexity, Wiley, New York, NY, 1987.
6. J. Borwein & P. Borwein, On the complexity of familiar functions and numbers, SIAM Review
30 (1988), 589–601.
7. J. Borwein, P. Borwein & D. H. Bailey, Ramanujan, modular equations and approximations
to pi, Amer. Math. Monthly 96 (1989), 201–219.
8. R. Brent, The parallel evaluation of general arithmetic expressions, J. Assoc. Comput. Mach.
21 (1974), 201–206.
9. S. Cook, A taxonomy of problems with fast parallel algorithms, Information and Control 64(1985), 2–22.
10. R. Crandall, K. Dilcher, and C. Pomerance, A search for Wieferich and Wilson primes(preprint).
11. R. Crandall and J. Buhler, On the evaluation of Euler sums, Experimental Mathematics 3,(1995), 275–285.
12. H. R. P. Ferguson & D. H. Bailey, Analysis of PSLQ, an integer relation algorithm (preprint).
13. E. R. Hansen, A Table of Series and Products, Prentice-Hall, Englewood Cliﬀs, NJ, 1975.
14. D. E. Knuth, The Art of Computer Programming. Vol. 2: Seminumerical Algorithms, AddisonWesley, Reading, MA, 1981.
15. L. Lewin, Polylogarithms and Associated Functions, North Holland, New York, 1981.
16. L. Lewin, Structural Properties of Polylogarithms, Amer. Math. Soc., RI., 1991.
17. N. Nielsen, Der Eulersche Dilogarithmus, Halle, Leipzig, 1909.
18. S. D. Rabinowitz and S. Wagon, A spigot algorithm for pi, Amer. Math. Monthly 103 (1995),
195-203.
19. A. Schönhage, Asymptotically fast algorithms for the numerical multiplication and divisionof polynomials with complex coeﬃcients, in: EUROCAM (1982) Marseille, Springer Lecture
Notes in Computer Science, vol. 144, 1982, pp. 3–15.
20. J. Todd, A problem on arc tangent relations, MAA Monthly 56 (1940), 517–528.
21. H. S. Wilf, Algorithms and Complexity, Prentice Hall, Englewood Cliﬀs, NJ, 1986.

13