在掌握了多项式的加法和乘法，并且通过 FFT 和 NTT 将时间复杂度降到了可以接受的 O(nlog⁡n)\mathcal{O}(n \log n)，我们就可以完成许多代数可以完成的基本运算了。

多项式乘法

P3803 【模板】多项式乘法（FFT） - 洛谷

前面已经详细讲过了，也是最基础的，这里放下完整代码。

FFT

记录详情 - 洛谷 | 计算机科学教育新生态

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
struct Complex
{
    double a, b;
    Complex(double _a = 0, double _b = 0) : a(_a), b(_b) {}
    
    Complex operator + (const Complex &rhs) const { return Complex(a + rhs.a, b + rhs.b); }
    Complex operator - (const Complex &rhs) const { return Complex(a - rhs.a, b - rhs.b); }
    Complex operator * (const Complex &rhs) const { return Complex(a * rhs.a - b * rhs.b, a * rhs.b + b * rhs.a); }
};

template <typename Container>
void FFT(int T, Container &A, int type)
{
    static int Rev[N];
    for (int i = 0; i < T; i ++) Rev[i] = (Rev[i >> 1] >> 1) | (i & 1 ? T >> 1 : 0);
    for (int i = 0; i < T; i ++) if (i < Rev[i]) swap(A[i], A[Rev[i]]);
    
    for (int cur = 1; cur < T; cur <<= 1)
    {
        Complex base(cos(PI / cur), sin(PI / cur) * type);

        static Complex comp[N];
        comp[0] = Complex(1, 0);
        for (int i = 1; i < cur; i ++) comp[i] = comp[i - 1] * base;

        for (int L = 0; L < T; L += cur << 1)
            for (int j = 0; j < cur; j ++)
            {
                Complex tempx = A[L | j], tempy = comp[j] * A[L | j | cur];
                A[L | j] = tempx + tempy, A[L | j | cur] = tempx - tempy;
            }
    }

    if (type == 1) return;
    for (int i = 0; i < T; i ++) A[i].a /= T;
}
// void FFT(int T, Complex *A, int type)
// {
//     vector<Complex> temp(A, A + T);
//     FFT(T, temp, type);
//     copy(temp.begin(), temp.end(), A);
// }

NTT

记录详情 - 洛谷 | 计算机科学教育新生态

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
const int mod = 998244353, ord = 3, ordinv = 332748118;

inline int ksm(int base, int k)
{
    int res = 1;
    while (k)
    {
        if (k & 1) res = res * base % mod;
        base = base * base % mod, k >>= 1;
    }
    return res;
}

template <typename Container>
void NTT(int T, Container &A, int type)
{
    static int Rev[N];
    for (int i = 0; i < T; i ++) Rev[i] = (Rev[i >> 1] >> 1) | (i & 1 ? T >> 1 : 0);
    for (int i = 0; i < T; i ++) if (i < Rev[i]) swap(A[i], A[Rev[i]]);

    for (int cur = 1; cur < T; cur <<= 1)
    {
        int base = ksm(type == 1 ? ord : ordinv, (mod - 1) / (cur << 1));

        static int comp[N];
        for (int i = comp[0] = 1; i < cur; i ++) comp[i] = comp[i - 1] * base % mod;

        for (int L = 0; L < T; L += cur << 1)
            for (int j = 0; j < cur; j ++)
            {
                int tempx = A[L | j], tempy = comp[j] * A[L | j | cur] % mod;
                A[L | j] = (tempx + tempy) % mod, A[L | j | cur] = (tempx - tempy + mod) % mod;
            }
    }

    if (type == 1) return;
    for (int i = 0, inv = ksm(T, mod - 2); i < T; i ++) A[i] = A[i] * inv % mod; 
}
void NTT(int T, int *A, int type)
{
    vector<int> temp(A, A + T);
    NTT(T, temp, type);
    copy(temp.begin(), temp.end(), A);
}

乘法函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
template<typename Container>
Container multi(int T, Container &F, Container &G)
{
    Container A(F), B(G), H(T);
    NTT(T, A, 1), NTT(T, B, 1);
    for (int i = 0; i < T; i ++) H[i] = A[i] * B[i] % mod;
    NTT(T, H, -1);
    return H;
}
vector<int> multi(int T, int *F, int *G)
{
    vector<int> A(F, F + T), B(G, G + T);
    return multi(T, A, B);
}

多项式乘法逆

P4238 【模板】多项式乘法逆 - 洛谷

给定一个多项式 F(x)F(x) ，请求出一个多项式 G(x)G(x)， 满足 F(x)×G(x)≡1(modxn)F(x) \times G(x) \equiv 1 \pmod{x^n}。系数对 998244353998244353 取模。

算法思路

做法很难想，考虑如果只有一项，此时 [x0]G(x)=[x0]F(x)−1[x^0]G(x) = [x^0]F(x)^{-1}，直接求出逆元即可。

假设我们现在已经求出了 F(x)F(x) 模 x⌈n2⌉x^{\lceil \frac{n}{2} \rceil} 的逆元 G′(x)G'(x)，现在对于模 xnx ^ n 求出 F(x)F(x) 的逆元 G(x)G(x)，数学表达为：

F×G′≡1(modx⌈n2⌉)F×G≡1(modxn)⇒F×G≡1(modx⌈n2⌉)F \times G' \equiv 1 \pmod {x ^ {\lceil \frac{n}{2} \rceil}} \\ F \times G \equiv 1 \pmod {x^n} \Rightarrow F \times G \equiv 1 \pmod {x^{\lceil \frac{n}{2} \rceil}}

两式相减再平方，有：

F×(G′−G)≡0(modx⌈n2⌉)(G′−G)≡G′2−2G′×G+G2≡0(modxn)F \times (G' - G) \equiv 0 \pmod {x ^ {\lceil \frac{n}{2} \rceil}} \\ (G' - G) \equiv G'^2 - 2G'\times G + G^2 \equiv 0 \pmod {x^n}

由于 F×G≡1(modxn)F \times G \equiv 1 \pmod {x^n}，两边同时乘 FF，有：

F×G′2−2G′+G≡0(modxn)⇒G≡2G′−F×G′2(modxn)F \times G'^2 - 2G' + G \equiv 0 \pmod {x ^ n} \Rightarrow \\ G \equiv 2G' - F \times G'^2 \pmod {x^n}

至此，我们就求出了在模 xnx ^ n 情况下的 GG 了，具体从 11 次项从下往上递推即可。

代码

记录详情 - 洛谷 | 计算机科学教育新生态

1
2
3
4
5
6
7
8
9
10
11
12
13
14
void getinv(int T, vector<int> &F, vector<int> &G)
{
    if (T == 1) return void(G[0] = ksm(F[0], mod - 2));

    getinv(T >> 1, F, G);
    vector<int> A(T << 1), B(T << 1), H(T << 1);
    for (int i = 0; i < T; i ++) A[i] = F[i], B[i] = G[i];
    NTT(T << 1, A, 1), NTT(T << 1, B, 1);

    for (int i = 0; i < T << 1; i ++) H[i] = (2 - A[i] * B[i] % mod + mod) % mod * B[i] % mod;
    NTT(T << 1, H, -1);

    for (int i = 0; i < T; i ++) G[i] = H[i];
}

时间复杂度分析：

T(n)=T(n2)+O(nlog⁡n)=∑k=0log⁡nO(2klog⁡2k)=O(nlog⁡n)\begin{aligned} T(n) &= T(\frac{n}{2}) + \mathcal{O}(n \log n) \\ &= \sum_{k = 0}^{\log n} \mathcal{O}(2^k \log 2^k) = \mathcal{O}(n \log n) \end{aligned}

因此多项式求逆的时间复杂度在 NTT 优化至 O(nlog⁡n)\mathcal{O}(n \log n)。

多项式除法

P4512 【模板】多项式除法 - 洛谷

给定一个 nn 次多项式 F(x)F(x) 和一个 mm 次多项式 G(x)G(x) ，请求出多项式 Q(x)Q(x), R(x)R(x)，满足以下条件：

  • Q(x)Q(x) 次数为 n−mn-m，R(x)R(x) 次数小于 mm
  • F(x)=Q(x)∗G(x)+R(x)F(x) = Q(x) * G(x) + R(x)

所有的运算在模 998244353998244353 意义下进行。

算法思路

同样十分巧妙，我们定义 FR(x)F_R(x) 是 F(x)F(x) 的系数置换：

FR(x):=xnF(1x),(n:=deg⁡F)F_R(x) := x ^ {n} F(\frac{1}{x}),(n := \deg F)

在代码中可以当作将多项式的系数数组翻转了，时间复杂度为 O(n)\mathcal{O}(n)。

将 1x\frac{1}{x} 带入等式中，得到：

F(x)=Q(x)∗G(x)+R(x)F(1x)=Q(1x)∗G(1x)+R(1x)xnF(1x)=xmG(1x)×xn−mQ(1x)+xm−1R(x)×xn−m+1F(x) = Q(x) * G(x) + R(x) \\ F(\frac{1}{x}) = Q(\frac{1}{x}) * G(\frac{1}{x}) + R(\frac{1}{x}) \\ x ^ n F(\frac{1}{x}) = x ^ m G(\frac{1}{x}) \times x^{n - m} Q(\frac{1}{x}) + x^{m - 1}R(x) \times x^{n - m + 1}

根据之前的定义，得到：

FR(x)=GR(x)×QR(x)+RR(x)×xn−m+1F_R(x) = G_R(x) \times Q_R(x) + R_R(x) \times x^{n - m + 1}

在模 xn−m+1x^{n - m + 1} 的意义下，不难发现：

FR(x)≡GR(x)×QR(x)+RR(x)×xn−m+1(modxn−m+1)FR(x)≡GR(x)×QR(x)(modxn−m+1)QR(x)≡FR(x)×GR(x)−1(modxn−m+1)F_R(x) \equiv G_R(x) \times Q_R(x) + R_R(x) \times x ^ {n - m + 1} \pmod {x ^ {n - m + 1}} \\ F_R(x) \equiv G_R(x) \times Q_R(x) \pmod {x ^ {n - m + 1}} \\ Q_R(x) \equiv F_R(x) \times G_R(x)^{-1} \pmod {x ^ {n - m + 1}}

因此我们只需要求出在模 xn−m+1x^{n - m + 1} 下的 FR(x)F_R(x) 和 GR(x)−1G_R(x)^{-1} 即可求出多项式的商的系数置换 QR(x)Q_R(x)，对于余数：

R(x)=F(x)−G(x)×Q(x)R(x) = F(x) - G(x) \times Q(x)

也可以求出。

总的时间复杂度是十分大常数的 O(nlog⁡n)\mathcal{O}(n \log n)，总共进行了六次 NTT 还有一次的多项式求逆。

代码

记录详情 - 洛谷 | 计算机科学教育新生态

注意到这里的 NTT 范围最高会到 2n2n，因此边界要足够大。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// while (T <= n << 1) T <<= 1;
void divide(int T, vector<int> &F, vector<int> &G, vector<int> &Q, vector<int> &R)
{
    vector<int> A(F), B(G);
    reverse(A.begin(), A.begin() + n + 1), reverse(B.begin(), B.begin() + m + 1);
    getinv(T, B, Q), B = Q;
    for (int i = n + 1; i < T; i ++) B[i] = 0;
    
    NTT(T, A, 1), NTT(T, B, 1);
    for (int i = 0; i < T; i ++) Q[i] = A[i] * B[i] % mod;
    NTT(T, Q, -1);
    for (int i = n - m + 1; i < T; i ++) Q[i] = 0;
    reverse(Q.begin(), Q.begin() + n - m + 1);

    B = G, R = Q;
    NTT(T, B, 1), NTT(T, R, 1);
    for (int i = 0; i < T; i ++) R[i] = R[i] * B[i] % mod;
    NTT(T, R, -1);
    for (int i = 0; i < m; i ++) R[i] = (F[i] - R[i] + mod) % mod;
}