Bostan-Mori, Tellegen's Principle, Power series composition

Bostan-Mori algorithm

흔히 "선형 점화식의 $N$번째 항"을 구하는 알고리즘으로 알려진 방법이다. 사실은 보다 일반적으로, $P(x)/Q(x)$ 꼴의 rational power series의 $N$차항 계수를 구하는 데 사용할 수 있다.

Theorem. $\deg P < \deg Q = d$, $Q(0) = 1$을 만족하는 다항식 $P, Q \in \mathbb{F}[x]$ 에 대해, $[x^N]P(x)/Q(x)$를 $\mathcal{M}(d)\log N$ 시간에 계산할 수 있다.

• 어떤 다항식 / power-series의 $N$차항 계수를 $[x^N]f(x)$로 쓴다.
• $\mathcal{M}(d)$는 $\mathbb{F}[x]$의 두 $d$차 다항식을 곱하는 데 걸리는 시간이다. Computer algebra에서는 sublogarithmic factor를 따지지만, PS 범위에서는 $\mathcal{M}(d) = O(d\log d)$로 생각해도 무방하다.
• $\mathbb{F}$는 덧셈, 뺄셈, 곱셈, 나눗셈이 충분히 빠른 어떤 체로 생각한다. $\mathbb{F} _ {p}$로 생각하면 제일 편하다.

Proof.

• $N=0$이라면 $P(0) / Q(0)$를 계산하면 된다.

• $N \ge 1$인 경우, 분자 분모에 $Q(-x)$를 곱해 분모를 even function으로 만들고, 분자를 odd-even으로 분리하자. 즉 $P(x)/Q(x) = P(x)Q(-x)/Q(x)Q(-x)$, $P(x)Q(-x) = U_{0}(x^2) + xU_{1}(x^2)$, $Q(x)Q(-x) = V(x^2)$으로 쓴다. 이 때 $$[x^N]P(x)/Q(x) = [x^{\lfloor N/2 \rfloor}] U_{N \bmod 2} (x) / V(x)$$가 성립하고, $\deg U < d = \deg V$이자 $V(0) = 1$이니 재귀적으로 해결할 수 있다.

따라서 $d$차식 곱셈 $O(\log N)$번으로 전체 문제를 해결할 수 있으니 시간 복잡도는 $O(\mathcal{M}(d) \log N)$이다. $\mathcal{F}$가 $\mathcal{Z} _ {998244353}$과 같이 DFT를 사용하기 좋은 환경인 경우, DFT Doubling을 활용하여 $\frac{2}{3}\mathcal{M}(d)\log N$ 까지 상수를 줄일 수 있음이 알려져 있다. 자세한 내용은 Aeren님의 글을 참조.

연습문제:

• BOJ 13725 RNG : "선형 점화식의 $N$차항"을 그대로 옮긴 문제이다. $P, Q$를 비교적 쉽게 계산할 수 있다.
• BOJ 13178 목공 : 꼴이 비선형 점화식이지만, 생성함수를 $x / (1-x)f(x)$ 꼴로 적을 수 있고 그 뒤는 Bostan-Mori algorithm을 적용하면 된다.

Tellegen's Principle

Degree가 각각 $n, m$인 두 다항식 $f, g \in \mathbb{F}[x]$를 곱하는 연산을 생각하자.
뜬금없지만, 이를 $(n + m + 1) \times (n+1)$ 행렬 $M _ {g} = ([x^{i-j}]g(x) ) _ {ij}$ 를 $\vec{f} = (f _ 0, \cdots, f _ n)$에 곱하는 행렬 연산으로 볼 수도 있다. ($f = \sum f _ {i} x^{i}, g = \sum g _ {i} x^{i}$로 편의상 표기)

다시 한 번 뜬금없지만, $M_{g}^{T}$는 어떤 연산자일까? 임의의 크기 $n+m-1$인 벡터 $w = (w_{0}, \cdots, w_{n + m})$를 생각하자. $({M _ {g}}^{T}w) _ {i} = \sum_{j=0}^{n+m} [x^{j-i}]g(x) w _ {j} = \sum_{j=0}^{m} g _ {j} w _ {i + j}$ 가 된다. 즉 $(\vec{g} \cdots w_{0:m}, \cdots \vec{g} \cdots w_{n:n+m})$ 과 같은 convolution 꼴들을 여럿 계산하게 된다.

사실 convolution과 다항식 곱셈이 비슷한 꼴을 갖는 것이 전혀 놀라운 일은 아니다. 그리고 둘을 계산하는 시간 복잡도가 동일하다는 사실도 잠시 생각해보면 그리 납득하기 어렵지 않다. Tellegen's principle은 이러한 "동치 관계"에 보다 일반적인 대응을 제시한다.

Thm. (Tellegen's principle) $m \times n$ 행렬 $M$에 대해 $f \mapsto Mf$를 계산하는 프로그램 $L$의 시간 복잡도가 $T_{1}$이라면, $g \mapsto M^{T}g$를 계산하는 $T_{1} - n + m + z_{0} - z_{1}$ 복잡도의 프로그램이 존재한다. 단, $z_0, z_1$은 각각 $M$의 zero row, zero column의 개수이다.

아무 프로그램이나 되는 것은 아니고, algebraic complexity theory를 하시는 분들의 정의인 "linear straight line program"이어야 한다.
간단히 말해, 모든 연산이 $x _ {i} \leftarrow a x _ {j} + b x _ {k}$ 처럼 모든 변수 ($x _ {\ast}$)에 대해 linear한 "단순한" 연산의 나열로 구성된 프로그램이어야 한다.

Linear Straight line program의 예시로는 Cooley-Tukey DFT (우리가 아는 그 DFT), Karatsuba algorithm 등이 있다.

당연히 두 변수 $x, y$에 대해 $x \leftarrow ax + by$ 같은 atomic한 연산은 $\begin{pmatrix} a & b \\ 0 & 1 \end{pmatrix}$로 나타날테니, 이 연산의 transpose는 $x \leftarrow ax, y \leftarrow bx + y$와 같이 간단히 쓸 수 있을 것이다. 그리고 $(f \circ g)^{t} = g^{t} \circ f^{t}$일 테니 모든 연산을 역순으로 조립하면 이론상 transposed program을 만들 수는 있다. 하지만 그 이상으로 Tellegen's principle이 명쾌한 conversion을 제공하진 않으며, 주로 몇 가지 building block들을 조합하여 transposed program을 만든다. Tellegen's principle into practice에서 polynomial division, polynomial multiplication 등 자주 쓰이는 연산에 대한 transposed program을 계산해두었다.

서로 transpose 관계인 문제들의 예시들은 아래와 같다. 각각은 input vector $f$ (혹은 polynomial)에 대해 linear함에 주의하자. 이외의 변수들에 대해서는 linear하지 않을 수 있다.

• (Primal) Multipoint Evaluation $f \mapsto f(a _ 0), \cdots f(a _ {q-1})$
• (Dual) Generalized Power Sum $(a _ {0}, \cdots, a _ {q-1}) \mapsto \sum_{i=0}^{q-1} a_{i}^{d}$ for $d = 0, \cdots, n$

이건 계산은 안해봤다.

• (Primal) Shift: $f(x) \mapsto f(x + 1)$

• (Dual) Falling factorial basis $\sum f_{j} x^{j} = \sum g_{j} x(x-1)\cdots (x-j+1)$일 때, $f \mapsto g$

• (Primal) Modular Composition: $f(x) \mapsto f(g(x)) \mod h(x)$. (단, $g(0) = 0$)

• (Dual) Power Projection: $(w_{0}, \cdots, w_{n-1})$와 다항식 $f$에 대해 $w(f) = \sum_{i} w_{i}f_{i}$로 정의할 때, $w( g(x)^{k} \mod h )$ for $k = 0, \cdots, \deg f$.

Power Series Composition in Near-Linear Time (FOCS '24)

주어진 다항식 $f, g, h$에 대해 modular composition $f(g(x)) \mod h(x)$를 계산하는 문제는 computer algebra 전반에 있어 관련성이 깊은 문제이다. 현재 이 문제에 대한 일반적인 near-linear time algorithm은 알려져 있지 않고, $O(n^{1.5})$ 에 동작하는 버전과 빠른 행렬곱에 의존하는 버전이 있는 것으로 알고 있다. (TODO: 참고자료 추가)

$g(x) = \log(1+x), e^{x} - 1$과 같이 특수한 경우에 대해서는 다른 문제로 formulate한 성과가 있었지만, 이번에 $h(x) = x^{n}$ 이라는, 상당히 일반화된 조건에서 이 문제를 해결한 논문이 발표되었다. 결과도 꽤 간단한 편인데, Bostan-mori algorithm과 Tellegen's principle을 배경지식으로 알고 있으면 좋다.

참고로 1저자인 Yasunori Kinoshita는 2024 ICPC WF에서 4위를 차지하였다.

논문에서는 Tellegen's principle에 의한 dual이 되는 Power Projection의 near-linear algorithm을 먼저 제시한다. 이로부터 자연스럽게 modular composition을 near linear time에 해결할 수 있는 셈이다.

$w^{R}(x) = \sum_{j=0}^{n-1} w_{n-1-j} x^{j}$라고 두자.

Power projection 문제의 답을 $p_{k} = w( g(x)^{k} \mod x^{n} )$라고 하면,

$$p_{k} = [x^{n-1}] w^{R}(x) (g(x)^{k} \mod x^{n}) = [x^{n-1}] w^{R}(x) g(x)^{k}$$
로 쓸 수 있다. 이제 generating function $\sum p_{k} y^{k}$를 생각하면

$$\sum_{k} p_{k}y^{k} = \sum_{k} [x^{n-1}] w^{R}(x) y^{k} g(x)^{k} = [x^{n-1}]\frac{ w^{R}(x) }{ 1 - yg(x) } \mod y^{m}$$

와 같이 쓸 수 있다. 여기서 $m = \deg g$.

Bostan-mori algorithm의 모든 부분은 division free라는 것을 생각하면, $\mathbb{F}[y][x]$에서도 대부분의 논리를 그냥 가져올 수 있다.

이변수 다항식의 분수꼴 $P(x, y) / Q(x, y)$의 $[x^{N}]$을 알고 싶다고 하자. 이는 $y$에 대한 다항식일 것이다. 똑같이 $Q(-x, y)$를 분자 분모에 곱해주고 $x$-차수의 기우성에 따라 분리해주면 된다. 이 때 $N$차항보다 더 큰 항은 앞으로 필요없으니 버리고, 재귀적으로 문제를 해결한다.

물론 이변수 다항식을 그렇게 속편하게 곱할 수는 없다. 다음의 정리를 사용하자.

Proposition. 두 이변수 다항식 $P(x, y), Q(x, y)$가 모두 $\deg_{x} < m, \deg_{y} < n$를 만족하면, $PQ$를 $O(\mathcal{M}(mn))$ 시간에 계산할 수 있다.

사실 $y = x^{2m-1}$ 정도로 치환해서 곱해주면 된다. 흔히 Kronecker's substitution이라고도 부른다.

이 과정의 시간 복잡도를 $P(x, y) = w^{R}(x)$, $Q(x, y) = 1 - yg(x)$, $N = n-1$에 대해 분석하면, $i$번째 iteration (즉, $N \leftarrow N / 2^{i}$인 시점) 에서는 $\deg_{x}Q \le \min(m-1, n/2^{i}), \deg_{y} Q \le \min(m-1, 2^{i+1})$이 성립한다. 마찬가지로 $\deg_{x} P \le n/2^{i}, \deg_{y} P \le 2^{i+1}$이 성립하므로 $2^{i} \le m$이라면 $i$번째 iteration에서 시간 복잡도는 $\mathcal{M}(n)$이다.

$2^{i} > m$인 경우, $x$에 대한 차수는 계속 줄어들지만 $y$에 대한 차수가 여전히 $m-1$까지 커질 수 있기 때문에 분석을 따로 해주어야 한다. $m < 2^{i} < n$에 대해 $\mathcal{M}(mn / 2^{i})$만큼의 시간이 들게 되는데 이는 $O(\mathcal{M}(n))$으로 bound 된다.

따라서 전체 시간 복잡도는 $O(\mathcal{M}(n)\log m)$이다.

당연히 Tellegen's principle만으로 이걸 뒤집어서 modular composition을 재현하기는 매우 어렵다. 본문의 4.2절에 직접 composition을 계산하는 부분이 있는데, 위 알고리즘과 대동소이하므로 궁금하신 분들은 일독을 권한다.

Reference

• Kinoshita, Yasunori, and Baitian Li. "Power Series Composition in Near-Linear Time." arXiv preprint arXiv:2404.05177 (2024). https://arxiv.org/pdf/2404.05177

• Bürgisser, Peter, Michael Clausen, and Mohammad A. Shokrollahi. Algebraic complexity theory. Vol. 315. Springer Science & Business Media, 2013. https://link.springer.com/book/10.1007/978-3-662-03338-8

• Bostan, Alin, Grégoire Lecerf, and Éric Schost. "Tellegen's principle into practice." Proceedings of the 2003 international symposium on Symbolic and algebraic computation. 2003.https://dl.acm.org/doi/pdf/10.1145/860854.860870