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(n2) 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.

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

x1 = 1, x2 = 2x1, x3 = 3x2, ¼, xn = nxn-1.

At iteration number k, the value xk contains O(klog(k)) digits, thus
the computation of xk+1 = kxk has cost O(klog(k)).
Finally, the total cost with this basic approach is
O(2log(2)+¼+n log(n)) = O(n2log(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

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.

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

Q(a,b) = (a+1)(a+2)¼b,

P(a,b) = b(b-1)¼(a+2) + b(b-1)¼(a+3) +¼+ (b-1)b + b + 1.

P(a,b) and Q(a,b) are integers which satisfy

P(a,b)Q(a,b)

=

1a+1

+

1(a+1)(a+2)

+ ¼+

1(a+1)(a+2)¼b

.

In particular, 1+P(0,K)/Q(0,K) are the first K terms of the series (*).

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,

m =

éê
ë

a+b2

ùú
û

(integer part).

Then we have

P(a,b) = P(a,m)Q(m,b)+P(m,b), Q(a,b) = Q(a,m)Q(m,b).

These operations are performed recursively.

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.

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, see The constant e).

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 10n < q2K+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.

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