Click here for a pdf version of this page.
Is it possible to compute directly the n-th digit of p without computing all its first n digits ? At first sight, the problem seems of the same complexity. Recently [1], a formula which permits to find directly the n-th bit of p with very little memory and in quasi-linear time was exhibited. In other words, the question has a positive answer in base 2.
1 Computing the n-th bit of p |
In [1], David Bailey, Peter Borwein and Simon Plouffe give the following formula
| (1) |
and use that to compute the n-th bit of p without computing all its first n bits. More precisely, their method permits to obtain the n-th bit of p in time O(nlog3n) and space O(log(n)).
For convenience, we call the k-th bit of a number the k-th bit of its fractional part.
We make use of the notation
|
The N+n-th bit of a real number a is obtained by computing the n-th bit of the fractional part of 2Na.Proof : This property is obvious since the binary expansion of a
|
|
|
Thus from a sufficiently accurate value of {2N a}, this result permits to know the first term of its binary expansion. We illustrate this simple idea on p. One has
|
|
We are thus lead to compute a sufficiently accurate value of the fractional part {2Np} of 2N p. Formula (1) entails that {2Np} is the fractional part of the series
|
|
|
|
The computation of 2n modulo m is obtained thanks to a binary powering method, using only O(log(n)) operations modulo m. It consists in writing n in base 2
|
|
Finally, formula (1) permits to obtain bits of p between position N+1 and N+K by computing the K first bits of the number
|
|
|
Other formulas of this type have been found for p. The fastest known is due to Fabrice Bellard and is approximately 43% faster than (1) to compute the n-th bit of p :
|
2 Computing of the n-th bit of other constants |
In [1], David Bailey, Peter Borwein and Simon Plouffe give similar kind of series to compute directy the n-th bit of the constants log(2), p2 and log2(2). Later [2], D. J. Broadhurst found other similar series for a wider class of constants, like z(3), z(5), the Catalan constant C, p3, p4, log3(2), log4(2), log5(2). More recently, Gery Huvent [4] found other formula of this kind.
Here are a sample of those formulas :
|
|
|
As for the Catalan constant G = åk ³ 0 (-1)k/(2k+1)2, we have the concise form, shown in [4]
|
The direct computation of the n-th bit of the classical constants e, Ö2 and g with a similar complexity remains an open problem.
3 Computing the n-th decimal digit |
Unhappily, no formula of the same kind to compute decimal digits of p have been found yet (or more generally, in any base that is not a power of two). In other words, no series of the form
|
The first progress in this direction is due to Simon Plouffe in 1997, who found an algorithm with time O(n3 log3n) and very little memory to compute the n-th digit of p (see the Plouffe page for a description of the method). Unhappily, the complexity is so high that his technique does not reasonably permit to reach decimal digits at position n higher than several thousands.
The Simon Plouffe method has been improved by Fabrice Bellard with a O(n2) algorithm (complexity in terms of elementary operations on numbers of size O(logn)). The Fabrice Bellard page describes the corresponding algorithm, and a C program is also available. The millionth decimal digit can be reached with that technique. Nethertheless, even using parallelization, the complexity remains too high to give a practical alternative to techniques in quasi-linear time which computes all the decimal digits.
Recently, Xavier Gourdon found a better algorithm for the n-th decimal digit computation of p, that runs in time O(n2loglogn/log2 n) in terms of elementary operations. The improvement is nearly a factor log2 n compared to Fabrice Bellard algorithm. The corresponding algorithm is described in an unpublished paper [3] that can be dowloaded here in pdf format.
Timings comparisons have been made between this program and the program of Fabrice Bellard.
Index of | pidec time | F. Bellard |
decimal digit | program time | |
5,000 | 0.96 sec | 4.85 sec |
10,000 | 3.13 sec | 18.10 sec |
20,000 | 10.34 sec | 68.29 sec |
40,000 | 35.96 sec | 259.8 sec |
100,000 | 185.1 sec | 1520 sec |
200,000 | 628.7 sec | 5703 sec |
500,000 | 3525 sec | 34730 sec |
1,000,000 | 15869 sec | not ran |
2,000,000 | 42316 sec | not ran |
4,000,000 | 168191 sec | not ran |
| (2) |
|
|
The problem is thus essentially reduced to computing powers and sum of binomial coefficients modulo small numbers. Details are fulfilled in [3].
On the one hand, we have an algorithm to compute directly the n-th decimal digit of p in nearly quadratic time and using very little memory; on the other hand, computing all the first n digits of p is possible in quasi-linear time and with memory O(n). It is natural to ask if it were possible to find intermediate algorithms, that use a limited amount of memory O(m) (for example memory O(Ön)) and obtain an intermediate cost between linear and quadratic. In fact, the question has a positive answer by the use of formula () together with the so-called binary splitting technique on numbers of size O(m). More precisely, it is possible to obtain an algorithm which uses O(m) memory with running time
|
Details of this algorithm will added soon in [3].
Colin Percival PiHex project : this project currently holds the world record for the bit computation of p with the forty trillionth bit (10 trillionth hexit) of pi using Bellard's formula.
Fabrice Bellard page on the n-th digit computation of p.
Simon Plouffe page on the n-th decimal digit computation of p.
numbers, constants and computation