computation of the n'th decimal digit
by Simon Plouffe
- We outline a method for computing the n'th decimal (or any
other base) digit of Pi in C*n^3*log(n)^3 time and with very
little memory. The computation is based on the recently discovered
Bailey-Borwein-Plouffe algorithm and the use of a new algorithm
that simply splits an ordinary fraction into its components. The
algorithm can be used to compute other numbers like Zeta(3),
Pi*sqrt(3), Pi**2 and 2/sqrt(5)*log(tau) where tau is the golden
ratio. The computation can be achieved without having to compute
the preceeding digits. We claim that the algorithm has a more
theoretical rather than practical interest, we have not found a
faster algorithm, nor have we proven that one does not exist.
The formula for Pi used is
- Key observation and the Splitting Algorithm.
- Other numbers.
The computation of the n'th digit of irrational or transcendental
numbers was considered either impossible or as difficult to compute
as the number itself. Last year (1995) [BBP], have found a
new way of computing the n'th binary digit of various constants like
Pi and log(2). An intensive computer search was then carried out to
find if that algorithm could be used to compute a number in an
arbitrary base. We present here a way of computing the n'th decimal
digit of Pi (or any other base) by using more time than the
[BBP] algorithm but still with very little memory.
Key observation and
- The observation is that a fraction 1/(a*b) can be split into
k1/a + k2/b by using the continued fraction algorithm of a/b. Here
a and b are two prime powers. This is equivalent to having to
solve a diophantine equation for k1 and k2 - it is always possible
to do so if (a,b) = 1 , see [HW] if they have no common
factor. If we have more than 2 prime factors then it can be done
by doing 2 at the time and then using the result to combine with
the third element. This way an arbitrary big integer M can be
split into small elements. If we impose the conditions on M of
having only small factors (meaning that the biggest prime power
size is of the order of a computer word), then an arbitrary M can
be represented. If this is true then a number of known series and
numbers can then be evaluated. For example the expression
1/C(2*n,n), the central binomials satisfies that : the prime
powers of this number are small when n is big.
- 1/binomial(100,50) = 1/100891344545564193334812497256 =
2 * 3 *11*13*17*19*29*31*53*59*61*67*71*73*79*83*89*97
- Now if we take 2 elements at the time and solve the simple
diophantine equation and proceed this way
- 1) 1/(a*b) = k1/a + k2/b
- 2) (k1/a + k2/b)/c = m1/a + m2/b + m3/c
- 3) we proceed with the next element.
- At each step the constants k1 and k2 are determined by simply
expanding a/b into a continued fraction and keeping the 'before
last' continuant, later m1, m2 and m3 are determined the same
- having finished with that number we quickly arrive at a number
which is (modulo 1) the same number but represented as a sum of
only small fractions.
- So, 1/100891344545564193334812497256 =
61 11 25 26 23 41 29
-3/8 - ---- -1/11 - ---- - 4/17 -9/19 - ---- - ---- + ---- + ---- + ----
81 13 29 31 53 59 61
37 33 19 36 13 88
+ ---- + ---- + ---- + ---- + 7/83 + ---- + ---- .
67 71 73 79 89 97
- The time taken to compute this expression is log(n)*n**2,
log(n) being the time spent to compute with the euclidian
algorithm on each number. We did not take into account the time
spent on finding what is the next prime in the expression simply
because we can consider (at least for the moment) that the
applicability of the algorithm is a few thousands digits and so
the time to compute a prime is really a matter of a few seconds in
that range for the whole process. Since we know by advance what is
the maximal prime there could be in C(2*n,n) then we can do it
with a greedy algorithm that pulls out the factors until we reach
2*n, and this can be done without having to compute the actual
number which would obviously not fit into a small space. It can be
part of the loop without having to store any number apart from the
current n. For any p in binomial(2*n,n) the maximal exponent is
(as Robert Israel pointed out).
- Equivalently, for p = 2 it gives the number of '1' in the
binary expansion of n, for p = 3 there is another clue with the
ternary expansion of the number and the number of times the
pattern '12' appears.
- Now looking at the sum(1/C(2*n,n),n=1..infinity), we can say
that the series is 'essentially' Pi*sqrt(3) since it differs only
by 4/9*Pi*sqrt(3) + 1/3, since these are 2 small rationals we can
use BBP algorithm to carry the computation to an arbitrary
position in almost no time. Having n/C(2*n,n) instead of (1) only
simplifies the process.
To compute the final result of each term we need only few
- 1 for the partial sum so far. (evaluated later with the BBP
- 4 for the current fractions k1/a and k2/b.
- 2 for the next element to be evaluated : 1/c.
- 1 for n itself.
- So with as little as 8 memory elements the sum for each term
of (1) can be carried out a without having to store any number
greater than a computer word in log(n) time, adding this for each
element the total cost for (1) is then n^3*log(n).
- The next thing we have to consider is that , if we have an
arbitrary large M and if M has only small factors then 1/M can be
computed. First, we need to represent 1/M as
Sum(a(i)/p(i)**(j),i=1..k) (2) where p(i)^(j) is a prime power and
a(i) is smaller than p(i)^(j).
- If we have 2^n/M then by using the binary method on
each element of the representation of 1/M with (2) is
possible in log(n) time. Again if we don't want to store the
elements of (2) in memory we can do it as we do the computation of
the first part at each step. In this algorithm we can either store
the powers of 2 to do the binary method or not. There is variety
of ways to do it, we refer to [Knuth vol. 2] for
- This step is important, essentially once we can represent
1/(a*b) by splitting them then to multiply by 2^n only adds log(n)
steps for each element and it can be done in arbitrary base since
we have the actual fraction for each element of (2). It only
pushes the decimal (or the the decimal point of the base chosen),
further. At any moment only one element in the expansion of 1/M is
considered with the current fraction, that same fraction can be
represented in base 10 at anytime if we want the decimal expansion
at that point. For this reason multiplying the current fraction by
2^n involves only small numbers and fractions.
- Once this is done, the total cost becomes n^3*log(n)^2. This
cost is for the computation of the k'th partial sum of
- If we want at each step to compute (the final n'th digit) then
we need log(n) steps to do it. It can be done in any base chosen
in advance, in BBP the computation could be done in base 2 but
here we have the actual explicit fraction which is independant of
the base. This is where we actually compute the decimal expansion
of the final fraction of the process.
- So finally the n'th digit of Pi can be computed in
- By looking at the plethora of formulas of the same type as (1)
or (3) we see that [PIagm], [RamI and IV] the
numbers Pi*sqrt(3) Pi, Pi**3, Zeta(3) and even powers of Pi can be
computed as well. The condition we need to enseure is :
- if any term of a series can be split into small fractions of
size no greater than that of a computer word, then it is part of
that class. This includes series of the type :
- where c is an integer, P(n) a polynomial and C(mn,n) is a
near central binomial coefficient. This class of series
contains many numbers that are not yet identified in terms of
known constants and conversely the known constants that are of
similar nature like Zeta(5) have not yet been identified as
members of the class. The process of identifying a series as being
expressed in terms of known constants and the exact reverse
process is what the Inverse
Symbolic Calculator tries to do.
- The number e or exp(1) which is sum(1/n!,n=0..infinity) does
not satisfies our condition because 1/n! eventually contains high
powers of 2, therefore cannot be computed to the n'th digit using
our algorithm. The factorisation of 1/n! has high powers of small
primes, the highest is 2**k and k is nearly the size of n. For
this particular number only a very few series are known and appear
to be only a variation on that first one.
- Others like gamma or Catalan do not seem to have a proper
series representation and computer search using Bailey's PSLQ or
LLL with MapleV and Pari-Gp gave no answer to this. Algebraic
numbers like sqrt(2) have not been yet been fully investigated and
we still do not know if those would fall into this category.
- There are many, but first and foremost we cannot resist
thinking at William Shanks who did the computation of Pi by hand
in 1853 - if he would had known this algorithm, he would have
certainly tried it before spending 20 years of his life computing
Pi (half of it on a mistake).
- Secondly, the algorithm shown here is theoretical and not
practical. We do not know if there is a way to improve it, and if
so then it is reasonable to think that it could then be used to
check long computations like the one that Yasumasa Kanada
conducted last year for the computation of Pi to 6.4 billion
digits. There could be a way to speed the algorithm to make it an
- Thirdly, so far there are 2 classes of numbers that can be
computed to the n'th digit : 1) the SC(2) class as in the
[BBP] algorithm which includes various polylogarithms. 2)
this new class of numbers. Now what's next ?, so
- far we do not know wheter, for example, series whose general
term is H(n)/2**n (where H(n) is the n'th harmonic number) whcich
fall into the first class, can be extended. We think that this new
approach is only the tip of the iceberg.
- Finally, it is interesting to observe that we can then compute
Pi to the 10000'th digit without having to store (hardly) any
array or matrix, so it can be computed using a small pocket
calculator. We also note that, in some way we have a way to
produce the digits of Pi without using memory, this means that the
number is compressible , if we consider that we could use
the algorithm to produce a few thousands digits of the number. We
think that other numbers are yet to come and that there is a
possibility (?) of having a direct formula for the n'th digit (in
any base) of a naturally occuring constant like log(2).
- We wish to thanks, Robert Israel (Univ. British Columbia) ,
David H. Bailey (NASA Ames) and Peter B. Borwein (Simon Fraser
Univ.) for their useful discussions.
- David H. Bailey, Peter B. Borwein and Simon Plouffe, On The
Rapid Computation of Various Polylogarithmic Constants, april
1997 in Mathematics of Computation.
- Bruce C. Berndt, Ramanujan Notebooks vols. I to IV. Springer
Verlag, New York.
- Pi and the AGM, by Jonathan M. Borwein and Peter B. Borwein, A
Study in Analytic Number Theory and Computational Complexity,
Wiley, N.Y., 1987.
- David H. Bailey and Simon Plouffe, Recognizing Numerical
Constants, preprint 1995.
- Robert Israel at University of British Columbia, personal
- M. Abramowitz and I. Stegun, Handbook of Mathematical
Functions, Dover, New York, 1964.
- G. H. Hardy and E. M. Wright, An Introduction to the Theory
of Numbers 5e, Oxford University Press, 1979.
- W. Shanks, Contributions to Mathematics Comprising Chiefly of
the Rectification of the Circle to 607 Places of Decimals, G.
Bell, London, 1853.
- D.E. Knuth, The Art of Computer Programming, Vol. 2:
Seminumerical Algorithms, Addison-Wesley, Reading, MA, 1981.