On the computation of the n'th decimal digit

of various transcendental numbers.

by Simon Plouffe

November 30, 1996

Revised in March 2003

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, we have not found a faster

algorithm, nor have we proven that one does not exist.

The formula for Pi used is

infinity

----- n 2

\ n 2 (n!)

) ---------- = Pi + 3 (1)

/ (2 n)!

-----

n = 1

* 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, 1/100891344545564193334812497256 =

-3/8 - 61/81 - 1/11 - 11/13 - 4/17 - 9/19 - 25/29 - 26/31 +

23/53 + 41/59 + 29/61 + 37/67 + 33/71 + 19/73 + 36/79 + 7/83

+ 13/89 + 88/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).

2 n

-----

\ / 2 n n \

) | [---] - [---]|

/ | k k |

----- \ p p /

k = 1

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

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

cost is for the computation of the k'th partial sum of

infinity

----- n

\ 2

) ----------- (2)

/ n C(2 n, n)

-----

n = 1

where C(2*n,n) is the central binomial coefficient. 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 :

infinity

----- n

\ c

) --------------- (3)

/ r

----- P(n) C(m n, n)

n = 1

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 Inverter[Inv] 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 highestis 2**k and k is nearly the size of n. For this

particular number only avery 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.

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

Acknowledgments

We wish to thanks, Robert Israel (Univ. British Columbia) and David H.

Bailey [NASA] for their helpful comments.

Bibliography

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

* Plouffe's Inverter at http://pi.lacim.uqam.ca/eng

Keywords : Pi, complexity, algorithm, n'th digit computation, Plouffe's Inverter.