多项式小寄
多项式小记
目录
序
十分的基础，还有很多东西 没有学习，准备 鸽了 学了 \(DS\) 之后再来补充
这篇博客主讲 \(FFT\)，\(NTT\) 的 基本原理 以及 多项式简单运算，而 题目涉及较少（没做
文中可能有不少 \(\LaTeX\) 书写错误，还有各种 \(n, N\)，\(m, M\) 混用 的情况，感谢包容
（应该不会 \(n, m\)，\(N, M\) 不分）
强烈建议 阅读每个 小专题 对应 引用的资料，都是在 对应板块 讲解的十分详细的
更新日志
$2024.04.28 ~ UPD: $ 修了一些锅，规范化一些语言，保留一些废话
\(2024.05.03 ~ UPD:\) 修了一些锅，补充了 多项式求逆 - 优化 - 循环卷积 的证明 关键
\(2024.05.06 ~ UPD:\) 修了一些锅，注意到 多项式开根 当次数为 \(N\) 时，补项 不取等 并没有正确性
什么是多项式
即形如 \(F(x) = a_1 x ^ {k_1} + a_2 x ^ {k_2} + ... + a_{k_1} x ^ 1 + a_{k_1 + 1} x ^ 0\) 的 式子，也可以理解为一个 函数
\(OI\) 中一般讨论的 多项式 默认 只有一元，也就是 一个变量 \(x\)
上式可以直接写作 左式 \(F(x)\) 的形式，也可以直接称作 多项式 \(F\)
这都不影响其 实际指代（也就是 右式）
（上述 右式 也可以写为 \(\sum a_k x^{k}\)）
\(eg. F(x) = x ^ 2 + 5x + 2\) 故有 \(F(5) = 52\)
多项式四则运算
加法
存在多项式 \(F(x) = \sum a_k x ^ k, G(x) = \sum b_k x ^ k\)
则有 \(F(x) + G(x) = \sum (a_k + b_k) x ^ k\)
\(eg. F(x) = 2 x ^ 2 + 5x + 4\)
\(G (x) = 4x ^ 3 + 3x + 6\)
则 \(F(x) + G(x) = (0 + 4) x ^ 3 + (2 + 0) x ^ 2 + (5 + 3) x ^ 1 + (4 + 6) x ^ 0 = 4 x ^ 3 + 2 x ^ 2 + 8 x + 10\)
减法
和加法 显然互逆，不多言
乘法
存在多项式 \(F(x) = \sum_i a_i x ^ i, G(x) = \sum_j b_j x ^ j\)
则有
除法
和乘法 显然互逆，注意 余式
\(FFT\)
以下内容多引用自 [1] 对应的文章
多项式的四则运算 中，加减法都可以在 线性时间 内很快的解决
而如果按照定义暴力计算，多项式乘法（加法卷积）则需要 \(O(N^2)\) 的时间来解决
题内话：多项式乘法 为啥是 加法卷积
我们一直在说 多项式乘法 是 加法卷积，但是它们的关系到底是什么？
考虑卷积的基本形式 长成这样
\[\sum _ {i + j = k} a_i b_j \]其中 \(k\) 是一个给定常数，\(i, j \in [0,k]，i, j \in N\)（有些时候 \(i, j \in [1, k]\)，无伤大雅）
然后我们注意到
它就是 两多项式 \(F(x) = \sum a _ i x ^ i\) 和 \(G(x) = \sum b _ i x ^ i\) 的 积多项式 的 \(k\) 次项系数
所以我们做 \(F(x) * G (y)\)，就相当于对于 所有的 \(k \in [0, N + M]\) 求得了 上述卷积式
（这里假设 \(F\) 是 \(N\) 次多项式，\(G\) 是 \(M\) 次多项式）
并且由于 多项式项 之间用 \(+\) 连接，故称 加法卷积
非常的慢！
而 \(FFT~(Fast ~ Fourier ~ Transformation)\) ，即 快速傅里叶离散变换，
就是在 \(O(N \log N)\) 时间下解决 多项式乘法 的一种算法
这玩意儿十分的精妙，他利用到一个非常重要的性质，即
多项式的点值表达
就像 二次函数 有 两点式（加上二次项系数确定这个二次函数） 一样
显然，一个 \(N\) 次多项式 可以由其式子 对应函数 上的 \(N + 1\) 个点来确认
所以我们可以把一个 \(N\) 次多项式从 系数表达 转化成 \(N + 1\) 个点来表达，这一步也叫 求值
这有什么好处呢？
可以发现，多项式的四则运算 可以 直接映射到 点坐标的四则运算 上
也就是若 \(F(x)\) 过点 \((x_1, y_1)\) ，而 \(G (x)\) 过点 \((x_2, y_2)\)
则 \(F(x_1) + G(x_2)\) 一定等于 \(y_1 + y_2\)，多项式的乘法 也满足，即 \(F(x_1) * G(x_2)\) 等于 \(y_1 y_2\)
而对于 \(N\) 次多项式 \(F(x)\)，\(M\) 次多项式 \(G(x)\)，它们的积最高只有 \(N + M\) 次
也就是我们至多在两个函数上 各找 \(N + M + 1\) 个点运算后就可以确定 \(F(x) * G (x)\)
显然 知道点值表达 后的运算只需要 \(O(N + M)\) 的时间 （\(N, M\) 同阶时认为是 \(O(N)\)）
于是重点就变成了 用系数表达求点值表达 和 用点值表达反推系数表达
亦即 求值 和 插值
\(DFT\)
单位根
按理说 求值 代入 任意点 都没有问题，可是 求取任意点 就只能 代入多项式 做到
这玩意儿求一个点就是 \(O(N)\) 的时间，而我们需要 \(N + M\) 个点，完蛋
于是我们需要一些 有性质的数 来带入求值，而 单位根 就是一种好东西
单位根 / \(n\) 次单位根 ：即 \(n\) 次幂 为 \(1\) 的 复数，用 \(w_n\) 表示
显然 \(n\) 次单位根 共 \(n\) 个，写作 \(w_n ^ k\) 的形式
根据 复数乘法 模长相乘，幅角相加 的法则，发现 单位根 只存在于 半径为 \(1\) 的 单位圆 上
且 \(n\) 次单位根 将此 单位圆 \(n\) 等分
故可以得到 单位根 的如下性质（之后会用到的）
-
\(w_n ^ 0 = 1 + 0i = 1\)
-
\(w_n ^ k = (w_n ^ 1) ^ k\)
-
\(w_n ^ k = w_n ^ {k \% n}\)
-
\((w_n ^ k) ^ p = w_n ^ {kp} = w_n ^ {kp \% n}\)
-
\(w_n ^ i * w_n ^ j = w_n ^ {i + j}\)
-
\(w_{an} ^ {ak} = w_n ^ k \to w_{2n} ^ {2k} = w_n ^ k\)
-
\(w_n ^ {k + n/2} = - w_n ^ {k} ~ (2 | n)\)
（通过 复数乘法 法则都可以简单证明）
算法原理
现在我们可以回到 \(DFT\) 的原理了
考虑一个 \(n - 1\) 次（\(n\) 项）多项式 \(F(x) = \sum_{k = 0} ^ {n - 1} a_k x^k\)
我们按 次数奇偶 分成两部分，这里保证 \(N\) 是 \(2\) 的整数次幂，从而可以 平均分开
得到
于是设
显然我们可以发现
我们代入 \(x = w_n ^ k\)，因为
故上式变为
显然 \(k < n / 2\)，故我们再代入 \(x = w_n ^ {k + n / 2}\)，有
又 \((w_n ^ k) ^ 2 = w_n ^ {2k}\)，有
又 \(w_n ^ k = w_n ^ {k \% n}\)，有
又 \(w_{2n} ^ {2k} = w_n ^ k\)，有
又 \(w_n ^ {k + n/2} = - w_n ^ {k} ~ (2 | n)\)，有
核心流程
把 刚刚的玩意儿 \(②\) 和 上面推出的一个式子 \(①\) 放在一起
艹？长得差不多？只有符号的区别？
非常的好，于是现在我们只需要知道 \(FL, FR\) 在 \(w_{n / 2} ^ {0}, w_{n / 2} ^ {1}, ..., w_{n / 2} ^ {n / 2 - 1}\) 处的 点值表达 即可
而 \(FL, FR\) 的性质与原式 \(F\) 完全相同（右式 形式一致）仅有最高次的区别
可以把它们看作次数 减小为一半 的一般多项式，我们很容易可以想到 分治
通过上式将 原问题 不断划分为规模只有 一半 的 子问题，递归求解，项数为 \(1\) 时 直接返回
最后按上述两式子，即可分别得出 \(F\) 在 \(w_n ^ 0 ... w_n ^ {n / 2 - 1}\) 与 \(w_n ^ {n / 2} ... w_n ^ {n - 1}\) 处的 点值表达 了
时间复杂度 \(T(N) = 2 T(N / 2) + O(N) = O(N \log N)\)
注意到我们需要 多项式项数 为 \(2\) 的 整次幂 从而保证 划分始终均匀
所以当 项数不足 时直接在 最高次补足项 即可，显然 补的项数不超过原项数，复杂度不变
代码
我们可以 重载复数运算 方便后续编写，这里提供一个 重载好的复数结构体
struct Complex {
double R, I;
inline Complex (double r = 0, double i = 0) {
R = r, I = i;
}
inline Complex operator + (const Complex &a) const {
return Complex (R + a.R, I + a.I);
}
inline Complex operator - (const Complex &a) const {
return Complex (R - a.R, I - a.I);
}
inline Complex operator * (const Complex &a) const {
return Complex (R * a.R - I * a.I, I * a.R + R * a.I);
}
// 除法是不必要的
};
这里给出一个 按照上述流程 的 \(DFT\) 递归实现，后面会讲优化
注意到由于存在 \(0\) 次项，以及为了后续分治均匀，下标最好从 \(0\) 开始
const double PI = acos (-1);
Complex S[MAXN];
inline void DFT (Complex* p, const int Len) {
if (Len == 1) return ;
Complex *l = p, *r = (p + (Len >> 1));
for (int i = 0; i < Len; ++ i) S[i] = p[i];
for (int i = 0; i < (Len >> 1); ++ i)
l[i] = S[i << 1], r[i] = S[i << 1 | 1];
DFT (l, Len >> 1), DFT (r, Len >> 1);
Complex U = Complex (cos (2 * PI / Len), sin (2 * PI / Len)); // 第一单位根
Complex Now = Complex (1, 0); // 当前单位根
for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
S[k] = l[k] + Now * r[k], S[k + (Len >> 1)] = l[k] - Now * r[k];
for (int i = 0; i < Len; ++ i) p[i] = S[i];
}
这个函数将多项式 系数数列 转化为 点值数列，初始传入 系数数列 \(p\)
做完后的 \(p\) 数组就是 多项式 分别在 \(w_n ^ 0 ... w_n ^ {n - 1}\) 处的 点值
\(IDFT\)
求值 是 求完了，我们得到了 一堆 表示 积多项式 的 点
但是这玩意儿它没用啊，我们要的是 系数表达，所以还需要 插值
也就是还得要个东西 把刚刚的流程反过来，得到 最终的多项式（系数表达）
结论
这东西流程上是很简单的，可以先来个 结论
把 \(DFT\) 中的 \(w_n ^ 1\) 换成 \(w_n ^ {-1}\)，做完之后 除以 \(n\) 即可
（可以做题了
证明
以下给出证明，设有 \(G = DFT (F)\)，也就是 \(G\) 数列 是 \(F\) 数列 \(DFT\) 变换后的 结果
模拟暴力代入，（设将 \(w_n ^ k\) 代入以求得第 \(k\) 个点 的 点值表达），则有
（\(G_k\) 表示 \(G\) 多项式的第 \(k\) 项，同形式同理）
我们需要证得 的结论是
于是将 \(①\) 代入 \(②\) 有
分两部分贡献讨论，（最后 两部分相加）若 \(j = k\)，则（有 \(w_n ^ 0 = 1\)）
若 \(j \neq k\)，设 \(p = j - k\)，则
即需要证得 \(F_k = R_1 + R_2\)，而我们注意到
等比数列求和有
所以
于是最终有 \(R_1 + R_2 = F_k + 0 = F_k\)，即 结论得证
代码
就，和 \(DFT\) 不能说 一模一样，只能说 长得很像
const double PI = acos (-1);
Complex S[MAXN];
inline void DFT (Complex* p, const int Len) {
if (Len == 1) return ;
Complex *l = p, *r = (p + (Len >> 1));
for (int i = 0; i < Len; ++ i) S[i] = p[i];
for (int i = 0; i < (Len >> 1); ++ i)
l[i] = S[i << 1], r[i] = S[i << 1 | 1];
DFT (l, Len >> 1), DFT (r, Len >> 1);
Complex U = Complex (cos (2 * PI / Len), - sin (2 * PI / Len)); // 第一单位根 (注意取负)
Complex Now = Complex (1, 0); // 当前单位根
for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
S[k] = l[k] + Now * r[k], S[k + (Len >> 1)] = l[k] - Now * r[k];
for (int i = 0; i < Len; ++ i) p[i] = S[i];
}
最后的实现
显然，有了 \(DFT / IDFT\)，最后的乘法是简单的
我们只需要 补齐项数，对最开始的两个多项式分别做 \(DFT\)
后面 深究了一下 补项 这个问题，还有点 神秘
如果对常数没要求的话，封装每个运算（乘法，求逆，开根...）的时候 暴力补满 即可
也就是
p <= N + M或如果 只有一个多项式 时p <= N
但如果 对常数有一些要求，那么这样凭空多出的 \(1\) 倍时间就 不可忽略
那具体而言什么时候 需要取等，什么时候 可以不取等 呢？
比如下面的代码，注意到实质上 \(N\) 次多项式有 \(N + 1\) 项
\(N\) 次 多项式 和 \(M\) 次 多项式 相乘 最多会有 \(N + M + 1\) 项，自然要 取等
而 后面的封装 传入多项式后 传入的是 项数 而非 次数
\(N\) 项 多项式 和 \(M\) 项 多项式 相乘 最多只会有 \(N + M - 1\) 项，不用 取等
然后将求出的 点值数列 相乘，对乘积做 \(IDFT\)
最后输出结果前 \(N + M\) 项即可，这部分代码如下
inline Mul () {
cin >> N >> M;
for (P = 1; P <= N + M; P <<= 1);
for (int i = 0; i <= N; ++ i) cin >> F[i].R;
for (int i = 0; i <= M; ++ i) cin >> G[i].R;
FFT (F, P, 1), FFT (G, P, 1);
for (int i = 0; i < P; ++ i) F[i] = F[i] * G[i];
FFT (F, P, - 1);
for (int i = 0; i <= N + M; ++ i)
cout << (int) (F[i].R / P + 0.49) << ' ';
}
融融
用 \(EPS\) 个 \(zhicheng\) 的注意力 可以发现
\(DFT\) 和 \(IDFT\) 的 代码差距 和 它们的 名字差距 一样巨大
据完全统计，有整整 一个字符
所以为了节省码量，我们可以把它们 绒绒 起来
const double PI = acos (-1);
Complex S[MAXN];
inline void FFT (Complex* p, const int Len, const int Typ) {
if (Len == 1) return ;
Complex *l = p, *r = (p + (Len >> 1));
for (int i = 0; i < Len; ++ i) S[i] = p[i];
for (int i = 0; i < (Len >> 1); ++ i)
l[i] = S[i << 1], r[i] = S[i << 1 | 1];
FFT (l, Len >> 1, Typ), FFT (r, Len >> 1, Typ);
Complex U = Complex (cos (2 * PI / Len), Typ * sin (2 * PI / Len));
Complex Now = Complex (1, 0);
for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
S[k] = l[k] + Now * r[k], S[k + (Len >> 1)] = l[k] - Now * r[k];
for (int i = 0; i < Len; ++ i) p[i] = S[i];
}
// Use DFT -> FFT (p, Len, 1), Use IDFT -> FFT (p, Len, -1)
感觉可以了？交一波 P3803 【模板】多项式乘法（FFT）
可以看到因为常数过大而 \(TLE\) 了,看来我们需要一个常数更小的写法 —— \(command\_block\)
（显然，\(command\_block\) 老师 没有体验到 洛 谷 4 . 0 的 强 大）
（以上代码可以 轻 松 通 过 ）
于是我们继续来做
没用的 常数优化
观察 如下代码
for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
S[k] = l[k] + Now * r[k], S[k + (Len >> 1)] = l[k] - Now * r[k];
显然 Now * r[k] 这种 复数乘法 常数极大，可以 尝试省去，如下
Complex tmp = (1, 0);
for (int k = 0; k < (Len >> 1); ++ k, Now = Now * U)
tmp = Now * r[k], S[k] = l[k] + tmp, S[k + (Len >> 1)] = l[k] - tmp;
实测 会快 \(EPS\) ？笔者更倾向于 编译器比较聪明
还是得来点 有用的常数优化
蝴蝶变换
\(S\) 这个临时数组 是 大 常 数 的 一 大 助 力，减少数组复制总能 显著减小常数
我们可以发现，递归下去的过程中 \(p\) 数组 值的位置在变化，具体来讲
0 1 2 3 4 5 6 7 - Flr1
0 2 4 6 | 1 3 5 7 - Flr2
0 4 | 2 6 | 1 5 | 3 7 - Flr3
0 | 4 | 2 | 6 | 1 | 5 | 3 | 7 - Flr4
可以发现，最终位置 是 原位置 的 二进制反转 的值，而这个值可以 \(O(N)\) 递推得到
for (int i = 0; i < N; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? N >> 1 : 0);
数组复制，没了，效 果 显 著！
递归变递推（迭代）
显然，递归 令人不喜，如果有能够 递推 的形式那再好不过
我们可以直接尝试把 递归的过程反过来
也就是 最开始 就把数放到 递归的时候对应数最终的位置上，然后 往上合并
inline void FFT (Complex *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1;
Complex Now = Complex (1, 0), tmp, U;
for (int k = 2; k <= Len; len = k, k <<= 1) {
U = Complex (cos (2 * PI / k), Typ * sin (2 * PI / k));
for (int i = 0; i < Len; i += k) {
Now = Complex (1, 0);
for (int l = i; l < i + len; ++ l, Now = Now * U)
tmp = Now * p[l + len], p[l + len] = p[l] - tmp, p[l] = p[l] + tmp;
}
}
}
// k 即是当前 合并后 的区间长, i 是区间起始处，我们每次将 [l, l + len) 与 [l + len, l + k) 这两个区间合并
进 步 十 分 巨 大，这个东西得配合上面的 蝴蝶变换 来用才行
三步变两步
简单来说，这基于 下面的式子，设存在多项式 \(F, G\)，则
这个东西不很显然，但是可以用 多项式四则运算的性质 与 实数 基本相同 来 感性理解
也就是 多项式乘法 本身具有一般意义上的 分配律
所以我们可以把多项式 \(F, G\) 分别 放到 \(p\) 数组的 实部，虚部
做一遍 \(DFT\)，然后自身乘方后做 \(IDFT\)，我们需要的是 \(F*G\)
于是最后输出 虚部 的 \(\dfrac {1}{2}\) 即可，这东西就是 不严格的 \(\dfrac {1}{3}\) 巨量提升
inline void Mul () {
cin >> N >> M;
for (P = 1; P < N + M; P <<= 1);
for (int i = 0; i < P; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? P >> 1 : 0);
for (int i = 0; i <= N; ++ i) cin >> F[i].R;
for (int i = 0; i <= M; ++ i) cin >> F[i].I;
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) F[i] = F[i] * F[i];
FFT (F, P, - 1);
for (int i = 0; i <= N + M; ++ i)
cout << (int) (F[i].I / P / 2 + 0.49) << ' ';
}
破除封建迷信
谁说自己重载 \(Complex\) 会快的？标准库的轮子 还是 有实力的
所以我们可以把
struct Complex {
...
};
改成（当然直接修改 每个地方的定义 也行）
#define Complex complex <double>
然后修改对应的 输入输出（这玩意儿好像不支持直接 cin）
// Input Part
int x;
for (int i = 0; i <= N; ++ i) cin >> x, F[i] += x;
for (int i = 0; i <= M; ++ i) cin >> x, F[i] += Complex (0, x);
// Output Part
for (int i = 0; i <= N + M; ++ i) cout << (int) (F[i].imag () / P / 2 + 0.499) << ' '; // 注意 .imag ()
好了，现在再交交看，可以发现 std::complex 十 分 强 大
最终代码
2024-02-23 17:33 - C++20 O2 - 1.24 KB / 952 ms / 40.54 MB | Record
#include <bits/stdc++.h>
#define Complex complex <double>
const int MAXN = 2100005;
const double PI = acos (-1);
using namespace std;
int U[MAXN];
inline void FFT (Complex *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1;
Complex Now = Complex (1, 0), tmp, U;
for (int k = 2; k <= Len; len = k, k <<= 1) {
U = Complex (cos (2 * PI / k), Typ * sin (2 * PI / k));
for (int i = 0; i < Len; i += k) {
Now = Complex (1, 0);
for (int l = i; l < i + len; ++ l, Now = Now * U)
tmp = Now * p[l + len], p[l + len] = p[l] - tmp, p[l] = p[l] + tmp;
}
}
}
Complex F[MAXN];
int N, M, P, x;
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> N >> M;
for (P = 1; P <= N + M; P <<= 1);
for (int i = 0; i < P; ++ i)
U[i] = (U[i >> 1] >> 1) | ((i & 1) ? P >> 1 : 0);
for (int i = 0; i <= N; ++ i) cin >> x, F[i] += x;
for (int i = 0; i <= M; ++ i) cin >> x, F[i] += Complex (0, x);
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) F[i] = F[i] * F[i];
FFT (F, P, - 1);
for (int i = 0; i <= N + M; ++ i)
cout << (int) (F[i].imag () / P / 2 + 0.499) << ' ';
return 0;
}
\(MTT\)
\(Matthew ~ Transform\)（\(Most ~ TLE ~ Transform\)），毛啸变换，即 拆系数 \(FFT\)
由 毛啸（\(matthew99\)）在 国家集训队 \(2016\) 年论文集 《再探快速傅里叶变换》中提出
其核心是将 给定多项式系数 \(F_i\) 拆成 \(A_ik + B_i\)，其中 \(k\) 为常数 的形式
从而解决 值域过大 或 精度要求高 的问题
原理是非常简单的，我们先考虑一下 \(P4245\) 这个题，就是板子
它这东西的值域最多可以上到 \(\max (A_i) ^ 2 * N = 10 ^ {23}\)，令人难受
我们考虑把给定多项式的 每个系数 按 上述方法拆开
于是一个 \(F(x)\) 可以变换得到 \(A_0(x), A_1(x)\) 两个多项式
同理 \(G(x)\) 变换得到 \(B_0(x), B_1(x)\) 两个多项式，简单推一下式子
于是我们只需要算出 \(A_0(x)B_0(x)\)，\(A_0(x)B_1(x) + A_1(x)B_0(x)\)，\(A_1(x)B_1(x)\) 即可
显然我们知道，当 \(k = \sqrt {\max (A_i)}\) 时，上式 值域最小，只有 \(\sqrt {\max (A_i)} ^ 2 * N = 10 ^ {14}\)
用用 long double 直接就过了
而一般情况下，我们可以直接使用 \(7\) 次 \(FFT\) 实现（\(4DFT+3IDFT\)），
也可以优化到 \(4\) 次，但是由于实现难度高一些，难调，而且优化力度有限，所以一般不推荐
给一个 \(7\) 次的代码，\(4\) 次的直接取 板子题题解区 学习即可
以下代码较为 抽象，可能需要使用 不低于 \(12.2\) 版本的 \(GCC\) 才能 编译通过
#include <bits/stdc++.h>
#define Complex complex <long double>
const int MAXN = 262150;
const long double PI = acos (-1);
using namespace std;
int U[MAXN];
inline void Init (const int P) {
memset (U, 0, sizeof U);
for (int i = 0; i < P; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? P >> 1 : 0);
}
inline void FFT (Complex *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1;
Complex Now = Complex (1, 0), tmp, U;
for (int k = 2; k <= Len; len = k, k <<= 1) {
U = Complex (std::cos (2 * PI / k), Typ * std::sin (2 * PI / k));
for (int i = 0; i < Len; i += k) {
Now = Complex (1, 0);
for (int l = i; l < i + len; ++ l, Now = Now * U)
tmp = Now * p[l + len], p[l + len] = p[l] - tmp, p[l] = p[l] + tmp;
}
}
}
Complex F[MAXN], A0[MAXN], A1[MAXN], B0[MAXN], B1[MAXN];
Complex 🐒[MAXN], 🐧[MAXN], 🏄[MAXN];
unsigned long long A[MAXN], B[MAXN], C[MAXN], Ans[MAXN];
const int 🐨 = (1 << 15);
int N, M, P, Mod, x;
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> N >> M >> Mod, ++ N, ++ M;
for (P = 1; P < N + M; P <<= 1);
Init (P);
for (int i = 0; i < N; ++ i) cin >> x, A0[i] = x / 🐨, A1[i] = x % 🐨;
for (int i = 0; i < M; ++ i) cin >> x, B0[i] = x / 🐨, B1[i] = x % 🐨;
FFT (A0, P, 1), FFT (A1, P, 1), FFT (B0, P, 1), FFT (B1, P, 1);
for (int i = 0; i < P; ++ i)
🐒[i] = A0[i] * B0[i], 🐧[i] = A0[i] * B1[i] + B0[i] * A1[i], 🏄[i] = A1[i] * B1[i];
FFT (🐒, P, -1), FFT (🐧, P, -1), FFT (🏄, P, -1);
for (int i = 0; i < P; ++ i)
A[i] = round (🐒[i].real () / P), B[i] = round (🐧[i].real () / P), C[i] = round (🏄[i].real () / P);
for (int i = 0; i < P; ++ i)
Ans[i] = ((A[i] % Mod * 🐨 * 🐨 % Mod) + (B[i] % Mod * 🐨 % Mod) + C[i] % Mod) % Mod;
for (int i = 0; i < N + M - 1; ++ i) cout << Ans[i] % Mod << ' ';
return 0;
}
\(NTT\)
全称是
日本电报电话公司\(Number-Theoretic ~ Transform\)，快速数论变换（也可以叫 \(Nityacke ~ Toptree ~ Timelimitexcceed\) ）
这东西还真只能接着 \(FFT\) 讲，没看 \(FFT\) 先看 \(FFT\)
咋来的
\(FFT\) 这玩意儿好吧？但是要用到 \(double\) 这 破东西
这就不可避免地带来了 精度 的问题，数据范围一大就 令人恼火
但是 单位根 这东西吧，它不好替代，至少
数学家们已经证明，在 复数域 中，单位根是 唯一一类 满足要求的数 —— \(command\_block\)
但是吧，大部分 \(OI\) 中的计数题都是在 模意义下的，于是咱对 ”这类数“ 的限制放宽了
在 模意义同余系 中，我们找到了 单位根 的 替代品，原根 \(g\)
原根
得讲讲这玩意儿是啥...
以下内容多引用自 [3] 对应的文章
对原根本身 比较了解 或 不想了解 的可以跳转至标题 "回到正题"
定义
如果 \(a\), \(m\) 互质，那么我们把 \(a ^ k \equiv 1 \pmod m\) 的 最小正整数 \(k\) 称为 \(a\) 的模 \(m\) 的 阶
我们可以把上面的东西写作 \(ord_m (a) = k\)
在 原文 中老师将前提 写成了 \(a, m\) 不互质
虽然这篇博客对笔者以下内容编写有很大作用，但...火大
于是显然，如果 \(a ^ n \equiv 1 \pmod m\)，则 \(k | n\)（余数的可乘性）
由 欧拉定理 可知，当 \(a, m\) 互质时，$a ^ {\varphi (m)} \equiv 1 \pmod m $
故 \(a\) 模 \(m\) 的阶 \(k\) 为 \(\varphi (m)\) 的因数，而
我们把满足 \(k = \varphi (m)\) 对应的 \(a\) 称为 模 \(m\) 的原根
\(eg. m = 3\)，有 \(\varphi (3) = 2\)，且 \(2\) 模 \(3\) 的阶 \(k = 2\)，故 \(2\) 是 \(3\) 的 原根
称之为 根，因为其是 \(x ^ {\varphi (m)} \equiv 1 \pmod m\) 的 一个解
称之为 原，因为 原神 因为
性质
先来点用不到的
-
（上面提到过） \(a\) 模 \(m\) 的阶 \(k\) 为 \(\varphi (m)\) 的因数
-
\(g\) 是模 \(m\) 的 原根 当且仅当 对于任何与 \(m\) 互质的数 \(a\) 都能找到 \(k\)，使得 \(a \equiv g ^ k \pmod m\)
显然由于 1. 如果 \(g\) 不是 模 \(m\) 的 原根，则 \(g\) 的 阶 \(k\) 小于 \(\varphi (m)\)
又根据阶的定义，则 \(g ^ k \equiv 1 \pmod m\)，故 \(g ^ {t} \equiv g ^ {t + ik} \pmod m\)
即由 \(g ^ i\) 生成的 模 \(m\) 意义下互不同余 的数至多 \(k\) 个，小于 \(\varphi (m)\) 个，必要性得证
而若 \(g\) 为 模 \(m\) 的 原根，则 \(g, g ^ 2, ..., g ^ {\varphi (m)}\) 除以 \(m\) 的余数 互不相同
反之若存在 \(g ^ i \equiv g ^ j \pmod m\)，不妨设 \(i < j\)，由于 \(g, m\) 互质，\(g ^ i\) 与 \(m\) 必然互质
根据 同余的除法性质，即 同余式两边可以同时除以 \(a\) 的条件是 \(a, m\) 互质
我们得到 \(g ^ {j - i} \equiv 1 \pmod m\)，与 原根 的 阶 \(k = \varphi (m)\) 矛盾
故 \(g, g ^ 2, ..., g ^ {\varphi (m)}\) 除以 \(m\) 的余数 互不相同 为 真命题，同上可知
\(g, g ^ 2, ... , g ^ {\varphi (m)}\) 均与 \(m\) 互质，而与 \(m\) 互质且小于 \(m\) 的数一共 \(\varphi (m)\) 个
又知 \(g, g ^ 2, ..., g ^ {\varphi (m)}\) 互不相同，故其一定为 所有小于 \(m\) 且与 \(m\) 互质的数 的 一种排列
-
设 \(g\) 是模 \(m\) 的一个原根，则 \(g ^ k\) 是原根 当且仅当 \(k, \varphi (m)\) 互质
这也是 原根的判定定理，可由此得出 原根的个数 以及延伸出 原根的存在性定理
对 必要性 考虑 反证法，因为 \(g ^ k\) 为 原根
则根据 2. \((g ^ k) ^ 1, (g ^ k) ^ 2, ... , (g ^ k) ^ {\varphi (m)}\) 模 \(m\) 互不相同，而若 \(gcd (k, \varphi (m)) = t ~ (t \neq 1)\)
则 \((g ^ {k}) ^ {\varphi (m) / t} \equiv (g ^ {\varphi (m)}) ^ {k / t} \equiv 1 \pmod m\)，显然 \(\varphi (m) / t < \varphi (m)\)
此处与 原根 的 阶 为 \(\varphi (m)\) 矛盾，故 必要性得证
充分性则考虑到若 \(k, \varphi (m)\) 互质，则总能找到 \(s, t\) 使得 \(sk + t\varphi (m) = 1\)
故可以有 \(g ^ i \equiv g ^ {i (sk + t \varphi (m))} \pmod m\)
则 \(g ^ i \equiv g ^ {isk} g ^ {it\varphi (m)} \pmod m\)，显然 \(g ^ {it \varphi (m)} \equiv g ^ {\varphi (m)} \equiv 1 \pmod m\)
故 \(g ^ i \equiv g ^ {isk} \equiv (g ^ k) ^ {is}\)，也就是 \(g ^ i\) 可以表示成 \(g ^ k\) 的 \(is\) 次方
而根据 2. \(g ^ i\) 可以生成所有与 \(m\) 互质的数，故 \(g ^ k\) 亦可生成，故 \(g ^ k\) 也为原根
故 如果模数 \(m\) 有原根存在，则模数 \(m\) 有 \(\varphi (\varphi (m))\) 个原根
（对于原根 \(g\)，满足 \(i, \varphi (m)\) 互质的 \(g ^ i\) 都是原根）
-
对于质数 \(p\)，模 \(p\) 的 原根 \(g\) 存在，且有 \(\varphi (\varphi (p)) = \varphi (p - 1)\) 个
这是 3. 的推广，同时会用到下面性质
4.1. 如果整数 \(a\) 模 \(p\) 的 阶 是 \(t\)，则 当且仅当 \(k, t\) 互质时，\(a ^ k\) 的 阶 也是 \(t\)
对 必要性 同样考虑 反证法，因为 \(a\) 模 \(p\) 的 阶 为 \(t\)，故 \(a ^ t \equiv 1 \pmod p\)
显然 \(a ^ {it} \equiv 1 \pmod p\)，又设若 \(k, t\) 不互质，且存在 公因数 \(c > 1\)
则有 \(t | kt / c\)，也就是 存在 \(i\) 使得 \(it = kt / c\)，故 \(a^{ {kt} / {c}} \equiv 1 \pmod p\)
即 \((a ^ k) ^ {t/c} \equiv 1 \pmod p\)，显然 \(t/c < t\)，又 \(t\) 为 \(a ^ k\) 的 阶，矛盾，故 必要性得证
充分性 也可以利用上述结论 证得，考虑因为 \(a\) 模 \(p\) 的 阶 是 \(t\)
则 根据阶的定义，显然，对于一切 \(a ^ i \equiv 1 \pmod p\)，存在 \(t | i\)
我们设 \(a ^ k\) 存在 阶 \(c\)，即有 \((a ^ k) ^ c \equiv a ^ {kc} \equiv 1 \pmod p\)
\(k, t\) 互质，即 \(lcm (k, t) = kt\)，即不存在数 \(c\) 使得 \(c < t\) 且 \(t | ck\)
故 \(c \ge t\)，又显然 \(c = t\) 时，\((a ^ k) ^ c \equiv 1 \pmod p\)，故此时 \(c ~ (t)\) 即是 \(a ^ k\) 模 \(p\) 的 阶
充分性得证
4.2. 对于 \(n\) 次 同余方程 \(f(x) \equiv 0 \pmod p\)，根 至多 有 \(n\) 个
考虑 数学归纳法，一次同余方程 \(ax + b \equiv 0 \pmod p\) 至多一个根
注意到这里 \(p\) 是 质数，故满足 一次同余方程唯一根 的前提，即 \(a, p\) 互质
考虑 \(n\) 次同余方程 \(f(x) \equiv 0 \pmod p\)
若 方程无根，命题成立，反之则 不妨设一根 \(t\)，代入则有
\(f (t) = \sum_i ^ n a_i t ^ i\)，又 原式 \(\sum _ i ^ n a_i x ^ i\)，两式相减有
\(f(x) = \sum _ i ^ n a _ i (x ^ i - t ^ i)\)，显然有 公因式 \((x - t)\)，故原式可化为
\(f (x) = (x - t) g(x)\)，其中 \(g(x)\) 为一个一般的 \(n - 1\) 次多项式
即 \(n\) 次多项式 \(f\) 至多 比 \(n - 1\) 次多项式 \(g\) 多出一个根 \(t\)，由于 一次同余方程 仅 一根
故 原命题得证
4.3. 如果 \(d | p - 1\)，则 同余方程 \(x ^ d \equiv 1 \pmod p\) 恰有 \(d\) 个根
考虑 \(p\) 为质数，故 费马小定理 可知 \(x ^ {p - 1} \equiv 1 \pmod p\)
又由于 \(4.2.\)，该 \(d\) 次 同余方程 \(x ^ d - 1 \equiv 0 \pmod p\) 至多 \(d\) 个根
只需证明其 下界为 \(d\)
又考虑到 同余方程 \(x ^ {p - 1} \equiv 1 \pmod p\) 显然存在 \(p - 1\) 个根
即是 正整数域 内所有 \(p\) 的 可能余数
设 \(p - 1 = qd\)，显然 \(q\) 为 正整数，根据 代数恒等式（等比数列求和逆向）
\(x ^ {p - 1} - 1 = (x ^ d) ^ q - 1 = (x ^ d - 1) ((x ^ d) ^ {q - 1} + (x ^ d) ^ {q - 2} + ... + (x ^ d) ^ 1 + 1)\)
注意到这里 原文老师 又 笔误 将 \(x ^ d - 1\) 写成 \(x ^ q - 1\) 了
即 \(x ^ {p - 1} - 1 = (x ^ d - 1) g (x)\)
而 \(g (x)\) 为 \(dq - d\) 次 多项式，根据 4.2.，其 至多贡献 \(dq - d\) 个根
又 \(x ^ {p - 1} - 1 \equiv 0 \pmod p\) 共有 \(p - 1 = dq\) 个根
故 \(x ^ d - 1 \equiv 0 \pmod p\) 至少贡献 \(d\) 个根，即其 根的数量下界为 \(d\)，原命题得证
4.4. 若 \(a\) 的 阶 为 \(k\)，则 同余方程 \(x ^ k \equiv 1 \pmod p\) 的 所有根 恰为 \(a, a ^ 2, ... a ^ k\)
\(p\) 为质数，显然 \(a, p\) 互质，根据 4.3.，同余方程 \(x ^ k \equiv 1 \pmod p\) 恰有 \(k\) 个根
且显然 \(a ^ k \equiv (a ^ 2) ^ k \equiv ... \equiv (a ^ k) ^ k \equiv 1 \pmod p\)
故只需证明 \(a \% p, (a ^ 2) \% p, ... , (a ^ k) \% p\) 互不相等 即可
考虑 反证法，若存在 \(a ^ i \equiv a ^ j \pmod p，1 \le i < j \le k\)
由于显然 \(a ^ i\) 与 \(p\) 互质，故 根据同余的除法性质 \(a ^ {j - i} \equiv 1 \pmod p\)
又考虑 \(j - i < k\)，而 \(a\) 的 阶 为 \(k\)，矛盾，故原命题得证
由以上四个引理，我们可以证明 4.
考虑我们记 \(n (t)\) 表示 满足 模 \(p\) 的 阶 为 \(t\) 且与 \(p\) 互质 的 \(a\) 的 个数
（显然与 \(p\) 互质只是摆设）而显然，对于所有与 \(p\) 互质的数 \(a\)，\(a\) 模 \(p\) 的 阶 存在且唯一，故 \(\sum n(t) = \varphi (p) = p - 1\)
唯一性 根据 阶的定义可知，存在性则，由于 \(p\) 是质数，根据 费马小定理 可构造
由于 \(p\) 是质数，\(\varphi(p) = p - 1\)，故根据 1.，对于所有 \(t\)，存在 \(t | \varphi (p) = p - 1\)
则显然 \(\sum_{t | (p - 1)} n(t) = \sum n(t)\ = varphi (p) = p - 1\)
而 满足 \(gcd (p - 1, a) = t\) 的 \(a\) 一定共有 \(\varphi ((p - 1) / t)\) 个（\(gcd ((p - 1)/t, a/t) = 1\)）
故 \(\sum_{t | p - 1} \varphi ((p - 1) / t) = \sum_{t | p - 1} \varphi (t) = p - 1\)
即 \(\sum_{t | p - 1} n (t) = \sum_{t | p - 1} \varphi (t)\)
而对于 \(t | p - 1\) 的任意 \(t\)
① 若 \(n(t) = 0\)，则显然 \(n (t) < \varphi (t)\)
② 反之则存在 模 \(p\) 的 阶 为 \(t\) 且与 \(p\) 互质的 \(a\)（此处 \(a\) 代指其一而非整体）
根据 4.4.，\(x ^ t \equiv 1 \pmod p\) 恰有 \(t\) 个解，且为 \(a, a ^ 2, ..., a ^ t\)
又考虑 模 \(p\) 的 阶 为 \(t\) 且与 \(p\) 互质的数 均为同余方程的根，也就是一定在 \(a, a ^ 2, ..., a ^ t\) 中
而根据 4.1.，这样的数应当为 满足 \(i, t\) 互质的 \(a ^ i\)，显然共有 \(\varphi (t)\) 个
而 模 \(p\) 的 阶 为 \(t\) 且与 \(p\) 互质的数 即是 \(n(t)\)，则表明这种情况下 \(n(t) = \varphi(t)\)
故有 当 \(n(t) = 0\) 时，\(n(t) < \varphi(t)\)，而其余情况 \(n(t) = \varphi(t)\)
又知 \(\sum_{t | p - 1} n (t) = \sum_{t | p - 1} \varphi (t)\)
故显然 \(n(t) = 0\) 的情况 不存在（否则 \(\sum_{t | p - 1} n (t) < \sum_{t | p - 1} \varphi (t)\)）
显然我们知道，原根的个数 就是 模 \(p\) 的 阶 为 \(\varphi (p) = p - 1\) 的数的个数，即 \(n(p - 1)\)
由于 \(n(p - 1) \neq 0\) 故此情况下 原根一定存在（这是重点）
又 \(n(p - 1) = \varphi (p - 1) = \varphi (\varphi (p))\) 故此情况下 原根个数为 \(\varphi (\varphi (p))\)，命题得证
回到正题
以下内容部分引用自 [2] 对应文章
来看看在 \(FFT\) 中用到了 单位根 \(w\) 的哪些性质？
- \(w_n ^ k = (w_n ^ 1) ^ k\)
- \(w_n ^ k = w_n ^ {k \% n}\)
- \((w_n ^ k) ^ p = w_n ^ {kp} = w_n ^ {kp \% n}\)
- \(w_n ^ i * w_n ^ j = w_n ^ {i + j}\)
- \(w_{an} ^ {ak} = w_n ^ k \to w_{2n} ^ {2k} = w_n ^ k\)
- \(w_n ^ {k + n/2} = - w_n ^ {k} ~ (2 | n)\)
显然，原根可以在 \(n = p - 1\) 时满足以上性质（通过 余数的性质 易得）
考虑 \(n\) 次单位根 就是把 单位圆 平均分 \(n\) 分
原根 这玩意儿吧就是 在模意义下 把 “单位圆” 平均分 \(p - 1\) 份
然后 模数 也同样具有相似的 可乘性 啥的，就，感性理解是容易的
简单推广得到当 \(n | p - 1\) 的时候也可以满足（几块并成一块）
但是啊，这个 \(p\) 这玩意儿不是 少得可怜？咋做？
仔细想想，发现 由于要 补项，\(n\) 通常情况下可以表示为 \(2 ^ i\) 的形式，也少得可怜
好家伙，于是我们只需要找质数 \(p\) 使得 \(p = 2 ^ k * m + 1\)，然后 \(k\) 巨大就行了
这里的 \(2 ^ k\) 其实决定了 可以做的长度的最大值，然后 我褐了 \(0.125\) 个表（来源 [4]）
| \(p = 2 ^ k * m + 1\) | \(k\) | \(g\) |
|---|---|---|
| \(998244353\) | \(23\) | \(3\) |
| \(1004535809\) | \(21\) | \(3\) |
| \(2013265921\) | \(27\) | \(31\) |
| \(1945555039024054273\) | \(56\) | \(5\) |
| \(4179340454199820289\) | \(57\) | \(3\) |
对于一般题，找到 一个 合适的模数 的 原根 \(g\)，那么 \(g ^ {(p - 1) / n}\) 就是一个 单位根 的完美替代
（可以感性理解成 \(g _ n ^ {(p - 1) / n} = w _ n ^ 1\)）
不要原根
这是从 [2] 直接褐的方法
\(from\) \(UOJ : hly1204\) \(\to\) \(command\_block\)，拜谢两位老师
我们先找到一个质数 \(p = 2 ^ r * s + 1\)，且使得 \(s\) 为 奇数（没有 \(2\) 了）
这时候 随机一个二次非剩余 \(v\)（期望次数 $ = 2$）
这时候如果存在原根 \(g\) 使得 \(v = g ^ k\)，则显然 \(k\) 为 奇数（否则 \(v\) 为 二次剩余）
有 \(v ^ {is} = g ^ {isk}\)，由于 \(g\) 为 原根，\(p\) 为 质数，故有 \(g ^ {p - 1} \equiv 1 \pmod p\)
显然如果 \(v ^ {is} \equiv g ^ {isk} \equiv g ^ {p - 1} \equiv 1 \pmod p\)，则存在 \(p - 1 | isk\)
注意到 \(p - 1 = 2 ^ r * s\)，而 \(s, k\) 均为 奇数，故 \(2 ^ r | i\)，即有 \(v ^ s\) 模 \(p\) 的 阶 为 \(2 ^ r\)
用 \(v ^ s\) 替代 单位根 即可，就是 \((v ^ s) ^ {2 ^ r / n}\) 代替 \(w _ n ^ 1\) 的作用
实现
简单替换一下，直接开搞就行
我们先拉出刚刚 \(FFT\) 的 最终版代码...的 \(FFT\) 函数
inline void FFT (Complex *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1;
Complex Now = Complex (1, 0), tmp, U;
for (int k = 2; k <= Len; len = k, k <<= 1) {
U = Complex (cos (2 * PI / k), Typ * sin (2 * PI / k));
for (int i = 0; i < Len; i += k) {
Now = Complex (1, 0);
for (int l = i; l < i + len; ++ l, Now = Now * U)
tmp = Now * p[l + len], p[l + len] = p[l] - tmp, p[l] = p[l] + tmp;
}
}
}
在 \(FFT\) 的基础上稍做修改
const int MOD = 998244353, g = 3, ig = 332748118;
inline void NTT (int *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1, Now = 1, tmp, u;
for (int k = 2; k <= Len; len = k, k <<= 1) {
u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
for (int i = 0; i < Len; i += k) {
Now = 1;
for (int l = i; l < i + len; ++ l, Now = 1ll * Now * u % MOD)
tmp = 1ll * Now * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD) % MOD, p[l] = (p[l] + tmp) % MOD;
}
}
}
最后 除以 \(n\) 的时候注意 求一下 \(n\) 的逆元 即可，直 接 就 过 了
当然 这玩意儿 是 不优的，显然在 取模 上面 优化空间很大
我们可以用 \(unsigned ~ long ~ long\) 扩大存数范围
减少取模次数
inline void NTT (unsigned long long *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1, Now = 1, tmp, u;
for (int k = 2; k <= Len; len = k, k <<= 1) {
u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
for (int i = 0; i < Len; i += k) {
Now = 1;
for (int l = i; l < i + len; ++ l, Now = 1ll * Now * u % MOD)
tmp = 1ll * Now * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);
}
if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
预处理单位根
也是 十分不错的方法
inline void NTT (unsigned long long *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1, tmp, u;
W[0] = 1;
for (int k = 2; k <= Len; len = k, k <<= 1) {
u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
for (int i = 1; i <= len; ++ i) W[i] = W[i - 1] * u % MOD;
for (int i = 0; i < Len; i += k)
for (int l = i; l < i + len; ++ l)
tmp = 1ll * W[l - i] * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);
if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
最后也放一份比较完整的
板子
2024-02-24 14:39 - C++20 O2 - 1.63 KB / 981 ms / 50.55 MB | Record
#include <bits/stdc++.h>
#define Complex complex <double>
const int MAXN = 2100005;
const int MOD = 998244353;
using namespace std;
int U[MAXN];
inline int Qpow (const int a, const int b = MOD - 2) {
if (b == 1) return a;
int Ret = Qpow (a, b >> 1);
if (b & 1) return 1ll * Ret * Ret % MOD * a % MOD;
return 1ll * Ret * Ret % MOD;
}
const int g = 3, ig = Qpow (g);
unsigned long long W[MAXN];
inline void NTT (unsigned long long *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1, tmp, u;
W[0] = 1;
for (int k = 2; k <= Len; len = k, k <<= 1) {
u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
for (int i = 1; i <= len; ++ i) W[i] = W[i - 1] * u % MOD;
for (int i = 0; i < Len; i += k)
for (int l = i; l < i + len; ++ l)
tmp = 1ll * W[l - i] * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);
if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
unsigned long long F[MAXN], G[MAXN];
int N, M, P, x, in;
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> N >> M;
for (P = 1; P <= N + M; P <<= 1);
for (int i = 0; i < P; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? P >> 1 : 0);
for (int i = 0; i <= N; ++ i) cin >> F[i];
for (int i = 0; i <= M; ++ i) cin >> G[i];
NTT (F, P, 1), NTT (G, P, 1);
for (int i = 0; i < P; ++ i) F[i] = 1ll * F[i] * G[i] % MOD;
NTT (F, P, - 1), in = Qpow (P);
for (int i = 0; i <= N + M; ++ i) cout << (1ll * F[i] * in) % MOD << ' ';
return 0;
}
多模 \(NTT\)
用于处理 巨大值域 的情况，当使用 __int128 类型时可以做到 \(10 ^ {27}\) 的值域
保留 \(NTT\) 的优良传统，没有误差，但是由于需要 跑多次，常数就很大
值域一般是 多项式长度 乘上 （单项系数最大值 的 平方）
思路是非常简单的，就是用多个 较小模数 分别做一次 \(NTT\)
一般使用 \(3\) 个 \(10 ^ 9\) 级别的模数
然后直接使用 \(CRT\) 合并出每项系数即可，支持的值域就是模数的积
下面给一个 任意模数多项式乘法 这个题的板子，常数很大，但完全能过
就是多模 \(NTT\) 求出精确值之后再按题目要求取模即可
#include <bits/stdc++.h>
const int MAXN = 300005;
// 这里的 模数 可以直接改改，比如上点啥 998244353，1004535809，469762049 的
const unsigned long long MOD1 = 104857601;
const unsigned long long MOD2 = 23068673;
const unsigned long long MOD3 = 167772161;
using namespace std;
static unsigned long long MOD = 1;
inline unsigned long long qpow (const unsigned long long a, const unsigned long long b = MOD - 2, const __int128 mod = MOD) {
if (b == 0) return 1;
if (b == 1) return a;
unsigned long long Ret = qpow (a, b >> 1, mod);
if (b & 1) return (unsigned __int128) Ret * Ret % mod * a % mod;
return (unsigned __int128) Ret * Ret % mod;
}
int U[MAXN];
unsigned long long W[MAXN];
unsigned long long g = 3, ig;
inline void Init (const int N) {
memset (U, 0, sizeof U);
for (int i = 0; i < N; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? N >> 1 : 0);
}
inline void NTT (unsigned long long *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1, tmp, u;
W[0] = 1;
for (int k = 2; k <= Len; len = k, k <<= 1) {
u = qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
for (int i = 1; i <= len; ++ i) W[i] = W[i - 1] * u % MOD;
for (int i = 0; i < Len; i += k)
for (int l = i; l < i + len; ++ l)
tmp = 1ll * W[l - i] * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);
if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
unsigned long long MI[5], Ans[MAXN];
__int128 T[MAXN], MP = 1;
inline void Mul (const unsigned long long * A, const unsigned long long * B, const int N, const int M, const int Mod) {
unsigned long long F[(N + M) << 1] = {0}, G[(N + M) << 1] = {0}, MI = MP / Mod;
int I, P;
memcpy (F, A, N << 3), memcpy (G, B, M << 3);
for (int i = 0; i < N; ++ i) F[i] = A[i] % Mod;
for (int i = 0; i < M; ++ i) G[i] = B[i] % Mod;
MOD = Mod, g = 3, ig = qpow (g);
for (P = 1; P < (N + M); P <<= 1);
Init (P), NTT (F, P, 1), NTT (G, P, 1);
for (int i = 0; i < P; ++ i) F[i] = F[i] * G[i] % Mod;
NTT (F, P, -1), I = qpow (P);
for (int i = 0; i < P; ++ i) Ans[i] = F[i] * I % Mod;
for (int i = 0; i < P; ++ i) T[i] += (__int128) Ans[i] * MI % MP * qpow (MI, Mod - 2, Mod) % MP, T[i] = T[i] % MP;
}
int N, M, P, I;
unsigned long long mod;
unsigned long long A[MAXN], B[MAXN];
unsigned long long F[MAXN], G[MAXN];
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> N >> M >> mod, ++ N, ++ M;
for (int i = 0; i < N; ++ i) cin >> A[i];
for (int i = 0; i < M; ++ i) cin >> B[i];
MP = (__int128) MOD1 * MOD2 * MOD3;
Mul (A, B, N, M, MOD1), Mul (A, B, N, M, MOD2), Mul (A, B, N, M, MOD3);
for (int i = 0; i < N + M - 1; ++ i) cout << (long long) (T[i] % mod) << ' ';
return 0;
}
多项式基本运算
以下内容多引用至 [2] 对应的文章
总算写完 多项式乘法 的部分了...
注意笔者这里将（暂时？）鸽掉 任意模数 \(NTT\)，\(MTT\)，这是因为 笔者还不会
为啥上面的 乘法 要花 那么多篇幅来讲？
因为 下面的 所有基本运算 都将以 多项式乘法 为基础
界
\(nxbj\)，在讲运算之前，先来了解一下 界 这个概念
显然，多数 多项式题 有效项 是有限的，为了保证 正确的复杂度
我们人为给 所有运算 设定了一个 次数上界，即 \(\pmod {x ^ n}\)，这就是 界
同 整数的模 相似，我们可以简单证明到 多项式加法 / 乘法 在 界限制 下满足如下性质
由于 以下所有基本运算 都将以 多项式加法 / 乘法 为基础定义的
故我们的运算可以 时刻保持 这个次数上界，以达到 复杂度正确 的目的
一些基本函数
\(nxbj\)，我们需要 封装一些 基础的操作 来保证后面的 代码复杂度不会爆炸
清空区间 和 两数组相乘
inline void Clr (unsigned long long *f, int L, int R) {
for (int i = L; i < R; ++ i) f[i] = 0;
}
inline unsigned long long * Cross (unsigned long long *F, unsigned long long *G, int N) {
for (int i = 0; i < N; ++ i) F[i] = F[i] * G[i] % MOD;
return F;
}
多项式加法
inline unsigned long long * Add (unsigned long long *F, unsigned long long *G, int N, int M) {
for (int i = 0; i <= max (N, M); ++ i) F[i] = (F[i] * G[i]) % MOD;
return F;
}
多项式乘法
inline unsigned long long * Mul (unsigned long long *F, unsigned long long *G, int N, int M) {
int P = 1, I;
for (P = 1; P < N + M; P <<= 1);
Init (P), NTT (F, P, 1);
if (F != G) NTT (G, P, 1);
Cross (F, G, P), NTT (F, P, -1), I = Qpow (P);
for (int i = 0; i <= N + M; ++ i) F[i] = F[i] * I % MOD;
return F;
}
多项式求逆
注意到这部分下面讲的 \(R'(x)\) 不是 \(R(x)\) 的 导数
先把四则运算给补齐了... 这里的除法 不带余数
来考虑一个 倍增 的方法
设我们现在有 \(R'(x) = F(x) ^ {-1} \pmod {x ^ n}\)，而我们想求 \(R(x) = F(x) ^ {-1} \pmod {x ^ n}\)
也就是我们想要把 \(\pmod {x ^ {n / 2}}\) 倍增到 \(\pmod {x ^ n}\)
最直接能想到的方法就是 直接平方，可以注意到 一定有 \(R(x) = F(x) ^ {-1} \pmod {x ^ n}\)
故而构造 \(R(x) - R'(x) = 0 ~ (mod x ^ {n / 2})\)，然后 两边平方
于是 \((R(x) - R'(x)) ^ 2 = 0 \pmod {x ^ n}\)，展开有 \(R(x) ^ 2 - 2 R(x) R'(x) + R'(x) ^ 2 = 0 \pmod {x ^ n}\)
我们是知道 \(R'(x)\) 的，所以 \(R'(x), R'(x) ^ 2\)都是好求的
而我们需要 \(R(x)\)，这东西从 \(R(x) ^ 2\) 中显然 不易获得，所以我们需要 消掉一个 \(R(x)\)
注意到 \(R(x) = F (x) ^ {-1} \pmod {x ^ n}\)，显然有 \(F(x)R(x) = 1 \pmod {x ^ n}\)
故给 展开式 左右同时 乘上一个 \(F(x)\)，有 \(R(x) - 2R'(x) + F(x)R'(x) ^ 2 = 0 \pmod {x ^ n}\)
移项即可，最终即得 \(R(x) = 2R'(x) - F(x)R'(x) ^ 2\)，倍增即可
上面这个式子 直接用确实没有问题，但是其实 暗藏玄机
我们发现由于 \(R (x) \equiv R’ (x) \pmod {x ^ {n / 2}}\)
那么其实 \(R(x)\) 的前面 \(n / 2\) 项并 不会在迭代的时候更新
（就是下面代码中当枚举到 长度为 \(k\)，最后统计答案时其实前 \(k / 2\) 项是 不变的）
而枚举到 界 大于 \(x ^ {n / 2}\) 时，由于 \(R' (x)\) 只有 \(n / 2\) 项，故其实它对应项系数为 \(0\)
也就是在程序中，从系数上看，我们实际可以得到 分两段的式子
这里设 \(R(x) = \sum r _ i x ^ i\)，\(R’ (x) = \sum r’ _i x ^ i\)，\(F (x) R'(x) ^ 2= \sum fr _i x ^ i\)
当 \(0 \le i < n / 2\) 时，有 \(r _ i = r'_ i\)
当 \(n / 2 \le i < n\) 时，有 \(r _ i = fr _ i\)
时间复杂度 \(T(n) = T(n / 2) + O(n \log n) = O(n \log n)\)（这里 \(O(n \log n)\) 在于 \(F(x)R'(x) ^ 2\)）
这一部分的代码，先有一个简单易懂的版本（没有玄机）
unsigned long long TmpA[MAXN], TmpB[MAXN], TmpC[MAXN];
inline unsigned long long * Inv (unsigned long long *F, int N) {
int P = 1, I;
for (P = 1; P < N; P <<= 1);
unsigned long long *A = TmpA, *B = TmpB, *C = TmpC;
A[0] = Qpow (F[0]);
for (int k = 2; k <= P; k <<= 1) {
memcpy (B, F, k << 3), memcpy (C, A, k << 2);
Clr (A, k >> 1, k), Init (k << 1);
NTT (A, k << 1, 1), NTT (B, k << 1, 1);
Cross (B, A, k << 1), Cross (B, A, k << 1);
NTT (B, k << 1, -1), I = Qpow (k << 1);
for (int i = 0; i < k; ++ i) B[i] = B[i] * I % MOD;
for (int i = 0; i < k; ++ i) A[i] = ((C[i] << 1) - B[i] + MOD) % MOD;
}
for (int i = 0; i <= N; ++ i) F[i] = A[i];
return F;
}
上面这个代码就是 完完全全按照 式子来的，常 数 稍 大，但其实可以接受
inline unsigned long long * Inv (unsigned long long *F, int N) {
int P = 1, I;
for (P = 1; P < N; P <<= 1);
unsigned long long *A = TmpA, *B = TmpB, *C = TmpC;
A[0] = Qpow (F[0]);
for (int k = 2; k <= P; k <<= 1) {
memcpy (B, F, k << 3), memcpy (C, A, k << 2);
Clr (A, k >> 1, k << 2), Mul (B, Mul (A, A, k, k), k << 1, k << 1);
for (int i = 0; i < k; ++ i) A[i] = ((C[i] << 1) - B[i] + MOD) % MOD;
}
for (int i = 0; i <= N; ++ i) F[i] = A[i];
return F;
}
你也可以利用 上面写的乘法，得到一个 常 数 更 大 的 写 法
如果对常数有一定要求，那么也可以这样写（利用上面讲的性质）
inline unsigned long long * Inv (unsigned long long *F, int N) {
int P = 1, I;
for (P = 1; P < N; P <<= 1);
unsigned long long A[MAXN], B[MAXN], C[MAXN];
A[0] = Qpow (F[0]);
for (int k = 2; k <= P; k <<= 1) {
memcpy (B, F, k << 3);
memcpy (C, A, k << 2);
Clr (A, k >> 1, k), Init (k), NTT (A, k, 1), NTT (B, k, 1);
Cross (B, A, k), NTT (B, k, -1), I = Qpow (k);
for (int i = 0; i < k; ++ i) B[i] = B[i] * I % MOD;
Clr (B, 0, k >> 1), NTT (B, k, 1), Cross (B, A, k), NTT (B, k, -1);
for (int i = 0; i < k; ++ i) B[i] = B[i] * I % MOD;
for (int i = 0; i < (k >> 1); ++ i) A[i] = C[i];
for (int i = (k >> 1); i < k; ++ i) A[i] = (MOD - B[i]) % MOD;
}
for (int i = 0; i <= N; ++ i) F[i] = A[i];
return F;
}
会快一点，但挺多的，主要是 实现的很好，但是 Inv::C[MAXN] 还是可以省掉的
考虑到 \(R(x) - R'(x) = 0 \pmod {x ^ \frac {k}{2}}\)
这意味着每次倍增时 只有 \((k / 2,k]\) 区间内的项 真正被改变了，前面 \(k / 2\) 项是保持不变的
于是对于 \(F(x) R'(x) ^ 2\) 这个 需要乘法的部分， 我们可以像上面那样 分两次算
而非直接 Cross (B, A, k << 1), Cross (B, A, k << 1) 一次解决
我们发现 两次计算 都只对前 \(k\) 次（而非 \(2k\)）做 \(NTT\)
这样虽然做了 \(5\) 次 \(NTT\)，但是比之 \(3\) 次，每次长为 \(2k\) 的做法，常数更小
而每次 前面 \(k / 2\) 项的答案保持不变，故我们 只关心后面部分
所以 \(k\) 次多项式 和 \(k / 2\) 次多项式 相乘，最终取 \(k\) 次，
由于 循环卷积 的特性，我们只可能 破坏前面 \(k / 2\) 项，也就是我们不需要的部分
\(UPD:\) 循环卷积 即 真实积的第 \(t\) 次项 在 以 \(k\) 为长度 的 变换后 会贡献到 \(t \% k\) 项上
证明考虑 单位根的性质 即可，即 \(w_n ^ t = w_n ^ {t \% k}\)
但是我们这里不能直接像 第一种写法 一样连续
Cross两次因为 \(NTT\) 的时候其实 前 \(k / 2\) 项已经被破坏了
也就是说相当于我第一次
Cross完，得到 \(R'(x)F(x)\) 其实 只有后 \(k / 2\) 项是对的前面 \(k / 2\) 项是 错误的，所以我们得先 转化为系数，把前面的 \(k / 2\) 项 清掉 再做
否则在错误的结果上再做一次
Cross之后再 \(NTT\)，就会直接寄掉
而不会影响答案，所以 正确性是有保证 的
多项式牛顿迭代
结论部分 请直接转到标题 “在多项式上应用”
以下内容多引用自 [5] 对应的文章
\(Newton's ~ Method\) 是 牛顿 提出的一种将 非线性方程线性化 的 近似方法。
它也可以 运用在多项式 中，求关于多项式的非线性方程在 模意义下的解。
泰勒级数
我们得先了解 泰勒级数 这个东西
它用一个 无限项连加式 表示（拟合）光滑函数 \(f(x)\)，有
其中 \(f ^ {(n)}\) 表示 \(f\) 的 \(n\) 阶 导数
而上式被称作 \(f(x)\) 在 \(a\) 处的 泰勒展开式，而 等号右边 的式子就是 \(f(x)\) 在 \(a\) 处的 泰勒级数
特别的，在 \(0\) 处的 泰勒展开式 / 泰勒级数 被称为 麦克劳林展开式 / 麦克劳林级数
牛顿迭代
然后来看 \(Newton's ~ Method\)（牛顿迭代）
它一般被用于（程式化的）求解非线性方程，以下是 求解 \(f(x) = 0\) 的根的 算法流程
-
选取一个（合适的）数 \(x_0\)
-
将 \(f(x)\) 在 \(x_0\) 处 泰勒展开为 \(f(x) = \sum_{i = 0} ^ {+∞} \dfrac {f ^ {(i)} (x_0)} {i!} (x - x_0) ^ i\)
-
只取 \(0,1\) 次项的 系数，令之为 \(0\)，即 \(\dfrac {f(x)} {0!} (x - x_0) ^ 0 + \dfrac {f'(x)} {1!} (x - x_0) ^ 1 = 0\)
-
化简为 \(f(x_0) + f'(x_0) (x - x_0) = 0\) 并求解得到根 \(x_1\)
-
回到 \(1.\)，并令 \(x_0 = x_1\)，继续迭代，可以证明 每次得到的解都更接近 \(f(x) = 0\) 的精确根
我们把 \(4.\) 中的式子化简后就可以得到一个好用的 递推式 \(x_{i + 1} = x_i - \dfrac {f (x _i)} {f'(x_i)}\)
已经证明，\(Newton's ~ Method\) 具有 平方收敛 的性能，也就是其 每次迭代，有效数字翻倍
这里提供一种 来源于 [6] 的 感性理解，个人感觉那篇博客写的很好，可以去原文看
可以注意到 取展开后 \(0, 1\) 次项系数得到的直线 就相当于 \(f(x)\) 在 \(a\) 处的 切线
于是 算法流程 在几何意义上就是
找在（任意一个）\(a\) 处的切线
求切线与 \(y = 0\) 交点
在交点 \(x_0\) 处作垂线 \(x = x_0\)，与原函数交于 \((x_0, f(x_0))\)
回到 \(1.\)，令 \(a = (x_0, f(x_0))\)
画图模拟 上述过程将会发现
切线交点 \((x_0, 0)\) 在 原函数是一般光滑函数 的时候会 不断接近原函数 \(f(x) = 0\) 的根
代数上的 严谨证明 似乎是 困难的，这里不做进一步说明（
不会，鸽了，会的快给我讲)
根据上面的 递推过程，我们容易得到下面的 代码框架
inline double f (double x) {} // 函数单点求值
inline double fd (double x) {} // 导函数单点求值
inline double Solve (double x, int dep) {
while (-- dep) x = x - f(x) / fd(x);
return x;
}
多项式积分 / 求导
主要是定义，下面计算中会用到这些东西，算法没有难度，时间都是 \(O(N)\) 的
多项式求导
\(F'(x) = \sum_{i = 1} i F_i x ^ {i - 1}\)
inline unsigned long long * Dlt (unsigned long long *F, int N) {
for (int i = 1; i < N; ++ i) F[i - 1] = F[i] * i % MOD;
F[N - 1] = 0;
return F;
}
多项式积分（需要处理 整数逆元）
\(F'(x) = Const + \sum_{i = 0} \dfrac {F_i x ^ {i + 1}} {i}\)
inline unsigned long long * Sum (unsigned long long *F, int N) {
for (int i = N; i; -- i) F[i] = F[i - 1] * inv[i] % MOD;
F[0] = 0;
return F;
}
在多项式上的应用
已知 函数 \(G\) 使得 \(G(F(x)) = 0\)，求 多项式 \(F \pmod {x ^ n}\)
公式就是把上面的 递推式 的 \(x\) 换成 多项式 \(F(x)\)（\(F ^ * (x) = F(x) \pmod {x ^ n}\)）
虽然 \(Newton's ~ Method\) 是一种 求近似解 的方法，
但是我们本身需要的也只是 模意义下的 “精确解”，所以 不用担心正确性问题
引用 [7] 中的证明如下
把 \(G(F(x))\) 在 \(F ^ * (x)\) 处 泰勒展开 有
\(G(F(x))\) = \(\sum_{i = 0} ^ {+ ∞} \dfrac {G(F ^ *(x))} {i!} (F(x) - F ^ * (x)) ^ i\)
又注意到 \(F(x) - F ^ *(x) = 0 \pmod {x ^ n}\)
于是 \(F(x) - F ^ *(x) \pmod {x ^ n}\) 的 最低次项 至少为 \(x ^ {n / 2}\)
故当展开式 \(i \ge 2\) 时，系数均为 \(0\)（\((F(x) - F ^ * (x)) ^ i = 0 \pmod {x ^ n}\)）
即有 \(G(F(x)) = G(F ^ *(x)) + G' (F ^ * (x)) (F(x) - F ^ *(x))\)，整理可得上式
然后 倍增 即可解决问题，时间复杂度 \(T(N) = T (N / 2) + O(N \log N) = O(N \log N)\)
然后以下 很多运算的推导中 都可能会用到这个方法
多项式开根
牛顿迭代，启动！
设 给定 多项式 \(F(x)\)，存在多项式 \(W (x)\)，使得 \(W (x) ^ 2 = F(x)\)
则 令 \(G(W(x)) = F(x) - W(x) ^ 2 = 0\)，可得 \(G'(W(x)) = - 2 W (x)\)
\(W(x)\) 为未知元，\(F(x)\) 为常数项
根据 牛顿迭代公式，有 \(W(x) = W ^ * (x) - \dfrac {F(x) - W ^ * (x) ^ 2} {- 2 W ^ * (x)}\)
当多项式 只剩常数项时，开方结果为 二次剩余，而 板子题 保证 常数为 \(1\) ，\(W(0) = 1\) 即可
inline unsigned long long * Sqrt (unsigned long long *F, int N) {
int P = 1;
for (P = 1; P <= N; P <<= 1);
unsigned long long B[MAXN], C[MAXN];
B[0] = 1;
for (int k = 2; k <= P; k <<= 1) {
for (int i = 0; i < (k >> 1); ++ i) C[i] = (B[i] << 1) % MOD;
Clr (C, k >> 1, k), Inv (C, k), Clr (C, k, k << 1), Mul (B, B, k, k);
for (int i = 0; i < k; ++ i) B[i] = (F[i] + B[i]) % MOD;
Mul (B, C, k, k), Clr (B, k, k << 2);
}
for (int i = 0; i <= N; ++ i) F[i] = B[i];
return F;
}
常 数 还 行，不想写 只用 \(NTT\) 的版本...
\(2024.05.06 ~ UPD:\) 当前的 提交记录 是 未取等版本，有正确性问题，但是上面 代码片段 无误
$Hack : $
17 1 41 81 83 65 16 75 61 98 33 67 72 18 85 63 92 62
如果这个时候 不保证常数为 \(1\)，那么就对 常数项 找一次 二次剩余
考虑一共只需要 找一次，所以时间不卡，可以使用 \(BSGS\) 在 \(O (\sqrt {p})\) 的时间做掉
甚至可以 开 \(k\) 次根
具体而言，如果存在 原根 \(g\)，常数 \(c\)，我们要找 \(c\) 的二次剩余
那么先找到 \(i\) 使得 \(g ^ i \equiv c \pmod p\)，于是 \(g ^ {i / 2} \pmod p\) 就是所求
\(k\) 次根 就是 \(g ^ {i / k}\)
而 \(BSGS\) 就是用来优化找 \(i\) 过程的，不然暴力枚举显然是 \(O(p)\) 的复杂度，爆炸
我们考虑设 \(g ^ i = g ^ {At - B}\)，于是有 \(g ^ {At} \equiv g ^ B c \pmod p\)，取 \(t = \lceil \sqrt {p} \rceil\)
枚举可能的 \(t\) 种 \(B\)，求出 \(g ^ B c \pmod p\)，将对应的 \(B\) 放入 std::unordered_map
然后枚举 \(t\) 中 \(A\)，求出 \(g ^ {At}\)，查询是否在 std::unordered_map 中存在
存在时的 \(At - B\) 即为所求，而总体时间复杂度 \(O(\sqrt {p}) + O(N \log N) = O (N \log N)\)
具体代码
inline unsigned long long Sqrt_BSGS (unsigned long long x) {
unordered_map <unsigned long long, int> mp;
int t = sqrt (MOD) + 1;
unsigned long long Now = 1, gg = g;
for (int b = 1; b <= t; ++ b) Now = (Now * gg) % MOD, mp[(Now * x) % MOD] = b;
gg = Now, Now = 1;
for (int a = 1; a <= t; ++ a) {
Now = (Now * gg) % MOD;
if (mp[Now]) return a * t - mp[Now];
}
return -1;
}
unsigned long long Sqrt_Mod (unsigned long long a, unsigned long long k = 2) {
int r = Qpow (g, Sqrt_BSGS (a) / k);
return r = min (r, MOD - r);
}
inline unsigned long long * Sqrt_2 (unsigned long long *F, int N) { //DNF
int P = 1, I = Qpow (2);
for (P = 1; P < N; P <<= 1);
unsigned long long A[MAXN], B[MAXN], C[MAXN];
B[0] = Sqrt_Mod (F[0], 2);
for (int k = 2; k <= P; k <<= 1) {
for (int i = 0; i < (k >> 1); ++ i) C[i] = B[i] % MOD;
Clr (C, k >> 1, k), Inv (C, k), Clr (C, k, k << 1);
memcpy (A, F, k << 3), Mul (C, A, k, k);
for (int i = 0; i < (k >> 1); ++ i) C[i] = (C[i] + B[i]) % MOD;
for (int i = 0; i < k; ++ i) B[i] = C[i] * I % MOD;
}
for (int i = 0; i <= N; ++ i) F[i] = B[i];
return F;
}
多项式带余除法
考虑 乘法逆元 实际上求的是 模意义下的精确解（因为无论如何，逆元乘原式 均等于 \(1\)）
就像 实数除法里面得到的分数结果一样，是准确的，没有余数
所以想要得到 整除后的 余式，得想个新的方法
我们能想到 先求出整除后的结果式，乘上除式，用被除式减去这个积，从而得到 余式
所以我们得先把 余式 给搞掉，使得原式能被整除
又考虑到 余式次数很低（\([0, x ^ {n - m + 1}]\)），容易在 界 内被保留下来
于是，来个 神秘操作，我们 反转一下 这个多项式，即 \(F ^ T (x) = x ^ n F(x ^ {-1})\)
从程序上讲，就是 反转系数数组，然后我们来推一下，消掉余数
设存在 \(F(x) = G(x) * T(x) + R(x)\)，我们已知 被除式 \(F(x)\)，除式 \(G(x)\)，分别为 \(N, M\) 次
于是考虑等式两边同时乘上 \(x ^ N\)，有 \(x ^ N F (x) = x ^ N G(x) * T (x) + x ^ N R (x)\)
注意到 \(G(x) * T(x)\) 的 最高次项 应当为 \(x ^ N\)，即其 次数和 为 \(x ^ N\)
而 余式 \(R(x)\) 次数不超过 除式次数 \(M\)，\(R(x)\) 应当为 \(M - 1\) 次，故可以 反转多项式 得
\(F ^ T (x) = G ^ T (x) * T ^ T (x) + x ^ {N - M + 1} R ^ T (x)\)，此时 \(\bmod x ^ {N - M + 1}\) 就有
\(F ^ T (x) = G ^ T (x) * T ^ T (x) \pmod {x ^ {N - M + 1}}\)，然后 直接做求逆 + 乘法 就可以算出 \(T ^ T (x)\)
时间复杂度 \(O (N \log N)\)
反转得到 \(T (x)\) ，乘 \(G(x)\)，再用 \(F(x)\) 减去 这个积，就可以得到余数 \(R(x)\)，代码
inline unsigned long long * Mof (unsigned long long * F, unsigned long long * G, int N, int M) {
unsigned long long A[MAXN], B[MAXN];
int P = N - M + 1;
Rev (F, 0, N), memcpy (A, F, P << 3), Rev (F, 0, N);
Rev (G, 0, M), memcpy (B, G, P << 3), Rev (G, 0, M);
Inv (B, P), Mul (B, A, P, P), Rev (B, 0, P), Clr (B, P, P << 1)
memcpy (A, B, P << 3), Mul (G, B, N, N), Clr (G, N, N << 1);
for (int i = 0; i < M - 1; ++ i) G[i] = (F[i] - G[i] + MOD) % MOD;
for (int i = 0; i < P; ++ i) F[i] = A[i];
Clr (G, M - 1, M << 1), Clr (F, P, P << 1);
return F;
}
多项式对数函数（\(\ln\)）
以下部分引用自 [8] 对应的文章
用 麦克劳林级数 来定义，因为 多项式的自然对数 无法在模意义下定义...
不知道 麦克劳林级数 是什么的 转到 多项式牛顿迭代 - 泰勒级数
看着就很恶心... 但这玩意儿 只是看个乐子，真正要求解 还得用 求导
考虑 \(\ln ' (x) = \dfrac {1} {x}\)，根据 链式法则 \(f(g(x))' = g'(x)f'(g(x))\)（一层一层 导出来）
于是 \(\ln (F(x))' = F'(x) \ln' (F(x)) = \dfrac {F'(x)} {F(x)}\)，直接做就行，时间 \(O (N \log N)\)
板子题 保证了 \(F_0 = 1\)，故 \(\ln (F(x))_0 = 0\) 即可，否则 \(\ln (F(x))_0 = \ln (F_0)\) 是 困难的
inline unsigned long long * In (unsigned long long *F, int N) {
unsigned long long A[MAXN];
memcpy (A, F, N << 3);
Inv (A, N), Dlt (F, N), Mul (F, A, N, N);
Sum (F, N - 1), Clr (F, N, N << 1);
return F;
}
多项式指数函数（\(\exp\)）
以下部分引用自 [8] 对应的文章
同样用 麦克劳林级数 定义，如下
然后来考虑 怎么求，设 \(T(x) = \exp (F(x))\)，要求 \(T(x)\)，两边同时 \(\ln\) 一下
于是 \(\ln (T(x)) = F(x)\)，还是不好搞，但是我们可以 设 \(G(T(x)) = \ln (T(x)) - F(x) = 0\)
这下就能用 牛顿迭代 了，于是有 \(T(x) = T ^ *(x) - \dfrac {G (T ^ * (x))} {G' (T ^ * (x))}\)
然后 \(G'(T ^ * (x)) = \ln'(T ^ * (x)) = \dfrac {1} {T ^ *(x)}\)
代入一下有 \(T(x) = T ^ * (x) - \dfrac {\ln (T ^ *(x)) - F (x)} {\dfrac {1} {T ^ *(x)}} = T ^ * (x) (1 - \ln(T ^ * (x)) + F (x))\)
然后又是 直接做时间
时间 \(T(N) = T(N / 2) + O (N \log N) = O (N \log N)\)
板子题 保证了 \(F_0 = 0\)，故 \(\exp (F(x))_0 = 1\) 即可，否则 \(\exp (F(x))_0 = e ^ {F_0}\) 是 困难的
inline unsigned long long * Exp_1 (unsigned long long *F, int N) {
int P = 1;
for (P = 1; P < N; P <<= 1);
unsigned long long A[MAXN], B[MAXN];
A[0] = 1;
for (int k = 2; k <= P; k <<= 1) {
memcpy (B, A, k << 2), Clr (B, k >> 1, k), Ln (B, k);
for (int i = 0; i < k; ++ i) B[i] = (- B[i] + F[i] + MOD) % MOD;
B[0] = (B[0] + 1) % MOD, Mul (A, B, k, k), Clr (A, k, k << 1);
}
for (int i = 0; i < N; ++ i) F[i] = A[i];
return F;
}
常 数 还 行，有没有好心老哥可以帮忙看看？
注意 用此处的 \(\exp\) 之前 切记预处理整数逆元
这道题还有一种 分治 \(FFT\) 的做法，但是笔者还不会，就先鸽了，可以参考 [9] 来学习
多项式三角函数
以下内容多引用自 [10] 对应的文章
三角恒等变形
考虑 欧拉公式（此处没有证明）
于是变形有
两式 和差
把 \(x\) 替换成 \(F(x)\) 即可，一点不影响
考虑这里有 \(i\)，令人难受，但我们注意到 \(i ^ 2 = -1\)，也就是 \(i ^ 2 \equiv p - 1 \pmod p\)
于是 做一次 \(p - 1\) 的 二次剩余 即可，用前面的 \(BSGS\)，时间不卡
最后用 \(\exp\) 做就行了，注意 预处理整数逆元，代码片段
inline unsigned long long * Sin (unsigned long long *F, int N) {
unsigned long long I = Sqrt_Mod (MOD - 1), A[MAXN];
for (int i = 0; i < N; ++ i) F[i] = F[i] * I % MOD;
init (N + 4), Exp_1 (F, N);
memcpy (A, F, N << 3), Inv (A, N), I = Qpow (I << 1);
for (int i = 0; i < N; ++ i) F[i] = I * (F[i] - A[i] + MOD) % MOD;
return F;
}
inline unsigned long long * Cos (unsigned long long *F, int N) {
unsigned long long I = Sqrt_Mod (MOD - 1), A[MAXN];
for (int i = 0; i < N; ++ i) F[i] = F[i] * I % MOD;
init (N + 4), Exp_1 (F, N);
memcpy (A, F, N << 3), Inv (A, N), I = Qpow (2);
for (int i = 0; i < N; ++ i) F[i] = I * (F[i] + A[i]) % MOD;
return F;
}
多项式快速幂
保证了 \(F_0 = 1\)，这是简单的，显然我们有 \(F(x) ^ k = e ^ {\ln (F (x)) * k}\) 这种 式子
于是一个 \(\exp\)，一个 \(\ln\) 就做掉了，时间复杂度 \(O(N \log N)\)，代码
inline unsigned long long * Qpow_1 (unsigned long long *F, int N, int K) {
Ln (F, N);
for (int i = 0; i < N; ++ i) F[i] = F[i] * K % MOD;
Exp_1 (F, N);
return F;
}
暴力 \(NTT\) 也是可以的！时间 \(O (N \log ^ 2 N)\)，代码
inline unsigned long long * Qpow_2 (unsigned long long *F, int N, int K) {
unsigned long long A[MAXN], B[MAXN], C[MAXN];
Clr (A, 0, N), Clr (B, 0, N), Clr (C, 0, N), A[0] = 1, memcpy (B, F, N << 3);
while (K > 0) {
if (K & 1) -- K, memcpy (C, B, N << 3), Mul (A, B, N, N), memcpy (B, C, N << 3), Clr (A, N, N << 1), Clr (B, N, N << 2);
Mul (B, B, N, N), Clr (B, N, N << 1), K >>= 1;
}
for (int i = 0; i < N; ++ i) F[i] = A[i];
return F;
}
还有当不保证 \(F_0 = 1\) 时，我们需要操作一下
如果 \(F_0\) 不为 \(0\)，其实直接 所有项除 \(F_0\)，正常快速幂，最后 乘个 \(F_0 ^ k\) 回来 就行
但是 \(F_0\) 为 \(0\) 呢？考虑我们去找 最小的 \(i\)，使得 \(F_i \neq 0\)，然后 除了再乘
相当于 平移了 \(i\) 项，注意到根据 费马小定理，我们的 \(k\) 次方 在 整数 时应当 \(\bmod p - 1\)
在多项式时应当 \(\bmod p\)，所以需要处理出 两种 \(k\)，具体见代码（这东西是真的恶心）
inline unsigned long long * Qpow_3 (unsigned long long *F, int N, long long K, long long K_, long long Siz) {
int I = 1, p = 0, c = 1;
while (F[p] == 0 && p < N) ++ p;
c = Qpow (F[p], K_), I = Qpow (F[p], MOD - 2);
if ((p && Siz > 6) || p * K > N) return Clr (F, 0, N), F;
for (int i = 0; i + p < N; ++ i) F[i] = F[i + p] * I % MOD;
Ln (F, N);
for (int i = 0; i < N; ++ i) F[i] = F[i] * K % MOD;
Exp_1 (F, N);
for (long long i = N - 1ll - 1ll * K * p; i >= 0; -- i) F[i + 1ll * p * K] = F[i] * c % MOD;
Clr (F, 0, min (1ll * p * K, 1ll * N));
return F;
}
in main :
for (int i = S.size () - 1, j = 1; i >= 0; -- i, j = 10ll * j % MOD)
K = (1ll * (S[i] - '0') * j + K) % MOD;
for (int i = S.size () - 1, j = 1; i >= 0; -- i, j = 10ll * j % (MOD - 1))
K_ = (1ll * (S[i] - '0') * j + K_) % (MOD - 1);
Qpow_3 (F, N, K, K_, S.size ());
\(OK, OK\) 基本运算搞定了，题后面再补
综述
给个表格，方便快速的 查阅 / 复习，或者找到对应的题
以下 \(n, N\) 均代表多项式项数，\(F ^ *\) 表示界为多项式 \(F\) 一半，其余相同的多项式
通常以 \(F(x)\) 表示 原多项式，或 原式的变换，\(R(x)\) 或 \(G(x)\) 表示 答案多项式
不理解 \(\to\) 回去看前文
| 操作 | 模板题 | 核心式子 | 时间复杂度 |
|---|---|---|---|
| 多项式乘法 | \(Luogu ~ P3803\) | \(F(w_n ^ k) = FL (w_{n/2} ^ k) + w_n ^ k FR(w_{n / 2} ^ k)\) \(F(w_n ^ {k + n / 2}) = FL (w_{n/2} ^ k) - w_n ^ k FR(w_{n / 2} ^ k)\) |
\(O(N \log N)\) |
| 多项式加法 | \(-\) | \(F(x) + G(x) = \sum_{i = 0} ^ n F_i + G_i\) | \(O(N)\) |
| 多项式求逆 | \(Luogu ~ P4238\) | \(R(x) = 2R ^ *(x) - F(x)R ^ *(x) ^ 2\) | \(O(N \log N)\) |
| 多项式牛顿迭代 | \(-\) | \(F(x) = F ^ * (x) - \dfrac {G (F ^ * (x))} {G' (F ^ * (x))}\) | \(O (N \log N)\) |
| 多项式求导 | \(-\) | \(F'(x) = \sum_ {i = 1} i F_i x ^ {i - 1}\) | \(O (N)\) |
| 多项式积分 | \(-\) | \(F'(x) = \sum _ {i = 0} \dfrac {F_i x ^ {i + 1}} {i}\) | \(O (N)\) |
| 多项式开根 \(EZ\) | \(Luogu ~ P5205\) | \(G(x) = \dfrac {G ^ *(x) ^ 2 + F(x)} {2 G ^ * (x)}\) | \(O (N \log N)\) |
| 多项式开根 \(HD\) | \(Luogu ~ P5277\) | \(x = g ^ {i / k} \pmod p\) | \(O (N \log N)\) |
| 多项式带余除法 | \(Luogu ~ P4512\) | \(F ^ T (x) = G ^ T (x)T ^ T (x) + x ^ {N - M + 1} R ^ T(x)\) | \(O (N \log N)\) |
| 多项式对数函数 | \(Luogu ~ P4725\) | \(\ln (F (x)) = \dfrac {F'(x)} {F (x)}\) | \(O (N \log N)\) |
| 多项式指数函数 | \(Luogu ~ P4726\) | \(G (x) = G ^ * (x) (1 + F (x) - \ln (G ^ * (x)))\) | \(O(N \log N)\) |
| 多项式三角函数 | \(Luogu ~ P5264\) | \(\cos x = \dfrac {e ^ {i x} + e ^ {-i x}} {2}\) \(\sin x = \dfrac {e ^ {i x} - e ^ {-i x}} {2i}\) |
\(O (N \log N)\) |
| 多项式快速幂 \(EZ\) | \(Luogu ~ P5245\) | \(Sol ~ 1. ~ F (x) ^ K = e ^ {\ln (F (x)) * K}\) \(Sol ~ 2. ~ (F(x) ^ i) ^ 2 = F (x) ^ {2i}\) |
\(Sol ~ 1. ~ O (N \log N)\) \(Sol ~ 2. ~ O(N \log ^ 2 N)\) |
| 多项式快速幂 \(HD\) | \(Luogu ~ P5273\) | \(F(x) ^ K = e ^ {\ln (F(x)) * K}\) | \(O (N \log N)\) |
完整板板
#include <bits/stdc++.h>
const int MAXN = 1000005;
const int MOD = 998244353;
using namespace std;
int U[MAXN];
inline int Qpow (const int a, const int b = MOD - 2) {
if (b == 0) return 1;
if (b == 1) return a;
int Ret = Qpow (a, b >> 1);
if (b & 1) return 1ll * Ret * Ret % MOD * a % MOD;
return 1ll * Ret * Ret % MOD;
}
const int g = 3, ig = Qpow (g);
unsigned long long W[MAXN];
inline void Init (const int N) {
memset (U, 0, sizeof U);
for (int i = 0; i < N; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? N >> 1 : 0);
}
inline unsigned long long * Cross (unsigned long long *F, unsigned long long *G, int N) {
for (int i = 0; i < N; ++ i) F[i] = F[i] * G[i] % MOD;
return F;
}
inline void Print (unsigned long long * f, int N) {
for (int i = 0; i < N; ++ i) cerr << f[i] << ' ';
cerr << endl;
}
inline void Clr (unsigned long long *f, int L, int R) {
for (int i = L; i < R; ++ i) f[i] = 0;
}
inline void Rev (unsigned long long *f, int L, int R) {
reverse (f + L, f + R);
}
inline void NTT (unsigned long long *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1, tmp, u;
W[0] = 1;
for (int k = 2; k <= Len; len = k, k <<= 1) {
u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
for (int i = 1; i <= len; ++ i) W[i] = W[i - 1] * u % MOD;
for (int i = 0; i < Len; i += k)
for (int l = i; l < i + len; ++ l)
tmp = 1ll * W[l - i] * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);
if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
inline unsigned long long * Add (unsigned long long *F, unsigned long long *G, int N, int M) {
for (int i = 0; i <= max (N, M); ++ i) F[i] = (F[i] * G[i]) % MOD;
return F;
}
inline unsigned long long * Mul (unsigned long long *F, unsigned long long *G, int N, int M) {
int P = 1, I;
for (P = 1; P < N + M; P <<= 1);
Init (P), NTT (F, P, 1);
if (F != G) NTT (G, P, 1);
Cross (F, G, P), NTT (F, P, -1), I = Qpow (P);
for (int i = 0; i <= N + M; ++ i) F[i] = F[i] * I % MOD;
return F;
}
inline unsigned long long * Inv (unsigned long long *F, int N) {
int P = 1, I;
for (P = 1; P < N; P <<= 1);
unsigned long long A[MAXN], B[MAXN], C[MAXN];
A[0] = Qpow (F[0]);
for (int k = 2; k <= P; k <<= 1) {
memcpy (B, F, k << 3);
memcpy (C, A, k << 2);
Clr (A, k >> 1, k), Init (k), NTT (A, k, 1), NTT (B, k, 1);
Cross (B, A, k), NTT (B, k, -1), I = Qpow (k);
for (int i = 0; i < k; ++ i) B[i] = B[i] * I % MOD;
Clr (B, 0, k >> 1), NTT (B, k, 1), Cross (B, A, k), NTT (B, k, -1);
for (int i = 0; i < k; ++ i) B[i] = B[i] * I % MOD;
for (int i = 0; i < (k >> 1); ++ i) A[i] = C[i];
for (int i = (k >> 1); i < k; ++ i) A[i] = (MOD - B[i]) % MOD;
}
for (int i = 0; i <= N; ++ i) F[i] = A[i];
return F;
}
int inv[MAXN];
inline void init (int N) {
inv[1] = 1;
for (int i = 2; i <= N; ++ i) inv[i] = 1ll * inv[MOD % i] * (MOD - MOD / i) % MOD;
}
inline unsigned long long * Dlt (unsigned long long *F, int N) {
for (int i = 1; i < N; ++ i) F[i - 1] = F[i] * i % MOD;
F[N - 1] = 0;
return F;
}
inline unsigned long long * Sum (unsigned long long *F, int N) {
for (int i = N; i; -- i) F[i] = F[i - 1] * inv[i] % MOD;
F[0] = 0;
return F;
}
inline unsigned long long Sqrt_BSGS (unsigned long long x) {
unordered_map <unsigned long long, int> mp;
int t = sqrt (MOD);
unsigned long long Now = 1, gg = g;
for (int b = 1; b <= t; ++ b) Now = (Now * gg) % MOD, mp[(Now * x) % MOD] = b;
gg = Now, Now = 1;
for (int a = 1; a <= t; ++ a) {
Now = (Now * gg) % MOD;
if (mp[Now]) return a * t - mp[Now];
}
return -1;
}
unsigned long long Sqrt_Mod (unsigned long long a, unsigned long long k = 2) {
int r = Qpow (g, Sqrt_BSGS (a) / k);
return r = min (r, MOD - r);
}
inline unsigned long long * Sqrt_1 (unsigned long long *F, int N) {
int P = 1;
for (P = 1; P < N; P <<= 1);
unsigned long long B[MAXN], C[MAXN];
B[0] = 1;
for (int k = 2; k <= P; k <<= 1) {
for (int i = 0; i < (k >> 1); ++ i) C[i] = (B[i] << 1) % MOD;
Clr (C, k >> 1, k), Inv (C, k), Clr (C, k, k << 1), Mul (B, B, k, k);
for (int i = 0; i < k; ++ i) B[i] = (F[i] + B[i]) % MOD;
Mul (B, C, k, k), Clr (B, k, k << 2);
}
for (int i = 0; i <= N; ++ i) F[i] = B[i];
return F;
}
inline unsigned long long * Sqrt_2 (unsigned long long *F, int N) {
int P = 1;
for (P = 1; P < N; P <<= 1);
unsigned long long B[MAXN], C[MAXN];
B[0] = Sqrt_Mod (F[0], 2);
cerr << B[0] << endl;
for (int k = 2; k <= P; k <<= 1) {
for (int i = 0; i < (k >> 1); ++ i) C[i] = (B[i] << 1) % MOD;
Clr (C, k >> 1, k), Inv (C, k), Clr (C, k, k << 1), Mul (B, B, k, k);
for (int i = 0; i < k; ++ i) B[i] = (F[i] + B[i]) % MOD;
Mul (B, C, k, k), Clr (B, k, k << 2);
}
for (int i = 0; i <= N; ++ i) F[i] = B[i];
return F;
}
inline unsigned long long * Mof (unsigned long long * F, unsigned long long * G, int N, int M) {
unsigned long long A[MAXN], B[MAXN];
int P = N - M + 1;
Rev (F, 0, N), memcpy (A, F, P << 3), Rev (F, 0, N);
Rev (G, 0, M), memcpy (B, G, P << 3), Rev (G, 0, M);
Inv (B, P), Mul (B, A, P, P), Rev (B, 0, P), Clr (B, P, P << 1), memcpy (A, B, P << 3);
Mul (G, B, N, N), Clr (G, N, N << 1);
for (int i = 0; i < M - 1; ++ i) G[i] = (F[i] - G[i] + MOD) % MOD;
for (int i = 0; i < P; ++ i) F[i] = A[i];
Clr (G, M - 1, M << 1), Clr (F, P, P << 1);
return F;
}
inline unsigned long long * Ln (unsigned long long *F, int N) {
unsigned long long A[MAXN];
Clr (A, 0, N << 1), memcpy (A, F, N << 3);
Inv (A, N), Dlt (F, N), Mul (F, A, N, N);
Sum (F, N - 1), Clr (F, N, N << 1);
return F;
}
inline unsigned long long * Exp_1 (unsigned long long *F, int N) {
int P = 1;
for (P = 1; P < N; P <<= 1);
unsigned long long A[MAXN], B[MAXN];
A[0] = 1, Clr (A, 1, P << 1), Clr (B, 0, P << 1);
for (int k = 2; k <= P; k <<= 1) {
memcpy (B, A, k << 2), Clr (B, k >> 1, k << 1), Ln (B, k);
for (int i = 0; i < k; ++ i) B[i] = (- B[i] + F[i] + MOD) % MOD;
B[0] = (B[0] + 1) % MOD, Mul (A, B, k, k), Clr (A, k, k << 1);
}
for (int i = 0; i < N; ++ i) F[i] = A[i];
return F;
}
inline unsigned long long * Qpow_1 (unsigned long long *F, int N, int K) {
Ln (F, N);
for (int i = 0; i < N; ++ i) F[i] = F[i] * K % MOD;
Exp_1 (F, N);
return F;
}
inline unsigned long long * Qpow_2 (unsigned long long *F, int N, int K) {
unsigned long long A[MAXN], B[MAXN], C[MAXN];
Clr (A, 0, N), Clr (B, 0, N), Clr (C, 0, N);
A[0] = 1, memcpy (B, F, N << 3);
if (K == 0) A[0] = F[0];
while (K > 0) {
if (K & 1) -- K, memcpy (C, B, N << 3), Mul (A, B, N, N), memcpy (B, C, N << 3), Clr (A, N, N << 1), Clr (B, N, N << 2);
Mul (B, B, N, N), Clr (B, N, N << 2), K >>= 1;
}
for (int i = 0; i < N; ++ i) F[i] = A[i];
return F;
}
inline unsigned long long * Qpow_3 (unsigned long long *F, int N, long long K, long long K_, long long Siz) {
int I = 1, p = 0, c = 1;
while (F[p] == 0 && p < N) ++ p;
c = Qpow (F[p], K_), I = Qpow (F[p], MOD - 2);
if ((p && Siz > 6) || p * K > N) return Clr (F, 0, N), F;
for (int i = 0; i + p < N; ++ i) F[i] = F[i + p] * I % MOD;
Ln (F, N);
for (int i = 0; i < N; ++ i) F[i] = F[i] * K % MOD;
Exp_1 (F, N);
for (long long i = N - 1ll - 1ll * K * p; i >= 0; -- i) F[i + 1ll * p * K] = F[i] * c % MOD;
Clr (F, 0, min (1ll * p * K, 1ll * N));
return F;
}
inline unsigned long long * Sin (unsigned long long *F, int N) {
unsigned long long I = Sqrt_Mod (MOD - 1), A[MAXN];
for (int i = 0; i < N; ++ i) F[i] = F[i] * I % MOD;
init (N + 4), Exp_1 (F, N);
memcpy (A, F, N << 3), Inv (A, N), I = Qpow (I << 1);
for (int i = 0; i < N; ++ i) F[i] = I * (F[i] - A[i] + MOD) % MOD;
return F;
}
inline unsigned long long * Cos (unsigned long long *F, int N) {
unsigned long long I = Sqrt_Mod (MOD - 1), A[MAXN];
for (int i = 0; i < N; ++ i) F[i] = F[i] * I % MOD;
init (N + 4), Exp_1 (F, N);
memcpy (A, F, N << 3), Inv (A, N), I = Qpow (2);
for (int i = 0; i < N; ++ i) F[i] = I * (F[i] + A[i]) % MOD;
return F;
}
unsigned long long F[MAXN], G[MAXN];
int main () {
// ...
return 0;
}
一些题
都很基础，煎蛋题居多，各位巨佬可以跳过这部分
Luogu P3338 [ZJOI2014] 力
力起来！
看式子，可以想到把这玩意儿 化为卷积形式，然后 \(FFT\) 加速计算
首先 注意到 可以 提取公因数
于是代入 \(E_j\) 有
我们令 \(A_i = i ^ {-2}\)，显然 \(A _ {i - j} = A _ {j - i}\) 于是有
可以看到 前面已经是卷积形式 了，但后面是不对的，因为 \(i = j + 1\) 进行了限制
但注意到它 求和上界顶到了次数上界 \(n\)，故而我们 反转 \(q\) 的系数（像 带余除法 的操作）
即反转后多项式 \(Q ^ T = \sum q ^ T _ i x ^ i\)，有 \(q ^ T _ i = q _ {n - i + 1}\)（代码里面下标从 \(0\) 开始，是 \(n - i - 1\)）
于是上式就可以写成（我们发现第 \(j\) 项加上不会有任何影响，但是会让式子变得好看）
这下是 卷积 的形式了，然后套 \(FFT\) 做两遍乘法就行，代码
#include <bits/stdc++.h>
#define Complex complex<long double>
const int MAXN = 400005;
const long double ZC = 1e11;
const long double PI = acos (-1);
using namespace std;
int U[MAXN];
inline void Init (const int N) {
memset (U, 0, sizeof U);
for (int i = 0; i < N; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? N >> 1 : 0);
}
Complex W[MAXN];
inline Complex * FFT (Complex *F, int N, int Typ) {
for (int i = 0; i < N; ++ i) if (i < U[i]) swap (F[i], F[U[i]]);
W[0] = Complex (1, 0);
for (int k = 2, len = 1; k <= N; len = k, k <<= 1) {
Complex u = Complex (cos (2 * PI / k), Typ * sin (2 * PI / k)), tmp;
for (int i = 1; i < len; ++ i) W[i] = W[i - 1] * u;
for (int i = 0; i < N; i += k)
for (int l = i; l < i + len; ++ l)
tmp = W[l - i] * F[l + len], F[l + len] = F[l] - tmp, F[l] = F[l] + tmp;
}
return F;
}
inline void Clr (Complex * F, int L, int R) {
for (int i = L; i <= R; ++ i) F[i] = {0, 0};
}
inline void Mul (Complex * F, Complex * G, int N, int M) {
int P = 1;
for (P = 1; P < N + M; P <<= 1);
Complex f[P + 2]; Init (P); Clr (f, 0, P);
for (int i = 0; i < P; ++ i) f[i] = {F[i].real (), G[i].real ()};
FFT (f, P, 1);
for (int i = 0; i < P; ++ i) f[i] = f[i] * f[i];
FFT (f, P, -1);
for (int i = 0; i < N + M; ++ i) F[i] = {f[i].imag () / P / 2.0 + 0.4999, 0};
}
int N, M;
Complex Q[MAXN], G[MAXN], K[MAXN];
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> N;
for (int i = 0; i < N; ++ i) cin >> Q[i], K[i] = (i == 0) ? i : ((long double) 1.0 / i / i), K[i] *= ZC;
for (int i = 0; i < N; ++ i) G[i] = Q[N - i - 1];
Mul (Q, K, N, N), Mul (G, K, N, N);
for (int i = 0; i < N + N - 1; ++ i) Q[i] = Q[i] - G[N - i - 1];
for (int i = 0; i < N; ++ i)
cout << fixed << setprecision (3) << Q[i].real () / ZC << '\n';
return 0;
}
这题非常恶心，卡卡精度，而且 \(i ^ {-2}\) 这玩意儿太小了，你需要 先乘一下，最后在除回来
甚至乘的这个数还有要求，我这里试出来 \(10 ^ {11}\) 是可以的
Luogu P4173 残缺的字符串
\(FFT\) 字符串匹配的 板子？
考虑虽然 \(KMP\)，\(Hash\) 都能在 线性时间 内解决 标准的字符串匹配问题
但是遇到这种 棘手的通配符，上面的算法就有些 捉襟见肘 了
而这个时候，在 \(O(N \log N)\) 时间解决问题的 \(FFT\) 字符串匹配就出现了
这个算法基于这样一种原理，先从 标准的字符串匹配 来讲
我们可以假设 有字符串 \(A, B\)，长为 \(N, M\)，令 \(C(x, y) = A_x - B_y\)
显然当 \(C(x, y) = 0\) 时，\(A\) 的第 \(x\) 位 和 \(B\) 的第 \(y\) 位 匹配
进而我们发现，如果存在 \(C(i + j, i) = 0，i \in [0, M - 1]\)
则相当于 \(A\) 从 \(j\) 开始的 长为 \(M\) 的串 与 \(B\) 完全匹配
于是我们定义 \(P(x) = \sum_{i = 0} ^ M C(x - M + i, i)\)
想要使得 \(P(x) = 0\) 时，代表 \(A\) 从 \(x\) 结尾的 \(M\) 位与 \(B\) 完全匹配
为了后续 凑卷积式，这里只能用 \(x\) 结尾，而不能用 \(x\) 开头
但是显然这个时候 \(C (x - M + i, i)\) 有正有负，我们并 不能实现上述目标
也就是我可能某一位 \(C = 1\)，某一位 \(C = -1\)，最后和为 \(0\)，但显然不匹配
所以我们 再平方一下，令 \(P(x) = \sum_{i = 0} ^ M (C(x - M + i, i)) ^ 2\)
此时对于每个 \(i\)，\((C(x - M + i, i)) ^ 2\) 均为 非负整数，这就实现了 上述目标
那怎么求这个 \(P(x)\) 呢？我们 展开这个式子
\[P(x) = \sum _ {i = 1} ^ M (A_{x - M + i} - B_i) ^ 2 \]\[= \sum _ {i = 1} ^ M A_{x - M + i} ^ 2 + B_i ^ 2 - 2A_{x - M + i} B_i \]我们发现 \(\sum B_i ^ 2\) 这部分 显然可以直接预处理
而 \(\sum A_{x - M + i} ^ 2\) 也可以通过 预处理 \(A_i ^ 2\) 的前缀和 后 \(O(1)\) 求得
难点就在 \(\sum _ {i = 0} ^ M 2A_{x - M + i} B_i\) 上了，看到 \(-M + i\) 和 \(i\)，恰好 \(B\) 的长度又是 \(M\)
我们又可以想到 反转 \(B\) 串 这个方法，于是令 \(B'_i = B_{M - i}\)
\(\sum _ {i = 0} ^ M 2A_{x - M + i}B_i\) 就变成了 \(\sum _ {i = 0} ^ M 2A_{x - M + i}B'_{M - i}\)，符合卷积形式
于是 \(FFT\) 加速即可 \(O(N \log N)\) 处理出所有 \(P(x)\)
标准的匹配就是这样，回到这道题，由于 通配符可以和任意字符匹配
也就是当 \(A_x\) 或 \(B_y\) 中 存在至少一个通配符 时，\(C(x, y) = 0\)
所以我们可以 改一下完全匹配函数 \(P(x)\)，令 \(P(x) = \sum_{i = 0} ^ M (C(x - M + i, i)) ^ 2 A_{x - M + i} B_i\)
并使得 当 \(A\) 第 \(i\) 位为通配符时，\(A_i = 0\)（\(B\) 同理）
简单模拟可以发现，改动后的 完全匹配函数 可以实现题目要求的匹配
于是我们继续拆式子
我们同理反转 \(B\) 串，令 \(B'_i = B_{M - i}\)，得到下式
容易发现 三个的部分 均 分别符合卷积形式，于是做三次 \(FFT\) 加速卷积
最后系数和差即可，代码
#include <bits/stdc++.h>
#define Complex complex<double>
const int MAXN = 1050005;
const double PI = acos (-1);
const double EPS = 1;
using namespace std;
int U[MAXN];
inline void Init (const int N) {
memset (U, 0, sizeof U);
for (int i = 0; i < N; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? N >> 1 : 0);
}
Complex W[MAXN];
inline void FFT (Complex *F, int N, int Typ) {
for (int i = 0; i < N; ++ i) if (i < U[i]) swap (F[i], F[U[i]]);
Complex tmp; W[0] = 1;
for (int k = 2, len = 1; k <= N; len = k, k <<= 1) {
Complex u = {cos (2 * PI / k), Typ * sin (2 * PI / k)};
for (int i = 1; i < len; ++ i) W[i] = W[i - 1] * u;
for (int i = 0; i < N; i += k)
for (int l = i; l < i + len; ++ l)
tmp = W[l - i] * F[l + len], F[l + len] = F[l] - tmp, F[l] = F[l] + tmp;
}
}
double A[MAXN], B[MAXN];
Complex F[MAXN], G[MAXN], T[MAXN];
vector <int> Q;
inline int Match (string S1, string S2) {
int l1 = S1.size (), l2 = S2.size (), P = 1;
reverse (S1.begin (), S1.end ());
for (P = 1; P < l1 + l2; P <<= 1);
for (int i = 0; i < l1; ++ i) A[i] = (S1[i] == '*') ? 0 : (S1[i] - 'a' + 1);
for (int i = 0; i < l2; ++ i) B[i] = (S2[i] == '*') ? 0 : (S2[i] - 'a' + 1);
Init (P);
for (int i = 0; i < P; ++ i) F[i] = {A[i] * A[i] * A[i], 100.0 * B[i]};
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) T[i] += F[i] * F[i] / 200.0;
for (int i = 0; i < P; ++ i) F[i] = {A[i] * A[i], B[i] * B[i]};
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) T[i] -= F[i] * F[i];
for (int i = 0; i < P; ++ i) F[i] = {100.0 * A[i], B[i] * B[i] * B[i]};
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) T[i] += F[i] * F[i] / 200.0;
FFT (T, P, -1);
for (int i = 0; i < P; ++ i) T[i] = (T[i].imag () / (double) P);
for (int i = l1 - 1; i < l2; ++ i) if (T[i].real () < EPS) Q.push_back (i - l1 + 2);
return Q.size ();
}
int N, M;
string x, y;
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> N >> M;
cin >> x >> y;
cout << Match (x, y) << endl;
for (auto i : Q) cout << i << ' ';
return 0;
}
Luogu P3723 [AH2017/HNOI2017] 礼物
我居然最开始能想到 做 \(100\) 次 \(FFT\)，真是天才
题目是简单的，考虑设 亮度改变量 \(c\)，于是有 差异值为 \(\sum _ {i = 1} ^ n (x_i - y_i - c) ^ 2, c \in Z\)
展开有 \(\sum _ {i = 1} ^ n ((x _ i ^ 2 + y _ i ^ 2) + (c ^ 2 - 2c (x_i - y _i)) - 2 x_i y_i)\)
可以看到这分成 三个部分
第一部分 \(x_i ^ 2 + y_i ^ 2\) 显然可以 直接预处理，是个定值
第二部分 \(c ^ 2 - 2c (x _i - y _ i)\)，把求和符号放进来，变成 \(nc ^ 2 - 2c (\sum_{i = 1} ^ n (x_i - y_i))\)
显然 \(\sum_{i = 1} ^ n (x_i - y_i)\) 可以预处理，于是变成一个 定系数二次函数 求最小值问题
而由于 \(n > 0\)，二次函数开口向上，最小值在顶点处
即 $c = - \dfrac {2n} {2 \sum_{i = 1} ^ n (x_i - y_i)} = - \dfrac {n} { \sum_{i = 1} ^ n (x_i - y_i)} $ 时 第二部分有最小值，代入求值即可
但此时 \(c\) 可能非整数，需要把 左右两边的数代入求值后取最小值
最后考虑第三部分，\(2 x_i y_i\) 看上去不好做，套路的 反转一下第一串 变成 \(2 x_{n - i} y_i\)
求和符号放进来就是 \(\sum _ {i = 1} ^ {n} 2 x_{n - i} y_i\)，典型的卷积形式
而由于 手环可以旋转，所以可以想到 倍增一串，长度变为 \(2n\)
最后 \(FFT\) 加速卷积来求 第三部分最大值，统计答案时只应当在 \([n, 2n]\) 区间统计
下标从 \(0\) 开始就是 \([n - 1, 2n)\)
然后 加上前面常数 即可，代码难度较小，并且由于都是整数，可以使用 \(NTT\)
#include <bits/stdc++.h>
const int MAXN = 262150;
const int MOD = 998244353;
using namespace std;
int U[MAXN];
inline int Qpow (const int a, const int b = MOD - 2) {
if (b == 0) return 1;
if (b == 1) return a;
int Ret = Qpow (a, b >> 1);
if (b & 1) return 1ll * Ret * Ret % MOD * a % MOD;
return 1ll * Ret * Ret % MOD;
}
inline void Init (const int N) {
memset (U, 0, sizeof U);
for (int i = 0; i < N; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? N >> 1 : 0);
}
unsigned long long W[MAXN];
unsigned long long g = 3, ig = Qpow (g);
inline void NTT (unsigned long long *p, const int Len, const int Typ) {
for (int i = 0; i < Len; ++ i) if (i < U[i]) swap (p[i], p[U[i]]);
int len = 1, tmp, u;
W[0] = 1;
for (int k = 2; k <= Len; len = k, k <<= 1) {
u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
for (int i = 1; i <= len; ++ i) W[i] = W[i - 1] * u % MOD;
for (int i = 0; i < Len; i += k)
for (int l = i; l < i + len; ++ l)
tmp = 1ll * W[l - i] * p[l + len] % MOD, p[l + len] = (p[l] - tmp + MOD), p[l] = (p[l] + tmp);
if (k == (1 << 16)) for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
for (int i = 0; i < Len; ++ i) p[i] %= MOD;
}
int N, M;
unsigned long long A[MAXN], B[MAXN], Ans, Tmp;
int SumA = 0, SumB = 0, X, Sum = 0, P = 1, I;
inline int Solve () {
for (int i = 0; i < N; ++ i) Ans += A[i] * A[i], SumA += A[i];
for (int i = 0; i < N; ++ i) Ans += B[i] * B[i], SumB += B[i];
X = - (SumA - SumB) / N, Sum = min ({1ll * N * X * X + 2 * X * (SumA - SumB), 1ll * N * (X + 1) * (X + 1) + 2 * (X + 1) * (SumA - SumB), 1ll * N * (X - 1) * (X - 1) + 2 * (X - 1) * (SumA - SumB)}), Ans += Sum;
reverse (A, A + N), memcpy (A + N, A, N << 3);
for (P = 1; P < N * 3; P <<= 1);
Init (P), I = Qpow (P), NTT (A, P, 1), NTT (B, P, 1);
for (int i = 0; i < P; ++ i) A[i] = A[i] * B[i] % MOD;
NTT (A, P, -1);
for (int i = N - 1; i < (N << 1); ++ i) Tmp = max (Tmp, (A[i] * I) % MOD * 2);
return (int) (Ans - Tmp + MOD) % MOD;
}
int main () {
cin >> N >> M;
for (int i = 0; i < N; ++ i) cin >> A[i];
for (int i = 0; i < N; ++ i) cin >> B[i];
cout << Solve () << endl;
return 0;
}
还 挺 快 的，进行快读之后，获取 \(rk1\)
Luogu P3763 [TJOI2017] DNA
好题，牛牛的一个套路 —— \(\textsf H\)\(\textsf {anghang}\)
也是一个字符串匹配，但是和上面的 带通配符 的方法不一样
这个是 容许 \(3\) 位差距（但是 不指定是哪 \(3\) 位，通配符相当于指定了位置）
然后看到这个题 字符集大小 只有 \(4\)，可以想到在这上面做手脚
我们对四个字符 分别处理，例如当前处理 \(A\)，两字符串分别为 \(S1, S2\)，长为 \(N, M\)
我们用两个 \(01\) 序列 \(F, G\) 表示 两个字符串 \(S1, S2\) 对应位置是不是 \(A\)
这时候我们发现 匹配函数 \(C(x, y) = F_x G_y\) 即可
这将保证当且仅当 \(S1_x, S2_y\) 均为 \(A\) 时，\(C(x, y)\) 为 \(1\)，也相当于 这一位匹配上了
于是同样的，完美匹配函数 \(P(x) = \sum_{i = 0} ^ M C(x + i, i)\)，我们 套路地将 \(S2\) 反转
设反转后对应的 \(01\) 序列 为 \(G'\)，于是最终 \(P(x) = \sum_{i = 0} ^ M F_{x + i} G'_{M - x - i}\)
对反转搞不清楚的话看 上面几道题
直接跑 \(FFT\) 加速卷积，然后 四个字符都这么处理，最后 系数加起来
得到的 \((\sum)P(x)\) 就是 文本串从第 \(x\) 位开始后 \(M\) 位 能与模式串匹配的位数
于是满足 \(P(x) + 3 \ge M\) 的 \(x\) 就是 符合条件的位，统计答案即可
上面写的可能有点抽象，但是 代码是好理解的
#include <bits/stdc++.h>
#define Complex complex<double>
const int MAXN = 262150;
const double PI = acos (-1);
const double EPS = 1e-8;
using namespace std;
Complex W[MAXN];
int U[MAXN];
inline void Init (const int N) {
memset (U, 0, sizeof U);
for (int i = 0; i < N; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? N >> 1 : 0);
}
inline void FFT (Complex *F, int N, int Typ) {
for (int i = 0; i < N; ++ i) if (i < U[i]) swap (F[i], F[U[i]]);
Complex u, tmp; W[0] = 1;
for (int k = 2, len = 1; k <= N; len = k, k <<= 1) {
u = {cos (2 * PI / k), Typ * sin (2 * PI / k)};
for (int i = 1; i < len; ++ i) W[i] = W[i - 1] * u;
for (int i = 0; i < N; i += k)
for (int l = i; l < i + len; ++ l)
tmp = W[l - i] * F[l + len], F[l + len] = F[l] - tmp, F[l] = F[l] + tmp;
}
}
namespace Value {
string S1, S2;
Complex F[MAXN], T[MAXN];
int A[MAXN], B[MAXN];
int P = 1;
inline void Solve () {
cin >> S1 >> S2;
int l1 = S1.size (), l2 = S2.size (), Ans = 0;
for (P = 1; P < l1 + l2; P <<= 1);
reverse (S2.begin (), S2.end ()), Init (P);
for (int i = 0; i < l1; ++ i) A[i] = (S1[i] == 'A');
for (int i = 0; i < l2; ++ i) B[i] = (S2[i] == 'A');
for (int i = 0; i < P; ++ i) F[i] = {1.0 * A[i], 1.0 * B[i]};
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) F[i] *= F[i], T[i] += F[i];
for (int i = 0; i < l1; ++ i) A[i] = (S1[i] == 'T');
for (int i = 0; i < l2; ++ i) B[i] = (S2[i] == 'T');
for (int i = 0; i < P; ++ i) F[i] = {1.0 * A[i], 1.0 * B[i]};
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) F[i] *= F[i], T[i] += F[i];
for (int i = 0; i < l1; ++ i) A[i] = (S1[i] == 'C');
for (int i = 0; i < l2; ++ i) B[i] = (S2[i] == 'C');
for (int i = 0; i < P; ++ i) F[i] = {1.0 * A[i], 1.0 * B[i]};
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) F[i] *= F[i], T[i] += F[i];
for (int i = 0; i < l1; ++ i) A[i] = (S1[i] == 'G');
for (int i = 0; i < l2; ++ i) B[i] = (S2[i] == 'G');
for (int i = 0; i < P; ++ i) F[i] = {1.0 * A[i], 1.0 * B[i]};
FFT (F, P, 1);
for (int i = 0; i < P; ++ i) F[i] *= F[i], T[i] += F[i];
FFT (T, P, -1);
for (int i = l2 - 1; i < l1; ++ i)
if (T[i].imag () / P / 2 + 3.0 + EPS >= l2) ++ Ans;
for (int i = 0; i < P; ++ i) T[i] = 0;
cout << Ans << endl;
}
}
int T;
int main () {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> T;
while (T --) Value::Solve ();
return 0;
}
Luogu P5488 差分与前缀和
我们把这个序列啊，整到一个 多项式 \(F\) 上，表示成这样 \(F(x) = \sum_{i = 1} ^ N a_i x ^ i\)
我们发现 差分就是后面一项减去前面一项
相当于 多项式中 \(x + 1\) 项系数减去 \(x\) 项系数，于是可以想到 \(xF(x) - F(x)\)
就是 \((x - 1) F(x)\)，显然，做一次差分 就相当于 多项式乘一次 \(x - 1\)
\(k\) 次差分的式子易得为 \((x - 1) ^ k F(x)\)
我们注意到 前缀和是差分的逆运算，可以猜想到前缀和就是需要 乘一个 \((x - 1) ^ {-1}\)
简单验证这东西显然正确，这里不讲
于是 \(k\) 维前缀和的式子也容易知道是 \((x - 1) ^ {-k} F(x)\)
但是 \((x - 1) ^ {\pm k}\) 如果用 多项式快速幂 算 你就输了
显然我们可以直接用 二项式定理 展开，然后 预处理一下组合数
最后 \(FFT\) 加速卷积即可，具体看代码
#include <bits/stdc++.h>
const int MAXN = 262150;
const unsigned long long MOD = 1004535809;
const unsigned long long g = 3, ig = 334845270;
using namespace std;
inline unsigned long long gi() {
char cc = getchar() ; unsigned long long cn = 0, flus = 1 ;
while( cc < '0' || cc > '9' ) { if( cc == '-' ) flus = - flus ; cc = getchar() ; }
while( cc >= '0' && cc <= '9' ) cn = ( cn * 10 + cc - '0' ) % MOD, cc = getchar() ;
return cn * flus ;
}
inline void prt (unsigned long long k) {
char s[20]; int a = 0;
while (k) s[++ a] = k % 10 + '0', k /= 10;
while (a) putchar (s[a]), a --;
putchar (' ');
}
inline long long Qpow (const long long a, const long long b = MOD - 2) {
if (b == 0) return 1;
if (b == 1) return a;
long long Ret = Qpow (a, b >> 1);
if (b & 1) return Ret * Ret % MOD * a % MOD;
return Ret * Ret % MOD;
}
unsigned long long W[MAXN];
int U[MAXN];
inline void NTT (unsigned long long *F, int N, int Typ) {
for (int i = 0; i < N; ++ i) if (i < U[i]) swap (F[i], F[U[i]]);
int u = 1, tmp; W[0] = 1;
for (int k = 2, len = 1; k <= N; len = k, k <<= 1) {
u = Qpow (Typ == 1 ? g : ig, (MOD - 1) / k);
for (int i = 1; i <= len; ++ i) W[i] = W[i - 1] * u % MOD;
for (int i = 0; i < N; i += k)
for (int l = i; l < i + len; ++ l)
tmp = W[l - i] * F[l + len] % MOD, F[l + len] = F[l] - tmp + MOD, F[l] = F[l] + tmp;
if (k == (1 << 16)) for (int i = 0; i < N; ++ i) F[i] = F[i] % MOD;
}
for (int i = 0; i < N; ++ i) F[i] = F[i] % MOD;
cerr << endl;
}
int N, K, T, P, I;
unsigned long long F[MAXN], B[MAXN];
inline void Init () {
for (P = 1; P < (N << 1); P <<= 1);
if (T == 0) for (int i = 1; i < N; ++ i) B[i] = B[i - 1] * (K + i - 1) % MOD * Qpow (i) % MOD;
if (T == 1) for (int i = 1; i < N; ++ i) B[i] = (- ((B[i - 1] * (K - i + 1 + MOD) % MOD) * Qpow (i) % MOD) + MOD) % MOD;
for (int i = 0; i < P; ++ i) U[i] = (U[i >> 1] >> 1) | ((i & 1) ? P >> 1 : 0);
}
int main () {
N = gi (), K = gi (), T = gi ();
for (int i = 0; i < N; ++ i) F[i] = gi ();
B[0] = 1, Init ();
NTT (F, P, 1), NTT (B, P, 1);
for (int i = 0; i < P; ++ i) F[i] = F[i] * B[i] % MOD;
NTT (F, P, -1), I = Qpow (P);
for (int i = 0; i < N; ++ i) prt (F[i] * I % MOD);
return 0;
}
尾
记一下差的 坑
- 三角函数 欧拉公式 的 证明
- P5050 多项式多点求值
- P5265 多项式反三角函数
- P5158 多项式快速插值
- P5373 多项式复合函数
- P5809 多项式复合逆
引用 & 鸣谢
这些文章都使笔者学到 褐了 很多内容，而且 质量较高，十分推荐
[1] 傅里叶变换(FFT)学习笔记 - command_block 的博客 - 洛谷博客 (luogu.com.cn)
[2] NTT与多项式全家桶 - command_block 的博客 - 洛谷博客 (luogu.com.cn)
[3] 原根的概念、性质及其存在性 - 知乎 (zhihu.com) - Mathis Wang（笔 误 の 神）
[4] FFT用到的各种素数 | Miskcoo’s Space
[5] 多项式牛顿迭代 - feiko - 博客园 (cnblogs.com)
[6] 如何通俗易懂地讲解牛顿迭代法？_牛拉迭代-CSDN博客（配图 十分清晰，好理解）
[7] 多项式计数杂谈 - command_block 的博客 - 洛谷博客 (luogu.com.cn)
[8] 浅谈多项式ln和exp - 知乎 (zhihu.com)
[9] 半在线卷积小记 - command_block 的博客 - 洛谷博客 (luogu.com.cn)
[10] 【题解】多项式三角函数 - 洛谷专栏 (luogu.com.cn)
[11] 国家集训队 \(2016\) 年论文集 《再探快速傅里叶变换》（\(P123 \sim P140\)）—— 毛啸
浙公网安备 33010602011771号