多项式科技初步

时间 2019-12-11

标签多项式科技初步繁體版

原文原文链接

写在前面

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

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

若是 MathJax 加载不出来或加载有误，请您多刷新几回。函数

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

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

参考资料

Picks 的博客 Picks's Blogui

Miskcoo 的博客 Miskcoo's Spacespa

多项式乘法

问题描述

给定两个多项式 $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\times b_{i-k} \]code

解决方法

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

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

代码实现

使用 [ 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 c\pmod{x}$ ，那么 $A^{-1}(x)$ 即为 $c^{-1}$。

对于 $n>1$ 的时候，设 $B(x)=A^{-1}(x)$ ，有
\[ A(x)B(x)\equiv 1\pmod{x^n} \]
由于此时模 $x^n$ 至关于只保留多项式前 $n$ 项，因此该同余式在模 $x^k,0\le k\le n$ 时都成立
\[ A(x)B(x)\equiv 1\pmod{x^{\lceil\frac{n}{2}\rceil}} \]
假设在 $\pmod {x^{\lceil \frac n2\rceil}}$ 意义下 $A(x)$ 的逆元是 $B′(x)$ 而且咱们已经求出，那么
\[ A(x)B'(x) \equiv 1 \pmod {x^{\lceil \frac{n}{2} \rceil}} \]
两式相减，得
\[ B(x)-B'(x) \equiv 0 \pmod {x^{\lceil \frac{n}{2} \rceil}} \]
两边平方，得
\[ B^2(x)-2B(x)B'(x)+B'^2(x)\equiv 0\pmod{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)$，求一个在 $\bmod{x^{n+1}}$ 意义下的多项式 $B(x)$，使得

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

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

解决方法

一样采用倍增的思想。

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

对于 $n>1$ 的时候，一样根据上一题的结论，咱们能够把问题范围缩小到 $\bmod {x^{\lceil \frac{n}2\rceil}}$ ，有
\[ B^2(x)\equiv A(x)\pmod{x^n}\ \Rightarrow\ B^2(x)\equiv A(x)\pmod{x^{\lceil \frac{n}2\rceil}} \]
不妨设咱们已经求出来了 $\bmod{x^{\lceil \frac{n}2\rceil}} $ 意义下的根 $D(x)$，即
\[ D^2(x)\equiv A(x)\pmod{x^{\lceil \frac{n}2\rceil}} \]
所以 $B(x)$ 与 $D(x)$ 在模 $x^{\lceil \frac{n}2\rceil}$ 意义下同余，移项得
\[ B(x)-D(x)\equiv 0\pmod{x^{\lceil \frac{n}2\rceil}} \]
两侧平方，得
\[ B^2(x)+D^2(x)-2B(x)D(x)\equiv 0\pmod{x^n} \]
模数能平方的缘由与上一题相同。

咱们知道$\bmod{x^n}$ 时 $B^2(x)$ 即为 $A(x)$ ，所以
\[ A(x)+D^2(x)-2B(x)D(x)\equiv 0\pmod{x^n} \]
移项，得
\[ B(x)\equiv \frac{D^2(x) + A(x)}{2D(x)}\equiv \bigg(D(x) +\frac{A(x)}{D(x)}\bigg)\times 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^n\sum_{i=0}^na_i\frac{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$ 。

因此咱们能够把多项式运算在 $\bmod{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 ] 多项式除法做为测试题，再也不区分多项式求逆部分的实现方式。

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

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

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$ 的贡献一块儿作。

考虑对 $x\in[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^i\sum_{j}f[i-j]g[j]=F(x)-f[0] \]
后一个等式成立的缘由是，注意到后一个求和就是 $f[i]$ 的形式，因此只有 $f[0]$ 没有被计数

求 $f$ 数组能够看做是 $\bmod{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(n\log 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)$，求一个 $\bmod{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);
}