多项式学习-完全版
持续更新中...
start in 14 July 2026
卷积变换
快速求 $$c_k = \sum _{i \oplus j = k} a_i b_j$$
关于卷积…… 如何通俗易懂地解释卷积? - 知乎 https://www.zhihu.com/question/22298352/answer/228543288
FFT -> NTT -> FMT/FWT
https://www.luogu.com.cn/training/1008#information
https://blog.messywind.top/posts/cp生成函数多项式题单/
快速傅里叶变换 FFT
题目大意:对多项式 \(F(x), G(X)\),以 \(O(n \log n)\) 的时间复杂度求得 \(F(x) G(x)\) 的乘积。
https://www.luogu.com.cn/article/jm6yxjpx
https://www.cnblogs.com/pks-t/p/9251147.html
普通的做法,直接相乘,时间复杂度是 \(O(n^2)\) 的,难以接受。
将 \(F(X), G(X)\) 用点值表示,这样就可以用 \(O(n)\) 来相乘了。而将系数表达式与点值表达式的互相转换,时间复杂度为 \(O(n \log n)\)。
如下是非递归写法。
#include <stdio.h>
#include <algorithm>
#include <math.h>
const double PI = acos(-1.0);
struct CPX {
double x, y; // x+yi
CPX (double x = 0.0, double y = 0.0) : x(x), y(y) { }
CPX operator + (const CPX &t) { return { x + t.x, y + t.y }; }
CPX operator - (const CPX &t) { return { x - t.x, y - t.y }; }
CPX operator * (const CPX &t) { return { x * t.x - y * t.y , x * t.y + y * t.x }; }
} a[3000003], b[3000003], c[3000003]; // note that N >= 2^upper(logn)
int n, m, len, rev[3000003];
inline void fft (CPX a[], double md) { // md: 1: FFT, -1: IFFT
for (int i = 0; i < len; ++i) if (i < rev[i]) std:: swap(a[i], a[rev[i]]);
for (int l = 1; l < len; l <<= 1) { // l: half of the len of the interval now
CPX W = { cos(PI / l), md * sin(PI / l) };
for (int i = 0; i < len; i += (l << 1)) { // for every interval...
CPX w(1.0, 0.0);
for (int j = 0; j < l; ++j) { // for every node in the interval now...
CPX u = a[i + j];
CPX v = w * a[i + j + l];
a[i + j] = u + v;
a[i + j + l] = u - v;
w = w * W;
}
}
}
}
inline void kagari () {
scanf("%d %d", &n, &m);
for (int i = 0; i <= n; ++i) scanf("%lf", &a[i].x);
for (int i = 0; i <= m; ++i) scanf("%lf", &b[i].x);
len = 1;
while (len <= n + m) len <<= 1;
for (int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
fft(a, 1);
fft(b, 1);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i];
fft(c, -1);
for (int i = 0; i <= n + m; ++i)
printf("%d ", (int)(c[i].x / len + 0.5));
puts("");
return;
}
int main () {
kagari();
return 0;
}
此外还有一个有趣的优化(from https://www.luogu.com.cn/article/pvwh9hpp ):
将 \(G(x)\) 放入 \(F(x)\) 的虚部中,这时 \((F+Gi)^2=(F^2-G^2)+2FGi\),虚部除以 \(2\) 即为答案。如此,只需做两次傅里叶变换(先对 \(F+Gi\) 做 FFT,开方后,再对 \(FG\) 做 IFFT)而非三次,在常数级别大大优化了时间复杂度。
原根
https://blog.csdn.net/zhang20072844/article/details/11541133
假设一个数 \(G\) 对于质数 \(P\) 来说是原根,那么 \(G^i \mod P\)的结果两两不同,且有 \(1<G<P, 1<i<P\),那么 \(G\) 可以称为是 \(P\) 的一个原根。\(P\) 的原根 $g $ 的 \(1,2,\dots,P-1\) 次幂 \(\mod P\) 的结果一定互不相同。
要求最小原根,则可以从 \(2\) 开始枚举 \(G\),直到 \(G^{P-1} \equiv 1 \pmod P\)。
p.s. 不是每个素数都存在原根。\(p=k\times2^a+1\)。
快速数论变换 NTT
https://zhuanlan.zhihu.com/p/347726949
在 FFT 的基础上,用原根代替单位根,以解决取模与 double 的精度问题。
那么显然 NTT 就难以处理 double 类型的数据了。
#include <stdio.h>
#include <algorithm>
#define ll long long
const ll MODN = 998244353;
const ll G = 3; // primitive root of MODN
const ll GI = 332748118; // inv of MODN
inline ll power (ll x, ll k) {
ll res = 1;
while (k) {
if (k & 1) res = res * x % MODN;
x = x * x % MODN;
k >>= 1;
}
return res;
}
int n, m, len, rev[3000003];
ll a[3000003], b[3000003], c[3000003];
inline void ntt (ll a[], ll md) { // md: 1: NTT, -1: INTT
for (int i = 0; i < len; ++i) if (i < rev[i]) std:: swap(a[i], a[rev[i]]);
for (int l = 1; l < len; l <<= 1) {
ll W = power(md == 1 ? G : GI, (MODN - 1) / (l << 1));
for (int i = 0; i < len; i += (l << 1)) {
ll w = 1;
for (int j = 0; j < l; ++j) {
ll u = a[i + j];
ll v = w * a[i + j + l] % MODN;
a[i + j] = (u + v) % MODN;
a[i + j + l] = (u - v + MODN) % MODN;
w = w * W % MODN;
}
}
}
}
inline void kagari () {
scanf("%d %d", &n, &m);
for (int i = 0; i <= n; ++i) scanf("%lld", &a[i]);
for (int i = 0; i <= m; ++i) scanf("%lld", &b[i]);
len = 1;
while (len <= n + m) len <<= 1;
for (int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
ntt(a, 1);
ntt(b, 1);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i] % MODN;
ntt(c, -1);
ll invlen = power(len, MODN - 2);
for (int i = 0; i <= n + m; ++i)
printf("%lld ", c[i] * invlen % MODN);
puts("");
return;
}
int main () {
kagari();
return 0;
}
快速莫比乌斯变换 FMT | 快速沃尔什变换 FWT
https://www.luogu.com.cn/article/2pavj2pd
https://www.luogu.com.cn/article/idd90f6t
FWT 是用于解决对下标进行位运算卷积问题的方法,与 FFT 类似。即从求 \(C_i= \Sigma_{i=j+k} A_j B_k\) 变为求 \(C_i= \Sigma_{i=j\oplus k} A_j B_k\),其中 \(\oplus\) 是二元位运算 and, or, xor 中的一个。
以 P4717 【模板】快速莫比乌斯 / 沃尔什变换 (FMT / FWT) 为例:
题目大意:给定长度为 \(2^n\) 两个序列 \(A,B\),设 \(C_i=\sum_{j\oplus k = i}A_j \times B_k\),分别当 \(\oplus\) 是 or, and, xor 时求出 \(C\),并对 \(998244353\) 取模。
对 or, and, xor 三种操作,代码差别只在于对 \(u, v\) 的处理区别。实际应用上在套用模板时,也很少会改变模板中的代码。
#include <stdio.h>
#define ll long long
#define MODN 998244353
#define inv2 499122177
int n, len;
ll a[2000003], b[2000003], c[2000003];
inline void fwtor (ll a[], ll md) {
for (int l = 1; l < len; l <<= 1)
for (int i = 0; i < len; i += (l << 1))
for (int j = 0; j < l; ++j) {
ll u = a[i + j];
ll v = a[i + j + l];
if (md == 1)
a[i + j + l] = (u + v) % MODN;
else
a[i + j + l] = (v - u + MODN) % MODN;
}
}
inline void fwtand (ll a[], ll md) {
for (int l = 1; l < len; l <<= 1)
for (int i = 0; i < len; i += (l << 1))
for (int j = 0; j < l; ++j) {
ll u = a[i + j];
ll v = a[i + j + l];
if (md == 1)
a[i + j] = (u + v) % MODN;
else
a[i + j] = (u - v + MODN) % MODN;
}
}
inline void fwtxor (ll a[], ll md) {
for (int l = 1; l < len; l <<= 1)
for (int i = 0; i < len; i += (l << 1))
for (int j = 0; j < l; ++j) {
ll u = a[i + j];
ll v = a[i + j + l];
if (md == 1)
a[i + j] = (u + v) % MODN,
a[i + j + l] = (u - v + MODN) % MODN;
else
a[i + j] = (u + v) % MODN * inv2 % MODN,
a[i + j + l] = (u - v + MODN) % MODN * inv2 % MODN;
}
}
inline void kagari () {
scanf("%d", &n);
len = (1 << n);
for (int i = 0; i < len; ++i) scanf("%lld", &a[i]);
for (int i = 0; i < len; ++i) scanf("%lld", &b[i]);
fwtor(a, 1), fwtor(b, 1);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i] % MODN;
fwtor(a, -1), fwtor(b, -1), fwtor(c, -1);
for (int i = 0; i < len; ++i) printf("%lld ", c[i]);
puts("");
fwtand(a, 1), fwtand(b, 1);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i] % MODN;
fwtand(a, -1), fwtand(b, -1), fwtand(c, -1);
for (int i = 0; i < len; ++i) printf("%lld ", c[i]);
puts("");
fwtxor(a, 1), fwtxor(b, 1);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i] % MODN;
fwtxor(a, -1), fwtxor(b, -1), fwtxor(c, -1);
for (int i = 0; i < len; ++i) printf("%lld ", c[i]);
puts("");
return;
}
int main () {
kagari();
return 0;
}
P1919 【模板】高精度乘法 / A*B Problem 升级版
题目大意:求 \(a \times b\),其中 \(a, b \leq 10^{10^6}\)
将 \(a, b\) 拆成数组,直接利用 NTT 即可。是 FFT 与 NTT 的经典简单应用。
const ll MODN = 998244353;
const ll G = 3;
const ll GI = 332748118;
int n, m, len, rev[8000003];
ll a[8000003], b[8000003], c[8000003];
inline void ntt (ll a[], ll md) {
for (int i = 0; i < len; ++i) if (i < rev[i]) std:: swap(a[i], a[rev[i]]);
for (int l = 1; l < len; l <<= 1) {
ll W = power(md == 1 ? G : GI, (MODN - 1) / (l << 1));
for (int i = 0; i < len; i += (l << 1)) {
ll w = 1;
for (int j = 0; j < l; ++j) {
ll u = a[i + j], v = w * a[i + j + l] % MODN;
a[i + j] = (u + v) % MODN;
a[i + j + l] = (u - v + MODN) % MODN;
w = w * W % MODN;
}
}
}
}
inline void nttmul (ll a[], ll b[], ll c[]) {
for (int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
ntt(a, 1), ntt(b, 1);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i] % MODN;
ntt(c, -1);
ll invlen = power(len, MODN - 2);
for (int i = 0; i < len; ++i) c[i] = c[i] * invlen % MODN;
}
inline void kagari () {
std:: string p, q; std:: cin >> p >> q; n = p.size(), m = q.size();
std:: reverse(p.begin(), p.end());
std:: reverse(q.begin(), q.end());
for (int i = 0; i < n; ++i) a[i] = p[i] - '0';
for (int i = 0; i < m; ++i) b[i] = q[i] - '0';
len = 1; while (len <= n + m) len <<= 1;
nttmul(a, b, c);
std:: string ans = "";
ll r = 0;
for (int i = 0; i < len; ++i) {
c[i] += r;
r = c[i] / 10;
c[i] %= 10;
ans.push_back(c[i] + '0');
}
while (ans.size() && ans.back() == '0') ans.pop_back();
std:: reverse(ans.begin(), ans.end());
std:: cout << ans << '\n';
return;
}
牛顿迭代法 | Newton's Method
在此介绍牛顿迭代法。
回到高等数学中,先考虑如何求一个非线性方程 \(f(x)=0\) 的解 \(x\)
我们任取一个点 \(x_0\),这时,\(x_1 = x_0 - \dfrac{f(x_0)}{f'(x_0)}\) 一定使得 \(f(x_1)\) 比 \(f(x_0)\) 更趋于 \(0\)。
以此类推,利用递推式 \(x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}\),无论 \(x_0\) 如何取值,重复足够多次后都可以取得近似解。
再将其应用于多项式中
在多项式运算里,我们往往是在模 \(x^n\) 意义下,求解一个未知的多项式 \(B(x)\)。这可以改写为:
假设我们已知 \(\mod{x^{\left\lfloor \frac{n}{2} \right\rfloor}}\) 时的 \(B_0(x)\)。
此时借用牛顿迭代公式,将 \(B(x)\) 当成变量 \(x_n\):
上式证明略。
而 \(B(x) \equiv 0 \pmod x^0\) 一般是极易求得的。
如此,只需倍增枚举取模次数,就可以 \(O(n \log n)\) 递推解决问题了。
应用
然而在求多项式 \(\ln\) 时是不需要使用牛顿迭代法的。在求多项式逆元 \(inv\),计算 \(exp\)等情况时则会用到。
以多项式求逆为例,我们要找 \(A(x)\) 的逆元 \(B(x)\),也就是满足 \(A \cdot B \equiv 1\),移项构造方程:
对未知多项式 \(B\) 求导,把 \(A\) 看作常数:
- 函数:\(F(B) = B^{-1} - A\)
- 导数:\(F'(B) = -B^{-2} = -\frac{1}{B^2}\)
现在,把 \(F(B_0)\) 和 \(F'(B_0)\) 塞进多项式牛顿迭代公式里:
化简即得:
注意!!!!!!!!!!!!!!!!!
时刻注意清零,不要将污染了的数据带入下一重计算。
同时注意数组大小,防止越界与 MLE
P4238 【模板】多项式乘法逆
题目大意:给定一个 \(n \leq 10^5\) 次多项式 \(F(x)\) ,求一个多项式 \(G(x)\), 使得 \(F(x) * G(x) \equiv 1 \pmod{x^n}\),系数对 \(998244353\) 取模。
通过数学关系推导,转化为使用 FFT/NTT 的方法。
已知:$$A B \equiv 1 \pmod{x^{\left\lfloor \frac{n}{2} \right\rfloor}}$$
而又有:$$A B' \equiv 1 \pmod{x^{\left\lfloor \frac{n}{2} \right\rfloor}}$$
相减得:$$B' - B \equiv 0 \pmod{x^{\left\lfloor \frac{n}{2} \right\rfloor}}$$
平方,得:$$B'^2 - 2BB' + B^2 \equiv 0 \pmod{x^n}$$
等式左右同乘 \(A\):$$A B'^2 - 2B' + B \equiv 0 \pmod{x^n}$$
移项得到递推式:$$B \equiv 2B' - A B'^2 \pmod{x^n}$$
也就是求 $$ f_n = 2f_{n-1}-A \cdot f_{n-1}^2$$
初始状态 \(f_0 = inv(A_0)\)。
由主定理知,此方法时间复杂度为 \(O(n \log n)\)。
注意在算 \(A B \equiv 1 \pmod{x^k}\) 的 \(B\) 时,\(B\) 所有 \(>k\) 次的系数均为 \(0\),否则会发生错误。
inline void ntt (ll a[], ll md) {
for (int i = 0; i < len; ++i) if (i < rev[i]) std:: swap(a[i], a[rev[i]]);
for (int l = 1; l < len; l <<= 1) {
ll W = power(md == 1 ? G : GI, (MODN - 1) / (l << 1));
for (int i = 0; i < len; i += (l << 1)) {
ll w = 1;
for (int j = 0; j < l; ++j) {
ll u = a[i + j], v = w * a[i + j + l] % MODN;
a[i + j] = (u + v) % MODN, a[i + + j + l] = (u - v + MODN) % MODN;
w = w * W % MODN;
}
}
}
if (md == -1) {
ll invlen = power(len, MODN - 2);
for (int i = 0; i < len; ++i) a[i] = a[i] * invlen % MODN;
}
}
inline void kagari () {
scanf("%d", &n);
for (int i = 0; i < n; ++i) scanf("%lld", &a[i]);
b[0] = inv(a[0]);
for (int k = 1; k < (n << 1); k <<= 1) {
for (int i = 0; i <= (k >> 1); ++i) d[i] = b[i];
len = 1;
while (len <= (k << 1)) len <<= 1;
for (int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
for (int i = (k >> 1) + 1; i < len; ++i) b[i] = 0;
for (int i = 0; i <= k; ++i) c[i] = a[i]; // note: regard ai = 0 when i > k
for (int i = k + 1; i < len; ++i) c[i] = 0;
ntt(c, 1), ntt(b, 1);
for (int i = 0; i < len; ++i) b[i] = c[i] * b[i] % MODN * b[i] % MODN;
ntt(b, -1);
for (int i = 0; i < len; ++i) b[i] = (d[i] * 2ll % MODN - b[i] + MODN) % MODN;
}
for (int i = 0; i < n; ++i) printf("%lld ", b[i]);
puts("");
return;
}
P4725 【模板】多项式对数函数(多项式 ln)
要求 $$ G = \ln F \pmod{998244353}$$
对两边求导,有
也就是求
其中对 \(F\) 求导与对右边积分的时间复杂度都是 \(O(n)\),求 \(F\) 的逆与卷积都是 \(O(n \log n)\)。
#include <stdio.h>
#include <algorithm>
#define ll long long
const ll MODN = 998244353;
const ll G = 3;
const ll GI = 332748118;
inline ll power (ll x, ll k) {
ll res = 1;
while (k) {
if (k & 1) res = res * x % MODN;
x = x * x % MODN;
k >>= 1;
}
return res;
}
inline ll inv (const ll &x) {
return power(x % MODN, MODN - 2);
}
int n, len, rev[800003];
ll a[800003], b[800003], c[800003], d[800003], e[800003], f[800003];
inline void ntt (ll a[], ll md) {
for (int i = 0; i < len; ++i) if (i < rev[i]) std:: swap(a[i], a[rev[i]]);
for (int l = 1; l < len; l <<= 1) {
ll W = power(md == 1 ? G : GI, (MODN - 1) / (l << 1));
for (int i = 0; i < len; i += (l << 1)) {
ll w = 1;
for (int j = 0; j < l; ++j) {
ll u = a[i + j], v = w * a[i + j + l] % MODN;
a[i + j] = (u + v) % MODN, a[i + + j + l] = (u - v + MODN) % MODN;
w = w * W % MODN;
}
}
}
if (md == -1) {
ll invlen = power(len, MODN - 2);
for (int i = 0; i < len; ++i) a[i] = a[i] * invlen % MODN;
}
}
inline void inv (ll a[], ll b[]) {
b[0] = inv(a[0]);
for (int k = 1; k < (n << 1); k <<= 1) {
for (int i = 0; i <= (k >> 1); ++i) d[i] = b[i];
len = 1;
while (len <= (k << 1)) len <<= 1;
for (int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
for (int i = (k >> 1) + 1; i < len; ++i) b[i] = 0;
for (int i = 0; i <= k; ++i) c[i] = a[i]; // note: regard ai = 0 when i > k
for (int i = k + 1; i < len; ++i) c[i] = 0;
ntt(c, 1), ntt(b, 1);
for (int i = 0; i < len; ++i) b[i] = c[i] * b[i] % MODN * b[i] % MODN;
ntt(b, -1);
for (int i = 0; i < len; ++i) b[i] = (d[i] * 2ll % MODN - b[i] + MODN) % MODN;
}
for (int i = n; i < len; ++i) b[i] = 0; // 注意 i>n 时 bi 应为 0
}
inline void qiudao (ll a[], ll b[]) {
for (int i = 1; i < n; ++i) b[i - 1] = a[i] * (ll)i % MODN;
}
inline void jifen (ll a[], ll b[]) {
for (int i = 0; i < n; ++i) b[i + 1] = a[i] * inv(i + 1) % MODN;
}
inline void nttmul (ll a[], ll b[], ll c[]) {
len = 1;
while (len <= (n << 1)) len <<= 1;
for (int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
ntt(a, 1), ntt(b, 1);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i] % MODN;
ntt(c, -1);
}
inline void kagari () {
scanf("%d", &n);
for (int i = 0; i < n; ++i) scanf("%lld", &a[i]);
qiudao(a, e);
inv(a, b);
nttmul(e, b, b);
jifen(b, f);
for (int i = 0; i < n; ++i) printf("%lld ", f[i]);
puts("");
return;
}
int main () {
kagari();
return 0;
}
P4726 【模板】多项式指数函数(多项式 exp)
题目大意:给定一个 \(n-1\) 次多项式 \(A(x)\),其中常数项 \(a_0 = 0\),求多项式 \(B(x)\),使得:\(B(x) \equiv e^{A(x)} \pmod{x^n}\)
使用牛顿迭代方法。
要求 $$B(x) \equiv e^{A(x)} \pmod{x^n}$$
两边同时取对数,移项,则构造出一个关于 \(B(x)\) 的函数方程:
由牛顿迭代公式:
于是递推即可。注意常数较大,有时需优化。
P5245 【模板】多项式快速幂
题目大意:给定一个 \(n-1\) 次多项式 \(A(x)\),其中 \(a_0=1\),求一个在 \(\bmod\ x^n\) 意义下的多项式 \(B(x)\),使得 \(B(x) \equiv (A(x))^k \ (\bmod\ x^n)\)。系数对 \(998244353\) 取模。
对 \(B(x) \equiv (A(x))^k \ (\bmod\ x^n)\) 两边取对数:
那么 \(e^{k \ln A(x)}\) 即为所求。
注意 \(k\) 可能很大,直接将其对 \(998244353\) 即可。
P5205 【模板】多项式开根
题目大意:给定一个 \(n-1\) 次多项式 \(A(x)\),其中 \(a_0=1\),求一个在 \(\bmod x^n\) 意义下的多项式 \(B(x)\),使得 \(B^2(x) \equiv A(x) \pmod{x^n}\)。系数对 \(998244353\) 取模。
令 \(F(x) = B^2-A\),则有
套上 NTT 与逆元即可。
多项式大模板
| 函数 | 空间大小 | |
|---|---|---|
| $$ntt(a, b, N)$$ | $$4N$$ | |
| \(ntt(a, b, N)\) | \(4N\) | |
例题
P3723 [AHOI2017/HNOI2017] 礼物
FFT/NTT 模板题。
将式子改写为:
左边两项是定值,中间两项通过枚举 \(c \in [-m, m] \in [-100, 100]\) (表示第一个手环加上的值)即可得最小值。于是,除了最右边的项全都是定值。
所以我们只需要求 k 取何值时, \(\sum_{i=1}^n x_i y_{i+k}\) 最大就行了。
翻转 \(x\),式子变为 \(\sum_{i=1}^n x_{n-i+1} y_{i+k}\),这是熟悉的卷积的形式,使用 NTT 进行求解。
const ll MODN = 998244353;
const ll G = 3;
const ll GI = 332748118;
int n, len, rev[140003];
ll m, a[140003], b[140003], c[140003];
inline void ntt (ll a[], ll md) {
for (int i = 0; i < len; ++i) if (i < rev[i]) std:: swap(a[i], a[rev[i]]);
for (int l = 1; l < len; l <<= 1) {
ll W = power(md == 1 ? G : GI, (MODN - 1) / (l << 1));
for (int i = 0; i < len; i += (l << 1)) {
ll w = 1;
for (int j = 0; j < l; ++j) {
ll u = a[i + j], v = w * a[i + j + l] % MODN;
a[i + j] = (u + v) % MODN;
a[i + j + l] = (u - v + MODN) % MODN;
w = w * W % MODN;
}
}
}
}
inline void nttmul () {
for (int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
ntt(a, 1), ntt(b, 1);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i] % MODN;
ntt(c, -1);
ll invlen = power(len, MODN - 2);
for (int i = 0; i < len; ++i) c[i] = c[i] * invlen % MODN;
}
inline void kagari () {
scanf("%d %lld", &n, &m);
for (int i = 1; i <= n; ++i) scanf("%lld", &a[i]);
for (int i = 1; i <= n; ++i) scanf("%lld", &b[i]), b[i + n] = b[i];
ll res = 0, ans = MODN, r = 0;
for (int i = 1; i <= n; ++i) r += a[i] - b[i];
for (ll i = -m; i <= m; ++i) ans = std:: min(ans, i * i * (ll)n + i * 2ll * r);
for (int i = 1; i <= n; ++i) ans += a[i] * a[i] + b[i] * b[i];
len = 1;
while (len <= (n << 1)) len <<= 1;
for (int i = 1; i <= (n >> 1); ++i) std:: swap(a[i], a[n - i + 1]);
nttmul();
for (int i = n + 1; i <= (n << 1); ++i) res = std:: max(res, c[i]);
ans = (ans - 2ll * res + MODN) % MODN;
printf("%lld\n", ans);
return;
}
P3338 [ZJOI2014] 力
较简单的题目。
题目大意:给出 \(n\) 个数 \(a_1,a_2, \dots a_n, n \leq 10^5\),对 \(1 \leq i \leq n\),求 \(E_i=\dfrac{F_i}{q_i}\),其中\(F_j~=~\sum_{i = 1}^{j - 1} \dfrac{a_i \times a_j}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \dfrac{a_i \times a_j}{(i - j)^2}\)。
即求
联想到 FFT 的形式 \(C_k= \sum A_i \times B_{k-i}\),令 \(b_i=\dfrac{1}{i^2}\),则有
注意到上式的左式虽易得,但右式并不好求,下标 \(i - j\) 有时为负,不妨扩大 \(b\) 的范围,将上式与 \(b_i\) 改写为:
答案就是 \(F_n+1\) 到 \(F_{2n}\)。
#include <stdio.h>
#include <algorithm>
#include <math.h>
const double PI = acos(-1.0);
struct CPX {
double x, y; // x+yi
CPX (double x = 0.0, double y = 0.0) : x(x), y(y) { }
CPX operator + (const CPX &t) { return { x + t.x, y + t.y }; }
CPX operator - (const CPX &t) { return { x - t.x, y - t.y }; }
CPX operator * (const CPX &t) { return { x * t.x - y * t.y , x * t.y + y * t.x }; }
} a[800003], b[800003], c[800003], d[800003], e[800003], f[800003]; // note that N >= 2^upper(logn)
int n, rev[800003];
inline void fft (CPX a[], double md, int len) {
for (int i = 0; i < len; ++i) if (i < rev[i]) std:: swap(a[i], a[rev[i]]);
for (int l = 1; l < len; l <<= 1) {
CPX W = { cos(PI / l), md * sin(PI / l) };
for (int i = 0; i < len; i += (l << 1)) {
CPX w(1.0, 0.0);
for (int j = 0; j < l; ++j) {
CPX u = a[i + j];
CPX v = w * a[i + j + l];
a[i + j] = u + v;
a[i + j + l] = u - v;
w = w * W;
}
}
}
}
inline void fftmul (CPX a[], CPX b[], CPX c[], int N) {
int len = 1; while (len <= (N << 1)) len <<= 1;
for (int i = 0; i < len; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (len >> 1) : 0);
fft(a, 1, len);
fft(b, 1, len);
for (int i = 0; i < len; ++i) c[i] = a[i] * b[i];
fft(c, -1, len);
for (int i = 0; i < len; ++i) c[i].x /= len, c[i].y = 0;
}
inline void kagari () {
scanf("%d", &n);
for (int i = 1; i <= n; ++i) scanf("%lf", &a[i].x);
for (int i = 1; i <= n; ++i) b[n + i].x = 1.0 / i / i;
for (int i = -1; i >= -n; --i) b[n + i].x = -1.0 / i / i;
b[n].x = 0.0;
fftmul(a, b, c, n << 1);
for (int i = n + 1; i <= n * 2; ++i) printf("%f\n", c[i].x);
return;
}
int main () {
kagari();
return 0;
}
P5488【模板】差分与前缀和
题目大意:给定一个长为 \(n\) 的序列 \(a\),求出其 \(k\) 次前缀和与差分(即前缀和/差分 \(k\) 次后的结果)。系数对 \(1004535809\) 取模。
在此以生成函数的方式来说明此题:
- 先考虑如何求前缀和 \(S(x)\)
令生成函数 \(A(x) = a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}\),进行一次前缀和的结果为 \(s_0, s_1, \dots, s_{n-1}\)。
对前缀和函数 \(S(x) = s_0+s_1x+s_2x^2+\dots+s_{n-1}x^{n-1}\)
两边乘以 \(x\),有 \(xS(x) = s_0x+s_1x^2+s_2x^3+\dots+s_{n-1}x^n\)
相减,有 \((1-x)S(x) = a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1} = A(X)\)。
于是 \(S(X)=\dfrac{A(X)}{1-x}\)。则 \(k\) 次前缀和就是 \(S_k(x) = \dfrac{A(X)}{(1-x)^k}=A(x)\times(1-x)^{-k}\)。
- 现在考虑差分函数 \(C(x) = c_0+c_1x+c_2x^2+\dots+c_{n-1}x^{n-1}\),其中 \(c_i = a_i - a_{i-1}\)
思考如何凑出 \(a_i - a_{i-1}\),考虑到对 \(A(x)\) 进行操作
那么要求 \(k\) 次差分,即求 \(C_k(x) = A(x)\times(1-x)^k\)。
综上,我们套用 NTT,多项式求逆与快速幂即可。
事实上,前缀和的式子也可以是 \(S(x) = A(x) \times (1+x+\dots+x^{n-1})\),而且感觉这个式子是更容易想到的。
组合数学方式的进一步优化,可提供一种时间优化方式 by gemini
直接去算多项式幂 \((1-x)^K\) 或 \((1-x)^{-K}\),常数极大,效率很低。但这道题里的多项式只有两项 \((1-x)\),我们可以直接利用二项式定理展开!
假设我们要求的多项式为 \(B(x) = \sum_{i=0}^{N-1} b_i x^i\):
- 情况一:\(K\) 阶前缀和 —— 计算 \((1-x)^{-K}\)
根据负整数次幂的广义二项式定理展开:
其第 \(i\) 项系数 \(b_i = \binom{K+i-1}{i}\)。
观察相邻两项的关系:
已知 \(b_0 = 1\),我们可以用 \(O(N)\) 的线性递推直接求出整个 \(B(x)\) 多项式:
- 情况二:\(K\) 阶差分 —— 计算 \((1-x)^K\)
根据二项式定理展开:
其第 \(i\) 项系数 \(b_i = (-1)^i \binom{K}{i}\)。
观察相邻两项的关系:
已知 \(b_0 = 1\),同样可以用 \(O(N)\) 线性递推:
inline void sum (const vell &a, vell &b, const std:: string &K, int N) { // K 次前缀
vell c(N << 3, 0), d(N << 3), e(N << 3);
for (int i = 0; i < N; ++i) d[i] = 1;
// or: c[0] = 1, c[1] = MODN - 1; inv(c, d, N);
power(d, e, K, N);
nttmul(a, e, b, N);
}
inline void fini (const vell &a, vell &b, const std:: string &K, int N) { // k 次差分
vell c(N << 3, 0), d(N << 3);
c[0] = 1, c[1] = MODN-1; // 注意不是 -1 而是 MODN - 1!!
power(c, d, K, N);
nttmul(a, d, b, N);
}
P5641 【CSGRound2】开拓者的卓识
题目大意:给定序列 \(a\),记一个序列 \([l,r]\) 的 \(k\) 阶子段和为 \(sum_{k,l,r}\),有\(sum_{k,l,r}=\begin{cases}\sum\limits_{i=l}^{r}a_i&,k=1\\\sum\limits_{i=l}^{r}\sum\limits_{j=i}^{r}sum_{k-1,i,j}&,k\geq 2\end{cases}\)。对所有 \(1 \leq i \leq n\),求 \(sum_{k,1,i}\)。答案对 \(998244353\) 取模。
from https://www.luogu.com.cn/article/bev8cd54 :
考虑 \(a_i\) 对答案的贡献。
\(\mathit{sum}_{1,l,r}\) 就是从 \(a_l\) 加到 \(a_r\)。
\(\mathit{sum}_{2,l,r}\) 根据 \(\mathit{sum}_{k,l,r}\) 的公式,可以看作是选择一个区间 \([l_1,r_1] \subseteq [l,r]\),再将 \(a_{l_1}\) 到 \(a_{r_1}\) 的和加入答案中。
\(\mathit{sum}_{3,l,r}\) 是选择一个区间 \([l_1,r_1] \subseteq [l,r]\),再对 \([l_1,r_1]\) 求一遍 \(\mathit{sum}_{2,l_1,r_1}\),也就是再选择一个区间 \([l_2,r_2] \in [l_1,r_1]\) 将 \(a_{l_2}\) 到 \(a_{r_2}\) 的和加入答案中。
以此类推,\(\mathit{sum}_{k,l,r}\) 就是选择 \(k-1\) 个区间 \([l_{k-1},r_{k-1}] \subseteq [l_{k-2},r_{k-2}] \subseteq \dots \subseteq [l_1,r_1] \subseteq [l,r]\),并将 \([l_k,r_k]\) 的和加入答案中。
考虑转换贡献体,枚举 \(a_i\) 会对 \(\mathit{sum}_{k,1,r}\) 产生多大的贡献,也就是有多少组区间 \([l_i,r_i]\) 满足 \(i \in [l_{k-1},r_{k-1}] \subseteq [l_{k-2},r_{k-2}] \subseteq \dots \subseteq [l_1,r_1] \subseteq [l_0=1,r_0=r]\),由于是相邻区间之间是包含关系,故 \(i \ge l_{k-1} \ge l_{k-2} \ge \dots \ge l_1 \ge 1\),记 \(d_k = i - l_{k-1},d_{k-1} = l_{k-1} - l_{k-2},\dots,d_2 = l_2 - l_1,d_1 = l_1 - 1\),那么有 \(d_1 + d_2 + d_3 + \dots + d_k = i - 1\),其中 \(d_i \ge 0\),采用隔板法可求得所有 \((l_1,l_2,\dots,l_{k-1})\) 的个数为\(\binom{i + k - 2}{k - 1}\)。同理可得所有 \((r_1,r_2,\dots,r_{k-1})\) 的个数为\(\binom{r - i + k - 1}{k - 1}\)。
故 \(a_i\) 对 \(r\) 产生的贡献为
也就是要对每个 \(1 \leq r \leq n\) 求
注意到其中类似卷积的形式,记
那么
计算组合数 \(B_i\) ,有递推式
其中 \(m = k-1\),\(B_1=\binom{k-1}{k-1}=1\)。
搞清楚了会发现代码很短。
inline void kagari () {
int n; ll k; scanf("%d %lld", &n, &k);
vell a(n << 3), b(n << 3), c(n << 3);
for (int i = 1; i <= n; ++i) scanf("%lld", &a[i]);
b[1] = 1;
for (int i = 1; i < n; ++i) b[i + 1] = b[i] * (ll)(k - 1 + i) % MODN * inv(i) % MODN;
for (int i = 1; i <= n; ++i) a[i] = a[i] * b[i] % MODN;
nttmul(a, b, c, n + 1);
for (int i = 1; i <= n; ++i) printf("%lld ", c[i + 1]);
puts("");
return;
}

浙公网安备 33010602011771号