# 多项式科技初步

## 参考资料

Picks 的博客 Picks's Blog

Miskcoo 的博客 Miskcoo's Space

## 多项式乘法

### 问题描述

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

$C(x)=\sum_{i=0}^{n} c_i x^i\ ,\ c_i=\sum_{k=0}^{i} a_k\times b_{i-k}$

### 代码实现

• 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;
}
}


## 多项式求逆

### 问题描述

$A(x) * B(x) \equiv 1 \pmod{ x^{n+1}}$

### 解决方法

$A(x)B(x)\equiv 1\pmod{x^n}$

$A(x)B(x)\equiv 1\pmod{x^{\lceil\frac{n}{2}\rceil}}$

$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}$

$B(x)\equiv2B'(x)-A(x)B'^2(x)\pmod{x^n}$

$B_i=2B'_i-A_iB'^2_i=B'_i(2-A_iB')$

$T(n) = T(\frac{n}{2}) + \mathcal O(n \log n) = \mathcal O(n \log n)$

### 代码实现

• 递归版本，使用 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;
}
}


## 多项式开根

### 问题描述

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

### 解决方法

$B^2(x)\equiv A(x)\pmod{x^n}\ \Rightarrow\ B^2(x)\equiv A(x)\pmod{x^{\lceil \frac{n}2\rceil}}$

$D^2(x)\equiv A(x)\pmod{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}$

$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)$

### 代码实现

• 递归版本，使用 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;
}
}


---

## 多项式除法和取模

### 问题描述

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

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

### 解决方法

$$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}$$

$$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)$$

$$A^R(x) \equiv D^R(x)B^R(x) \pmod{x^{n-m+1}}$$

$$\frac{A^R(x)}{B^R(x)} \equiv D^R(x) \pmod{x^{n-m+1}}$$

### 代码实现

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

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

### 问题描述

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

### 解决方法

#### 分治求解

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

$T(n) = 2T(\frac{n}{2}) + \mathcal O(n \log n) = \mathcal O(n \log^2 n)$

#### 多项式求逆

$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$ 数组可以看作是 $\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$

### 代码实现

• 分治版本，使用 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);
}


## 多项式求导和积分

### 问题描述

$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}$

### 代码实现

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

### 问题描述

$B(x) \equiv \ln A(x)\pmod{x^{n+1}}$

### 解决方法

$F(x)=\ln x$，则 $B(x)=F(A(x))$

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

$B'(x)=F'(A(x))A'(x)=\frac{A'(x)}{A(x)}$

### 代码实现

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  阅读(235)  评论(0编辑  收藏