未完成的 · 初三多项式的运算练习 题解
（标题是中二属性溢出导致的，来源于《未完成的·乐章》。）
美好的下午时光要拿来写题解呜呜呜，一篇一篇地鸽得了，大概写完了捏。
有些题要用到 GF 的知识，或许我可以找时间讲一下？
说真的，不理解多项式卷积的组合意义的是真的没法做 2，5 题，真别以为 GF 是什么高深的玩意，你把这些元素塞到次数里面那就是一个经典的用 OGF 计数背包方案的做法，别再说什么我没用 GF 就做出来这题了，你凭你直觉打出来的那玩意就是个 OGF。
贴一份我的 FFT 和 NTT 的板子。
FFT：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int n, m, limit, f[1 << 22], g[1 << 22], h[1 << 22];
namespace poly {
typedef complex<double> cplx;
const double pi = acos(-1.0);
void FFT(cplx *f, int limit, double type) {
static int r[1 << 22];
for(int i = 1, l = __lg(limit); i < limit; ++i) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
if(i < r[i]) swap(f[i], f[r[i]]);
}
cplx w, wn, f1, f2;
for(int mid = 1; mid < limit; mid <<= 1) {
wn = cos(pi / mid) + type * 1i * sin(pi / mid);
for(int i = 0; i < limit; i += mid << 1) {
w = 1.0;
for(int j = 0; j < mid; ++j, w *= wn) {
f1 = f[i + j], f2 = w * f[i + mid + j];
f[i + j] = f1 + f2;
f[i + mid + j] = f1 - f2;
}
}
}
}
void mul(int* F, int n, int *G, int m, int *H) {
static cplx f[1 << 22], g[1 << 22];
int l = __lg(n + m) + 1;
int limit = 1 << l;
fill(f, f + limit, 0);
fill(g, g + limit, 0);
for(int i = 0; i <= n; ++i) f[i].real(F[i]);
for(int i = 0; i <= m; ++i) g[i].real(G[i]);
FFT(f, limit, 1.0);
FFT(g, limit, 1.0);
for(int i = 0; i < limit; ++i) f[i] *= g[i];
FFT(f, limit, -1.0);
for(int i = 0; i <= n + m; ++i) H[i] = int(f[i].real() / limit + 0.5);
}
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n >> m;
for(int i = 0; i <= n; ++i) cin >> f[i];
for(int i = 0; i <= m; ++i) cin >> g[i];
poly::mul(f, n, g, m, h);
for(int i = 0; i <= n + m; ++i) cout << h[i] << " ";
return 0;
}
NTT：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
namespace poly {
constexpr int mod = 998244353, G = 3, gi = 332748118, img = 86583718, inv2 = 499122177, inv2i = 954952494;
struct modint {
int x;
modint(ll v = 0): x((v % mod + mod) % mod) {}
friend modint ksm(modint x, ll y) {
modint ret = 1;
while(y) {if(y & 1) ret *= x; x *= x, y >>= 1;}
return ret;
}
friend modint operator + (const modint& x, const modint& y) {return x.x + y.x >= mod ? x.x + y.x - mod : x.x + y.x;}
friend modint operator - (const modint& x, const modint& y) {return x.x < y.x ? x.x - y.x + mod : x.x - y.x;}
friend modint operator * (const modint& x, const modint& y) {return (ll)x.x * y.x % mod;}
friend modint operator / (const modint& x, const modint& y) {return x * ksm(y, mod - 2);}
friend modint operator - (const modint& x) {return mod - x.x;}
modint operator = (const modint& _) {x = _.x; return *this;}
modint operator += (const modint& _) {*this = *this + _; return *this;}
modint operator -= (const modint& _) {*this = *this - _; return *this;}
modint operator *= (const modint& _) {*this = *this * _; return *this;}
modint operator /= (const modint& _) {*this = *this / _; return *this;}
bool operator == (const modint& _) const {return x == _.x;}
bool operator != (const modint& _) const {return x != _.x;}
friend istream& operator >> (istream& in, modint& _) {static ll tmp; in >> tmp; _.x = tmp % mod; return in;}
friend ostream& operator << (ostream& out, const modint& _) {return out << _.x;}
};
modint ksm(ll x, ll y) {return ksm(modint(x), y);}
int r[1 << 22];
modint power[2][23];
void init(int len) {
int l = __lg(len) + 1;
int limit = 1 << l;
for(int i = 1; i < limit; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
power[0][l] = ksm(gi, (mod - 1) / limit);
power[1][l] = ksm(G, (mod - 1) / limit);
for(int i = l; i >= 1; --i) {
power[0][i - 1] = power[0][i] * power[0][i];
power[1][i - 1] = power[1][i] * power[1][i];
}
}
void NTT(modint* f, int n, bool type) {
for(int i = 1; i < n; ++i) {
if(i < r[i]) swap(f[i], f[r[i]]);
}
modint w, wn, h1, h2;
for(int mid = 1, lg = 0; mid < n; mid <<= 1, ++lg) {
wn = power[type][lg + 1];
for(int i = 0; i < n; i += mid << 1) {
w = 1;
for(int j = 0; j < mid; ++j, w *= wn) {
h1 = f[i + j], h2 = w * f[i + mid + j];
f[i + j] = h1 + h2, f[i + mid + j] = h1 - h2;
}
}
}
if(type) return;
modint inv = ksm(n, mod - 2);
for(int i = 0; i < n; ++i) f[i] *= inv;
}
void mul(modint* F, int n, modint *G, int m, modint *H) {
static modint f[1 << 22], g[1 << 22];
int l = __lg(n + m) + 1;
int limit = 1 << l;
copy(F, F + limit, f);
copy(G, G + limit, g);
init(n + m);
NTT(f, limit, true);
NTT(g, limit, true);
for(int i = 0; i < limit; ++i) H[i] = f[i] * g[i];
NTT(H, limit, false);
}
}
using poly::modint;
int n, m;
modint f[1 << 22], g[1 << 22], h[1 << 22];
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n >> m;
for(int i = 0; i <= n; ++i) cin >> f[i];
for(int i = 0; i <= m; ++i) cin >> g[i];
poly::mul(f, n, g, m, h);
for(int i = 0; i <= n + m; ++i) cout << h[i] << " ";
return 0;
}
力老师
先把 \(q_{j}\) 从 \(F_{j}\) 中消掉。
再把 \(F_{j}\) 分成两部分算，一个 \(\sum\limits_{i = 1}^{j - 1}\dfrac{q_{i}}{(i - j)^{2}}\) 一个 \(-\sum\limits_{i = j + 1}^{n}\dfrac{q_{i}}{(i - j)^{2}}\)。
令多项式 \(Q(x) = \sum\limits_{i = 1}^{n}q_{i}x^{i}\)，也就是 \([x^{i}]Q(x) = q_{i}\)，特别的，\([x^{0}]Q(x) = 0\)。
令多项式 \(G(x) = \sum\limits_{i = 1}^{n}\dfrac{x^{i}}{i^{2}}\)，也就是 \([x^{i}]G(x) = \dfrac{1}{i^{2}}\)，特别的，\([x^{0}]G(x) = 0\)。
然后把上面的式子替换成多项式的系数的运算，有：
\[\begin{aligned}
\sum\limits_{i = 1}^{j - 1}\dfrac{q_{i}}{(i - j)^{2}} &= \sum\limits_{i = 1}^{j - 1}[x^{i}]Q(x) \times [x^{j - i}]G(x)\\
&= \sum\limits_{i = 0}^{j}[x^{i}]Q(x) \times [x^{j - i}]G(x)\\
\sum\limits_{i = j + 1}^{n}\dfrac{q_{i}}{(i - j)^{2}} &= \sum\limits_{i = j + 1}^{n}[x^{i}]Q(x) \times [x^{i - j}]G(x)\\
&= \sum\limits_{i = j}^{n}[x^{i}]Q(x) \times [x^{i - j}]G(x)
\end{aligned}
\]
可以发现上面的部分就是一个卷积的形式，下面的部分是一个「差卷积」，把其中一个多项式的次数反转一下即可（即令 \([x^{i}]H(x) = [x^{n - i + 1}]G(x)\)）。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
constexpr int lim = 1 << 18;
namespace poly {
typedef complex<double> cplx;
const double pi = acos(-1.0);
void FFT(cplx *f, int limit, double type) {
static int r[lim];
for(int i = 1, l = __lg(limit); i < limit; ++i) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
if(i < r[i]) swap(f[i], f[r[i]]);
}
cplx w, wn, f1, f2;
for(int mid = 1; mid < limit; mid <<= 1) {
wn = cos(pi / mid) + type * 1i * sin(pi / mid);
for(int i = 0; i < limit; i += mid << 1) {
w = 1.0;
for(int j = 0; j < mid; ++j, w *= wn) {
f1 = f[i + j], f2 = w * f[i + mid + j];
f[i + j] = f1 + f2;
f[i + mid + j] = f1 - f2;
}
}
}
}
void mul(double* F, int n, double *G, int m, double *H) {
static cplx f[lim], g[lim];
int l = __lg(n + m) + 1;
int limit = 1 << l;
fill(f, f + limit, 0);
fill(g, g + limit, 0);
for(int i = 0; i <= n; ++i) f[i].real(F[i]);
for(int i = 0; i <= m; ++i) g[i].real(G[i]);
FFT(f, limit, 1.0);
FFT(g, limit, 1.0);
for(int i = 0; i < limit; ++i) f[i] *= g[i];
FFT(f, limit, -1.0);
for(int i = 0; i <= n + m; ++i) H[i] = f[i].real() / limit;
}
}
int n;
double q[lim], i2[lim], res1[lim], res2[lim];
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n;
for(int i = 1; i <= n; ++i) {
cin >> q[i];
i2[i] = 1.0 / i / i;
}
poly::mul(q, n, i2, n, res1);
reverse(q + 1, q + 1 + n);
poly::mul(q, n, i2, n, res2);
for(int i = 1; i <= n; ++i) {
cout << fixed << setprecision(3) << res1[i] - res2[n - i + 1] << '\n';
}
return 0;
}
Triple
芙芙不讲武德捏。
用到了一点点的多项式卷积的组合意义（我不说是 GF 了好吧），然后就是容斥。
具体来说，把每把斧头的价值塞到多项式的次数里面，发现让其相乘就可以使价值相加，并且其系数就是得到这个价值的方案数。
假设两把斧头价值为 \(v1\) 和 \(v2\)，那么对于多项式 \(F(x)\) 有 \([x^{v1}]F(x) = [x^{v2}]F(x) = 1\)，那么让它自己卷一下得到的结果中必然会有 \([x^{v1}]F(x) \times [x^{v2}]F(x) = 1\) 这一贡献，并且贡献在 \([x^{v1 + v2}]F(x)\) 这里，其它选择方法同样会造成贡献。
但是这样会算重，直接大力/暴力容斥一下即可。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
constexpr int lim = 1 << 20;
namespace poly {
constexpr ll mod = 79164837199873, G = 5, gi = 47498902319924;
struct modint {
ll x;
modint(ll v = 0): x((v % mod + mod) % mod) {}
friend modint ksm(modint x, ll y) {
modint ret = 1;
while(y) {if(y & 1) ret *= x; x *= x, y >>= 1;}
return ret;
}
friend modint operator + (const modint& x, const modint& y) {return x.x + y.x >= mod ? x.x + y.x - mod : x.x + y.x;}
friend modint operator - (const modint& x, const modint& y) {return x.x < y.x ? x.x - y.x + mod : x.x - y.x;}
friend modint operator * (const modint& x, const modint& y) {return (__int128)x.x * y.x % mod;}
friend modint operator / (const modint& x, const modint& y) {return x * ksm(y, mod - 2);}
friend modint operator - (const modint& x) {return mod - x.x;}
modint operator = (const modint& _) {x = _.x; return *this;}
modint operator += (const modint& _) {*this = *this + _; return *this;}
modint operator -= (const modint& _) {*this = *this - _; return *this;}
modint operator *= (const modint& _) {*this = *this * _; return *this;}
modint operator /= (const modint& _) {*this = *this / _; return *this;}
bool operator == (const modint& _) const {return x == _.x;}
bool operator != (const modint& _) const {return x != _.x;}
friend istream& operator >> (istream& in, modint& _) {static ll tmp; in >> tmp; _.x = tmp % mod; return in;}
friend ostream& operator << (ostream& out, const modint& _) {return out << _.x;}
};
modint ksm(ll x, ll y) {return ksm(modint(x), y);}
int r[lim];
modint power[2][23];
void init(int len) {
int l = __lg(len) + 1;
int limit = 1 << l;
for(int i = 1; i < limit; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
power[0][l] = ksm(gi, (mod - 1) / limit);
power[1][l] = ksm(G, (mod - 1) / limit);
for(int i = l; i >= 1; --i) {
power[0][i - 1] = power[0][i] * power[0][i];
power[1][i - 1] = power[1][i] * power[1][i];
}
}
void NTT(modint* f, int n, bool type) {
static modint w, wn, h1, h2;
for(int i = 1; i < n; ++i) {
if(i < r[i]) swap(f[i], f[r[i]]);
}
for(int mid = 1, lg = 1; mid < n; mid <<= 1, ++lg) {
wn = power[type][lg];
for(int i = 0; i < n; i += mid << 1) {
w = 1;
for(int j = 0; j < mid; ++j, w *= wn) {
h1 = f[i + j], h2 = w * f[i + mid + j];
f[i + j] = h1 + h2, f[i + mid + j] = h1 - h2;
}
}
}
if(type) return;
modint inv = ksm(n, mod - 2);
for(int i = 0; i < n; ++i) f[i] *= inv;
}
void mul(modint* F, int n, modint *G, int m, modint *H) {
static modint f[lim], g[lim];
int l = __lg(n + m) + 1;
int limit = 1 << l;
copy(F, F + limit, f);
copy(G, G + limit, g);
init(n + m);
NTT(f, limit, true);
NTT(g, limit, true);
for(int i = 0; i < limit; ++i) H[i] = f[i] * g[i];
NTT(H, limit, false);
if(n + m < limit) fill(H + n + m + 1, H + limit, 0);
}
void polyinv(modint* F, int n, modint* H) {
static modint f[lim], g[lim];
fill(f, f + (n << 1), 0);
copy(F, F + n + 1, f);
fill(H, H + (n << 1), 0);
H[0] = 1 / f[0];
for(int len = 2, limit; len <= (n << 1); len <<= 1) {
limit = len << 1;
copy(f, f + len, g);
if(len > n + 1) fill(g + n + 1, g + len, 0);
fill(H + len, H + limit, 0);
init(len);
NTT(g, limit, true);
NTT(H, limit, true);
for(int i = 0; i < limit; ++i) H[i] *= 2 - g[i] * H[i];
NTT(H, limit, false);
fill(H + len, H + limit, 0);
}
fill(H + n + 1, H + (1 <<__lg(n << 1)), 0);
}
}
using poly::modint;
ll n, a, ans[120005];
modint f1[lim], f2[lim], f3[lim], res[lim], sub[lim];
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n;
for(int i = 1; i <= n; ++i) {
cin >> a;
f1[a] = 1;
f2[a * 2] = 1;
f3[a * 3] = 1;
}
poly::mul(f1, 40000, f2, 80000, sub);
for(int i = 1; i <= 40000; ++i) ans[i] += f1[i].x;
poly::mul(f1, 40000, f1, 40000, res);
for(int i = 1; i <= 80000; ++i) ans[i] += (res[i].x - f2[i].x) >> 1;
poly::mul(f1, 40000, res, 80000, res);
for(int i = 1; i <= 120000; ++i) ans[i] += (res[i].x - 3 * (sub[i].x - f3[i].x) - f3[i].x) / 6;
for(int i = 1; i <= 120000; ++i) {
if(ans[i]) cout << i << " " << ans[i] << '\n';
}
return 0;
}
万径人踪灭
首先不管「不能是连续子序列」这个条件。
考虑枚举对称轴（类似于 manacher，在两两字符之间加一个 #），实际上就是求 \(2^{k} - 1\)，\(k\) 为两侧对称的字母数量。
令枚举的对称轴位置为 \(i\)，则 \(k = \sum\limits_{1 \leqslant i - j \leqslant i + j \leqslant |S|}[S_{i - j} = S_{i + j}]\)。
发现 \((i - j) + (i + j) = 2i\)，这是一个卷积的形式！
注意到 \(0 \times 0 = 0, 0 \times 1 = 0, 1 \times 1 = 1\)。
所以可以令 \([x^{i}]F(x) = [S_{i} == \texttt{a}]\)。
让 \(F(x)\) 自己卷一下，就可以统计出关于对称轴 \(i\) 对称且为 \(\texttt{a}\) 的位置数量，\(\texttt{b}\) 的情况同理。
那怎么去掉连续子序列的情况？manacher 即可。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
constexpr int lim = 1 << 20;
namespace poly {
constexpr int mod = 998244353, G = 3, gi = 332748118, img = 86583718, inv2 = 499122177, inv2i = 954952494;
struct modint {
int x;
modint(ll v = 0): x((v % mod + mod) % mod) {}
friend modint ksm(modint x, ll y) {
modint ret = 1;
while(y) {if(y & 1) ret *= x; x *= x, y >>= 1;}
return ret;
}
friend modint operator + (const modint& x, const modint& y) {return x.x + y.x >= mod ? x.x + y.x - mod : x.x + y.x;}
friend modint operator - (const modint& x, const modint& y) {return x.x < y.x ? x.x - y.x + mod : x.x - y.x;}
friend modint operator * (const modint& x, const modint& y) {return (ll)x.x * y.x % mod;}
friend modint operator / (const modint& x, const modint& y) {return x * ksm(y, mod - 2);}
friend modint operator - (const modint& x) {return mod - x.x;}
modint operator = (const modint& _) {x = _.x; return *this;}
modint operator += (const modint& _) {*this = *this + _; return *this;}
modint operator -= (const modint& _) {*this = *this - _; return *this;}
modint operator *= (const modint& _) {*this = *this * _; return *this;}
modint operator /= (const modint& _) {*this = *this / _; return *this;}
bool operator == (const modint& _) const {return x == _.x;}
bool operator != (const modint& _) const {return x != _.x;}
friend istream& operator >> (istream& in, modint& _) {static ll tmp; in >> tmp; _.x = tmp % mod; return in;}
friend ostream& operator << (ostream& out, const modint& _) {return out << _.x;}
};
modint ksm(ll x, ll y) {return ksm(modint(x), y);}
int r[lim];
modint power[2][23];
void init(int len) {
int l = __lg(len) + 1;
int limit = 1 << l;
for(int i = 1; i < limit; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
power[0][l] = ksm(gi, (mod - 1) / limit);
power[1][l] = ksm(G, (mod - 1) / limit);
for(int i = l; i >= 1; --i) {
power[0][i - 1] = power[0][i] * power[0][i];
power[1][i - 1] = power[1][i] * power[1][i];
}
}
void NTT(modint* f, int n, bool type) {
static modint w, wn, h1, h2;
for(int i = 1; i < n; ++i) {
if(i < r[i]) swap(f[i], f[r[i]]);
}
for(int mid = 1, lg = 1; mid < n; mid <<= 1, ++lg) {
wn = power[type][lg];
for(int i = 0; i < n; i += mid << 1) {
w = 1;
for(int j = 0; j < mid; ++j, w *= wn) {
h1 = f[i + j], h2 = w * f[i + mid + j];
f[i + j] = h1 + h2, f[i + mid + j] = h1 - h2;
}
}
}
if(type) return;
modint inv = ksm(n, mod - 2);
for(int i = 0; i < n; ++i) f[i] *= inv;
}
void mul(modint* F, int n, modint *G, int m, modint *H) {
static modint f[lim], g[lim];
int l = __lg(n + m) + 1;
int limit = 1 << l;
copy(F, F + limit, f);
copy(G, G + limit, g);
init(n + m);
NTT(f, limit, true);
NTT(g, limit, true);
for(int i = 0; i < limit; ++i) H[i] = f[i] * g[i];
NTT(H, limit, false);
if(n + m < limit) fill(H + n + m + 1, H + limit, 0);
}
}
using poly::modint;
constexpr int md = 1000000007;
string str, s = "V";
int len, l = 0, r = -1, p[200005], res[lim];
ll ans;
modint f[lim], g[lim];
ll ksm(ll x, ll y) {
ll ret = 1;
while(y) {
if(y & 1) ret = ret * x % md;
x = x * x % md, y >>= 1;
}
return ret;
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> str;
for(const auto& i : str) s.push_back(i), s.push_back('V');
len = s.length();
for(int i = 0; i < len; ++i) {
if(i <= r) p[i] = min(p[l + r - i], r - i + 1);
while(i - p[i] >= 0 && i + p[i] < len && s[i - p[i]] == s[i + p[i]]) ++p[i];
if(i + p[i] - 1 > r) l = i - p[i] + 1, r = i + p[i] - 1;
ans -= (p[i] >> 1);
// cerr << "> " << i << " - " << p[i] << '\n';
}
ans = ans % md + md;
for(int i = 0; i < len; ++i) f[i] = (s[i] == 'a'), g[i] = (s[i] == 'b');
poly::mul(f, len - 1, f, len - 1, f);
poly::mul(g, len - 1, g, len - 1, g);
for(int i = 0; i < (len << 1); ++i) res[i] = (f[i].x + g[i].x) % md;
for(int i = 0; i < len; ++i) ans = (ans + ksm(2, (res[i << 1] + 1) >> 1) - 1 + md) % md;
cout << ans;
return 0;
}
两个串
这玩意是正则表达式嘛？
沿用上一题的套路！直接枚举字符集的每个字符，让 \(S(x)\) 和 \(T(x)\) 卷起来求出相同位置的数量，最后枚举每一个位置看是不是和 \(T\) 中非 \(\texttt{?}\) 字符数量相同。
让后你就喜提一份 \(\mathcal{O}(\Sigma n \log n)\) 的代码！实现得漂亮应该能过，我是卡着时限过的。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
constexpr int lim = 1 << 20;
namespace poly {
constexpr int mod = 998244353, G = 3, gi = 332748118, img = 86583718, inv2 = 499122177, inv2i = 954952494;
struct modint {
int x;
modint(ll v = 0): x((v % mod + mod) % mod) {}
friend modint ksm(modint x, ll y) {
modint ret = 1;
while(y) {if(y & 1) ret *= x; x *= x, y >>= 1;}
return ret;
}
friend modint operator + (const modint& x, const modint& y) {return x.x + y.x >= mod ? x.x + y.x - mod : x.x + y.x;}
friend modint operator - (const modint& x, const modint& y) {return x.x < y.x ? x.x - y.x + mod : x.x - y.x;}
friend modint operator * (const modint& x, const modint& y) {return (ll)x.x * y.x % mod;}
friend modint operator / (const modint& x, const modint& y) {return x * ksm(y, mod - 2);}
friend modint operator - (const modint& x) {return mod - x.x;}
modint operator = (const modint& _) {x = _.x; return *this;}
modint operator += (const modint& _) {*this = *this + _; return *this;}
modint operator -= (const modint& _) {*this = *this - _; return *this;}
modint operator *= (const modint& _) {*this = *this * _; return *this;}
modint operator /= (const modint& _) {*this = *this / _; return *this;}
bool operator == (const modint& _) const {return x == _.x;}
bool operator != (const modint& _) const {return x != _.x;}
friend istream& operator >> (istream& in, modint& _) {static ll tmp; in >> tmp; _.x = tmp % mod; return in;}
friend ostream& operator << (ostream& out, const modint& _) {return out << _.x;}
};
modint ksm(ll x, ll y) {return ksm(modint(x), y);}
int r[lim];
modint power[2][23];
void init(int len) {
int l = __lg(len) + 1;
int limit = 1 << l;
for(int i = 1; i < limit; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
power[0][l] = ksm(gi, (mod - 1) / limit);
power[1][l] = ksm(G, (mod - 1) / limit);
for(int i = l; i >= 1; --i) {
power[0][i - 1] = power[0][i] * power[0][i];
power[1][i - 1] = power[1][i] * power[1][i];
}
}
void NTT(modint* f, int n, bool type) {
static modint w, wn, h1, h2;
for(int i = 1; i < n; ++i) {
if(i < r[i]) swap(f[i], f[r[i]]);
}
for(int mid = 1, lg = 1; mid < n; mid <<= 1, ++lg) {
wn = power[type][lg];
for(int i = 0; i < n; i += mid << 1) {
w = 1;
for(int j = 0; j < mid; ++j, w *= wn) {
h1 = f[i + j], h2 = w * f[i + mid + j];
f[i + j] = h1 + h2, f[i + mid + j] = h1 - h2;
}
}
}
if(type) return;
modint inv = ksm(n, mod - 2);
for(int i = 0; i < n; ++i) f[i] *= inv;
}
void mul(modint* F, int n, modint *G, int m, modint *H) {
static modint f[lim], g[lim];
int l = __lg(n + m) + 1;
int limit = 1 << l;
copy(F, F + limit, f);
copy(G, G + limit, g);
init(n + m);
NTT(f, limit, true);
NTT(g, limit, true);
for(int i = 0; i < limit; ++i) H[i] += f[i] * g[i];
// NTT(H, limit, false);
// if(n + m < limit) fill(H + n + m + 1, H + limit, 0);
}
}
using poly::modint;
vector<int> ans;
ll n, m, sum;
string s, t;
modint f[lim], g[lim], h[lim];
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> s >> t;
reverse(t.begin(), t.end());
n = s.length() - 1, m = t.length() - 1;
for(char ch = 'a'; ch <= 'z'; ++ch) {
for(int i = 0; i <= n; ++i) f[i] = (s[i] == ch);
for(int i = 0; i <= m; ++i) g[i] = (t[i] == ch), sum += (t[i] == ch);
poly::mul(f, n, g, m, h);
}
int l = __lg(n + m) + 1;
int limit = 1 << l;
poly::NTT(h, limit, false);
if(n + m < limit) fill(h + n + m + 1, h + limit, 0);
for(int i = m; i <= n; ++i) {
if(h[i].x == sum) ans.push_back(i - m);
}
cout << ans.size() << '\n';
for(const auto& i : ans) cout << i << '\n';
return 0;
}
3-idiots
用到了一点点的多项式卷积的组合意义（我不说是 GF 了好吧），然后就是容斥。
具体地，类似于 Triple 那题，把三角形的边长丢到次数，方案数丢进系数，然后卷一下就可以得到组合之后的方案数了。
再利用 \(a + b > c(a \leqslant b \leqslant c)\) 的判定条件，从小到大枚举每条边作为 \(c\)，再分选一条/选两条相同的/选三条相同的分别计算后加起来一下就好了，小小地容斥一下也可以。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
constexpr int lim = 1 << 20;
namespace poly {
constexpr int mod = 998244353, G = 3, gi = 332748118, img = 86583718, inv2 = 499122177, inv2i = 954952494;
struct modint {
int x;
modint(ll v = 0): x((v % mod + mod) % mod) {}
friend modint ksm(modint x, ll y) {
modint ret = 1;
while(y) {if(y & 1) ret *= x; x *= x, y >>= 1;}
return ret;
}
friend modint operator + (const modint& x, const modint& y) {return x.x + y.x >= mod ? x.x + y.x - mod : x.x + y.x;}
friend modint operator - (const modint& x, const modint& y) {return x.x < y.x ? x.x - y.x + mod : x.x - y.x;}
friend modint operator * (const modint& x, const modint& y) {return (ll)x.x * y.x % mod;}
friend modint operator / (const modint& x, const modint& y) {return x * ksm(y, mod - 2);}
friend modint operator - (const modint& x) {return mod - x.x;}
modint operator = (const modint& _) {x = _.x; return *this;}
modint operator += (const modint& _) {*this = *this + _; return *this;}
modint operator -= (const modint& _) {*this = *this - _; return *this;}
modint operator *= (const modint& _) {*this = *this * _; return *this;}
modint operator /= (const modint& _) {*this = *this / _; return *this;}
bool operator == (const modint& _) const {return x == _.x;}
bool operator != (const modint& _) const {return x != _.x;}
friend istream& operator >> (istream& in, modint& _) {static ll tmp; in >> tmp; _.x = tmp % mod; return in;}
friend ostream& operator << (ostream& out, const modint& _) {return out << _.x;}
};
modint ksm(ll x, ll y) {return ksm(modint(x), y);}
int r[lim];
modint power[2][23];
void init(int len) {
int l = __lg(len) + 1;
int limit = 1 << l;
for(int i = 1; i < limit; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
power[0][l] = ksm(gi, (mod - 1) / limit);
power[1][l] = ksm(G, (mod - 1) / limit);
for(int i = l; i >= 1; --i) {
power[0][i - 1] = power[0][i] * power[0][i];
power[1][i - 1] = power[1][i] * power[1][i];
}
}
void NTT(modint* f, int n, bool type) {
static modint w, wn, h1, h2;
for(int i = 1; i < n; ++i) {
if(i < r[i]) swap(f[i], f[r[i]]);
}
for(int mid = 1, lg = 1; mid < n; mid <<= 1, ++lg) {
wn = power[type][lg];
for(int i = 0; i < n; i += mid << 1) {
w = 1;
for(int j = 0; j < mid; ++j, w *= wn) {
h1 = f[i + j], h2 = w * f[i + mid + j];
f[i + j] = h1 + h2, f[i + mid + j] = h1 - h2;
}
}
}
if(type) return;
modint inv = ksm(n, mod - 2);
for(int i = 0; i < n; ++i) f[i] *= inv;
}
void mul(modint* F, int n, modint *G, int m, modint *H) {
static modint f[lim], g[lim];
int l = __lg(n + m) + 1;
int limit = 1 << l;
copy(F, F + limit, f);
copy(G, G + limit, g);
init(n + m);
NTT(f, limit, true);
NTT(g, limit, true);
for(int i = 0; i < limit; ++i) H[i] = f[i] * g[i];
NTT(H, limit, false);
if(n + m < limit) fill(H + n + m + 1, H + limit, 0);
}
}
using poly::modint;
int t, a;
ll n, num, s, sum[lim];
vector<int> visited;
modint f[lim], g[lim];
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> t;
while(t--) {
cin >> n;
visited.clear();
for(int i = 1; i <= n; ++i) {
cin >> a;
visited.push_back(a);
f[a] += 1;
}
visited.erase((stable_sort(visited.begin(), visited.end()), unique(visited.begin(), visited.end())), visited.end());
poly::mul(f, 100000, f, 100000, g);
for(int i = 1; i <= 100000; ++i) g[i << 1] -= f[i];
for(int i = 1; i <= 200000; ++i) {
sum[i] = sum[i - 1] + (g[i].x >> 1);
}
num = s = 0;
for(const auto& i : visited) {
num += f[i].x * (s * (s - 1) / 2 - sum[i]) + f[i].x * (f[i].x - 1) / 2 * s + f[i].x * (f[i].x - 1) * (f[i].x - 2) / 6;
s += f[i].x;
}
for(const auto& i : visited) f[i] = 0;
cout << fixed << setprecision(7) << num / (double)(n * (n - 1) * (n - 2) / 6) << '\n';
}
return 0;
}
分治 FFT/NTT
我觉得我讲不明白，请移步至洛谷题解区。
草草地讲一下，直接套一个 CDQ 的壳子，计算左半部分对右半部分的贡献时，让左半部分的 \(F(x)\) 和右半部分的 \(G(x)\) 卷一下就好。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 998244353, G = 3, gi = 332748118;
int n, limit, r[1 << 19];
ll f[1 << 19], g[1 << 19], power[2][23];
ll ksm(ll x, ll y) {
ll ret = 1;
while(y) {
if(y & 1) ret = ret * x % mod;
x = x * x % mod;
y >>= 1;
}
return ret;
}
void get_rev(int len) {
int l = __lg(len) + 1;
limit = 1 << l;
for(int i = 1; i < limit; ++i) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}
power[0][l] = ksm(gi, (mod - 1) / limit);
power[1][l] = ksm(G, (mod - 1) / limit);
for(int i = l; i >= 1; --i) {
power[0][i - 1] = power[0][i] * power[0][i] % mod;
power[1][i - 1] = power[1][i] * power[1][i] % mod;
}
}
void NTT(ll *h, int n, bool type) {
for(int i = 0; i < n; ++i) {
if(i < r[i]) {
swap(h[i], h[r[i]]);
}
}
ll w, wn, h1, h2;
for(int mid = 1, lg = 0; mid < n; mid <<= 1, ++lg) {
wn = power[type][lg + 1];
for(int i = 0; i < n; i += (mid << 1)) {
w = 1;
for(int j = 0; j < mid; ++j, w = w * wn % mod) {
h1 = h[i + j], h2 = w * h[i + j + mid] % mod;
h[i + j] = (h1 + h2) % mod;
h[i + j + mid] = (h1 - h2 + mod) % mod;
}
}
}
if(type) return;
ll inv = ksm(n, mod - 2);
for(int i = 0; i < n; ++i) {
h[i] = h[i] * inv % mod;
}
}
void CDQ(int l, int r) {
static ll F[1 << 19], G[1 << 19];
if(l == r) {
if(!l) f[l] = 1;
return;
}
const int mid = (l + r) >> 1;
CDQ(l, mid);
copy(f + l, f + mid + 1, F);
copy(g, g + r - l + 1, G);
get_rev(mid - l + r - l);
fill(F + mid - l + 1, F + limit, 0);
fill(G + r - l + 1, G + limit, 0);
NTT(F, limit, true);
NTT(G, limit, true);
for(int i = 0; i < limit; ++i) {
F[i] = F[i] * G[i] % mod;
}
NTT(F, limit, false);
for(int i = mid + 1; i <= r; ++i) {
f[i] = (f[i] + F[i - l]) % mod;
}
CDQ(mid + 1, r);
}
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n; --n;
for(int i = 1; i <= n; ++i) {
cin >> g[i];
}
CDQ(0, n);
for(int i = 0; i <= n; ++i) {
cout << f[i] << " ";
}
return 0;
}
差分与前缀和
这是真的要用 OGF 的知识了啊喂！
简单提一嘴，\(\dfrac{1}{1 - x}\) 被称为前缀和算子，一个 OGF 乘上它可以达到做前缀和的效果，而 \(1 - x\) 被称为差分算子，一个 OGF 乘上它可以达到做差分的效果。
学了 GF 就会明白了，真的。
要做 \(k\) 次前缀和/差分的话用一个多项式快速幂就好。可以先 \(\ln\) 一下，乘 \(k\) 后在 \(\exp\) 回来，用牛顿迭代的话就是单 \(\log\)（不过牛顿迭代的常数好像大得比双 \(\log\) 还慢？），或者就套一个普通快速幂上去也行。
当然你用牛顿二项式展开一样可行，但是我不会，并且复杂度也差不到哪儿去，就直接无脑快速幂就是了。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1004535809, G = 3, gi = 334845270;
int limit, r[1 << 21];
ll power[2][23];
ll ksm(ll x, ll y) {
ll ret = 1;
while(y) {
if(y & 1) ret = ret * x % mod;
x = x * x % mod;
y >>= 1;
}
return ret;
}
void get_rev(int len) {
int l = __lg(len) + 1;
limit = 1 << l;
for(int i = 1; i < limit; ++i) {
r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}
power[0][l] = ksm(gi, (mod - 1) / limit);
power[1][l] = ksm(G, (mod - 1) / limit);
for(int i = l; i >= 1; --i) {
power[0][i - 1] = power[0][i] * power[0][i] % mod;
power[1][i - 1] = power[1][i] * power[1][i] % mod;
}
}
void polyderivate(ll *h, int n, ll *f) {
for(int i = 1; i <= n; ++i) {
f[i - 1] = i * h[i] % mod;
}
f[n] = 0;
}
void polyintegrate(ll *h, int n, ll *f) {
for(int i = n; i >= 0; --i) {
f[i + 1] = h[i] * ksm(i + 1, mod - 2) % mod;
}
f[0] = 0;
}
void NTT(ll *h, int n, bool type) {
for(int i = 0; i < n; ++i) {
if(i < r[i]) {
swap(h[i], h[r[i]]);
}
}
ll w, wn, h1, h2;
for(int mid = 1, lg = 0; mid < n; mid <<= 1, ++lg) {
wn = power[type][lg + 1];
for(int i = 0; i < n; i += (mid << 1)) {
w = 1;
for(int j = 0; j < mid; ++j, w = w * wn % mod) {
h1 = h[i + j], h2 = w * h[i + j + mid] % mod;
h[i + j] = (h1 + h2) % mod;
h[i + j + mid] = (h1 - h2 + mod) % mod;
}
}
}
if(type) return;
ll inv = ksm(n, mod - 2);
for(int i = 0; i < n; ++i) {
h[i] = h[i] * inv % mod;
}
}
void polyinv(ll *h, int n, ll *f) {
static ll g[1 << 21];
fill(f, f + n + n, 0);
f[0] = ksm(h[0], mod - 2);
for(int len = 2, limit; len <= (n << 1); len <<= 1) {
limit = len << 1;
copy(h, h + len, g);
if(len > n + 1) fill(g + n + 1, g + len, 0);
fill(g + len, g + limit, 0);
get_rev(len);
NTT(f, limit, true);
NTT(g, limit, true);
for(int i = 0; i < limit; ++i) {
f[i] = f[i] * (2ll - f[i] * g[i] % mod + mod) % mod;
}
NTT(f, limit, false);
fill(f + len, f + limit, 0);
}
fill(f + n + 1, f + (1 << 21), 0);
}
void polyln(ll *h, int n, ll *f) {
static ll g[1 << 21];
polyinv(h, n, g);
polyderivate(h, n, f);
get_rev(n + n - 1);
NTT(g, limit, true);
NTT(f, limit, true);
for(int i = 0; i < limit; ++i) {
f[i] = f[i] * g[i] % mod;
}
NTT(f, limit, false);
polyintegrate(f, n - 1, f);
}
void CDQ(ll *f, ll *g, int l, int r) {
static ll F[1 << 21], G[1 << 21];
if(l + 1 == r) {
if(!l) f[l] = 1;
else f[l] = f[l] * ksm(l, mod - 2) % mod;
return;
}
const int mid = (l + r) >> 1;
CDQ(f, g, l, mid);
copy(f + l, f + mid, F);
copy(g, g + r - l, G);
get_rev(r - l);
fill(F + mid - l, F + limit, 0);
fill(G + r - l, G + limit, 0);
NTT(F, limit, true);
NTT(G, limit, true);
for(int i = 0; i < limit; ++i) {
F[i] = F[i] * G[i] % mod;
}
NTT(F, limit, false);
for(int i = mid; i < r; ++i) {
f[i] = (f[i] + F[i - l - 1] % mod) % mod;
}
CDQ(f, g, mid, r);
}
void polyexp(ll *h, int n, ll *f) {
static ll g[1 << 21];
polyderivate(h, n, g);
CDQ(f, g, 0, n + 1);
}
void polypower(ll *h, int n, ll x, ll *f) {
static ll g[1 << 21];
polyln(h, n, g);
for(int i = 0; i <= n; ++i) g[i] = g[i] * x % mod;
polyexp(g, n, f);
}
string str;
ll n, k, t, f[1 << 21], g[1 << 21], h[1 << 21];
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n >> str >> t;
for(const auto& i : str) k = (k * 10 + i - '0') % mod;
--n;
for(int i = 0; i <= n; ++i) cin >> f[i];
if(!n) {
cout << f[0];
return 0;
}
g[0] = 1, g[1] = mod - 1;
polypower(g, n, k, h);
if(!t) polyinv(h, n, g);
else copy(h, h + n + 1, g);
get_rev(n + n);
NTT(f, limit, true);
NTT(g, limit, true);
for(int i = 0; i < limit; ++i) h[i] = f[i] * g[i] % mod;
NTT(h, limit, false);
for(int i = 0; i <= n; ++i) cout << h[i] << " ";
return 0;
}
多项式乘法
板题，不讲。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
namespace poly {
constexpr int mod = 998244353, G = 3, gi = 332748118, img = 86583718, inv2 = 499122177, inv2i = 954952494;
struct modint {
int x;
modint(ll v = 0): x((v % mod + mod) % mod) {}
friend modint ksm(modint x, ll y) {
modint ret = 1;
while(y) {if(y & 1) ret *= x; x *= x, y >>= 1;}
return ret;
}
friend modint operator + (const modint& x, const modint& y) {return x.x + y.x >= mod ? x.x + y.x - mod : x.x + y.x;}
friend modint operator - (const modint& x, const modint& y) {return x.x < y.x ? x.x - y.x + mod : x.x - y.x;}
friend modint operator * (const modint& x, const modint& y) {return (ll)x.x * y.x % mod;}
friend modint operator / (const modint& x, const modint& y) {return x * ksm(y, mod - 2);}
friend modint operator - (const modint& x) {return mod - x.x;}
modint operator = (const modint& _) {x = _.x; return *this;}
modint operator += (const modint& _) {*this = *this + _; return *this;}
modint operator -= (const modint& _) {*this = *this - _; return *this;}
modint operator *= (const modint& _) {*this = *this * _; return *this;}
modint operator /= (const modint& _) {*this = *this / _; return *this;}
bool operator == (const modint& _) const {return x == _.x;}
bool operator != (const modint& _) const {return x != _.x;}
friend istream& operator >> (istream& in, modint& _) {static ll tmp; in >> tmp; _.x = tmp % mod; return in;}
friend ostream& operator << (ostream& out, const modint& _) {return out << _.x;}
};
modint ksm(ll x, ll y) {return ksm(modint(x), y);}
int r[1 << 22];
modint power[2][23];
void init(int len) {
int l = __lg(len) + 1;
int limit = 1 << l;
for(int i = 1; i < limit; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
power[0][l] = ksm(gi, (mod - 1) / limit);
power[1][l] = ksm(G, (mod - 1) / limit);
for(int i = l; i >= 1; --i) {
power[0][i - 1] = power[0][i] * power[0][i];
power[1][i - 1] = power[1][i] * power[1][i];
}
}
void NTT(modint* f, int n, bool type) {
for(int i = 1; i < n; ++i) {
if(i < r[i]) swap(f[i], f[r[i]]);
}
modint w, wn, h1, h2;
for(int mid = 1, lg = 0; mid < n; mid <<= 1, ++lg) {
wn = power[type][lg + 1];
for(int i = 0; i < n; i += mid << 1) {
w = 1;
for(int j = 0; j < mid; ++j, w *= wn) {
h1 = f[i + j], h2 = w * f[i + mid + j];
f[i + j] = h1 + h2, f[i + mid + j] = h1 - h2;
}
}
}
if(type) return;
modint inv = ksm(n, mod - 2);
for(int i = 0; i < n; ++i) f[i] *= inv;
}
void mul(modint* F, int n, modint *G, int m, modint *H) {
static modint f[1 << 22], g[1 << 22];
int l = __lg(n + m) + 1;
int limit = 1 << l;
copy(F, F + limit, f);
copy(G, G + limit, g);
init(n + m);
NTT(f, limit, true);
NTT(g, limit, true);
for(int i = 0; i < limit; ++i) H[i] = f[i] * g[i];
NTT(H, limit, false);
}
}
using poly::modint;
int n, m;
modint f[1 << 22], g[1 << 22], h[1 << 22];
int main() {
ios::sync_with_stdio(0);
cin.tie(0), cout.tie(0);
cin >> n >> m;
for(int i = 0; i <= n; ++i) cin >> f[i];
for(int i = 0; i <= m; ++i) cin >> g[i];
poly::mul(f, n, g, m, h);
for(int i = 0; i <= n + m; ++i) cout << h[i] << " ";
return 0;
}
咕咕咕捏，大概写完了，有什么不懂可以问我，让我直接上台讲也不是不行。
本文来自博客园，作者：A_box_of_yogurt，转载请注明原文链接：https://www.cnblogs.com/A-box-of-yogurt/p/18083310
【推荐】国内首个AI IDE，深度理解中文开发场景，立即下载体验Trae
【推荐】编程新体验，更懂你的AI，立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包，你的智能百科全书，全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell：AI 加持，快人一步
· 地球OL攻略 —— 某应届生求职总结
· 周边上新：园子的第一款马克杯温暖上架
· Open-Sora 2.0 重磅开源！
· 提示词工程——AI应用必不可少的技术
· .NET周刊【3月第1期 2025-03-02】