Journal of Symbolic Computation 26 (1998), 339–341.

COMPOSING POWER SERIES OVER A FINITE RING IN ESSENTIALLY LINEAR TIME

DANIEL J. BERNSTEIN

Abstract. Fix a finite commutative ring R. Let u and v be power series over R, with v(0) = 0. This paper presents an algorithm that computes the first n terms of the composition u(v), given the first n terms of u and v, in n^{1+o(1)} ring operations. The algorithm is very fast in practice when R has small characteristic.

1. Introduction

Let f be a polynomial over a commutative ring R, and let g be an element of R[x]/x^{n}. The point of this paper is a simple new algorithm to compute f(g) when R has finite characteristic.

For prime characteristic, see section 2. For prime-power characteristic more generally, see section 3. For other characteristics, use the Chinese Remainder Theorem, as in [5, equation 4.3.2–9].

Applications. The problem of computing f(g), under the restrictions deg f < n and g(0) = 0, is known as order-n power series composition. Here’s the point: given power series u and v over R, with v(0) = 0, define f = u(z) mod z^{n} and g = v(x) mod x^{n}; then f(g) = u(v(x)) mod x^{n}.

Power series composition is the bottleneck in reversion and iteration of power series. See [2] and [5, section 4.7].

Previous work. Brent and Kung describe two power series composition algorithms in [2]. The first algorithm computes f(g) in n^{α} ring operations for some α > 1.5, depending on the speed of matrix multiplication. This algorithm can be applied to the more general problem of modular composition; see [4]. The second algorithm computes f(g) in about n^{1.5} ring operations, provided that g′ is invertible and that all primes up to about \sqrt{n} are cancellable in R.

My algorithm takes n^{1+o(1)} ring operations if R is fixed. It is the method of choice for power series composition over rings whose characteristic is a product of small primes; in particular, fields of small prime characteristic. The second Brent-Kung algorithm remains the fastest method for fields of large prime characteristic.

2. Prime characteristic

Fix a ring R of prime characteristic p. I will reduce the problem of computing f(g), with deg f < d and g ∈ R[x]/x^{n}, to a sequence of p subproblems with d and n both reduced by a factor of p. Write m = ⌈n/p⌉.

Date: 19971007.

1991 Mathematics Subject Classification. Primary 13P99.

The author was supported by the National Science Foundation under grant DMS–9600083.

Notice that g^{p} is a polynomial in x^{p}: there is a polynomial h ∈ R[y]/y^{m} with g^{p} = h(x^{p}). The coefficient of y^{j} in h is the pth power of the coefficient of x^{j} in g.

Split f as f(z) = f_{0}(z^{p}) + zf_{1}(z^{p}) + ⋯ + z^{p−1}f_{p−1}(z^{p}). Compute each f_{j}(h) recursively by the same procedure; substitute y ↦ x^{p} into f_{j}(h) to obtain f_{j}(g^{p}); finally apply Horner’s rule to evaluate f(g) = f_{0}(g^{p}) + ⋯ + g^{p−1}f_{p−1}(g^{p}).

The recursion stops when d is sufficiently small. For example, f(g) is simply f(0) when d = 1.

3. Prime-power characteristic

Fix a ring R of characteristic p^{k}, with p prime and k ≥ 1.

Write A = R[x]/x^{n}. Also write m = ⌈n/p⌉ and B = R[y]/y^{m}. Embed B into A by y ↦ x^{p}; this embedding, which amounts to some copying inside the computer, is not stated explicitly in the following algorithm.

Algorithm C. Given f ∈ R[z] and g ∈ A, to compute f(g + ϵ) in the ring A[ϵ]/(ϵ^{k}, pϵ^{k−1}, …, p^{k−1}ϵ):

1. If deg f < 1: Print f(0) and stop.

2. Find h ∈ B with g^{p} − h ∈ pA. Set β ← (g + ϵ)^{p} − h.

3. Set j ← p − 1 and s ← 0.

4. Compute f_{j}(h + δ) in B[δ]/(δ^{k}, pδ^{k−1}, …, p^{k−1}δ) by Algorithm C recursively, where f(z) = f_{0}(z^{p}) + zf_{1}(z^{p}) + ⋯ + z^{p−1}f_{p−1}(z^{p}).

5. Set s ← (g + ϵ)s + ∑ b_{i}β^{i}, where f_{j}(h + δ) = ∑ b_{i}δ^{i}.

6. If j = 0: Print s and stop.

7. Decrease j by 1 and return to step 4.

The idea of Algorithm C is as follows. Consider f(g + pt) in the polynomial ring A[t]. It equals f_{0}((g+pt)^{p})+⋯+(g+pt)^{p−1}f_{p−1}((g+pt)^{p}). To compute f_{j}((g+pt)^{p}), find h ∈ B with g^{p} − h ∈ pA, and find v ∈ A[t] satisfying h + pv = (g + pt)^{p}; recursively compute f_{j}(h + pu) in the polynomial ring B[u]; then substitute u ↦ v to obtain f_{j}((g + pt)^{p}).

To avoid multiplications and divisions by p, Algorithm C works with polynomials in ϵ = pt and δ = pu. Thus f(g+pt) becomes f(g+ϵ), f_{j}(h+pu) becomes f_{j}(h+δ), and u ↦ v becomes δ ↦ (g + ϵ)^{p} − h.

Algorithm C computes f(g + ϵ), not merely f(g). One can save some time by eliminating ϵ at the top level of the recursion, if the goal is to compute f(g). The recursive call in step 4 still needs f_{j}(h + δ), not merely f_{j}(h).

Details and improvements. I represent the ring A[ϵ]/(ϵ^{k}, pϵ^{k−1}, …, p^{k−1}ϵ) by A[ϵ]/ϵ^{k}. To multiply in A[ϵ]/ϵ^{k}, I do k(k + 1)/2 multiplications in A. For large k there are faster algorithms; see, e.g., [3].

To compute (g + ϵ)^{p} in step 2, I perform p − 1 multiplications by g + ϵ in A[ϵ]/ϵ^{k}. Each multiplication by g + ϵ is implemented with k multiplications in A. There are many faster algorithms; see, e.g., [5, section 4.6.3].

I choose h in step 2 as ∑ h_{i}y^{i} where h_{i} is the coefficient of x^{pi} in g^{p}. Then h_{i} can be extracted from (g + ϵ)^{p} with no extra arithmetic.

The sum ∑ b_{i}β^{i} = f_{j}(h + β) in step 5 can be viewed as an order-k composition over A. Horner’s rule, which uses k − 1 multiplications in A[ϵ]/ϵ^{k}, suffices for small k. The first Brent-Kung algorithm is better for large k.

COMPOSING POWER SERIES OVER A FINITE RING IN ESSENTIALLY LINEAR TIME 3

One can skip the multiplication of g+ϵ by s in step 5 when j = p−1, since s = 0. One can speed up some of the remaining multiplications by taking advantage of the sparseness of B inside A.

Speed. Let μ be a nondecreasing function such that elements of R[x]/x^{n} can be multiplied in time nμ(n). I usually assume fast multiplication, meaning that μ(n) ∈ n^{o(1)} for n → ∞. See [3] for a fast multiplication method. See [1] for a survey of multiplication methods.

Algorithm C’s run time is dominated by multiplications. Its multiplication time is at most

c (e(n − 1) + \frac{p^{e} − 1}{p − 1}) μ(n)

if deg f < p^{e}, where c = (2p − 1)k + p(k − 1)k(k + 1)/2. (Here c is the number of multiplications in A performed in steps 2 and 5 of Algorithm C. It can be reduced in several ways, as discussed above.) Indeed, for deg f < 1, Algorithm C uses no multiplications. Otherwise it performs p recursive calls, taking time at most

pc ((e − 1)(m − 1) + \frac{p^{e−1} − 1}{p − 1}) μ(m) ≤ c ((e − 1)(n − 1) + \frac{p^{e} − p}{p − 1}) μ(n)

by induction, and c multiplications in A, taking time at most cnμ(n).

In particular, order-n power series composition takes time O(nμ(n) log n) for fixed characteristic.

References

[1] Daniel J. Bernstein, Multidigit multiplication for mathematicians, preprint available from

   http://pobox.com/~djb/papers/m3.dvi.

[2] Richard P. Brent, H. T. Kung, Fast algorithms for manipulating formal power series, Journal of the ACM 25 (1978), 581–595.

[3] David G. Cantor, Erich Kaltofen, On fast multiplication of polynomials over arbitrary algebras, Acta Informatica 28 (1991), 693–701.

[4] Erich Kaltofen, Victor Shoup, Subquadratic-time factoring of polynomials over finite fields, preprint.

[5] Donald E. Knuth, The art of computer programming, volume 2: seminumerical algorithms, 2nd edition, Addison-Wesley, Reading, Massachusetts, 1981.

Department of Mathematics, Statistics, and Computer Science, The University of Illinois at Chicago, Chicago, IL 60607–7045

E-mail address: djb@pobox.com