多项式求逆,顾名思义就是求一个多项式的逆元。
其实数学中很多元素都有逆元素,比如对于一个非零实数 $r$,那么它的逆就是 $\dfrac 1r$,方便起见,记作 $r^{-1}$。在线性代数中,对于一个矩阵 $\mathbf A$,如果它的行列式 $\det \mathbf A = \left| \mathbf A \right| \neq 0$,那么就会有逆矩阵 $\mathbf A^{-1} = \dfrac {\mathbf A^*} {|\mathbf A|}$。对于一个函数 $f(x)$,它也有它自己的反函数 $f^{-1} (x)$。在群论中,任何一个群都要求有逆元素,等等。
那么,对于多项式,它的逆又是什么捏?它有逆元素的充要条件是什么?
显然,对于一个多项式 $f(x)$,如果要找到一个多项式 $g(x)$,使得 $f(x) \cdot g(x) = 1$,那么 $f(x), g(x)$ 都只能为 (非零的) 常数多项式,否则 $g(x)$ 就是无穷级数 (比如 $(1-x) (1+x+x^2+x^3 + \cdots) = 1$),那也太扫兴了。因此我们可以稍稍减弱条件 (像模意义逆一样),只需在 $\bmod x^m$ 意义下为 $1$ 即可 (即一次项到 $m$ 次项均为 $0$)。
对于多项式 $f(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0$,如果存在低于 $m$ 次的多项式 $g(x) = b_{m-1} x^{m-1} + b_{m-2} x^{m-2} + \cdots + b_2 x^2 + b_1 x + b_0$,使得 $f(x) \cdot g(x) \equiv 1 \pmod {x^m}$ (即 $f(x) \cdot g(x) = 1 + \sum\limits_{i > m} c_i x^i$),则称 $g(x)$ 为 $f(x)$ 在模 $x^m$ 意义下的逆元,记作 $f^{-1} (x) \pmod {x^m}$。
那是不是所有的多项式都有逆元吗?如果多项式 $f(x)$ 的常数项 $a_0 = 0$ 的话,那么对于任意的 $g(x)$,都有 $c_0 = a_0 \cdot b_0 = 0$,因此 $f(x)$ 不存在逆元。
如果 $a_0 \neq 0$,那么是否一定存在逆元呢?答案是肯定的。我们按照上述设 $f(x)$ 和 $g(x)$,则
$$ f(x) \cdot g(x) = \left( \sum_{i=0}^{n-1} a_i x^i \right) \cdot \left( \sum_{j=0}^{n-1} b_j x^j \right) = \sum_i \sum_{j=0}^i a_{i-j} b_j x^i = \sum_i c_i x^i = 1 + \sum_{i > m} c_i x^i $$
故可得如下方程组:
\begin{cases} a_0 \cdot b_0 = c_0 = 1 \\ a_0 \cdot b_1 + a_1 \cdot b_0 = c_1 = 0 \\ a_0 \cdot b_2 + a_1 \cdot b_1 + a_2 \cdot b_0 = 0 \\ \cdots \\ a_0 \cdot b_{m-1} + a_1 \cdot b_{m-2} + \cdots + a_{m-1} \cdot b_0 = 0 \end{cases}可以轻松解得 $b_0 = \dfrac 1 {a_0}, b_1 = - \dfrac {a_1 \cdot b_0} {a_0}, \cdots, b_k = - \dfrac {\sum\limits_{j=0}^{k-1} a_{k-j} \cdot b_j} {a_0}$ ($k = 1, 2, \cdots, m-1$),因此我们就求得了它们的系数,$m$ 个未知数恰好 $m$ 个方程,因此解也是唯一的。
因此得到结论:一个多项式有逆元当且仅当它的常数项存在逆元,且一个多项式的逆元是唯一的。
刚刚的求解 (证明) 过程中我们已经得到了一个暴力算法,就是按照上面的公式去计算,时间复杂度 $O(m^2)$。
伪代码如下:
b0 ← 1 / a0
for i ← 1 to m - 1
bi ← 0
for j ← 0 to i - 1
bi ← bi - bj * ai - j
bi ← bi / a0
显然这个算法还不太优秀。可以想一想,多项式乘法可以从暴力的 $O(n^2)$ 优化到快速傅里叶变换的 $O(n \log n)$,那么多项式求逆是否也能优化呢?显然可以,使用倍增思想即可。
前面提到了倍增,那么我们来实现一下。如果知道了 $b_0, b_1, b_2, \cdots, b_{k-1}$ ($k > 0$),我们希望通过 (不是暴力的) 方法求得 $b_k, b_{k+1}, b_{k+2}, \cdots, b_{2k-1}$。
此时,我们可以抛开模意义去谈这件事情了。如果一个无穷级数是一个多项式的逆,则它的系数就被无穷多个方程组成的方程组约束。前面也看到过,如果模的是 $m$ 次,就相当于取前 $m$ 个方程,最终能解出来的,当然也就是前 $m$ 项喽。
现在已知 $k$ 项,那么它们就应该满足前 $k$ 个方程,即这个 $k-1$ 次多项式为 $g_k (x)$,因此有 $f(x) g_k (x) \equiv 1 \pmod {x^k}$,根据前述,有 $g_{2k} (x) \equiv g_k (x) \pmod {x^k}$,因此有
$$ \left( g_{2k} (x) - g_k (x) \right) \equiv 0 \pmod {x^k} $$
两边同时平方 (注意模也可以平方),就得到了 $$ \left( g_{2k} (x) - g_k (x) \right)^2 = g_{2k}^2 (x) - 2 g_{2k} (x) g_k (x) + g_k^2 (x) \equiv 0 \pmod {x^{2k}} $$
只需两边同乘 $f(x)$,就得到了:$$ g_{2k} (x) - 2 g_k (x) + f(x) g_k^2 (x) \equiv 0 \pmod {x^{2k}} $$
也就是说,$$ g_{2k} (x) = g_k (x) \left( 2 - f(x) g_k (x) \right) \pmod {x^{2k}} $$
话说你们看到这个式子有没有感到很熟悉?我们在使用牛顿迭代求 $a$ 的倒数 (即 $f(x) = ax - 1 = 0$ 的根) 的时候,所得到的公式也是 $b_{n+1} = b_n \left( 2 - a_n b_n \right)$,有着异曲同工之妙,因此这种方法亦称为多项式的牛顿迭代。
好的,回到正题。那么它其实只是几次多项式乘法而已。显然可以通过 FFT/FNTT 来优化到 $O(k \log k)$ (只需求出 $f$ 和 $g_k$ 的点值形式即可)。
总复杂度即由下列递推式确定:$T(n) = T(n/2) + O(n \log n)$,不管怎么搞都能求出 $T(n) = O(n \log n)$,经实(li)际(lun)测(fen)试(xi),常数大概是单次多项式乘法的两倍。
#include <bits/stdc++.h>
#define N 1109771
using namespace std;
typedef long long ll;
const ll mod = 998244353, half_mod = 499122177, root = 31;
int l, n, D;
int rev[N];
ll x[N], y[N], t[N];
ll p[N], buf[N], ans[N];
ll PowerMod(ll a, int n, ll c = 1) {for(; n; n >>= 1, (a *= a) %= mod) if(n & 1) (c *= a) %= mod; return c;}
void NTT_init(int length){
n = 1 << (l = length);
ll g = PowerMod(root, 1 << 23 - l);
x[0] = 1; rev[0] = 0;
for(int i = 1; i < n; ++i){
x[i] = x[i - 1] * g % mod;
rev[i] = (i & 1 ? rev[i - 1] | n >> 1 : rev[i >> 1] >> 1);
}
}
void DNTT(ll *d){
int i, j, k, len = 1, delta = n; ll r;
for(i = 0; i < n; ++i) t[i] = d[rev[i]];
for(i = 0; i < l; ++i){
delta >>= 1;
for(k = j = 0; j < len; ++j, k += delta) y[j] = x[k];
for(j = 0; j < n; j += len << 1)
for(k = j; k < j + len; ++k){
r = y[k - j] * t[k + len] % mod;
(t[k + len] = t[k] - r) < 0 ? t[k + len] += mod : 0;
(t[k] += r) >= mod ? t[k] -= mod : 0;
}
len <<= 1;
}
// memcpy(d, t, n << 3);
}
void PolyInv(ll *a, ll *b, int degree){
int len, i; ll inv = half_mod;
a[0] = p[0]; a[1] = p[1];
// B_{i+1} = B_i (2 - A B_i)
for(len = 0; 1 << len < degree; ++len){
NTT_init(len + 2);
DNTT(b); memcpy(b, t, n << 3); DNTT(a);
for(i = 0; i < n; ++i)
((b[i] *= 2 - t[i] * b[i] % mod) %= mod) < 0 ? b[i] += mod : 0;
DNTT(b); reverse(t + 1, t + n);
inv = (inv >> 1) + (inv & 1 ? half_mod : 0);
for(i = 0; i < 1 << len + 1; ++i) b[i] = t[i] * inv % mod;
for(; i < n; ++i) {b[i] = 0; a[i] = p[i];}
}
}
int main(){
int deg, Modeg, i;
scanf("%d%d", °, &Modeg);
for(i = deg; i; --i) scanf("%lld", p + i);
ans[0] = p[0] = 1;
PolyInv(buf, ans, Modeg);
for(i = Modeg - 1; i; --i) printf("%lld%c", ans[i], i == 1 ? 10 : 32);
return 0;
}