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

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²) 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)³) (or O(n log(n)²) for e), which is much better.

1.1  The example of the factorial

Suppose you would like to compute the digits of the factorial n!. A basic approach would consist in computing successively the following values

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

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

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

1.2  Generalization

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.

We want to compute the first n decimal digits of the series

We should first stop the summation as soon as k! > 10ⁿ, which occurs for k = K @ d/log(d). We make use of the notations

To compute P(a,b) and Q(a,b) by the binary splitting method, we proceed as follows. We set m to the half of a+b,

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

We now present the binary splitting method to compute digits of the arctan series

Notice that

The splitting process consists in recursively performing the operations (m denotes the half of a+b)

Now, to compute n decimal digits of arctan(1/q), one must stop the summation at k = K as soon as 10ⁿ < 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

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

[1]

J. M. and P. B. Borwein Pi and the AGM John Wiley and sons, (1987)

[2]

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)

[3]

R. P. Brent, Fast multiple-precision evaluation of elementary functions, Journal of the ACM, 23 (1976) 242-251.

[4]

B. Haible and T. Papanikolaou, Fast multiprecision evaluation of series of rational numbers, report TI-97-7.binsplit, TH Darmstadt. Available here.

[5]

E. A. Karatsuba, Fast evaluation of transcendental functions, Problems of Information Transmission 27 (1991) 339-360.

[6]

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.

[7]

E. A. Karatsuba, Fast evaluation of the Hurwitz Zeta function and Dirichlet L-series, Problems of Information Transmission 34 (1998) 342-353.

[8]

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.