多项式全文背诵

不要把技能树点歪了。

math-数学前置

复数

$x+yi$ 的数称为复数，其中 $i = \sqrt{-1}$。

一般把复数放到复平面上，那么复数 $x+yi$ 就可以看做一个平面向量 $\large{x \brack y}$，用两个 double 存一下就好了。

乘法运算法则：$(x+yi)(x'+y'i)=(xx'-yy')+(xy'+x'y)i$

复数类

struct Complex { // 本质就是向量
    double x, y;
    Complex(double x_ = 0, double y_ = 0) {x = x_, y = y_;} // 不要写成 int
    inline double len() {return sqrt(x * x + y * y);}
    friend Complex operator + (Complex a, Complex b) {return Complex(a.x + b.x, a.y + b.y);}
    friend Complex operator - (Complex a, Complex b) {return Complex(a.x - b.x, a.y - b.y);}
    friend Complex operator * (Complex a, Complex b) {return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);} // 相乘，用 i^2 = -1 化简
};

单位根

$n$ 次单位根：就是把复平面上单位圆 $n$ 等分（从 $(1,0)$ 开始），从圆心到这些点构成的向量中幅角为正且最小的向量对应的复数。

单位根公式，证明看图就行：

\[\omega_{n}^{k} = \cos(k\frac{2 \pi}{n}) + \sin(k\frac{2 \pi}{n})i \]

pic

推论 1：$\omega_{n}^{0}=\omega_n^n=1$。
推论 2：若 $k < \frac{n}{2}, 2 | n$，则有 $\omega_{n}^{k}=-\omega_{n}^{k+\frac{n}{2}}$。

证明：本质等价于向量旋转 180 度。

推论三：$\omega_{n}^{k}=\omega_{2n}^{2k}=\omega_{n}^{k+n}$。

poly-多项式

明确：FFT, NTT 本质都是把多项式在点值表示法和系数表示法之间转化，也就是快速求出多项式在一些点（横坐标是单位根的点）上的取值。

FFT

为了方便，若输入多项式不是 $2^k$ 项，在后面补 $0$ 使其变为 $2^k$ 项。

对于 $n-1=2^k-1$ 次多项式 $A(x) = a_0 + a_1x+a_2x^2 + \dots + a_{n-1}x^{n-1}$，首先对下标奇偶分类：

\[\left\{ \begin{aligned} A_l(x) &= a_0 + a_2x^1 + \dots + a_{n-2} x^{\frac{n}{2}-1}\\ A_r(x) &= a_1 + a_3x^1 + \dots + a_{n-1} x^{\frac{n}{2}-1} \end{aligned} \right. \]

不难观察到：

\[A(x) = A_l(x^2)+xA_r(x^2) \]

现在，如果我们取 $n$ 次单位根作为这个 $n-1$ 次多项式点值表达式的 $n$ 个关键点的横坐标带入：

\[\begin{aligned} A(\omega_{n}^{k}) &= A_l(\omega_{n}^{2k})+\omega_{n}^{k}A_r(\omega_{n}^{2k})\\ &= A_l(\omega_{n/2}^{k})+\omega_{n}^{k}A_r(\omega_{n/2}^{k})\\ \\ A(\omega_{n}^{k+n/2}) &= A_l(\omega_{n}^{2k+n})-\omega_{n}^{k}A_r(\omega_{n}^{2k+n})\\ &= A_l(\omega_{n/2}^{k})-\omega_{n}^{k}A_r(\omega_{n/2}^{k})\\ \end{aligned} \]

奇妙的事情发生了：$A_l(\omega_{n/2}^{k})$ 和 $A_r(\omega_{n/2}^{k})$ 转化为大小为 $\frac{n}{2}$ 的子问题，因此直接递归求解就行。

但这只是将系数->点值，点值->系数又该怎么办呢？很简单，死记硬背一下，把传单位根的时候 $\sin$ 前面加个负号就行。

不要忘了，最后答案要除 $n$！

递归实现

inline void SFT(int len, vector<Complex> &a, int op) {
    if (len == 1) return;
    vector<Complex> la, ra; la.resize(len >> 1), ra.resize(len >> 1);
    for (int i = 0; i < len; i += 2) la[i >> 1] = a[i], ra[i >> 1] = a[i + 1]; // 奇偶分类。
    SFT(len >> 1, la, op), SFT(len >> 1, ra, op); // 此时不难发现有 a = la(x^2) + x * ra(x^2)
    Complex omega0 = Complex(cos(PI * 2.0 / len), sin(PI * 2.0 / len) * op), omega = Complex(1, 0); // len 次单位根公式
    for (int i = 0; i < (len >> 1); ++i, omega = omega * omega0) {
        Complex butterfly = omega * ra[i]; // 蝶形优化
        a[i] = la[i] + butterfly, a[i + (len >> 1)] = la[i] - butterfly; // 通过单位根性质节约复杂度。可以推式子，在多项式 a 中带入 w(i) 与 w(i + (len / 2)) 化简得到。
    } 
}

但这样太慢了，有没有更快的写法？

从下往上枚举 SFT 递归的每一层，每一个分治区间处理一下，就可以非递归了。注意为了使我们每一项在最开始就在分治结束后的位置，我们需要一个奇妙东西叫位逆序置换 $rev$，将 $i$ 与 $rev_i$ 交换就会得到分治结束后的 $a$。（代码为了更快，把 $rev$ 预处理放外面了）

而位逆序置换的递推公式是 $rev_i = \frac{rev_{i/2}}{2}+(i \bmod 2 = 1) \times \frac{n}{2}$。

非递归实现

inline void FFT(int LEN, vector<Complex> &a, vector<int> &rev, int op) {
    REP(i, 0, LEN - 1) if (i < rev[i]) swap(a[i], a[rev[i]]); // 一步到位的奇偶分类
    for (int len = 2; len <= LEN; len <<= 1) { // 从下往上枚举 SFT 中分治层
        Complex omega0 = Complex(cos(PI * 2.0 / len), sin(PI * 2.0 / len) * op); // 此时 len 次单位根可求。
        for (int l = 0; l < LEN; l += len) { // 枚举当前分治层分治区间（左端点）
            Complex omega = Complex(1, 0);
            for (int i = l; i < l + (len >> 1); ++i, omega = omega * omega0) { // 和 SFT 中一模一样
                Complex butterfly = omega * a[i + (len >> 1)];
                a[i + (len >> 1)] = a[i] - butterfly, a[i] = a[i] + butterfly;
            }
        }
    }
}

多项式乘法

就是把 $A,B$ 转成点值表示法，对位相乘，最后用逆变换转回来。

有个小优化：将 $A$ 放到实部，$B$ 放到虚部，这样做一遍 FFT，平方，最后出来的答案的虚部就是两倍的 $A \times B$ 的值了，最后只需要两遍 FFT，常数更小！

证明就是写出式子用平方和公式暴力展开，取中间项即可得证。

两遍 FFT 多项式乘法实现

inline vector<int> Mul(vector<int> &a, vector<int> &b) {
    int len = 1, Max = (int)a.size() + (int)b.size() - 2;
    while (len <= Max) len <<= 1;
    vector<int> c(Max + 1), rev(len);
    vector<Complex> fa(len);
    REP(i, 0, len - 1) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0); // 位逆序置换.
    REP(i, 0, len - 1) fa[i] = Complex(i < a.size() ? a[i] : 0, i < b.size() ? b[i] : 0);
    FFT(len, fa, rev, 1);
    REP(i, 0, len - 1) fa[i] = fa[i] * fa[i];
    FFT(len, fa, rev, -1);
    REP(i, 0, Max) c[i] = ((int)(fa[i].y / len + 0.5) / 2);
    return c;
}

NTT

适用于模数有原根的情况，如 $998244353$，它的原根为 $3$，模数为 $10^9+7$ 就用不了一点。

把 FFT 中的 Complex 类换成整形，点值转系数的时候将单位根换成原根，系数转点值的时候将单位根换成原根的逆元就可以避免精度误差。因为原根具有单位根在这个算法中所用到的性质。

和 FFT 同理，最后要乘 $n$ 的逆元。

非递归 NTT 实现（$P=998244353$）

const ll g = 3;
const ll gi = 332748118;
const ll P = 998244353;
inline ll add(ll x) {return (x >= P) ? x - P : x;}
inline ll sub(ll x) {return (x < 0) ? x + P : x;}
inline void NTT(int LEN, vector<ll> &a, vector<int> rev, int op) {
    REP(i, 0, LEN - 1) if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int len = 2; len <= LEN; len <<= 1) {
        ll omega0 = qpow(op == 1 ? g : gi, (P - 1) / len);
        for (int l = 0; l < LEN; l += len) {
            ll omega = 1;
            for (int i = l; i < l + (len >> 1); ++i, omega = omega * omega0 % P) {
                ll butterfly = omega * a[i + (len >> 1)] % P;
                a[i + (len >> 1)] = sub(a[i] - butterfly), a[i] = add(a[i] + butterfly);
            }
        }
    }
}

MTT

任意模数多项式乘法。

设 $F(x),G(x)$ 的系数分别为 $f_{0\dots n-1},g_{0 \dots n-1}$。

FFT 的问题在于值域过大精度会爆，所以我们拆系数，不妨令 $f_i=(f_i \bmod B) + (\frac{f_i}{B})B$，$g_i=(g_i \bmod B) + (\frac{g_i}{B})B$，当 $B$ 取 $\sqrt{V}$ 时，值域缩小为 $O(\sqrt{V})$。

不妨令 $A_i = (f_i \bmod B), B_i = (\frac{f_i}{B})$, $C_i = (g_i \bmod B), D_i = (\frac{g_i}{B})$, $ T_i = C_i + D_ii$。

将 $A,B,T$ 转成点值形式，乘起来得到 $A \times T, B \times T$，将这两个多项式转回系数形式，发现它们的实、虚部中正好存储着 $A\times C,A \times D,B \times C, B \times D$。

而我们最终的乘积系数：

\[\begin{aligned} h_i &= (A_i + B_iB) \times (C_i + D_iB)\\ &= A_iC_i + B(A_iD_i + B_iC_i) + B^2(B_iD_i) \end{aligned} \]

可以快速计算。

五次 FFT 实现

inline vector<int> MTT(int n, int m, vector<int> &F, vector<int> &G) {
    vector<int> H(n + m + 1);
    int len = 1;
    while (len <= n + m) len <<= 1;
    vector<Complex> A(len), B(len), T(len), AT(len), BT(len); rev.resize(len);
    REP(i, 0, len - 1) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
    REP(i, 0, n) A[i] = Complex(F[i] & Mask, 0), B[i] = Complex(F[i] >> 15, 0);
    REP(i, 0, m) T[i] = Complex(G[i] & Mask, G[i] >> 15);
    FFT(len, A, 1), FFT(len, B, 1), FFT(len, T, 1);
    REP(i, 0, len - 1) AT[i] = A[i] * T[i], BT[i] = B[i] * T[i];
    FFT(len, AT, -1), FFT(len, BT, -1);
    REP(i, 0, n + m) H[i] = (((ll)(AT[i].x / len + 0.5) % P) + (((ll)(AT[i].y / len + 0.5) % P) << 15ll) % P + (((ll)(BT[i].x / len + 0.5) % P) << 15ll) % P + (((ll)(BT[i].y / len + 0.5) % P) << 30ll) % P) % P;
    return H;
}

多项式求逆

最坑的一个。

设我们求出了 $f_0^{-1}(x) \equiv f^{-1}(x) \bmod x^{n/2}$ 考虑如何求出 $f_1^{-1}(x)\equiv f^{-1}(x) \bmod x^{n}$。

\[\begin{aligned} &\left\{\begin{aligned} f(x)f_1^{-1}(x) \equiv f(x)f_0^{-1}(x) \equiv 0 \bmod x^{n/2}\\ f_1^{-1}(x)-f_0^{-1}(x) \equiv 0 \bmod x^{n/2}\\ \end{aligned}\right.\\ \\ &\begin{aligned} (f_1^{-1}(x)-f_0^{-1}(x))^2 &\equiv 0 \bmod x^{n/2}\\ f_1^{-2}(x)+f_0^{-2}(x)-2f_0^{-1}(x)f_1^{-1}(x)&\equiv 0 \bmod x^{n/2}\\ f_1^{-1}(x)+f(x)f_0^{-2}(x)-2f_0^{-1}(x)&\equiv 0 \bmod x^{n/2} \end{aligned} \\ \\ &f_1^{-1}(x) = 2f_0^{-1}(x) - f(x)(f_0^{-1}(x))^2 \end{aligned} \]

直接倍增递推即可，注意 $f(x)$ 项数的问题！！！

NTT 实现多项式乘法逆

inline long long qpow(long long x, long long y) {
    long long Res = 1;
    while (y) {
        if (y & 1) (Res *= x) %= P;
        (x *= x) %= P, y >>= 1;
    } return Res;
} 
inline vector<int> PolyInv(vector<int> a) {
    vector<int> Inv[2]; int op = 0; Inv[0].emplace_back(qpow(a[0], P - 2));
    for (int len = 2; len >= 0; ((len <= a.size()) ? (len <<= 1) : (len = -1))) {
        op ^= 1, Inv[op] = Inv[op ^ 1], Inv[op].resize(len);
        for (int i = 0; i < len; ++i) Inv[op][i] = add(Inv[op][i] << 1);
        Inv[op ^ 1] = Mul(Inv[op ^ 1], Inv[op ^ 1]), Inv[op ^ 1].resize(len), Inv[op ^ 1] = Mul(Inv[op ^ 1], a), Inv[op ^ 1].resize(len);
        for (int i = 0; i < len; ++i) Inv[op][i] = sub(Inv[op][i] - Inv[op ^ 1][i]);
    }
    return Inv[op];
}

任意模数多项式乘法逆

inline long long qpow(long long x, long long y) {
    long long Res = 1;
    while (y) {
        if (y & 1) (Res *= x) %= P;
        (x *= x) %= P, y >>= 1;
    } return Res;
} 
inline int add(int x) {return (x >= P) ? (x - P) : x;}
inline int sub(int x) {return (x <  0) ? (x + P) : x;}
inline vector<int> polyinv(vector<int> &F) {
    vector<int> Inv[2]; int op = 0; Inv[0].emplace_back(qpow(F[0], P - 2));
    for (int len = 2; len >= 0; (len <= F.size() ? len <<= 1 : len = -1)) {
        op ^= 1, Inv[op] = Inv[op ^ 1], Inv[op].resize(len);
        for (int i = 0; i < len; ++i) Inv[op][i] = add(Inv[op][i] << 1);
        Inv[op ^ 1] = MTT(Inv[op ^ 1], Inv[op ^ 1]), Inv[op ^ 1].resize(len), Inv[op ^ 1] = MTT(Inv[op ^ 1], F), Inv[op ^ 1].resize(len);
        for (int i = 0; i < len; ++i) Inv[op][i] = sub(Inv[op][i] - Inv[op ^ 1][i]);
    }
    return Inv[op];
}

$O(n \log n)$。

分治 FFT

就是分治，先算左区间，再算右区间，区间信息合并用个 FFT 就行了。

$O(n \log^2 n)$。

不如推式子+多项式求逆。

多项式 exp

咕咕咕。

多项式 ln