Если Вы найдете статью полезной и интересной  не сочтите за труд,
переведите материал или хотя бы часть его и отправьте на адрес
algolist@manual.ru.
Задать вопросы или просто написать письмо можно
также на странице контактов.
Is it possible to compute directly the nth 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 nth
bit of p with very little memory and in quasilinear 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
ц ч
ш
116^{k}
.
(1)
and use that to compute the nth bit of p without computing
all its first n bits. More precisely, their method permits to obtain
the nth bit of p in time O(nlog^{3}n) and space O(log(n)).
For convenience, we call the kth bit of a number the kth bit of
its fractional part.
Computation of the Nth 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+nth bit of a real number a is obtained by computing
the nth bit of the fractional part of 2^{N}a.
Proof :
This property is obvious since the binary expansion of a
a =
е
k
a_{k}2^{k}
, (a_{k} О {0,1})
entails
2^{N} a =
е
k
a_{k}2^{kN}
= A+B,
where
A =
е
k Ј N
a_{k} 2^{Nk} and B =
е
k > N
a_{k}2^{kN}
=
е
k > 0
a_{N+k}2^{k}
.
Thus A is an integer and we have 0 Ј B < 1.
In other words, B is the fractional part of 2^{N} a
(which we denote by {2^{N}a}), and the kth bit of B is the
N+kth bit of a.
·
Thus from a sufficiently accurate value of {2^{N} a},
this result permits to know the first term of its
binary expansion.
We illustrate this simple idea on p. One has
2^{12} p = 12867.9635ј,
so that {2^{12} p} = 0.9635ј. The 4digits 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 13th bit of p is 1, the 14th bit of p is 1,
the 17th bit of p is 0.
Powering
We are thus lead to compute a sufficiently accurate value of the
fractional part {2^{N}p} of 2^{N} p.
Formula (1) entails that {2^{N}p} is the
fractional part of the series
p_{N} =
Ґ е
n = 0
ж з
и
м н
о
4 ·2^{N4n}8n+1
ь э
ю

м н
о
2·2^{N4n}8n+4
ь э
ю

м н
о
2^{N4n}8n+5
ь э
ю

м н
о
2^{N4n}8n+6
ь э
ю
ц ч
ш
= A_{N} + B_{N},
where A_{N} and B_{N} are the summations
A_{N} =
е
0 Ј n Ј N/4
, B_{N} =
е
n > N/4
.
The first few digits of B_{N} are obtained easily. We thus
concentrate on A_{N}.
The general term of the series in A_{N} has the form
м н
о
p·2^{n}m
ь э
ю
where p, n and m are non negative integers. If p·2^{n} = qm+r is the
euclidian division
of p·2^{n} by m, we have {p·2^{n}/m} = r/m. In other words
м н
о
p·2^{n}m
ь э
ю
=
p·2^{n} mod mm
.
The computation of 2^{n} 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
n_{k} 2^{k},
so that
2^{n} mod m =
Х
0 Ј k Ј K, n_{k} = 1
2^{2k} mod m.
The values P_{k} = (2^{2k} mod m) are computed by successive squaring
modulo m.
Since K @ log_{2}(n), only O(log(n)) operations modulo m are
finally needed to compute 2^{n} 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
p_{N} = A_{N}+B_{N}
where
A_{N} =
е
0 Ј n Ј N/4
ж з
и
4·2^{N4n} mod 8n+18n+1

2·2^{N4n} mod 8n+48n+4

2^{N4n} mod 8n+58n+5

2^{N4n} mod 8n+68n+6
ц ч
ш
and
B_{N} =
е
N/4 < n Ј 2+(N+K)/4
ж з
и
48n+1

28n+4

18n+5

18n+6
ц ч
ш
12^{4nN}
.
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 nth bit of p :
In [1], David Bailey, Peter Borwein and Simon Plouffe give similar
kind of series to compute directy the nth bit of the constants
log(2), p^{2} and log^{2}(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, p^{3}, p^{4}, log^{3}(2), log^{4}(2),
log^{5}(2).
Here are a sample of those formulas :
log(2) =
е
k > 0
1k2^{k}
p^{2} =
е
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}
ц ч
ш
116^{k}
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}
ц ч
ш
116^{k}
+
4
е
k і 0
ж з
и
12^{3}(8k+1)^{3}
+
12^{4}(8k+2)^{3}

12^{6}(8k+3)^{3}

22^{7}(8k+4)^{3}

12^{9}(8k+5)^{3}
+
12^{10}(8k+6)^{3}
+
12^{12}(8k+7)^{3}
ц ч
ш
14096^{k}
.
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}
ц ч
ш
116^{k}
+
2
е
k і 0
ж з
и
12^{3}(8k+1)^{2} </td>
+</td>
12^{4}(8k+2)^{2}
+
12^{6}(8k+3)^{2}

12^{9}(8k+5)^{2}

12^{10}(8k+6)^{2}

12^{12}(8k+7)^{2}
ц ч
ш
14096^{k}
The direct computation of the nth 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)
110^{n}
, P, Q polynomials
is known for p. No such series either is known for other classical
constant like log(2), p^{2}, ј.
The first progress in this direction is due to Simon Plouffe in 1997, who
found an algorithm with time O(n^{3} log^{3}n) and very little memory to
compute the nth 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(n^{2}) 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
quasilinear time which computes all the decimal digits.
Related links on the nth digit computation
The miraculous BaileyBorweinPlouffe Pi algorithm :
the Steven Finch page on the nth 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 nth digit computation of p.
Simon Plouffe
page on the nth 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. 903913