闲话 22.12.28

闲话

这篇写成闲话，主要是最近想推歌了（
《シャボン(肥皂泡)》 by 蜂屋ななし feat. 初音ミク
《イカサマダンス(欺诈舞蹈)》 by まふまふ feat. 鏡音リン
《Shamer》 by Chinozo feat. flower
就今天给我推的歌（
我记性留不到第二天（

线性规划先等一天（
咕了咕了！

关于 \(\text{Lonesum Matrix}\) 的讨论

又名 W2B。

\(\text{Lonesum Matrix}\) 计数

定义一个 \(n\times m\) 的 \(01\) 矩阵是 \(\text{Lonesum Matrix (LM)}\)，当且仅当其能由行列和向量唯一地确定。行和向量是对每一行分别求得 \(1\) 的个数后按原顺序排列得到的 \(n\) 维向量，列和向量同理。

给定 \(n, m\)，计数 \(n\times m\) 的 \(\text{LM}\) 个数。

我们将给出两种方案，其最终能得到相同的结果。

1. 组合意义

请和我默念：组合意义天地灭

我们考虑使用一些列向量来拼合一个 \(\text{LM}\) 非全 \(0\) 全 \(1\) 的部分。我们称拼合的矩阵为部分矩阵。
定义不合法分量为 \(\left[\ \begin{aligned} & 1 & 0 \\ & 0 & 1 \end{aligned}\ \right]\) 和 \(\left[\ \begin{aligned} & 0 & 1 \\ & 1 & 0 \end{aligned}\ \right]\)。
我们定义不合法的部分矩阵是存在不合法分量为子矩阵的部分矩阵。可以发现存在这样子矩阵的矩阵能够将对应子矩阵交换两行/列得到行列和向量相同的不同矩阵。

我们定义行/列的权值是本行/列中 \(1\) 的个数。容易发现在合法的部分矩阵中不会出现相同权的列向量。需要注意这里没有使用全 \(0\) 全 \(1\) 的列向量。
可以发现对行列进行重排不会改变其是否是 \(\text{LM}\)，因为不合法分量是彼此镜像的。

随后就能计数了。
我们取得一个部分矩阵。假设这个部分矩阵由 \(k - 1\) 个列向量按权值降序组成，则每行的权值在 \([0, k)\) 之间，这就构成了 \(k\) 个等价类。对 \(k\) 个彼此不同的等价类放入 \(n\) 个行计数得到 \(k! {n \brace k}\)。
\(k - 1\) 个列的部分矩阵肯定是 \(k\) 个列的部分矩阵的子集，因此做容斥，系数是 \((-1)^{n + m}\)。
最后还需要确定 \(k - 1\) 个列的部分矩阵所能生成的矩阵数量。可以在最后加入一行全 \(1\)，因此行权值在 \([0, k]\) 之间。有 \((k + 1)^m\) 种可能。

得到计数

\[\sum_{k=1}^n (-1)^{n + k} k! {n \brace k} (k + 1)^m \]

感觉不是那么好理解？看看下面？

2. 代数推导

请和我默念：组合意义天地灭，代数推导保平安。

当然需要声明的是：本方法来自 \(\text{A}\color{red}{\text{PJifengc}}\)。
此外，本文部分内容完全来自于她解题时的想法与评价。

首先需要考虑的是不合法情况，即子矩阵出现不合法分量。然后设 \(f(n, m)\) 为 \(n \times m\text{ LM}\) 的数量。

首先钦定第一行选什么。然后发现顺序不重要，有 \(01\) 个数就行了。枚举 \(1\) 的数量。随后考虑如果第二行有一个数和第一行的数相同，那其他第二行对应第一行位置上数和第一行的数相同的数就必须和第二行这个数相同。这就等价于 \(0\) 或 \(1\) 必须有一个数和第一行完全一样，直接计数，减掉都一样的情况就行了。如果下面的每一行都满足这个条件，那这些行之间就不用再考虑了，直接拆成了子问题。

也就是说：

\[f(n, m) = \sum_{i=0}^{m}\sum_{j=0}^n \sum_{k = 0}^{n-j} (-1)^j \binom{m}{i}\binom{n-1}{j}\binom{n-1-j}{k} f(k, i) f(n - 1 - j - k, m - i) \]

到这一步就是 \(O(n^5)\) 的了，下面我们要把它用代数工业优化成 \(O(n\log n)\)。

其实还挺顺的。
考虑生成这个 \(f\)，看到二项式系数后采用 EGF。设

\[F(x, y) = \sum_{i=0}^{\infty} \sum_{j=0}^{\infty} f(i, j)\frac{x^i}{i!}\frac{y^j}{j!} \]

暴力带入后拆解各项系数可以得到一长串式子。首先这应当是 \(F^2\) 的形式，但你发现系数在 \(x\) 方向有偏移，因此需要对 \(x\) 求积分，发现系数也对了。但是 \(x = 0\) 的地方缺了一项，加入 \(e^y\) 即可。

也就是说

\[F(x, y) = \int F^2(x, y) \text d x + e^y \]

解微分方程可得

\[F(x, y) = (e^{-x} + e^{- y} - 1)^{-1} \]

我们需要的就是提取 \(F\) 的 \([x^ny^m / n!m!]\) 次项。这时候你可能会想 “不会！”、”我要是会做我就会做了！”、“麻了，这个式子我拆拆就拆回最开始写的式子了”、”感觉封闭形式推了个寂寞“……。
但是不要担心，万事皆有可能！我们最需要的其实是良好的心态——

让我们从最简单的手段开始：

\[\begin{aligned} &\left[\frac{x^ny^m}{n!m!}\right] (e^{-x} + e^{- y} - 1)^{-1} \\ = & \ \left[\frac{x^ny^m}{n!m!}\right] \sum_{i=0}^{\infty} \binom{-1}{i} (e^{-x} - 1)^i e^{-y( - 1 - i)} \qquad & (广义二项式定理) \\ = & \ \left[\frac{x^ny^m}{n!m!}\right] \sum_{i=0}^{\infty} (-1)^i (e^{-x} - 1)^i e^{y(i + 1)} \qquad & (上指标反转) \\ = & \ \left[\frac{x^ny^m}{n!m!}\right] \sum_{i=0}^{\infty} (-1)^i (e^{-x} - 1)^i \sum_{j=0}^{\infty} \frac{y^j(i + 1)^j}{j!} \qquad& (e^x\ 的 \text{ EGF } 形式) \\ = & \ \left[\frac{x^ny^m}{n!m!}\right] \sum_{j=0}^{\infty} \left(\sum_{i=0}^{\infty} (-1)^i (e^{-x} - 1)^i (i + 1)^j\right) \frac{y^j}{j!} \\ = & \ \left[\frac{x^n}{n!}\right] \sum_{i=0}^{\infty} (-1)^i (e^{-x} - 1)^i (i + 1)^m \end{aligned}\]

我们已经去掉了一个变量了！只剩一个了！

我们发现，\((e^{-x} - 1)^i\) 这样的形式有些眼熟，令人想起第二类斯特林数的 EGF。考虑配凑出形式。

\[\begin{aligned} &\left[\frac{x^n}{n!}\right] \sum_{i=0}^{\infty} (-1)^i (i + 1)^m (e^{-x} - 1)^i \\ = & \ \left[\frac{x^n}{n!}\right] \sum_{i=0}^{\infty} (-1)^i (i + 1)^m i! \frac{(e^{-x} - 1)^i}{i!} \qquad & (如上方法配凑) \\ = & \ \left[\frac{x^n}{n!}\right] \sum_{i=0}^{\infty} (-1)^i (i + 1)^m i! \sum_{j=i}^{\infty} {j \brace i} \frac{(-x)^j}{j!} \qquad & (展开为 \ \text{EGF}) \\ = & \ \left[\frac{x^n}{n!}\right] \sum_{j=0}^{\infty} \sum_{i=0}^{j} (-1)^i (i + 1)^m i! {j \brace i} \frac{(-x)^j}{j!} \\ = & \ \left[\frac{x^n}{n!}\right] \sum_{j=0}^{\infty} \left( \sum_{i=0}^{j} (-1)^{i+j} i! {j \brace i}(i + 1)^m \right) \frac{x^j}{j!} \\ = & \ \sum_{i=0}^{n} (-1)^{i+n} i! {n \brace i}(i + 1)^m \end{aligned}\]

我们成功了！形式和组合方法得到的完全相同！让我们听听 APJifengc 先生是如何想的吧：

真是太强了！

好到这里两种做法就都告一段落了。可以使用生成函数法求得第二类斯特林数，随后模拟即可。总时间复杂度为 \(O(n\log n)\)。

最后需要说明的是，\(f(n, m)\) 即为 \(B_n^{(-m)}\)，多项式伯努利数，A099594。

code

#include <bits/stdc++.h>
using namespace std; using pii = pair<int,int>; using vi = vector<int>; using vp = vector<pii>; using ll = long long; 
using ull = unsigned long long; using db = double; using ld = long double; using lll = __int128_t;
mt19937 rnd(chrono::steady_clock::now().time_since_epoch().count());
template <typename T> T rand(T l, T r) { return uniform_int_distribution<T>(l, r)(rnd); }
template <typename T> void get(T & x) {
	x = 0; char ch = getchar(); bool f = false; while (ch < '0' or ch > '9') f = f or ch == '-', ch = getchar();
	while ('0' <= ch and ch <= '9') x = (x << 1) + (x << 3) + ch - '0', ch = getchar(); f && (x = -x); 
} template <typename T, typename ... Args> void get(T & a, Args & ... b) { get(a); get(b...); }
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
#define timer cerr << 1. * clock() / CLOCKS_PER_SEC << '\n';
template <typename T1, typename T2> T1 min(T1 a, T2 b) { return a < b ? a : b; }
template <typename T1, typename T2> T1 max(T1 a, T2 b) { return a > b ? a : b; }
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)
#define ufile(x) 
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
#define Aster(i, s) for (int i = head[s], v; i; i = e[i].next)
#define all(s) s.begin(), s.end()
#define eb emplace_back
#define pb pop_back
#define em emplace
const int N = 3e6 + 10, mod = 998244353, g = 3;
using poly = vector<int>;
int n, m, ans;
poly F, G; 

int qp(int a, int b) {
    int ret = 1;
    while (b) {
        if (b & 1) ret = 1ll * a * ret % mod;
        a = 1ll * a * a % mod;
        b >>= 1;
    } return ret;
} int invg = qp(g, mod - 2);

int btrs[N];
int initrs(int mx) {
    int lim = 1;
    while (lim <= mx) lim <<= 1;
    for (int i = 0; i < lim; ++ i)
        btrs[i] = (btrs[i >> 1] >> 1) | ((i & 1) ? lim >> 1 : 0);
    return lim;
}

void NTT(poly& a, int len, const int typ) {
    a.resize(len);
    for (int i = 0; i < len; ++ i) if (i < btrs[i]) swap(a[i], a[btrs[i]]);
    for (int mid = 1, t, wn; mid < len; mid <<= 1) {
        wn = qp(typ == 1 ? g : invg, (mod - 1) / (mid << 1));
        for (int j = 0, w = 1; j < len; j += (mid << 1)) {
            w = 1;
            for (int k = 0; k < mid; ++ k, w = 1ll * w * wn % mod) {
                t = 1ll * w * a[j + k + mid] % mod;
                a[j + k + mid] = a[j + k] - t + mod; if (a[j + k + mid] >= mod) a[j + k + mid] -= mod;
                a[j + k] += t; if (a[j + k] >= mod) a[j + k] -= mod;
            }
        }
    } if (typ == 1) return;
    int inv = qp(len, mod - 2);
    for (int i = 0; i < len; ++ i) a[i] = 1ll * a[i] * inv % mod;
}

poly A, B;
poly mul(const poly &a, const poly &b) {
    int o = a.size() + b.size() - 1, len = initrs(a.size() + b.size() - 1);
    A = a, B = b;
    NTT(A, len, 1), NTT(B, len, 1);
    for (int i = 0; i < len; ++ i) A[i] = 1ll * A[i] * B[i] % mod;
    NTT(A, len, 0);
    A.resize(o);
    return A;
}

int fac[N], inv[N];
void init(int M) {
    fac[0] = fac[1] = inv[0] = inv[1] = 1;
    rep(i,2,M) inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
    rep(i,1,M) fac[i] = 1ll * fac[i - 1] * i % mod, 
        inv[i] = 1ll * inv[i - 1] * inv[i] % mod;
}

signed main() {
    file(w2b);
    cin >> n >> m;
    init(max(n, m)); 
    F.resize(n + 1), G.resize(n + 1);
    rep(i,0,n) F[i] = 1ll * qp(i, n) * inv[i] % mod, G[i] = ((i & 1) ? mod - inv[i] : inv[i]);
    F = mul(F, G);
    rep(i,1,n) {
        int t = 1ll * fac[i] * F[i] % mod * qp(i + 1, m) % mod;
        ans += ((n + i) & 1 ? mod - t : t);
        if (ans >= mod) ans -= mod;
    } cout << ans << '\n';
}

upd(8 months later): 看到 qwaszx 老师写过这个东西 link