重学多项式
目录更好的阅读体验戳此进入
写在前面
关于多项式曾经也写过一篇 Blog，同时亦是我写过的第一篇 Blog，当然可能实际上也没有很理解多项式，所以写的或许比较乱，也可以看一下 FFT & NTT - 快速傅里叶变换 & 快速数论变换。
关于这篇Blog
有些地方可能是一些个人理解，不保证正确性。
包括但不限于：单位根、原根、FFT、NTT、多项式取模，多项式求逆，多项式。。。
单位根
定义
对于 $ \omega^n = 1 (\omega \neq 1)$，称 $ \omega $ 为 $ n $ 次单位根，其可以为模意义下的或复数意义下的。
复数意义下的的求法
考虑将单位圆 $ n $ 等分，并取这 $ n $ 个点对应的素数，从 $ x $ 轴，即 $ (1, 0) $ 开始取，逆时针从 $ 0 $ 开始编号，对于第 $ k $ 个复数记作 $ \omega_n^k $。所以以下结论显然：
对于幅角为 $ 0 $ 的最小复数 $ \omega_n^1 $，由复数乘法规则可知 $ (\omega_n1)k = \omega_n^k $，所以 $ (\omega_n1)n = \omega_n^n = \omega_n^0 = 1 $，则显然 $ \omega_n^1 $ 为复数意义下的 $ n $ 次单位根。
则显然对于 $ n $ 次单位根的 $ k $ 次方，用高中知识转换为三角函数即可，即：
\[\omega_n^k = \cos(\dfrac{2\pi}{n} \times k) + \sin(\dfrac{2\pi}{n} \times k)i \]性质
复数意义下的单位根，易证有如下性质：
性质一：
\[\omega_n^k = \omega_{an}^{ak} \]性质二：
\[\omega_n^{k + \tfrac{n}{2}} = -\omega_n^k \]等比数列求和公式
没什么可说的，公比为 $ q $，则有下式。详细证明。
\[S_n = a_1 \dfrac{1-q^n}{1-q} \]FFT
本质
在我的理解里，FFT 的本质就是将多项式的系数表示法转换为点值表示法，然后进行快速处理，再转换为一般的系数表示法。
因为显然朴素的多项式乘法为 $ O(n^2) $ 的，而点值表示法的多项式乘法是 $ O(n) $ 的，所以我们需要将朴素的 $ O(n^2) $ 的系数与点值转换过程优化为 $ O(n \log n) $ 的。
目的
显然目的就是求两个多项式的卷积，我们定义 $ A(x), B(x), C(x) $ 为多项式，$ a(x), b(x), c(x) $ 为多项式的对应系数。即 $ C(x) = A(x) \ast B(x) $，或者表示为：
\[c(i) = \sum_{j = 0}^ia(j) \times b(i - j) \]离散傅里叶变换
考虑我们前文提过的单位根，首先我们将 $ \omega_n^0, \omega_n1,\omega_n2, \cdots, \omega_n^{n - 1} $ 代入多项式的点值表示称作离散傅里叶变换。
关于为什么一定要选择这些值，在我以前写的里面的推式子有些繁琐，当然也可以证明。这里参考一下网上的另一种更简洁的方法。
首先我们记 $ (y_0, y_1, y_2, \cdots, y_{n - 1}) $ 为多项式 $ A(x) = \sum_{i = 0}^{n - 1}a_ix^i $ 的离散傅里叶变换。然后令多项式 $ B(x) = \sum_{i = 0}^{n - 1}y_ix^i $，然后这里我们用所有单位根 $ k $ 次方的倒数，也就是对应的共轭复数，带入 $ B(x) $，即带入 $ \omega_n^0, \omega_n{-1},\omega_n{-2}, \cdots, \omega_n^{-n + 1} $，得到一个新的 “离散傅里叶变换”，记作 $ (z_0, z_1, z_2,.\cdots, z_{n - 1}) $，则有：
\[\begin{aligned} z_k &= \sum_{i = 0}^{n - 1}y_i(\omega_n^{-k})^i \\ &= \sum_{i = 0}^{n - 1} (\sum_{j = 0}^{n - 1}a_j(\omega_n^i)^j) (\omega_n^{-k})^i \\ &= \sum_{j = 0}^{n - 1} a_j (\sum_{i = 0}^{n - 1}(\omega_n^{j - k})^i) \end{aligned} \]然后对于 $ \sum_{i = 0}^{n - 1}(\omega_n^{j - k})^i $，显然 $ j = k $ 的时候该式为 $ n $，反之通过等比数列求和公式易证其为 $ 0 $，则显然有 $ z_k = n \times a_k $，所以有：
\[a_i = \dfrac{z_i}{n} \]所以不难发现，我们转为系数表示法之后，即 DFT 过程之后，对于得到的离散傅里叶变换结果，再次当作一个新的多项式的系数以单位根倒数再跑一遍转换系数表示法的过程，即 IDFT，然后将每个数的实数部分除以 $ n $，最终得到的就是我们想要的结果
所以现在我们就需要优化 $ O(n^2) $ 的转化过程了。
快速傅里叶变换
首先仍考虑一个一般的多项式：
\[A(x) = \sum_{i = 0}^{n - 1}a_ix^i \]我们按照下标的奇偶性分类，则有如下两个多项式：
\[\begin{aligned} A_1(x) = a_0 + a_2x + \cdots + a_{n - 2}x^{\tfrac{n}{2} - 1} \\ A_2(x) = a_1 + a_3x + \cdots + a_{n - 1}x^{\tfrac{n}{2} - 1} \end{aligned} \]然后不难发现有：
\[A(x) = A_1(x^2) + xA_2(x^2) \]然后考虑对于 $ k \lt \dfrac{n}{2} $ 的部分，带入 $ \omega_n^k $，有：
\[\begin{aligned} A(\omega_n^k) &= A_1(\omega_n^{2k}) + \omega_n^kA_2(\omega_n^{2k}) \\ &= A_1(\omega_{\tfrac{n}{2}}^{k}) + \omega_n^kA_2(\omega_{\tfrac{n}{2}}^{k}) \end{aligned} \]然后对于 $ \ge \dfrac{n}{2} $ 的部分，考虑带入 $ \omega_n^{k + \tfrac{n}{2}} $：
\[\begin{aligned} A(\omega_n^{k + \tfrac{n}{2}}) &= A_1(\omega_n^{2k + n}) + \omega_n^{k + \tfrac{n}{2}}A_2(\omega_n^{2k + n}) \\ &= A_1(\omega_{\tfrac{n}{2}}^{k} \omega_{n}^{n}) + \omega_n^{k + \tfrac{n}{2}}A_2(\omega_{\tfrac{n}{2}}^{k} \omega_{n}^{n}) \\ &= A_1(\omega_{\tfrac{n}{2}}^{k}) - \omega_n^{k}A_2(\omega_{\tfrac{n}{2}}^{k}) \end{aligned} \]此时发现问题规模减半，合并是 $ O(n) $ 的，即：
\[T(n) = T(\dfrac{n}{2}) + O(n) \]则由主定理可知复杂度为 $ O(n \log n) $。
Tips：同时需要注意，上文即下文所用的长度 $ n $ 若无特殊限制则均应保证满足 $ n = 2^t $，若不满足可在最后添加 $ 0 $ 直至满足。
优化
这样虽然就可以实现 FFT 了，但是它是不优雅的，我们可以将递归改为非递归实现。
考虑 Cooley-Tukey 算法，首先有以下递归的例子：（图片来自 一小时学会快速傅里叶变换（Fast Fourier Transform））
不难发现，每个位置变化的数即为将其二进制下的所有数反转，如 $ 1 $ 从 $ 001 $ 变为 $ 110 $。于是可以据此进行非递归优化。
首先朴素的模拟这个反转的过程应该是 $ O(n \log n) $ 的，不够优秀，仍存在如下写法：
int pos[len + 10];
memset(pos, 0, sizeof(pos));
for(int i = 0; i < len; ++i){
pos[i] = pos[i >> 1] >> 1;
if(i & 1)pos[i] |= len >> 1;
}
for(int i = 0; i < len; ++i)if(i < pos[i])swap(pol[i], pol[pos[i]]);
不难发现这个的复杂度是线性的，严谨证明略，可以尝试通过举个例子来理解：
假设我们有一个二进制数 $ 0101110 $ 我们想要对其进行 Reverse，因为我们要进行递推，所以需要分解为子问题，可以考虑将其右移一位，即变为 $ 0010111 $，然后 Reverse，变为 $ 1110100 $，在对比我们想要求得的 $ 0111010 $，发现前者去掉最后一位，后者忽略第一位是完全相同的，那么将前者右移一位后在考虑最前面一位是 $ 1 $ 还是 $ 0 $ 即可。
此时则可以在过程中不进行递归，而从下至上地取合并，不难想到枚举长度后再枚举每个起点即可。
对于朴素的 FFT 提供如下实现：（以前写的代码，可能实现不是很精细）
对应题目：P3803 【模板】多项式乘法（FFT）。
Code
#define _USE_MATH_DEFINES
#include <bits/stdc++.h>
#include <mmintrin.h>
#define PI M_PI
#define E M_E
#define DFT true
#define IDFT false
#define eps 1e-6
#define comp complex < double >
/******************************
abbr
pat -> pattern
pol/poly -> polynomial
omg -> omega
******************************/
using namespace std;
mt19937 rnd(random_device{}());
int rndd(int l, int r){return rnd() % (r - l + 1) + l;}
typedef unsigned int uint;
typedef unsigned long long unll;
typedef long long ll;
class Polynomial{
private:
int lena, lenb;
int len;
comp A[2100000], B[2100000];
public:
comp Omega(int, int, bool);
void Init(void);
void FFT(comp*, int, bool);
void Reverse(comp*);
void MakeFFT(void);
}poly;
template<typename T = int>
inline T read(void);
int main(){
poly.Init();
poly.MakeFFT();
fprintf(stderr, "Time: %.6lf\n", (double)clock() / CLOCKS_PER_SEC);
return 0;
}
void Polynomial::MakeFFT(void){
FFT(A, len, DFT), FFT(B, len, DFT);
for(int i = 0; i <= len; ++i)A[i] *= B[i];
FFT(A, len, IDFT);
for(int i = 0; i <= lena + lenb - 2; ++i)
printf("%d%c", int(A[i].real() / len + eps + 0.5), i == lena + lenb - 2 ? '\n' : ' ');
}
void Polynomial::Reverse(comp* pol){
int pos[len + 10];
memset(pos, 0, sizeof(pos));
for(int i = 0; i < len; ++i){
pos[i] = pos[i >> 1] >> 1;
if(i & 1)pos[i] |= len >> 1;
}
for(int i = 0; i < len; ++i)if(i < pos[i])swap(pol[i], pol[pos[i]]);
}
void Polynomial::FFT(comp* pol, int len, bool pat){
Reverse(pol);
for(int size = 2; size <= len; size <<= 1){
for(comp* p = pol; p != pol + len; p += size){
int mid(size >> 1);
for(int i = 0; i < mid; ++i){
auto tmp = Omega(size, i, pat) * p[i + mid];
p[i + mid] = p[i] - tmp;
p[i] = p[i] + tmp;
}
}
}
}
void Polynomial::Init(void){
lena = read(), lenb = read();
for(int i = 0; i <= lena; ++i)A[i].real((double)read());
for(int i = 0; i <= lenb; ++i)B[i].real((double)read());
len = 1;
lena++, lenb++;
while(len <= lena + lenb)len <<= 1;
}
comp Polynomial::Omega(int n, int k, bool pat){
if(pat == DFT)return comp(cos(2 * PI * k / n), sin(2 * PI * k / n));
return conj(comp(cos(2 * PI * k / n), sin(2 * PI * k / n)));
}
template<typename T>
inline T read(void){
T ret(0);
int flag(1);
char c = getchar();
while(c != '-' && !isdigit(c))c = getchar();
if(c == '-')flag = -1, c = getchar();
while(isdigit(c)){
ret *= 10;
ret += int(c - '0');
c = getchar();
}
ret *= flag;
return ret;
}
原根
详细定义可参考 知乎 或 OI-WIKI，简而言之就是，对于模 $ P $ 意义下的原根 $ g $，有 $ g^1, g^2, \cdots, g^{P - 2} \bmod m $ 的值各不相同。
NTT
中文全称快速数论变换，和快速傅里叶变换的区别就是前者在模意义下，后者在复数意义下。
而可以证明，模意义下的可以通过原根代替单位根，特别地，假设存在模 $ P $ 意义下的原根 $ g $，若满足 $ n \mid P - 1 $，则令：
\[g_n^k = (g^{\tfrac{P - 1}{n}})^k \]我们只需要验证 $ g_n^k $ 满足单位根需要的几个性质即可。
显然 $ g_n^n = g^{P - 1} $ ，且 $ g_{an}^{ak} = (g^\tfrac{P - 1}{an})^{ak} = g_n^k $，且 $ g_n^0 = 1 $ 显然成立。
然后由于 $ g \bot P $ 且若 $ P $ 为质数，则根据欧拉定理可知 $ g^{P - 1} \equiv 1 \pmod{P} $ 则显然 $ g_n^{k + \tfrac{n}{2}} = g_{2n}^{2k + n} = g_{2n}^{2k} \times g_n^n = g_{2n}^{2k} \times g^{P - 1} = g_{2n}^{2k} = g_{n}^{k} $。
故可以直接用 $ g_n^k $ 代替 $ \omega_n^k $。
对于单位根的倒数，在 FFT 中通过共轭复数实现，这里我们便通过乘法逆元实现。
本质上并无太大区别，这里同样提供一个实现：（以前写的代码，可能实现不够精细）
对应题目：P3803 【模板】多项式乘法（FFT）。
Code
#define _USE_MATH_DEFINES
#include <bits/stdc++.h>
#include <mmintrin.h>
#define PI M_PI
#define E M_E
#define DFT true
#define IDFT false
#define eps 1e-6
#define MOD 998244353
/******************************
abbr
pat -> pattern
pol/poly -> polynomial
omg -> omega
******************************/
using namespace std;
mt19937 rnd(random_device{}());
int rndd(int l, int r){return rnd() % (r - l + 1) + l;}
typedef unsigned int uint;
typedef unsigned long long unll;
typedef long long ll;
ll kpow(int a, int b){
ll ret(1ll), mul((ll)a);
while(b){
if(b & 1)ret = (ret * mul) % MOD;
b >>= 1;
mul = (mul * mul) % MOD;
}
return ret;
}
class Polynomial{
private:
int lena, lenb;
int len;
int g, inv_g;
int A[2100000], B[2100000];
public:
int Omega(int, int, bool);
void Init(void);
void NTT(int*, int, bool);
void Reverse(int*);
void MakeNTT(void);
}poly;
template<typename T = int>
inline T read(void);
int main(){
// freopen("P3803_4.in", "r", stdin);
poly.Init();
poly.MakeNTT();
fprintf(stderr, "Time: %.6lf\n", (double)clock() / CLOCKS_PER_SEC);
return 0;
}
void Polynomial::MakeNTT(void){
NTT(A, len, DFT), NTT(B, len, DFT);
for(int i = 0; i <= len; ++i)A[i] = ((ll)A[i] * B[i]) % MOD;
NTT(A, len, IDFT);
int mul_inv = kpow(len, MOD - 2);
for(int i = 0; i <= lena + lenb - 2; ++i)
printf("%d%c", (ll)A[i] * mul_inv % MOD, i == lena + lenb - 2 ? '\n' : ' ');
}
void Polynomial::Reverse(int* pol){
int pos[len + 10];
memset(pos, 0, sizeof(pos));
for(int i = 0; i < len; ++i){
pos[i] = pos[i >> 1] >> 1;
if(i & 1)pos[i] |= len >> 1;
}
for(int i = 0; i < len; ++i)if(i < pos[i])swap(pol[i], pol[pos[i]]);
}
void Polynomial::NTT(int* pol, int len, bool pat){
Reverse(pol);
for(int size = 2; size <= len; size <<= 1){
int gn = kpow(pat == DFT ? g : inv_g, (MOD - 1) / size);
for(int* p = pol; p != pol + len; p += size){
int mid(size >> 1);
int g(1);
for(int i = 0; i < mid; ++i, g = ((ll)g * gn) % MOD){
auto tmp = ((ll)g * p[i + mid]) % MOD;
p[i + mid] = (p[i] - tmp + MOD) % MOD;
p[i] = (p[i] + tmp) % MOD;
}
}
}
}
void Polynomial::Init(void){
lena = read(), lenb = read();
for(int i = 0; i <= lena; ++i)A[i] = read();
for(int i = 0; i <= lenb; ++i)B[i] = read();
len = 1;
lena++, lenb++;
while(len < lena + lenb)len <<= 1;
g = 3;
inv_g = kpow(g, MOD - 2);
}
template<typename T>
inline T read(void){
T ret(0);
short flag(1);
char c = getchar();
while(c != '-' && !isdigit(c))c = getchar();
if(c == '-')flag = -1, c = getchar();
while(isdigit(c)){
ret *= 10;
ret += int(c - '0');
c = getchar();
}
ret *= flag;
return ret;
}
多项式取模
这个没什么可说的，是个概念性问题，可以认为如对于 $ A(x) \bmod{x^n} $，即表示保留多项式 $ A(x) $ 的 $ [0, n - 1] $ 次方项，也就是保留前 $ n $ 项。
多项式求逆
显然是在模意义下的，简而言之就是给定 $ F(x) $，求 $ G(x) $ 满足 $ F(x) \ast G(x) \equiv 1 \pmod{x^n} $。
模板题：P4238 【模板】多项式乘法逆。
尝试将问题规模缩减。
考虑如果我们已知存在 $ H(x) $ 满足 $ F(x) \ast H(x) \equiv 1 \pmod{x^{\lceil \tfrac{n}{2} \rceil}} $。
又显然有 $ F(x) \ast G(x) \equiv 1 \pmod{x^{\lceil \tfrac{n}{2} \rceil}} $。
两式相减得 $ F(x) \ast (G(x) - H(x)) \equiv 0 \pmod{x^{\lceil \tfrac{n}{2} \rceil}} $。
显然 $ F(x) $ 已给定，则有 $ G(x) - H(x) \equiv 0 \pmod{x^{\lceil \tfrac{n}{2} \rceil}} $。
同余式两侧分别平方，令 $ P(x) = (G(x) - H(x))^2 $，则有 $ P_i = \sum_{j = 0}^{i} (G(x) - H(x))j(G(x) - H(x)){i - j} $，显然 $ j $ 和 $ i - j $ 中至少有一个不大于 $ \lceil \dfrac{n}{2} \rceil $，所以一定有 $ P(x) \equiv 0 \pmod{x^n} $。
展开即为：$ G(x)^2 + H(x)^2 - 2G(x)H(x) \equiv 0 \pmod{x^{n}} $。
同余式两边各乘一个 $ F(x) $ 则有：$ G(x) \equiv 2H(x) - F(x)H(x)^2 \pmod{x^n} $。
此时可以发现我们将问题规模缩小了一半，则可以如此递归下去。
显然边界为 $ n = 1 $ 的时候，多项式求逆退化成一般的乘法逆元，直接求即可。
通过主定理分析复杂度，不难发现合并可以通过 NTT 优化到 $ O(n \log n) \(，则有：\) T(n) = T(\dfrac{n}{2}) + O(n \log n) $。
故最终复杂度应为 $ O(n \log^2 n) $。
然后这个代码实现比较乱，因为空间卡的比较严格，然后我最开始想封装一下，但是递归导致空间爆炸，并且实现上在 DFT 之后可以直接合并而不需要这么多操作，不过懒得改了，最后卡过去了，下次一定。
Code
#define _USE_MATH_DEFINES
#include <bits/stdc++.h>
#define PI M_PI
#define E M_E
#define npt nullptr
#define SON i->to
#define OPNEW void* operator new(size_t)
#define ROPNEW(arr) void* Edge::operator new(size_t){static Edge* P = arr; return P++;}
using namespace std;
mt19937 rnd(random_device{}());
int rndd(int l, int r){return rnd() % (r - l + 1) + l;}
bool rnddd(int x){return rndd(1, 100) <= x;}
typedef unsigned int uint;
typedef unsigned long long unll;
typedef long long ll;
typedef long double ld;
#define MOD (998244353ll)
#define DFT (true)
#define IDFT (false)
template < typename T = int >
inline T read(void);
int N;
int lim;
ll g = 3, inv_g;
int pos[280000];
ll qpow(ll a, ll b){
ll ret(1), mul(a);
while(b){
if(b & 1)ret = ret * mul % MOD;
b >>= 1;
mul = mul * mul % MOD;
}return ret;
}
class Polynomial{
private:
public:
int len;
int P[280000];
void Reverse(void){
pos[0] = 0;
for(int i = 1; i <= len - 1; ++i){
pos[i] = pos[i >> 1] >> 1;
if(i & 1)pos[i] |= len >> 1;
}for(int i = 0; i <= len - 1; ++i)if(i < pos[i])swap(P[i], P[pos[i]]);
}
void NTT(bool pat){
Reverse();
for(int siz = 2; siz <= len; siz <<= 1){
ll gn = qpow(pat ? g : inv_g, (MOD - 1) / siz);
for(int* p = P; p != P + len; p += siz){
int mid = siz >> 1; ll g(1);
for(int i = 0; i < mid; ++i, (g *= gn) %= MOD){
ll tmp = (ll)g * p[i + mid] % MOD;
p[i + mid] = ((ll)p[i] - tmp) % MOD;
p[i] = ((ll)p[i] + tmp) % MOD;
}
}
}
if(!pat){
ll inv_len = qpow(len, MOD - 2);
for(int i = 0; i <= len - 1; ++i)P[i] = ((ll)P[i] * inv_len) % MOD;
}
}
void Print(void){
for(int i = 0; i <= len - 1; ++i)printf("%d%c", P[i], i == len - 1 ? '\n' : ' ');
}
Polynomial operator %= (const ll &mod){
for(int i = mod; i <= len - 1; ++i)P[i] = 0;
len = mod;
return *this;
}
Polynomial operator *= (Polynomial B){
int rlen(1);
while(rlen < len + B.len)rlen <<= 1;
len = B.len = rlen;
// printf("len = %d\n", rlen);
// ret.len = rlen;
// A.Print(), B.Print();
NTT(DFT), B.NTT(DFT);
// A.Print(), B.Print();
for(int i = 0; i <= len - 1; ++i)P[i] = (ll)P[i] * B.P[i] % MOD;
NTT(IDFT);
// ret.Print();
return (*this) %= lim;
}
Polynomial operator *= (const ll &mul){
for(int i = 0; i <= len - 1; ++i)P[i] = ((ll)P[i] * mul) % MOD;
return *this;
}
Polynomial operator -= (const Polynomial &B){
len = max(len, B.len);
for(int i = 0; i <= len - 1; ++i)P[i] = (((ll)P[i] - (ll)B.P[i]) % MOD + MOD) % MOD;
return *this;
}
}F, H, tmp;
Polynomial Make(int len = N){
// printf("making len = %d\n", len);
if(len == 1){
Polynomial ret;
ret.len = 1, ret.P[0] = qpow(F.P[0], MOD - 2);
return ret;
}
H = Make(int(ceil((double)len / 2.0)));
// H.Print();
lim = len;
// printf("len = %d, 2H:", len); (H * 2).Print();
// printf("F: "); F.Print();
// printf("F * H: "); (F * H).Print();
// printf("F * H * H: "); (F * H * H).Print();
tmp = H;
return (H *= 2) -= ((tmp *= H) *= F);
// return H * 2 - F * H * H;
}
int main(){
// freopen("P4238_16.in", "r", stdin);
inv_g = qpow(g, MOD - 2);
F.len = N = read();
for(int i = 0; i <= N - 1; ++i)F.P[i] = read();
// decltype(F) tmp; tmp.len = 1; tmp.P[0] = 1;
// F.Print();
// tmp.Print();
// (F * tmp).Print();
Make().Print();
// for(int i = 0; i <= N - 1; ++i)printf("%lld%c", G.P[i], i == N - 1 ? '\n' : ' ');
fprintf(stderr, "Time: %.6lf\n", (double)clock() / CLOCKS_PER_SEC);
return 0;
}
template < typename T >
inline T read(void){
T ret(0);
int flag(1);
char c = getchar();
while(c != '-' && !isdigit(c))c = getchar();
if(c == '-')flag = -1, c = getchar();
while(isdigit(c)){
ret *= 10;
ret += int(c - '0');
c = getchar();
}
ret *= flag;
return ret;
}
大概是一些比较常用的，后面还有很多高端操作，先去补补其它算法再回来继续搞多项式。
UPD
update-2022_12_22 FFT NTT 多项式求逆
标签：int,多项式,void,ret,omega,define From： https://www.cnblogs.com/tsawke/p/17032830.html