Если Вы найдете статью полезной и интересной - не сочтите за труд,
переведите материал или хотя бы часть его и отправьте на адрес
algolist@manual.ru.
Задать вопросы или просто написать письмо можно
также на странице контактов.
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.
In [1], 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)).
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).
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 :
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
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 </td>
-
18(8k+5)2
+
18(8k+6)2 </td>
-
116(8k+7)2
ц ч
ш
116k
+
2
е
k і 0
ж з
и
123(8k+1)2 </td>
+</td>
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.
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.
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
