Математика: Быстрое вычисление функций и констант.

N-th digit computation

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

ж
з
и

8k+1

8k+4

8k+5

8k+6

ц
ч
ш

16^k

(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(nlog³n) 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 2^Na.

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^k-N

= A+B,

where

A =

е
k Ј N

a_k 2^N-k and B =

е
k > N

a_k

2^k-N

е
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^Na}), and the k-th bit of B is the N+k-th 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¹² p = 12867.9635ј,

so that {2¹² p} = 0.9635ј. The 4-digits precision of this fractional part permits to have the first term of its binary expansion

0.9635ј =

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 {2^Np} of 2^N p. Formula (1) entails that {2^Np} is the fractional part of the series

p_N =

Ґ
е
n = 0

ж
з
и

м
н
о

4 ·2^N-4n

8n+1

ь
э
ю

м
н
о

2·2^N-4n

8n+4

ь
э
ю

м
н
о

2^N-4n

8n+5

ь
э
ю

м
н
о

2^N-4n

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ⁿ

ь
э
ю

where p, n and m are non negative integers. If p·2ⁿ = qm+r is the euclidian division of p·2ⁿ by m, we have {p·2ⁿ/m} = r/m. In other words

м
н
о

p·2ⁿ

ь
э
ю

p·2ⁿ mod m

The computation of 2ⁿ 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ⁿ mod m =

Х
0 Ј k Ј K, n_k = 1

2^{2^k} mod m.

The values P_k = (2^{2^k} mod m) are computed by successive squaring modulo m. Since K @ log₂(n), only O(log(n)) operations modulo m are finally needed to compute 2ⁿ 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^N-4n mod 8n+1

8n+1

2·2^N-4n mod 8n+4

8n+4

2^N-4n mod 8n+5

8n+5

2^N-4n mod 8n+6

8n+6

ц
ч
ш

and

B_N =

е
N/4 < n Ј 2+(N+K)/4

ж
з
и

8n+1

8n+4

8n+5

8n+6

ц
ч
ш

2^4n-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 =

2⁶

Ґ
е
n = 0

(-1)ⁿ

2¹⁰ⁿ

ж
з
и

2⁵

4n+1

4n+3

2⁸

10n+1

2⁶

10n+3

2²

10n+5

2²

10n+7

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), p² and log²(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³, p⁴, log³(2), log⁴(2), log⁵(2).

Here are a sample of those formulas :

log(2) =

е
k > 0

k2^k

p² =

е
k і 0

ж
з
и

(8k+1)²

(8k+2)²

(8k+3)²

(8k+4)²

(8k+5)²

(8k+6)²

(8k+7)²

ц
ч
ш

16^k

z(3)

е
k і 0

ж
з
и

2(8k+1)³

2(8k+2)³

4(8k+3)³

4(8k+4)³

8(8k+5)³

8(8k+6)³

16(8k+7)³

ц
ч
ш

16^k

е
k і 0

ж
з
и

2³(8k+1)³

2⁴(8k+2)³

2⁶(8k+3)³

2⁷(8k+4)³

2⁹(8k+5)³

2¹⁰(8k+6)³

2¹²(8k+7)³

ц
ч
ш

4096^k

As for the Catalan constant G = е_{k і 0} (-1)^k/(2k+1)², we have

е
k і 0

ж
з
и

2(8k+1)²

2(8k+2)²

4(8k+3)²

8(8k+5)²

8(8k+6)²

16(8k+7)²

ц
ч
ш

16^k

е
k і 0

ж
з
и

2³(8k+1)²

2⁴(8k+2)²

2⁶(8k+3)²

2⁹(8k+5)²

2¹⁰(8k+6)²

2¹²(8k+7)²

ц
ч
ш

4096^k

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)

10ⁿ

, P, Q polynomials

is known for p. No such series either is known for other classical constant like log(2), p², ј.

The first progress in this direction is due to Simon Plouffe in 1997, who found an algorithm with time O(n³ log³n) 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(n²) 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), preprint (1998).

Вверх по странице, к оглавлению и навигации.