N-th digit computation

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

p = ¥
å
k = 0 
æ
ç
è
4
8k+1
- 2
8k+4
- 1
8k+5
- 1
8k+6
ö
÷
ø
1
16k
.
(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)).

1.1  Principle of the algorithm

For convenience, we call the k-th bit of a number the k-th bit of its fractional part.

Computation of the N-th bit of a number

We make use of the notation

{x} = x-[x]        (fractional part of a real number x).
The basic idea depends on the following easy result :

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
a =
å
k 
ak
2k
,       (ak Î {0,1})
entails
2N a =
å
k 
ak
2k-N
= A+B,
where
A =
å
k £ N 
ak 2N-k       and      B =
å
k > N 
ak
2k-N
=
å
k > 0 
aN+k
2k
.
Thus A is an integer and we have 0 £ B < 1. In other words, B is the fractional part of 2N a (which we denote by {2Na}), and the k-th bit of B is the N+k-th bit 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

212 p = 12867.9635¼,
so that {212 p} = 0.9635¼. The 4-digits precision of this fractional part permits to have the first term of its binary expansion
0.9635¼ = 1
2
+ 1
4
+ 1
8
+ 1
16
+ 0
32
+ 1
64
+ 1
128
+ ¼ = [0.1111011¼]base 2
Thus the 13-th bit of p is 1, the 14-th bit of p is 1, the 17-th bit of p is 0.

Powering

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

pN = ¥
å
n = 0 
æ
ç
è
ì
í
î
4 ·2N-4n
8n+1
ü
ý
þ
- ì
í
î
2·2N-4n
8n+4
ü
ý
þ
- ì
í
î
2N-4n
8n+5
ü
ý
þ
- ì
í
î
2N-4n
8n+6
ü
ý
þ
ö
÷
ø
= AN + BN,
where AN and BN are the summations
AN =
å
0 £ n £ N/4 
,      BN =
å
n > N/4 
.
The first few digits of BN are obtained easily. We thus concentrate on AN. The general term of the series in AN has the form
ì
í
î
p·2n
m
ü
ý
þ
where p, n and m are non negative integers. If p·2n = qm+r is the euclidian division of p·2n by m, we have {p·2n/m} = r/m. In other words
ì
í
î
p·2n
m
ü
ý
þ
= p·2n mod m
m
.

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

n = K
å
k = 0 
nk 2k,
so that
2n mod m =
Õ
0 £ k £ K, nk = 1 
22k mod m.
The values Pk = (22k mod m) are computed by successive squaring modulo m. Since K @ log2(n), only O(log(n)) operations modulo m are finally needed to compute 2n mod m.

Final result

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

pN = AN+BN
where
AN =
å
0 £ n £ N/4 
æ
ç
è
4·2N-4n mod 8n+1
8n+1
- 2·2N-4n mod 8n+4
8n+4
- 2N-4n mod 8n+5
8n+5
- 2N-4n mod 8n+6
8n+6
ö
÷
ø
and
BN =
å
N/4 < n £ 2+(N+K)/4 
æ
ç
è
4
8n+1
- 2
8n+4
- 1
8n+5
- 1
8n+6
ö
÷
ø
1
24n-N
.
The computation must be carried with an absolute precision of at least 2-K. Using a binary powering method, this formula easily permits to obtain the given bits with O(Nlog(N)) operations on numbers of size log(N).

1.2  Other BBP formulas for p

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 :

p = 1
26
¥
å
n = 0 
(-1)n
210n
æ
ç
è
- 25
4n+1
- 1
4n+3
+ 28
10n+1
- 26
10n+3
- 22
10n+5
- 22
10n+7
+ 1
10n+9
ö
÷
ø
(see here for a proof).

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 :

log(2) =
å
k > 0 
1
k2k

p2 =
å
k ³ 0 
æ
ç
è
16
(8k+1)2
- 16
(8k+2)2
- 8
(8k+3)2
- 16
(8k+4)2
- 4
(8k+5)2
- 4
(8k+6)2
+ 2
(8k+7)2
ö
÷
ø
1
16k

7
8
z(3)
=
6
å
k ³ 0 
æ
ç
è
1
2(8k+1)3
- 7
2(8k+2)3
- 1
4(8k+3)3
+ 10
4(8k+4)3
- 1
8(8k+5)3
- 7
8(8k+6)3
+ 1
16(8k+7)3
ö
÷
ø
1
16k
+
4
å
k ³ 0 
æ
ç
è
1
23(8k+1)3
+ 1
24(8k+2)3
- 1
26(8k+3)3
- 2
27(8k+4)3
- 1
29(8k+5)3
+ 1
210(8k+6)3
+ 1
212(8k+7)3
ö
÷
ø
1
4096k
.

As for the Catalan constant G = åk ³ 0 (-1)k/(2k+1)2, we have the concise form, shown in [4]

G
=
3
4

å
k ³ 0 
æ
ç
è
2
(4k+1)2
- 2
(4k+2)2
+ 1
(4k+3)2
ö
÷
ø
(-1)k
4k
-
1
32

å
k ³ 0 
æ
ç
è
8
(4k+1)2
4
(4k+2)2
+ 1
(4k+3)2
ö
÷
ø
(-1)k
64k
.

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


å
n 
P(n)
Q(n)
  1
10n
,       P, Q polynomials
is known for p. No such series either is known for other classical constant like log(2), p2, ¼.

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.

A new algorithm for the n-th decimal digit computation

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.

The program pidec

An implementation of the algorithm has been made with the program pidec. The source code and the windows executable are available for download. Notice that this implementation has not been fully optimized for clarity.

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

Figure 1: Comparison of timings on the n-th decimal digit computation, between the pidec program by Xavier Gourdon and Fabrice Bellard program, on a Pentium III 900 Mhz. (Fabrice Bellard program has been compiled using the long long option, which is faster).

Principle of the algorithm

The algorithm relies on a different technique compared to the one of Simon Plouffe and Fabrice Bellard. Our starting point was the arctan formula
p
4
= arctan(1) = ¥
å
k = 0 
(-1)k
2k+1
.
(2)
In this form, the formula is well suited to n-th decimal digit computation, but its convergence is too slow. This problem is overcomed with the use of a general alternating series acceleration technique, from Cohen, Villegas and Zagier, described in the page Acceleration of the convergence of series. Choosing families of polynomials of the form Pm(x) = (xM (1-x))N permits to obtain the following result. Setting
M = 2 é
ê
ê
n
log3 n
ù
ú
ú
    and    N = é
ê
ê
(n+n0+1) log(10)
log(2eM)
ù
ú
ú
.
the fractional part of 10n p is approximated with an error < 10-n0 by the fractional part of B-C, where
B
=
(M+1)N-1
å
k = 0 
(-1)k   4 ×10n mod (2k+1)
2k+1
,
(3)
C
=
N-1
å
k = 0 
(-1)k  5N-2 10n-N+2  sk mod (2MN+2k+1)
2MN+2k+1
,        sk = k
å
j = 0 
æ
ç
è
N
j
ö
÷
ø
(4)

The problem is thus essentially reduced to computing powers and sum of binomial coefficients modulo small numbers. Details are fulfilled in [3].

Generalization for intermediate memory size

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

O æ
ç
è
n2log3 n loglogn
m log2(n/m)
ö
÷
ø
.
(the value of m should be not too close to 0). For example, it is possible to obtain the n-th decimal digit of p in time O(n3/2logn loglogn) and memory O(n1/2). This technique uses small memory and it is likely that in the followings years, distributed computations on home computers with algorithms of this kind will be used to increase the reachable decimal digit of p with home computers, even if no quasi-linear complexity technique is found. Thousands of home computers could go higher than super computers ?

Details of this algorithm will added soon in [3].


Related links on the n-th digit computation


References

[1]
D.H. Bailey, P.B. Borwein and S. Plouffe, On the Rapid Computation of Various Polylogarithmic Constants, Mathematics of Computation, (1997), vol. 66, p. 903-913

[2]
D. J. Broadhurst, Polylogarithmic ladders, hypergeometric series and the ten millionth digits of z(3) and z(5), preprint (1998).

[3]
X. Gourdon, Computation of the n-th decimal digit of p with low memory, preprint (2003) nthdecimaldigit.pdf

[4]
G. Huvent, Formules BBP, Preprint (2001)


Back to

numbers, constants and computation


File translated from TEX by TTH, version 2.32.
On 12 Feb 2003, 02:08.