         ╧ґҐ№: ╠рҐхьрҐшър » ┴√ёҐЁ√х т√ішёыхэш  » ┬√ішёыхэшх N-ую чэрър
═р яЁртрє Ёхъырь√
ёюфхЁцрэъш ш ёяюэёюЁ√.
┬√ішёыхэшх N-ую чэрър  ┼ёыш ┬√ эрщфхҐх ёҐрҐ№■ яюыхчэющ ш шэҐхЁхёэющ - эх ёюіҐшҐх чр ҐЁґф, яхЁхтхфшҐх ьрҐхЁшры шыш єюҐ  с√ ірёҐ№ хую ш юҐяЁрт№Ґх эр рфЁхё algolist@manual.ru. ╟рфрҐ№ тюяЁюё√ шыш яЁюёҐю эряшёрҐ№ яшё№ью ьюцэю Ґръцх эр ёҐЁрэшІх ъюэҐръҐют.

╤ҐрҐ№  ы■схчэю яЁхфюёҐртыхэр:
web: numbers.computation.free.fr

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 , 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 , David Bailey, Peter Borwein and Simon Plouffe give the following formula

 p = е х k = 0 цч ш 48k+1 - 28k+4 - 18k+5 - 18k+6 Іі Ї 116k .
(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 ak2k ,       (ak ╬ {0,1})
entails
 2N a = х k ak2k-N = A+B,
where
 A = х k г N ak 2N-k       and      B = х k > N ak2k-N = х k > 0 aN+k2k .
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╝ = 12 + 14 + 18 + 116 + 032 + 164 + 1128 + ╝ = [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-4n8n+1 №¤ ■ - ьэ ю 2·2N-4n8n+4 №¤ ■ - ьэ ю 2N-4n8n+5 №¤ ■ - ьэ ю 2N-4n8n+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·2nm №¤ ■
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·2nm №¤ ■ = p·2n mod mm .

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+18n+1 - 2·2N-4n mod 8n+48n+4 - 2N-4n mod 8n+58n+5 - 2N-4n mod 8n+68n+6 Іі Ї
and
 BN = х N/4 < n г 2+(N+K)/4 цч ш 48n+1 - 28n+4 - 18n+5 - 18n+6 Іі Ї 124n-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 = 126 е х n = 0 (-1)n210n цч ш - 254n+1 - 14n+3 + 2810n+1 - 2610n+3 - 2210n+5 - 2210n+7 + 110n+9 Іі Ї

2  Computing of the n-th bit of other constants  In , 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 , 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).

Here are a sample of those formulas :

 log(2) = х k > 0 1k2k

 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 Іі Ї 116k

 78 z(3)
 =
 6 х k │ 0 цч ш 12(8k+1)3 - 72(8k+2)3 - 14(8k+3)3 + 104(8k+4)3 - 18(8k+5)3 - 78(8k+6)3 + 116(8k+7)3 Іі Ї 116k
 +
 4 х k │ 0 цч ш 123(8k+1)3 + 124(8k+2)3 - 126(8k+3)3 - 227(8k+4)3 - 129(8k+5)3 + 1210(8k+6)3 + 1212(8k+7)3 Іі Ї 14096k .

As for the Catalan constant G = хk 0 (-1)k/(2k+1)2, we have

 G
 =
 3 х k │ 0 цч ш 12(8k+1)2 - 12(8k+2)2 + 14(8k+3)2 - 18(8k+5)2 + 18(8k+6)2 - 116(8k+7)2 Іі Ї 116k
 +
 2 х k │ 0 цч ш 123(8k+1)2 + 124(8k+2)2 + 126(8k+3)2 - 129(8k+5)2 - 1210(8k+6)2 - 1212(8k+7)2 Іі Ї 14096k

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) 110n ,       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. 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.

Related links on the n-th digit computation

• The miraculous Bailey-Borwein-Plouffe Pi algorithm : the Steven Finch page on the n-th digit computation of p.
• 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.

References  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


D. J. Broadhurst, Polylogarithmic ladders, hypergeometric series and the ten millionth digits of z(3) and z(5), preprin998).