:: алгоритмы  и методы :: :: олимпиадные задачи :: :: связь :: :: о сайте ::
Путь: Математика » Быстрые вычисления » Вычисление N-го знака
  Вычисление N-го знака



Если Вы найдете статью полезной и интересной - не сочтите за труд, переведите материал или хотя бы часть его и отправьте на адрес algolist@manual.ru. Задать вопросы или просто написать письмо можно также на странице контактов.

Статья любезно предоставлена:
© Xavier Gordon
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 [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
ц
ч
ш

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

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

G
=
3
е
k і 0 
ж
з
и
1
2(8k+1)2
- 1
2(8k+2)2
+ 1
4(8k+3)2
<
/td>
- 1
8(8k+5)2
+ 1
8(8k+6)2
<
/td>
- 1
16(8k+7)2
ц
ч
ш
1
16k
+
2
е
k і 0 
ж
з
и
1
23(8k+1)2
<
/td>
+</td> 1
24(8k+2)2
+ 1
26(8k+3)2
- 1
29(8k+5)2
- 1
210(8k+6)2
- 1
212(8k+7)2
ц
ч
ш
1
4096k

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



[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), preprin998).