Если Вы найдете статью полезной и интересной - не сочтите за труд,
переведите материал или хотя бы часть его и отправьте на адрес
algolist@manual.ru.
Задать вопросы или просто написать письмо можно
также на странице контактов.
Автор Xavier Gordon
Most classical series formulaes to compute constants (like the
exponential series for e, arctan formulaes for p, Chudnovsky
or Ramanujan formulaes for p, ...) have a time cost of O(n^{2}) to
compute n digits of the series using a classical approach.
Using FFT based multiplication, the binary splitting approach
(foreshadowed in [2], [3] ; see
also [1][p. 329]) usually permits to evaluate such series in
time O(n log(n)^{3}) (or O(n log(n)^{2}) for e), which is much
better.
At iteration number k, the value x_{k} contains O(klog(k)) digits, thus
the computation of x_{k+1} = kx_{k} has cost O(klog(k)).
Finally, the total cost with this basic approach is
O(2log(2)+ј+n log(n)) = O(n^{2}log(n)).
A better approach is the binary splitting : it just consists in
recursively cutting the product of m consecutive integers in half.
It leads to better results when
products on large integers are performed with a fast method.
More precisely, the computation of p(a,b), where
p(a,b) є (a+1)(a+2) ј(b-1) b =
b!a!
,
is done by performing the product
p(a,b) = p
ж з
и
a,
a+b2
ц ч
ш
p
ж з
и
a+b2
,b
ц ч
ш
,
where the two terms in the product are computed recursively in the same
way. This permits to compute n! = p(0,n).
We now concentrate on the estimation of the timing cost C(a,b) of
this process.
We denote by M(d) the cost of multiplying two integers of size
d.
The terms p(a,(a+b)/2) and p((a+b)/2,b) have size
O((b-a)log(b)), thus
(the superlinearity of M(d) should be used : 2M(d/2) Ј M(d)).
Finally, we have proved that computing n! can be done in time
C(0,n) = O(log(n)M(n log(n)).
This is always better than the classical approach when products on large
integers are done with a fast method. When FFT is used for large
multiplication (see FFT based multiplication of large
numbers), it leads to the bound
The example of the factorial can be generalized for many series.
Instead of presenting a very general and complicated approach,
the next sections are dedicated to illustrate the techniques on
classical series.
The cost of the process to compute P(0,K) and Q(0,K)
is easily proved to be O(log(K)M(Klog(K))). A final division of
P(0,K) by Q(0,K) is then required to compute the first K terms of
the series (its cost is O(M(n)), see Inverse and
n-th roots of large numbers).
Since K @ n/log(n), we have finally obtained a process to
compute e in time O(log(n) M(n)). An FFT based multiplication
satisfies M(n) = O(n log(n)) and leads to a total cost of
O(n log(n)^{2}) to compute n digits of e.
An easy code illustrating the binary splitting method to compute e
can be found in here. It uses libraries described here.
Practical implementation
The computation of the exponential series thanks to the binary
splitting approach is quite easy. In the practice, one should stop the
splitting recursion when b-a becomes small (the threshold should be
an optimized value) to use a classical approach. An slight improvement
can also be obtained by using a classical multiplication when P(a,b)
and Q(a,b) have a small number of digits.
Thanks to the binary splitting and using the FFT code of
PiFast, the first million digits of e was
computed in less than 10 seconds (pentium II 350). The binary splitting
method to compute e is better than any other approaches (much better
than the AGM based approach.
It must be pointed out that binary splitting needs more memory than the
classical method (this is nearly always the case : fast methods are
memory consuming). Netherveless, memory is quite easy to manage with
this process.
Now, to compute n decimal digits of arctan(1/q), one must stop the
summation at k = K as soon as 10^{n} < q^{2K+1}, which gives
K @ n log(10)/2. The total cost of the process to compute the
first n digits of the series is
O(log(K)M(Klog(K))) = O(log(n)M(n log(n))).
With FFT based multiplication, this leads to the bound O(n log(n)^{3}).
Notice that the process is more complicated and expensive than for the
exponential series for several reasons :
An auxilliary large integer R(a,b) is needed.
More products should be performed at each step.
The large integers P(0,K) and Q(0,K) used to
compute the first K terms of the series have size O(n log(n)) which
is much bigger than n.
The last point is important. Nethertheless, the problem can be
overpassed by working with a maximal precision of n decimal digits
during the process (this requires a treatment for large integers with
limited precision). This solution reduces the memory to
O(n) instead of O(n log(n)) but only reduces the cost of
O(n log(n)^{3}) by a constant factor (but appreciable).
It must be pointed out that this crucial remark for practical
implementations does not often appear on the litterature concerning
binary splitting.
The technique to compute arctanh series is very similar and can be
efficiently used to compute digits of log(2).
Binary splitting generalizes to all hypergeometric series which
involves only integer coefficients (see [1][p. 335]).
A succint and recursive general
approach is also presented in [4] (note: the authors were
unfortunately unaware of the original work of E. A. Karatsuba on FEE,
just presented below).
The general algorithm also includes the evaluation of such series for non
rational parameters.
A method of called FEE (Fast E-function Evaluation) has been pioneered
by Ekatherine A. Karatsuba, and while reminiscent of splitting, FEE has
a unique and general power of its own. (for details and
examples of FEE see Borwein, Bradley, Crandall
"Computational Strategies for the Riemann Zeta Function").
The FEE method of E. A. Karatsuba can also be used to compute other
constants which do not write in a geometrically convergent
hypergeometric series : examples include z(n) for integer values of
n [6], Hurwitz Zeta function and Dirichlet
L-Series [7].
R. P. Brent, The Complexity of Multiple-Precision
Arithmetic, Complexity of Computational Problem Solving,
R. S. Andressen and R. P. Brent, Eds (Univ. of Queensland Press,
Brisbane, 1976)
E. A. Karatsuba,
Fast Calculation of the Riemann Zeta function z(s) for Integer
Values of the Argument s,
Problems of Information Transmission 31 (1995) 353-362.
E. A. Karatsuba,
Fast evaluation of hypergeometric function by FEE,
Computational Method and Function Theory 1997, Proceeding of the third
CMFT Conference (N. Papamichael, St. Ruscheweyh, E.B. Saff, eds)
in Approximation and Decomposition. World Scientific (1999),-314.