多项式科技初步

写在前面

为了体现简洁，在每一部分只会放关键代码。

代码具有一定的通用性，保证代码之间函数的调用是合法的。

如果 MathJax 加载不出来或加载有误，请您多刷新几次。

总结可能会比较长，可以点右边的小火箭回到目录。

有什么问题请联系我，万分感谢。

---

参考资料

Picks 的博客 Picks's Blog

Miskcoo 的博客 Miskcoo's Space

---

多项式乘法

问题描述

给定两个多项式 $A(x) $ 和 (B(x)) ：

[A(x)=sum_{i=0}^{n}a_ix^i , B(x)=sum_{i=0}^{n} b_i x^i ]

求卷积 (C(x) =A(x) * B(x)) ，满足

[C(x)=sum_{i=0}^{n} c_i x^i , c_i=sum_{k=0}^{i} a_k imes b_{i-k} ]

解决方法

这部分上面提到的神仙们讲解的都十分详细。

我的小结限于篇幅，就直接挂 pdf 版本的连接了：快速傅里叶变换初步

代码实现

使用 [ UOJ 34 ] 多项式乘法 作为测试题。

• Fast Fourier Transform，使用的是 3 次变换的最基本写法，用时 362 ms

  inline void FFT(Complex *f, int len, int o) {
      for (int i = 0; i < len; ++i)
        if (rev[i] > i) swap(f[i], f[rev[i]]);
      for (int i = 1; i < len; i <<= 1) {
        Complex wn = Complex(cos(PI / i) , o * sin(PI / i));
        for (int j = 0; j < len; j += (i << 1)) {
          Complex w = Complex(1, 0), x, y;
          for (int k = 0; k < i; ++k, w = w * wn) {
            x = f[j + k]; y = w * f[i + j + k];
            f[i + j + k] = x - y; f[j + k] = x + y;
          }
        }
      }
      if (o == -1) for (int i = 0; i < len; ++i) f[i].x /= len;
  }
  

• Fast Number-Theoretic Transform，模数为 998244353，用时 403 ms

  inline void NTT(int *f, int len, int o) {
    for (int i = 1; i < len; ++i)
      if (i > rev[i]) swap(f[i], f[rev[i]]);
    for (int i = 1; i < len; i <<= 1) {
      int wn = qpow(3, (mod - 1) / (i << 1));
      if (o == -1) wn = qpow(wn, mod - 2);
      for (int j = 0; j < len; j += (i << 1)) {
        int w = 1, x, y;
        for (int k = 0; k < i; ++k, w = 1ll * w * wn % mod) {
          x = f[j + k]; y = 1ll * w * f[i + j + k] % mod;
          f[j + k] = mo(x + y); f[i + j + k] = mo(x - y + mod);
        }
      }
    }
    if (o == -1) {
      int invl = qpow(len, mod - 2);
      for (int i = 0; i < len; ++i) f[i] = 1ll * f[i] * invl % mod;
    }
  }
  

---

多项式求逆

问题描述

给定一个 (n) 次多项式 (A(x)) ，求出一个多项式 (B(x))， 满足

[A(x) * B(x) equiv 1 pmod{ x^{n+1}} ]

系数对 998244353 取模。

解决方法

采用倍增的思想。

考虑只有常数项的时候，(A(x)equiv cpmod{x}) ，那么 (A^{-1}(x)) 即为 (c^{-1})。

对于 (n>1) 的时候，设 (B(x)=A^{-1}(x)) ，有

[A(x)B(x)equiv 1pmod{x^n} ]

因为此时模 (x^n) 相当于只保留多项式前 (n) 项，所以该同余式在模 (x^k,0le kle n) 时都成立

[A(x)B(x)equiv 1pmod{x^{lceilfrac{n}{2} ceil}} ]

假设在 (pmod {x^{lceil frac n2 ceil}}) 意义下 (A(x)) 的逆元是 (B′(x)) 并且我们已经求出，那么

[A(x)B'(x) equiv 1 pmod {x^{lceil frac{n}{2} ceil}} ]

两式相减，得

[B(x)-B'(x) equiv 0 pmod {x^{lceil frac{n}{2} ceil}} ]

两边平方，得

[B^2(x)-2B(x)B'(x)+B'^2(x)equiv 0pmod{x^n} ]

模数平方的合法性在于，大于 (n) 的系数卷积中每组乘法必然有一项下标小于 (n) ，因此乘起来必然为 (0) 。

两侧同称 (A(x)) ，整理得

[B(x)equiv2B'(x)-A(x)B'^2(x)pmod{x^n} ]

因此只需将 (B'(x)) 和 (A(x)) 在模 (x^n) 意义下的插值求出，有

[B_i=2B'_i-A_iB'^2_i=B'_i(2-A_iB') ]

因此一遍 NTT 就可以由 (B'(x)) 求出 (B(x)) 了。

总的时间复杂度为

[T(n) = T(frac{n}{2}) + mathcal O(n log n) = mathcal O(n log n) ]

由此过程也可以得到一个结论：一个多项式有没有逆元完全取决于其常数项是否有逆元。

代码实现

使用 [ Luogu P4238 ] 多项式求逆 作为测试题。

• 递归版本，使用 O2 优化，用时 562 ms

  inline void Inv(int *a, int *b, int n) {
    if (n == 1) {b[0] = qpow(a[0], mod - 2); return;}
    Inv(a, b, (n + 1) >> 1);
    int len = Rev(n << 1);
    for (int i = 0; i < n; ++i) tmp[i] = a[i];
    for (int i = n; i < len; ++i) b[i] = tmp[i] = 0;
    NTT(b, len, 1); NTT(tmp, len, 1);
    for (int i = 0; i < len; ++i)
      b[i] = (2ll - 1ll * tmp[i] * b[i] % mod + mod) * b[i] % mod;
    NTT(b, len, -1);
    for (int i = 0; i < len; ++i) tmp[i] = 0;
    for (int i = n; i < len; ++i) b[i] = 0;
  }
  

• 迭代版本，使用 O2 优化，用时 570 ms

  inline void Inv(int *a, int *b, int n) {
    b[0] = qpow(a[0], mod - 2);
    int len;
    for (int l = 1; l < (n << 1); l <<= 1) {
      len = Rev(l << 1);
      for (int i = 0; i < l; ++i) tmp[i] = a[i];
      NTT(b, len, 1); NTT(tmp, len, 1);
      for (int i = 0; i < len; ++i)
        b[i] = (2ll - 1ll * tmp[i] * b[i] % mod + mod) * b[i] % mod;
      NTT(b, len, -1);
      for (int i = l; i < len; ++i) b[i] = 0;
    }
  }
  

---

多项式开根

问题描述

给定一个 (n) 次多项式 (A(x))，求一个在 (mod{x^{n+1}}) 意义下的多项式 (B(x))，使得

[B^2(x) equiv A(x) pmod{x^{n+1}} ]

多项式的系数在模 998244353 意义下进行运算，保证常数项 (a_0=1)。

解决方法

同样采用倍增的思想。

考虑只有常数项的时候，(A(x)equiv cpmod x) ，那么 (sqrt {A(x)}) 即为 (sqrt c equiv 1pmod{x}) （二次剩余）。

对于 (n>1​) 的时候，同样根据上一题的结论，我们可以把问题范围缩小到 (mod {x^{lceil frac{n}2 ceil}}​) ，有

[B^2(x)equiv A(x)pmod{x^n} Rightarrow B^2(x)equiv A(x)pmod{x^{lceil frac{n}2 ceil}} ]

不妨设我们已经求出来了 (mod{x^{lceil frac{n}2 ceil}} ​) 意义下的根 (D(x)​)，即

[D^2(x)equiv A(x)pmod{x^{lceil frac{n}2 ceil}} ]

因此 (B(x)​) 与 (D(x)​) 在模 (x^{lceil frac{n}2 ceil}​) 意义下同余，移项得

[B(x)-D(x)equiv 0pmod{x^{lceil frac{n}2 ceil}} ]

两侧平方，得

[B^2(x)+D^2(x)-2B(x)D(x)equiv 0pmod{x^n} ]

模数能平方的原因与上一题相同。

我们知道(mod{x^n}) 时 (B^2(x)) 即为 (A(x)) ，因此

[A(x)+D^2(x)-2B(x)D(x)equiv 0pmod{x^n} ]

移项，得

[B(x)equiv frac{D^2(x) + A(x)}{2D(x)}equiv igg(D(x) +frac{A(x)}{D(x)}igg) imes 2^{-1}pmod{x^n} ]

因此倍增时进行多项式求逆即可，总的时间复杂度为

[T(n) = T(frac{n}{2}) + mathcal O(n log n) = mathcal O(n log n) ]

代码实现

使用 [ Luogu P5205 ] 多项式开根 作为测试题，多项式求逆部分均采用递归版本。

• 递归版本，使用 O2 优化，用时 3081 ms

  inline void Sqrt(int *a, int *b, int n) {
    if (n == 1) {b[0] = 1; return;}
    Sqrt(a, b, (n + 1) >> 1);
    Inv(b, b0, n);
    int len = Rev(n << 1);
    for (int i = 0; i < n; ++i) a0[i] = a[i];
    for (int i = n; i < len; ++i) a0[i] = 0;
    NTT(a0, len, 1); NTT(b0, len, 1);
    for (int i = 0; i < len; ++i) a0[i] = 1ll * a0[i] * b0[i] % mod;
    NTT(a0, len, -1);
    for (int i = 0; i < n; ++i) b[i] = 1ll * (b[i] + a0[i]) % mod * inv2 % mod;
    for (int i = n; i < len; ++i) b[i] = 0;
  }
  

• 迭代版本，使用 O2 优化，用时 3104 ms

inline void Sqrt(int *a, int *b, int n) {
b[0] = 1;
int len;
for (int l = 1; l < (n << 1); l <<= 1) {
Inv(b, b0, l);
len = Rev(l << 1);
for (int i = 0; i < l; ++i) a0[i] = a[i];
for (int i = l; i < len; ++i) a0[i] = 0;
NTT(a0, len, 1); NTT(b0, len, 1);
for (int i = 0; i < len; ++i) a0[i] = 1ll * a0[i] * b0[i] % mod;
NTT(a0, len, -1);
for (int i = 0; i < l; ++i) b[i] = 1ll * (b[i] + a0[i]) % mod * inv2 % mod;
for (int i = l; i < len; ++i) b[i] = 0;
}
}

---

## 多项式除法和取模

### 问题描述

给定一个 $n$ 次多项式 $A(x)$ 和一个 $m$ 次多项式 $B(x)$，求出多项式 $D(x)$, $R(x)$，满足

$$
A(x) = D(x)B(x) + R(x)
$$

$D(x)$ 次数为 $n-m$，$R(x)$ 次数小于 $m$ ，所有的运算在模 998244353 意义下进行。

### 解决方法

注意到带着 $R(x)​$ 在这里很麻烦，前人们想到了一个神奇的解决办法。

设 $A^R(x)=x^nA(frac{1}{x})$ ，我们将右侧展开：
$$
A^R(x)=x^nA(frac{1}{x})=x^nsum_{i=0}^na_ifrac{1}{x^i}=sum_{i=0}^na_ix^{n-i}
$$
所以 $A^R(x)$ 就是将 $A(x)$ 的**系数反转** 。

我们将所求的等式中 $x$ 全部换成 $frac{1}{x}$ ，然后两侧同乘 $x^n$ ：
$$
x^n A(frac{1}{x}) = x^{n - m}D(frac{1}{x}) x^mB(frac{1}{x}) + x^{n - m + 1}x^{m - 1}R(frac{1}{x})
$$
$$
A^R(x) = D^R(x)B^R(x) + x^{n - m + 1}R^R(x)
$$

注意到 $D^R(x)$ 最高次反转后不变，依然为 $n-m$。

而右侧的 $R^R(x)$ 因为前面有 $x^{n-m+1}$ 所以最低次为 $n-m+1$ 。

所以我们可以把多项式运算在 $mod{x^{n-m+1}}$ 意义下进行，这样 $R(x)$ 就消失了：
$$
A^R(x) equiv D^R(x)B^R(x) pmod{x^{n-m+1}}
$$
因此就可以得到 $D^R(x)$ 的解法：
$$
frac{A^R(x)}{B^R(x)} equiv D^R(x) pmod{x^{n-m+1}}
$$
一个多项式求逆就可以求出 $D^R(x)$了，再将 $D^R(x)$ 进行反转就得到了答案。

将求出的 $D(x)$ 回代，再进行一次减法即可求出 $R(x)$。

复杂度与多项式求逆同阶，为 $mathcal O(n log n)​$ 。

### 代码实现

使用 [[ Luogu P4512 ] 多项式除法](https://www.luogu.org/problemnew/show/P4512) 作为测试题，不再区分多项式求逆部分的实现方式。

[多项式求逆递部分归版本](https://blog.gyx.me/code/template/polynomial/div.cpp)，使用 O2 优化，用时 710 ms

注意求 $D(x)$ 部分的那次卷积是在 $mod{x^{n-m+1}}$意义下的，所以 $A^R(x)$ 和 $B^R(x)$ 中高于 $n-m+1$ 的项需要清空（因为 FFT 卷积过程中高次系数也会对低次系数造成影响）



```C++
inline void Div(int *a, int *b, int n, int m) {
for (int i = 0; i <= n; ++i) ar[i] = a[n - i];
for (int i = 0; i <= m; ++i) br[i] = b[m - i];
for (int i = n - m + 2; i <= n; ++i) ar[i] = 0;
for (int i = n - m + 2; i <= m; ++i) br[i] = 0;
Inv(br, invb, n - m + 1);
int len = Rev((n - m + 1) << 1);
NTT(ar, len, 1); NTT(invb, len, 1);
for (int i = 0; i < len; ++i) ar[i] = 1ll * ar[i] * invb[i] % mod;
NTT(ar, len, -1);
for (int i = 0; i <= n - m; ++i) tmp[i] = d[i] = ar[n - m - i];
len = Rev(n << 1);
for (int i = n - m + 1; i <= len; ++i) tmp[i] = 0;
NTT(b, len, 1); NTT(tmp, len, 1);
for (int i = 0; i < len; ++i) b[i] = 1ll * b[i] * tmp[i] % mod;
NTT(b, len, -1);
for (int i = 0; i <= m; ++i) r[i] = mo(a[i] - b[i] + mod);
}

---

分治 FFT

问题描述

给定长度为 (n-1) 的数组 (g[1],dots,g[n-1])，求长度为 (n) 的数组 (f[0],dots,f[n-1])，其中

[f[i]=sum_{j=1}^if[i-j]g[j] ]

边界为 (f[0]=1)，运算在模 998244353 下进行。

解决方法

分治求解

暴力做是 (mathcal O(n^2)) 的，考虑将相同的转移一起做以达到优化的目的。

考虑使用类似 CDQ 分治的思想，每次我们求出 ([L, mid]) 范围内的 (f) 数组之后，把这部分 (f) 对 ([mid+1, R]) 范围内 (f) 的贡献一起做。

考虑对 (xin[mid + 1, R]) 的 (f[x]) 的贡献 (w_x)，有

[w_x=sum_{i=l}^{mid} f[i]g[x-i] ]

因此 (w) 数组可以卷积求了，注意求 (w_x) 时后半段的 (f) 需要认为是 (0)，否则就存在右区间内部的贡献了。

总的时间复杂度为

[T(n) = 2T(frac{n}{2}) + mathcal O(n log n) = mathcal O(n log^2 n) ]

多项式求逆

一阶分治 FFT 是可以看作卷积处理的。

不妨设将数组看成多项式，有

[F(x)=sum_{i=0}^{n-1}f[i]x^i , G(x)=sum_{i=0}^{n-1}g[i]x^i ]

将两个多项式卷积，有

[F(x)G(x)=sum_{i=0}^{infty}x^isum_{j}f[i-j]g[j]=F(x)-f[0] ]

后一个等式成立的原因是，注意到后一个求和就是 (f[i]) 的形式，所以只有 (f[0]) 没有被计数

求 (f) 数组可以看作是 (mod{x^n}) 意义下进行的，因此有

[F(x)G(x) equiv F(x)-f[0] pmod x^n Rightarrow F(x) equiv frac{f_0}{1-G(x)} equiv (1-G(x))^{-1} pmod x^n ]

于是一遍多项式求逆就可以求出来了，复杂度为 (mathcal O(nlog n))

代码实现

使用 [ Luogu P4721 ] 分治 FFT 作为测试题，多项式求逆采用递归版本。

• 分治版本，使用 O2 优化，用时 974 ms

• 多项式求逆版本，使用 O2 优化，用时 261 ms

  inline void solve(int *a, int n) {
    a[0] = 1;
    for (int i = 1; i < n; ++i) a[i] = mo(mod - a[i]);
    Inv(a, b, n);
    for (int i = 0; i < n; ++i) print(b[i], 0);
  }
  

---

多项式求导和积分

问题描述

给定一个 (n) 次多项式 (A(x)) ，求一个 (n-1) 次多项式 (B(x)) ，和一个 (n+1) 次多项式 (C(x))，满足

[B(x)=A'(x) , C(x)=int A(x) ]

解决方法

直接按照定义做即可，有

[B(x)=sum_{i=1}^{n}ia_ix^{i-1} , C(x)=sum_{i=1}^{n}frac{a_ix^{i+1}}{i+1} ]

我们一般认为 (C(x)) 的常数项为 (0) ，复杂度显然为 (mathcal O(n))。

代码实现

inline void Der(int *a, int n) {
  for (int i = 1; i < n; ++i) a[i - 1] = 1ll * i * a[i] % mod;
  a[n - 1] = 0;
}

inline void Int(int *a, int n) {
  for (int i = n; i; --i) a[i] = 1ll * a[i - 1] * qpow(i, mod - 2) % mod;
  a[0] = 0;
}

---

多项式 ln

问题描述

给出 (n) 次多项式 (A(x))，求一个 (mod{x^{n+1}}) 下的多项式 (B(x))，满足

[B(x) equiv ln A(x)pmod{x^{n+1}} ]

所有运算在模 998244353 下进行。

解决方法

设 (F(x)=ln x)，则 (B(x)=F(A(x))) 。

对 (B(x)) 求导，根据链式法则，有

[B'(x)=F'(A(x))A'(x)=frac{A'(x)}{A(x)} ]

因此对 (A(x)) 分别进行求导和求逆，卷积即可求出 (B'(x)) ，再对其进行积分即可。

复杂度与多项式求逆同阶，为 (mathcal O(n log n)) 。

代码实现

使用 [ Luogu P4725 ] 多项式对数函数 作为测试题，不再区分多项式求逆部分的实现方式。

多项式求逆递部分归版本，使用 O2 优化，用时 682 ms

inline void Ln(int *a, int *b, int n) {
  Inv(a, b, n); Der(a, n);
  int len = Rev(n << 1);
  NTT(a, len, 1); NTT(b, len, 1);
  for (int i = 0; i < len; ++i) a[i] = 1ll * a[i] * b[i] % mod;
  NTT(a, len, -1); Int(a, n);
}