Solution -「LOJ #6895」Yet Another NPC Problem
\(\mathscr{Description}\)
Link.
给定 \(l,m\),求当 \(k=m,m+1,\dots,m+l-1\) 时,所有 \(k\) 阶有标号简单无向图中,最大团大小为 \(m\) 的图的数量奇偶性。
\(l\le10^6\),\(m\le2^{10^6}\)。
\(\mathscr{Solution}\)
这道题的解题过程含有密集的神谕(oracle),希望大家能承受来自 祂 的神力。而且体悟神谕之前尽量自己想一想。(
首先自然是固定 \(n=m+k\),想一想如何求此时 \((n,m)\) 的答案。
- Motivation 奇偶性 \(\rightarrow\) 两两配对抵消。
不妨用 \(1,2,\dots,n\) 直接表示无向图 \(G\) 中的结点,\(S_u\) 表示结点 \(u\) 的邻接点集,\(M(G)\) 表示 \(G\) 的最大团大小。我们尝试在 \(G\) 上构造一些性质,若在这些性质下,\(G\) 的数量一定是偶数,我们就可以直接忽略具有这一性质的所有 \(G\)。
- Oracle 任取两结点,不妨取 \(1,2\),则 \(S_1\setminus\{2\}\neq S_2\setminus\{1\}\) 的 \(G\) 有偶数个。
Proof
这个 oracle 比较平凡。对于满足这一性质的 $G$,交换 $1,2$ 的邻接点集得到 $G'$,显然 $M(G)=M(G')$ 且 $G\neq G',G=G''$,所以这样的 $G$ 有偶数个。 $\square$注意到这个 oracle 带来的限制非常强,它使得 \(1,2\) 在参与团构成的过程中“几乎没有区别”。进一步,很自然地想到继续讨论 \((1,2)\in E\) 的情况。
- 若 \((1,2)\in E\),若 \(1\) 在某个团中,\(2\) 必然可以加入这个团,反之亦然。
- 若 \((1,2)\notin E\),若 \(1\) 在某个团中,\(2\) 必然可以将 \(1\) 替换,反之亦然。
我们把这一讨论抽象到图上,引入 缩点 操作,将诸如此类的点集用一个点 \(u\) 表示。当然,若 \(u\) 被囊括入某个团,\(u\) 对团大小的变化就不一定是 \(1\) 了。我们进一步定义 \(u\) 的大小 \(s(u)\) 来刻画 \(u\) 对团大小的贡献。延续上文的讨论,设两个点(原图的点或者缩过的点)\(1,2\) 将要缩成点 \(u\),那么:
- 若 \((1,2)\in E\),则 \(u\) 成团时,必然可以同时取 \(1,2\)。因此 \(s(u)=s(1)+s(2)\)。
- 若 \((1,2)\notin E\),则 \(u\) 成团时,一定选 \(1,2\) 中较大者,让团更大。因此 \(s(u)=\max\{s(1),s(2)\}\)。
事实上,两个点缩或者不缩是我们说了算,现在还应当规定缩点的规则。
- Motivation 避免分类讨论。
为了避免 \(\max\) 的比较过程,我们规定,当且仅当 \(s(1)=s(2)\),我们才对 \(1,2\) 缩点。此时有 \(s(u)=s(1)+s(2)\) 或者 \(s(u)=s(1)=s(2)\)。缩点的顺序?规定每次取标号最小且等大的点对 \(u<v\) 进行缩点,缩点得到新点的标号为 \(u\),这显然是良定义的。
当图 \(G\) 无法进行缩点时,必然有 \(\forall u<v,s(u)>s(v)\),且 \(s(u)=2^k,k\in\mathbb N\),\(\sum_u s(u)\le n\)。因此,一个整数 \(b\in(0,n]\) 的二进制表示足以描述图 \(G\) 的 \(\{s(u)\}_{u\in V}\)。此后,我们还需要确定这些无法再缩的点之间的连边。
形式地,定义 \(f(n,b)\) 表示把初始有 \(n\) 个结点的图缩点最终得到状态 \(b\) 的方案数奇偶性,\(g(b,m)\) 表示在状态 \(b\) 对应的点集内连边,使得最大团为 \(m\) 的方案数奇偶性。那么答案显然是
先研究 \(g\),神谕:
- Oracle \(g(b,m)=[m=b\lor m=b-\operatorname{low}(b)]\).
Proof
- 对于 \(g(b,b)\),有且仅有连完全图一种方案,所以方案数为奇数。
- 对于 \(g(b,b-\operatorname{low}(b))~(b\neq2^k,k\in\mathbb N)\),\(\operatorname{low}(b)\) 所对应的点有 \((2^{\operatorname{count}(b)-1}-1)\equiv 1\pmod 2\) 种向前连边的方案,其余点内部连完全图。
- 对于其他 \(g(b,m)\),令 \(x\) 为 \(\operatorname{low}(b)\) 对应的点,\(y\) 为标号最小(\(s\) 最大)的不在最大团中的点 \(y\),显然,\((x,y)\) 的存在性不可能不影响最大团。因此这样的图的方案数一定为偶数。
综上。 \(\square\)
可见 \(g\) 取 \(1\) 的位置非常稀疏。
再研究 \(f\),注意到缩点的过程很像二进制加法的进位过程,所以尝试用增量法 DP 求 \(f\)。具体地,考虑在已有 \(n-1\) 个点的不可缩的图上加入一个标号为 \(n\),大小为 \(1\) 的点,讨论:
- 新点参与的最后一次缩点是 "\(s(u)=s(1)+s(2)\)",那显然 \(f(n,b)\gets f(n,b)+f(n-1,b-1)\)。
- 新点参与的最后一次缩点是 "\(s(u)=s(1)=s(2)\)",在进位的角度,\(b\) 本来有一堆低位 bit,这个 \(1\) 把它们一路进位上去,但最后却没有加入 \(b\),全部扔掉了。扔了多少?一定是 \(\operatorname{low}(b)\)。因此,\(f(n,b)\gets f(n,b)+f(n-1,b+\operatorname{low}(b)-1)\)。
综上,
基于此,我们已有了 \(\mathcal O(m(m+l))\) 的算法。
祂 说:"GF."
令 \(F_b(x)=\sum_nf(n,b)x^n\),那么
这是基于上文递推的形式。另外,根据组合意义,设 \(b=\sum_{i}2^ib_i\),那么
因此,我们只需要研究 \(\mathcal F_k(x)=F_{2^k}(x)\) 的性质。
套用递推式并变形:
整理一下,
利用前缀积与 \(\mathcal F_{k-1}(x)\) 的相互转化,得到 \(\mathcal F\) 的递推式:
根据 祂 传授的经验,“这种东西上下齐次,因此分子分母分别设出来,把分式递推转化成整式递推”,因此,设 \(\mathcal F_k(x)=R_k(x)/Q_k(x)\),递推:
因此令
设答案序列为 \(\mathcal R=\lang r_i\rang_{i=0}^{l-1}\),\(m=\sum_{i=1}^h2^{m_i}\),其中 \(m_1>m_2>\cdots>m_h\),利用上文的结论转化:
模 \(2\) 意义下多项式求逆和多项式卷积,可以做到 \(\mathcal O(l\log l\log m)\)。
模 \(2\) 多项式卷积,直接 FFT 加持算出卷积结果每项模 \(2\) 即可。模 \(2\) 意义下多项式求逆,考虑倍增求逆的过程,本质上求逆只是若干次卷积的封装,每次卷积结果拿出来模 \(2\) 即可。注意这警示我们必须严格遵从“两个多项式 DFT 后点值相乘,乘完立马 IDFT 回去最后模 \(2\)”,不可以拿着 DFT 的结果进行连续的点值乘法。
回顾一下 \(Q\) 的式子:
发现一个有趣的事实,在模 \(2\) 意义下,任意多项式 \(G(x)\) 有 \(G^2(x)=G(x^2)\),因而
观察到,对于 \(k>0\),有
于是,答案式中的求和与求积的上界都能收到 \(\log l\) 的进一步约束,继而得到 \(\mathcal O(l\log^2l)\) 的算法。
精彩的来力!
Oracle \(Q_k(x)\) 是稀疏多项式。
这是 祂 在题解里留下的猜想,并未给出证明。于是验题的时候,兔子老老实实写 NTT 结果被 祂 的暴力杀穿了。兔子感叹于自己的渺小,但为了膜拜 祂,她抄起 __uint128_t
打 \(Q_k(x)\) 系数数量的表,经过长时间以神遇而不以目视的观察后,尝试续写神的旨意……
Oracle 令 \(f(k)\) 表示 \(Q_k(x)\) 模 \(2\) 意义下系数 \(1\) 的数量,那么
(可以想象到观察过程有多么奇妙。) 当然这个描述得太傻了,稍微化简一下:
据此能够说明,\(Q_k(x)\) 的系数 \(1\) 的数量是 \(\mathcal O(k)=\mathcal O(\log l)\) 级别!
她震撼之余再次感叹于自己的渺小,于是拿着结论求助 祂,最后得到了证明。
Proof
观察 $Q_k(x)$ 系数 $1$ 的指标位置,注意到除 $x^{2^k}$ 这一项,其他的 $x^t$ 中,$t$ 都有 $(11\dots100\dots0)_2$ 的形式。令 $T_k(x)=Q_k(x)-x^{2^k}$,带入递推式得到 $T_k(x)=T_k^2(x)+x^{2^{k-1}}T_{k-1}(x)$,注意刚刚提到的 trick,进一步得到 $T_k(x)=T_{k-1}(x^2)+x^{2^{k-1}}T_{k-1}(x)$。组合意义!有一个二进制数,初始为空,操作 $k$ 次,每次在数的左边写一个 $1$ 或者在数的右边写一个 $0$,$x^t$ 就是得到二进制数 $t$ 的方案数奇偶性!而显然,枚举 $t$ 后缀 $0$ 的个数 $c$,这一方案数可由组合数刻画: $$ Q_k(x)=x^{2^k}+\sum_{c=0}^k\binom{k}{c}x^{2^k-2^c}. $$ $\binom{k}{c}\equiv1\pmod2$,Lucas 啊!所以 $c$ 的数量就是 $2^{\operatorname{count}(k)}$。加上 $x^{2^k}$,得到 $f(k)$。 $\square$利用结论,实现时,暴力存 \(Q_k(x)\) 系数 \(1\) 的指标位置,每次除法时保证分母稀疏,对于分母有 \(a\) 个系数 \(1\),分子有 \(b\) 个系数 \(1\) 的除法,暴力 \(\mathcal O(ab)\) 地实现。考察复杂度,共 \(\log l\) 次除法,每次除法假设分子系数数量为 \(\mathcal O(l)\),设 \(2^{2^{k-1}}\le l<2^{2^k}\),那么
在诸多神谕的加持下,我们最终得到了 \(\mathcal O(l\log^{\log_23} l)\) 的算法。
\(\mathscr{Code}\)
- 多项式 ver.(注意这里将 \(m\) 直接作为
long long
输入了。)
/*+Rainybunny+*/
#include <bits/stdc++.h>
#define rep(i, l, r) for (int i = l, rep##i = r; i <= rep##i; ++i)
#define per(i, r, l) for (int i = r, per##i = l; i >= per##i; --i)
typedef long long LL;
const int MAXL = 1 << 18;
namespace PolyOper {
const int MOD = 998244353, MG = 3;
int omega[19][MAXL];
inline int mul(const int u, const int v) { return 1ll * u * v % MOD; }
inline void subeq(int& u, const int v) { (u -= v) < 0 && (u += MOD); }
inline int sub(int u, const int v) { return (u -= v) < 0 ? u + MOD : u; }
inline void addeq(int& u, const int v) { (u += v) >= MOD && (u -= MOD); }
inline int add(int u, const int v) { return (u += v) < MOD ? u : u - MOD; }
inline int mpow(int u, int v) {
int ret = 1;
for (; v; u = mul(u, u), v >>= 1) ret = mul(ret, v & 1 ? u : 1);
return ret;
}
inline void init(const int h) {
rep (i, 1, h) {
int* wi = omega[i]; wi[0] = 1, wi[1] = mpow(MG, MOD - 1 >> i);
rep (j, 2, (1 << i) - 1) wi[j] = mul(wi[j - 1], wi[1]);
}
}
inline void ntt(const int len, int* u, const int type) {
static int rev[MAXL], las = -1;
if (las != len) {
las = len;
rep (i, 1, len - 1) rev[i] = rev[i >> 1] >> 1 | (i & 1) * (len >> 1);
}
rep (i, 1, len - 1) if (i < rev[i]) u[i] ^= u[rev[i]] ^= u[i] ^= u[rev[i]];
for (int i = 1, stp = 1; stp < len; ++i, stp <<= 1) {
const int* wi = omega[i];
for (int j = 0; j < len; j += stp << 1) {
rep (k, j, j + stp - 1) {
int x = u[k], y = mul(wi[k - j], u[k + stp]);
u[k] = add(x, y), u[k + stp] = sub(x, y);
}
}
}
if (!~type) {
int inv = MOD - (MOD - 1) / len;
rep (i, 0, len - 1) u[i] = mul(inv, u[i]);
std::reverse(u + 1, u + len);
}
}
inline void polyMul(const int len, const int* u, const int* v, int* r) {
static const int BF_THRES = 1 << 5;
static int tmp[2][MAXL];
if (len <= BF_THRES) {
rep (i, 0, len - 1) tmp[0][i] = 0;
rep (i, 0, len - 1) rep (j, 0, repi - i) tmp[0][i + j] ^= u[i] & v[j];
rep (i, 0, len - 1) r[i] = tmp[0][i];
} else {
rep (i, 0, len - 1) tmp[0][i] = u[i], tmp[1][i] = v[i];
rep (i, len, (len << 1) - 1) tmp[0][i] = tmp[1][i] = 0;
ntt(len << 1, tmp[0], 1), ntt(len << 1, tmp[1], 1);
rep (i, 0, (len << 1) - 1) tmp[0][i] = mul(tmp[0][i], tmp[1][i]);
ntt(len << 1, tmp[0], -1);
rep (i, 0, len - 1) r[i] = tmp[0][i] & 1;
}
}
inline void polyInv(const int len, const int* u, int* r) {
if (len == 1) return assert(r[0] = u[0]), void();
polyInv(len >> 1, u, r);
static int tmp[2][MAXL];
rep (i, 0, (len >> 1) - 1) tmp[0][i] = r[i], tmp[1][i] = u[i];
rep (i, len >> 1, len - 1) tmp[0][i] = 0, tmp[1][i] = u[i];
rep (i, len, (len << 1) - 1) tmp[0][i] = tmp[1][i] = 0;
ntt(len << 1, tmp[0], 1), ntt(len << 1, tmp[1], 1);
rep (i, 0, (len << 1) - 1) tmp[0][i] = mul(tmp[0][i], tmp[0][i]);
ntt(len << 1, tmp[0], -1);
rep (i, 0, (len << 1) - 1) tmp[0][i] &= 1;
ntt(len << 1, tmp[0], 1);
rep (i, 0, (len << 1) - 1) tmp[0][i] = mul(tmp[0][i], tmp[1][i]);
ntt(len << 1, tmp[0], -1);
rep (i, 0, len - 1) r[i] = tmp[0][i] & 1;
}
} // namespace PolyOper.
const int MAXLG = 18;
int L, Q[MAXLG + 5][MAXL], R[MAXL], T[MAXL], ans[MAXL];
LL M;
int main() {
scanf("%d %lld", &L, &M);
PolyOper::init(33 - __builtin_clz(L));
int len = 1; while (len < L) len <<= 1;
Q[0][0] = R[0] = T[0] = 1, Q[0][1] = T[1] = 1;
rep (i, 0, 31 - __builtin_clz(int(std::min((L << 1) - 1ll, M)))) {
if (M >> i & 1) PolyOper::polyMul(len, Q[i], R, R);
if (i < repi) {
rep (j, 0, len - 1) Q[i + 1][j] = T[j];
if (1ll << i << 1 < len) Q[i + 1][1 << i << 1] ^= 1;
PolyOper::polyMul(len, Q[i + 1], T, T);
}
}
PolyOper::polyInv(len, R, T);
rep (i, 0, L - 1) ans[i] ^= (R[i] = T[i]);
rep (i, 0, 31 - __builtin_clz(int(std::min(L - 1ll, M)))) {
if (M >> i & 1) break;
PolyOper::polyInv(len, Q[i], T);
PolyOper::polyMul(len, R, T, T);
rep (j, 1 << i, L - 1) ans[j] ^= T[j - (1 << i)];
}
rep (i, 0, L - 1) putchar(ans[i] ^ '0');
putchar('\n');
return 0;
}
- 正解 ver.
/*+Rainybunny+*/
#include <bits/stdc++.h>
#define rep(i, l, r) for (int i = l, rep##i = r; i <= rep##i; ++i)
#define per(i, r, l) for (int i = r, per##i = l; i >= per##i; --i)
typedef long long LL;
typedef std::vector<int> Poly;
const int MAXLG = 23, MAXL = 1e7, MAXML = 1e6;
int L, len;
Poly Q[MAXLG + 5];
char M[MAXML + 5];
inline void adapt(Poly& u) {
Poly res;
std::sort(u.begin(), u.end());
for (int s = int(u.size()), l = 0, r; l < s; l = r) {
for (r = l; r < s && u[r] == u[l]; ++r);
if ((r - l) & 1) res.push_back(u[l]);
}
u.swap(res);
}
inline Poly operator * (const Poly& u, const Poly& v) {
Poly ret;
for (int x: u) for (int y: v) {
if (x + y >= L) break;
ret.push_back(x + y);
}
return adapt(ret), ret;
}
inline Poly operator + (const Poly& u, const Poly& v) {
Poly ret;
for (size_t i = 0, j = 0; i < u.size() || j < v.size();) {
if (i < u.size() && j < v.size() && u[i] == v[j]) ++i, ++j;
else if (i < u.size() && (j == v.size() || u[i] < v[j])) {
ret.push_back(u[i++]);
} else {
ret.push_back(v[j++]);
}
}
return ret;
}
/// Given that the polynomial v is very sparse.
inline Poly operator / (const Poly& u, const Poly& v) {
assert(!v[0]); // [x^0]v=1.
static bool res[MAXL + 5];
Poly ret;
for (int i = 0, j = 0; i < L; ++i) {
while (j < u.size() && u[j] < i) ++j;
bool tar = j < u.size() && i == u[j];
for (int x: v) if (x) {
if (x > i) break;
tar ^= res[i - x];
}
if ((res[i] = tar)) ret.push_back(i);
}
return ret;
}
int main() {
scanf("%d %s", &L, M), len = strlen(M);
std::reverse(M, M + len);
int h = 31 - __builtin_clz((L << 1) - 1);
Poly T({ 0, 1 }); Q[0] = T;
rep (i, 0, h - 1) {
Q[i + 1] = T;
if (1 << i << 1 < L) Q[i + 1].push_back(1 << i << 1), adapt(Q[i + 1]);
T = T * Q[i + 1];
}
Poly R({ 0 });
rep (i, 0, h) {
if (M[i] ^ '0') break;
R = R + Poly({ 1 << i }) / Q[i];
}
rep (i, 0, h) if (i < len && M[i] ^ '0') R = R / Q[i];
for (int i = 0, j = 0; i < L; ++i) {
while (j < R.size() && R[j] < i) ++j;
putchar('0' ^ (j < R.size() && R[j] == i));
}
putchar('\n');
return 0;
}