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$ 的方案数奇偶性。那么答案显然是

$$\sum_bf(n,b)\cdot g(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)$。

  综上，

$$f(n,b)=f(n-1,b-1)+f(n-1,b+\operatorname{low}(b)-1).$$

基于此，我们已有了 $\mathcal O(m(m+l))$ 的算法。

---

  祂 说："GF."

  令 $F_b(x)=\sum_nf(n,b)x^n$，那么

$$F_b(x)=x(F_{b-1}(x)+F_{b+\operatorname{low}(b)-1}(x)).$$

这是基于上文递推的形式。另外，根据组合意义，设 $b=\sum_{i}2^ib_i$，那么

$$F_b(x)=\prod_{b_i=1}F_{2^i}(x).$$

因此，我们只需要研究 $\mathcal F_k(x)=F_{2^k}(x)$ 的性质。

  套用递推式并变形：

$$\begin{aligned} \mathcal F_k(x) &= x(F_{2^k-1}(x)+F_{2^{k+1}-1}(x))\\ &= x\prod_{i=0}^{k-1}\mathcal F_i(x)+x\mathcal F_k(x)\prod_{i=0}^{k-1}\mathcal F_i(x). \end{aligned}$$

整理一下，

$$\mathcal F_k(x)=\frac{x\prod_{i=0}^{k-1}\mathcal F_i(x)}{1-x\prod_{i=0}^{k-1}\mathcal F_i(x)},\\ x\prod_{i=0}^{k-1}\mathcal F_i(x)=\frac{\mathcal F_k(x)}{1+\mathcal F_k(x)}.$$

利用前缀积与 $\mathcal F_{k-1}(x)$ 的相互转化，得到 $\mathcal F$ 的递推式：

$$\mathcal F_0(x)=\frac{x}{1-x},\\ \mathcal F_k(x)=\frac{\mathcal F_{k-1}^2(x)}{1+\mathcal F_{k-1}(x)-\mathcal F_{k-1}^2(x)}.$$

  根据 祂 传授的经验，“这种东西上下齐次，因此分子分母分别设出来，把分式递推转化成整式递推”，因此，设 $\mathcal F_k(x)=R_k(x)/Q_k(x)$，递推：

$$\begin{aligned} \frac{R_k(x)}{Q_k(x)} &= \frac{R_{k-1}^2(x)/Q_{k-1}^2(x)}{1+R_{k-1}(x)/Q_{k-1}(x)-R_{k-1}^2(x)/Q_{k-1}^2(x)}\\ &= \frac{R_{k-1}^2(x)}{Q_{k-1}^2(x)+R_{k-1}(x)Q_{k-1}(x)-R_{k-1}^2(x)}. \end{aligned}$$

因此令

$$R_k(x)=R_{k-1}^2(x)=x^{2^k},\\ Q_k(x)=Q_{k-1}^2(x)+R_{k-1}(x)Q_{k-1}(x)-R_{k-1}^2(x)=Q_{k-1}^2(x)+x^{2^{k-1}}Q_{k-1}(x)-x^{2^k}.$$

  设答案序列为 $\mathcal R=\lang r_i\rang_{i=0}^{l-1}$，$m=\sum_{i=1}^h2^{m_i}$，其中 $m_1>m_2>\cdots>m_h$，利用上文的结论转化：

$$\begin{aligned} \mathcal R &= [x^{m..m+l-1}]\left(F_m(x)+\sum_{i=0}^{m_l-1}F_{m+2^i}(x)\right)\\ &= [x^{m..m+l-1}]\left(1+\sum_{i=0}^{m_l-1}\mathcal F_i(x)\right)\prod_{i=1}^h\mathcal F_{m_i}(x)\\ &= [x^{0..l-1}]\frac{\left(1+\sum_{i=0}^{m_l-1}\frac{x^{2^i}}{Q_i(x)}\right)}{\prod_{i=1}^hQ_{m_i}(x)}. \end{aligned}$$

模 $2$ 意义下多项式求逆和多项式卷积，可以做到 $\mathcal O(l\log l\log m)$。

  模 $2$ 多项式卷积，直接 FFT 加持算出卷积结果每项模 $2$ 即可。模 $2$ 意义下多项式求逆，考虑倍增求逆的过程，本质上求逆只是若干次卷积的封装，每次卷积结果拿出来模 $2$ 即可。注意这警示我们必须严格遵从“两个多项式 DFT 后点值相乘，乘完立马 IDFT 回去最后模 $2$”，不可以拿着 DFT 的结果进行连续的点值乘法。

---

  回顾一下 $Q$ 的式子：

$$Q_k(x)=Q_{k-1}^2(x)+x^{2^{k-1}}Q_{k-1}(x)-x^{2^k}.$$

发现一个有趣的事实，在模 $2$ 意义下，任意多项式 $G(x)$ 有 $G^2(x)=G(x^2)$，因而

$$Q_k(x)=Q_{k-1}(x^2)+x^{2^{k-1}}Q_{k-1}(x)-x^{2^k}.$$

观察到，对于 $k>0$，有

$$Q_k(x)\equiv 1\pmod{x^{2^{k-1}}}.$$

于是，答案式中的求和与求积的上界都能收到 $\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$ 的数量，那么

$$f(0)=2,\\ f(k)=2^{g(k-\operatorname{high}(k)+1)}+1,\\ g(k)=\operatorname{count}(k)+\operatorname{ctz}(k).$$

（可以想象到观察过程有多么奇妙。） 当然这个描述得太傻了，稍微化简一下：

$$f(k)=2^{\operatorname{count}(k)}+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}$，那么

$$\begin{aligned} \sum_{i=0}^{2^k-1}l\cdot f(i) &= l\sum_{i=0}^{2^k-1}2^{\operatorname{count}(i)}+1\\ &= l\left(2^k+\sum_{i=0}^{k-1}\binom{k}{i}2^i\right)\\ &= \mathcal O(l3^k)\\ &= \mathcal O(l3^{\log_2\log_2l})\\ &= \mathcal O(l\log_2^{\log_23} l). \end{aligned}$$

  在诸多神谕的加持下，我们最终得到了 $\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;
}