多项式科技初步

写在前面

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

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

如果 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\times 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 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 ] 多项式除法](https://www.luogu.org/problemnew/show/P4512) 作为测试题，不再区分多项式求逆部分的实现方式。

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

注意求 $D(x)$ 部分的那次卷积是在 $\bmod{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$ 的贡献一起做。

考虑对 $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);
}

posted @ 2019-03-02 11:33 SGCollin 阅读(391) 评论(0) 编辑收藏举报

刷新页面返回顶部

SGColin's Space

Dreaming in the memory is not as good as waiting for the paradise in the hell.

多项式科技初步

写在前面

参考资料

多项式乘法

问题描述

解决方法

代码实现

多项式求逆

问题描述

解决方法

代码实现

多项式开根

问题描述

解决方法

代码实现

分治 FFT

问题描述

解决方法

分治求解

多项式求逆

代码实现

多项式求导和积分

问题描述

解决方法

代码实现

多项式 ln

问题描述

解决方法

代码实现

公告