Power Series Composition in Near-Linear Time

In this paper, $\mathbb{A}$ denotes an effective ring, i.e., a commutative ring with unity, and one can perform the basic ring operations $(+,-,\times)$ with unit cost in $\mathbb{A}$. The polynomial ring and the formal power series ring over $\mathbb{A}$ are denoted by $\mathbb{A}[x]$ and $\mathbb{A}[[x]]$, respectively. We use $\mathbb{A}[x]_{<n}$ to denote the set of polynomials of degree less than $n$. For any polynomial or power series $f(x)=\sum_{i}f_{i}x^{i}$, $[x^{i}]f(x)$ denotes the coefficient $f_{i}$ of $x^{i}$, and $f(x)\bmod x^{n}$ denotes the truncation $\sum_{i=0}^{n-1}f_{i}x^{i}$.

For a bivariate polynomial $f(x,y)\in\mathbb{A}[x,y]$, $\deg_{x}f(x,y)$ and $\deg_{y}f(x,y)$ denote its degree with respect to $x$ and $y$, respectively. The bidegree of $f(x,y)$ is the pair $\mathop{\mathrm{bideg}}f(x,y)=(\deg_{x}f(x,y),\deg_{y}f(x,y))$. Inequalities between bidegrees are component-wise. For example, $\mathop{\mathrm{bideg}}f(x,y)\leq(n,m)$ means $\deg_{x}f(x,y)\leq n$ and $\deg_{y}f(x,y)\leq m$.

Informally, we measure the complexity of an algorithm by counting the number of arithmetic operations in the ring $\mathbb{A}$. The underlying computation model is the straight-line program. Readers are referred to the textbook of Bürgisser, Clausen and Shokrollahi [BCS10, Sec. 4.1] for a detailed explanation.

Our algorithm calls polynomial multiplication as a crucial subroutine. Over different rings, the time complexity of our algorithm may vary due to the different polynomial multiplication algorithms.

For an effective ring $\mathbb{A}$, we call a function $\mathop{\mathsf{M}}_{\mathbb{A}}(n)$ a multiplication time for $\mathbb{A}[x]$ if polynomials in $\mathbb{A}[x]$ of degree less than $n$ can be multiplied in $\mathop{\mathrm{O}}(\mathop{\mathsf{M}}_{\mathbb{A}}(n))$ operations. When the ring $\mathbb{A}$ is clear from the context, we simply write $\mathop{\mathsf{M}}(n)$.

Throughout this paper, we assume a regularity condition on $\mathop{\mathsf{M}}(n)$, that for some constants $\alpha,\beta\in(0,1)$, the function $\mathop{\mathsf{M}}$ satisfies the condition $\mathop{\mathsf{M}}(\lceil\alpha n\rceil)\leq\beta\mathop{\mathsf{M}}(n)$ for all sufficiently large $n$. All the examples of $\mathop{\mathsf{M}}$ presented in this paper satisfy this condition.

The fast Fourier transform (FFT) [CT65] gives the bound when $\mathbb{A}$ is a field possessing roots of unity of sufficiently high order; for example, the complex numbers $\mathbb{C}$ satisfies $\mathop{\mathsf{M}}_{\mathbb{C}}(n)=\mathop{\mathrm{O}}(n\log n)$. For any $\mathbb{F}_{p}$-algebra $\mathbb{A}$ where $p$ is a prime, Harvey and van der Hoeven [HH19] proved that $\mathop{\mathsf{M}}_{\mathbb{A}}(n)=\mathop{\mathrm{O}}(n\log n\,4^{\log^{*}n})$³³3$\log^{*}$ is the iterated logarithm.. For general rings, Cantor and Kaltofen [CK91] showed a uniform bound $\mathop{\mathsf{M}}_{\mathbb{A}}(n)=\mathop{\mathrm{O}}(n\log n\log\log n)$. This implies that our algorithm has a uniform bound $\mathop{\mathrm{O}}(n\log n\log m\log\log n+m\log m\log\log m)$.

We may also obtain complexity bounds for various computation models, depending on the fast arithmetic of $\mathbb{A}$ and the polynomial multiplication algorithm on the corresponding model. For example, for the bit complexity of boolean circuits and multitape Turing machines, since the polynomial multiplication over $\mathbb{F}_{p}$ can be done in $\mathop{\mathsf{M}}_{\mathbb{F}_{p}}(n)=\mathop{\mathrm{O}}((n\log p)\log(n\log p)4^{\max(\log^{*}p,\log^{*}n)})$ bit operations [HH19], substituting such a bound into our algorithm yields a bound for the bit complexity of the composition problem.

We need the following operations on power series and polynomials that are used in our algorithm.

Using Kronecker’s substitution, one can reduce the problem of computing the multiplication of bivariate polynomials to univariate cases.

We also require the following lemma, which bounds the total time complexity of polynomial multiplications with exponentially decreasing degrees.

The transposition principle, a.k.a. Tellegen’s principle, is an algorithmic theorem that allows one to produce a new algorithm of the transposed problem of a given linear problem.

Roughly speaking, the transposition principle allows one to compute the transposed problem with the same time complexity as the original problem, provided that the algorithm is $\mathbb{A}$-linear, in the sense that only linear operations in the coefficients of $b$ are performed.

The transposed problem of power series composition is the power projection problem, which is defined as follows.

For a given polynomial $g(x)$, consider the $n\times m$ matrix $M$ given by $M_{ij}=[x^{i}]g(x)^{j}$. The power series composition problem is equivalent to computing $Mf$, where $f$ is the column vector of coefficients of $f(x)$, and the power projection problem is equivalent to computing $M^{\mathsf{T}}b$, where $b$ is the column vector of coefficients of the linear form $w$. It is important to note that the power projection and composition problems are linear only when the polynomial $g$ is fixed, as the output of these problems is nonlinear with respect to the coefficients of $g$.

Our objective is to compute $f_{i}=w(g(x)^{i}\bmod x^{n})$ for all $i=0,\dots,m-1$. We define a polynomial $w^{\mathrm{R}}$ as $w^{\mathrm{R}}(x)=\sum_{j=0}^{n-1}w_{n-1-j}x^{j}$. Let $P(x,y)\coloneqq w^{\mathrm{R}}(x)$ and $Q(x,y)\coloneqq 1-yg(x)$. Then, the generating function $f(y)\coloneqq\sum_{i=0}^{m-1}f_{i}y^{i}$ is rearranged as:

If $n-1=0$, then $u=P(0,y)/Q(0,y)$. If $n-1\geq 1$, we multiply $Q(-x,y)$ by the numerator and the denominator, resulting in

Let $A(x,y)\coloneqq Q(x,y)Q(-x,y)\bmod x^{n}\bmod y^{m}$. $A(x,y)$ satisfies the equation $A(x,y)=A(-x,y)$, which implies $[x^{i}]A(x,y)=0$ for all odd $i$; thus, there exists a unique polynomial $V(x,y)$ such that $V(x^{2},y)=A(x,y)$. Now, we have

for $U(x,y)\coloneqq P(x,y)Q(-x,y)\bmod x^{n}\bmod y^{m}$. Since $1/V(x^{2},y)$ is an even formal power series with respect to $x$, i.e., $[x^{i}](1/V(x,y))=0$ for all odd $i$, we can ignore the odd (or even) part of $U(x,y)$ if $n-1$ is even (or odd). Let $U_{e}(x,y)$ and $U_{o}(x,y)$ be the unique polynomials such that $U(x,y)=U_{e}(x^{2},y)+xU_{o}(x^{2},y)$. Then, we have

Now the problem is reduced to a case where $n$ is halved. We repeat this reduction until $n-1$ becomes $0$. The resulting algorithm is summarized in Algorithm 1.

The algorithm presented in this section is essentially the transposition of the algorithm presented in Section 3.1. However, to enhance comprehension, we present a derivation that does not rely on the transposition principle. To this end, we rewrite the objective $f(g(x))\bmod x^{n}$.

Let $P(y)\coloneqq y^{m-1}f(y^{-1}),Q(x,y)\coloneqq 1-yg(x)$. Then, we observe that

Let us define a operator $\mathop{\mathcal{F}}\nolimits$ by

Then our objective is to compute $\mathop{\mathcal{F}}\nolimits_{m-1,m}(P(y)/Q(x,y))\bmod x^{n}$. For algorithmic convenience, we generalize our aim to computing $\mathop{\mathcal{F}}\nolimits_{d,m}(P(y)/Q(x,y))\bmod x^{n}$ for a given $d$. If $n=1$, then $\mathop{\mathcal{F}}\nolimits_{d,m}(P(y)/Q(x,y))\bmod x^{n}=\mathop{\mathcal{F}}\nolimits_{d,m}(P(y)/Q(0,y))$. If $n>1$, similar to the transformation described in Section 3.1, we can derive that

where $V(x^{2},y)=Q(x,y)Q(-x,y)\bmod x^{n}\bmod y^{m}$. Let $e\coloneqq\max\{0,d-\deg_{y}Q(x,y)\}$. Then, we can discard terms of order less than $e$ with respect to $y$ among $P(y)/V(x^{2},y)$ to compute $u$, because they has no contribution to term with order $d$ or higher with respect to $y$ when multiplied by $Q(-x,y)$. Therefore, $u$ can be rewritten as

Now the problem is reduced to a case where $n$ is halved. We repeat this reduction until $n$ becomes $1$. The resulting algorithm is summarized in Algorithm 2.