On the computation of the n'th decimal digit
of various numbers.
by Simon Plouffe
November 30, 1996
Research Associate at the Centre for Experimental &
Constructive Mathematics
Simon Fraser University, 8888 University Drive, Burnaby BC, CANADA
e-mail : plouffe@cecm.sfu.ca
Note : For best results when printing this document, set the printer
at 80% size.
Dedicated to all the people that computed the number Pi over the past
2000 years.
Abstract
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.
The algorithm have been improved by Fabrice Bellard to O(n^2).
Also, the Catalan constant can be computed using the Ramanujan
series in O(n^2) time. See bottom for new classes of series.
* Introduction.
* Key observation and the Splitting Algorithm.
* Other numbers.
* Conclusion.
* Bibliography.
Introduction
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 formula
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.
Example:
1/C(100,50) = 1/100891344545564193334812497256 =
1
-------------------------------------------------------
3 4
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
way.
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, we have,
[INLINE]
[INLINE]
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 C(2*n,n) the maximal
exponent is (as Robert Israel pointed out).
[INLINE]
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
memory elements,
* 1 for the partial sum so far. (evaluated later with the BBP
algorithm).
* 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 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 explanations.
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.
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 C*n**3*log(n)**3
steps.
Other Numbers.
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 :
[INLINE]
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 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.
Conclusions
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
efficient algorithm.
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) which falls
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).
Acknowledgments.
We wish to thanks, Fabrice Bellard, Robert Israel (Univ. British
Columbia) , David H. Bailey (NASA Ames), Victor Adamchik from
Wolfram Research and Steve Finch of MathSoft for their useful
discussions.
Bibliography
* Lehmer, D. H. Interesting series involving the central binomial
coefficient. Amer. Math. Monthly 92 (1985), no. 7, 449--457.
* David H. Bailey, Peter B. Borwein and Simon Plouffe, "On The Rapid
Computation of Various Polylogarithmic Constants", to appear in
Mathematics of Computation, vol. 66, no. 218 (April 1997), pp.
903-913. Postscript.
* Bruce C. Berndt, Ramanujan Notebooks vols. I to IV. Springer
Verlag, New York.
* A Passion for Pi an article from Ivars Peterson, Science
News/Mathland chronicle
* The Miraculous Bailey-Borwein-Plouffe Pi algorithm.
* Pi: A 2000-Year Search Changes Direction
* 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
communication.
* 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.