闲话 25.3.(4.5)

闲话

大学生活？退休生活！

just cause 3 真好玩，大家都去玩 /cy
lean game 真好玩，大家都去玩 /cy
没钱买 prove u can win，大家可以玩 /cy

BTW 快来做 STARGAZERS

推歌：Explore by 晴一番p feat. 初音未来；倒瞰清白 by Creuzer et al.。

详细揭秘：什么是单项式多项式？

单项式多项式(monopoly)

当 \(\text{M}\) 在 \(\mathbb{F}_p^{n\times m}\) 上秩为 \(k\) 的矩阵中均匀随机取值时，计算 \(\text{M}\) 中非零元素个数的期望。

\(1\le n, m, k, p \le 10^{18}\)，保证 \(p\) 为质数。

jjdw 版题解，就是删了很多不重要的内容。

虽然有期望，做法还是线性代数 + q-analog，但这题的难度真是省选 T1.5 吧！你要信我！

当不需要考虑权重时，计数秩为 \(k\) 的矩阵是 well-known 的，我们可以通过 \(\text{q-analog}\) 的手法处理其中一些问题。可以知道，\(\mathbb{F}_q^{n\times m}\) 中秩为 \(r\) 的矩阵数量为

\[\frac{\prod_{i = 0}^{r - 1}(q^n - q^i)\prod_{i = 0}^{r - 1}(q^m - q^i)}{\prod_{i = 0}^{r - 1}(q^r - q^i)} \]

证明：我们断言，若矩阵 \(A\in \text{F}_q^{n\times m}\) 的秩为 \(r\)，则定存在列满秩矩阵 \(B\in \text{F}_q^{n\times r}\) 与行满秩矩阵 \(C \in \text{F}_q^{r\times m}\)，使得 \(A = BC\)。

对断言的证明

为简便，记一个 \(\in \text{F}_q^{n\times m}\) 的矩阵 \(A\) 为 \(A_{n\times m}\)，并记 \(n\times m\) 全零矩阵为 \(\mathbf 0 _{n\times m}\)，\(r\times r\) 单位方阵为 \(E_r\)。

由于初等变换不改变矩阵的秩，总存在可逆矩阵 \(P, Q\in \text{F}_q^{n\times m}\) 使得

\[PAQ = \begin{bmatrix} E_r & \mathbf 0_{r\times (n - r)} \\ \mathbf 0 _{(m - r)\times r} & \mathbf 0 _{(m - r)\times (n - r)} \end{bmatrix} \]

即

\[\begin{aligned} A &= P^{-1} \begin{bmatrix} E_r & \mathbf 0_{r\times (n - r)} \\ \mathbf 0 _{(m - r)\times r} & \mathbf 0 _{(m - r)\times (n - r)} \end{bmatrix} Q^{-1} \\&= P^{-1} \begin{bmatrix} E_r \\ \mathbf 0 _{(m - r)\times r} \end{bmatrix}\begin{bmatrix} E_r & \mathbf 0_{r\times (n - r)} \end{bmatrix} Q^{-1} \\&= B_{n\times r} C_{r\times m} \end{aligned}\]

下面需要证明 \(B, C\) 的秩为 \(r\)。

由秩的性质可知，\(r = \text{rk}(A) \le \min(\text{rk}(B), \text{rk}(C))\)。而矩阵的秩不可能大于矩阵行数、列数中更小的值，则有 \(\text{rk}(B)\le r, \text{rk}(C)\le r\)。综合两个结论可知 \(\text{rk}(B)= r, \text{rk}(C)= r\)。

此方法被称为满秩分解。

只知道存在 \(B, C\) 无法计数，我们还需要知道一个 \(A\) 对应了多少对 \((B, C)\)。自然可以知道 \(A = BD_{r\times r}D^{-1}_{r\times r}C\)，其中 \(D_{r\times r}\) 为一个满秩矩阵。这样可以知道，一个 \(A\) 对应了 \([r : r]_q\) 对彼此不同的 \((B, C)\)。分别计数 \(B, C\)，再除以 \([r:r]_q\) 即可。\(\square\)

上面的证明给予了我们一定的启发。下面将 \(A = (\alpha_1, \dots, \alpha_n)\) 视作一个元素位于模 \(p\) 剩余类域 \(\mathbb F_p\) 内的 \(n \times m\) 矩阵处理，或 \(A \in M_{n\times m} (\mathbb F_p)\)。要求的就是所有秩为 \(k\) 的矩阵 \(A\) 中非零元素的期望。

随机变量 \(M\) 也可以写成 \(\sum_{i, j} M_{i,j}\)，其中 \(M_{i,j}\) 对位于 \((i, j)\) 的元素定义：\(M_{i,j} = 1\) 当且仅当 \(A\) 位于 \((i, j)\) 的元素非零。根据期望的线性性，可以将计算 \(M\) 转换为计算每个 \(M_{i, j}\) 并求和。

一个重要的观察是，对每个 \((i, j)\)，\(M_{i, j}\) 的分布都是相同的。这是由于初等变换不改变矩阵的秩，因此我们可以考虑一个双射：将给定矩阵的第 \(1\) 行与第 \(i\) 行交换、第 \(1\) 列与第 \(j\) 列交换。这就将每个元素都与位于 \((1, 1)\) 位置的元素建立了等价关系，自然有上述观察。那么我们知道 \(M\) 的期望即为 \(nm\times M_{1, 1}\)。因此，不需要分析整个矩阵。

考虑一个矩阵的行简化阶梯形式，即进行消元后，所有零行都在最下方，剩余行最左侧元素均为 \(1\)，位置按阶梯形由左上角向右下角排布的形式。这形式可以用来产生行等价类；两个矩阵行等价说明二者可以通过一系列初等行变换互化，也表明二者的行向量空间相同。

容易知道，一个秩为 \(k\) 的 \(n\times m\) 矩阵 \(A\) 的行简化阶梯形式有 \(k\) 个非零行。令 \(R\) 为这 \(k\) 行组成的 \(k\times m\) 矩阵。可以知道 \(R\) 的行向量空间是 \(A\) 的向量空间的一组基，那么进行满秩分解，得到对应的矩阵 \(C\)。

通过上面的表示方法，我们可以将秩为 \(k\) 的矩阵组成的集合表示为秩为 \(k\) 的 \(n\times k\) 矩阵组成的集合与秩为 \(k\) 的 \(k\times m\) 行简化阶梯矩阵组成的集合的笛卡尔积（或 \(A\) 和 \((C, R)\) 间形成了双射）。这说明可以通过独立地选择 \(C\) 与 \(R\) 来生成 \(A\)，并且这样的生成是均匀随机的。

可以知道

\[A_{1, 1} = \sum_{i = 1}^k C_{1, i} R_{i, 1} \]

由于 \(R\) 是行简化阶梯矩阵，\(R_{1, 1}\) 可能是 \(1\) 或 \(0\)，而 \(R_{i, 1}\) 对 \(i > 1\) 都为 \(0\)。因此可以知道，\(A_{1, 1} = C_{1, 1} R_{1, 1}\)。从而

\[\Pr[A_{1, 1} \neq 0] = \Pr[C_{1, 1}\neq 0]\times \Pr[R_{1, 1}\neq 0] \]

\(C\) 是秩为 \(k\) 的 \(n\times k\) 矩阵，因此 \(C\) 的第一列是任意长为 \(n\) 的非零向量，数量为 \(p^n - 1\)。这些非零向量中有 \(p^{n - 1} - 1\) 个向量首元素为 \(0\)，因此

\[\Pr [C_{1, 1}\neq 0] = \frac{p^n - p^{n - 1}}{p^n - 1} \]

\(R\) 和 \(C\) 携带的信息不同，应用的计算方式就不同。选择 \(R\) 就是在选择 \(A\) 的行向量空间。只要行向量空间中存在一个向量首项非零，\(R_{1, 1}\neq 0\) 就成立。正难则反，先计算 \(R_{1, 1} = 0\) 的方案数。这指出所有向量的首项都为 \(0\)，因此可以认为这些向量都是 \(m - 1\) 维的，在这其中选择一个 \(k\) 维子空间的方案数为

\[{m - 1\brack k}_p = \frac{\prod_{i = 0}^{k - 1} (p^{m - 1} - p^i)}{\prod_{i = 0}^{k - 1} (p^k - p^i)} \]

不难得到

\[\Pr[R_{1, 1}\neq 0] = 1 - \frac{{m - 1 \brack k}_p}{{m\brack k}_p} = \frac{p^m - p^{m - k}}{p^m - 1} \]

因此可以知道

\[\Pr[A_{1, 1} \neq 0] = \left(\frac{p^n - p^{n - 1}}{p^n - 1}\right)\left(\frac{p^m - p^{m - k}}{p^m - 1}\right) = \frac{(1 - 1/p)(1 - 1/p^k)}{(1 - 1/p^n)(1 - 1/p^m)} \]

最终答案为 \(nm \Pr[A_{1, 1} \neq 0]\)，复杂度显然是 \(O(\log nmk)\) 的。

提供一份暴力代码以验证正确性

#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); }
#define iot ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
#define timer cerr << 1. * clock() / CLOCKS_PER_SEC << '\n';
#define multi int T; cin >> T; while ( T -- )
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].nxt)
#define all(s) s.begin(), s.end()
#define eb emplace_back
#define pb pop_back
#define em emplace
const int N = 3e6 + 10;
const int inf = 0x3f3f3f3f;
const ll infll = 0x3f3f3f3f3f3f3f3fll;

// int mod;
const int mod = 3;
// const int mod = 10007;
// const int mod = 469762049, g = 3;
// const int mod = 998244353; // const int g = 3;
// const int mod = 1004535809, g = 3;
// const int mod = 1e9 + 7;
// const int mod = 1e9 + 9;
// const int mod = 1e9 + 3579, bse = 131;
struct FastMod { int m; ll b; void init(int _m) { m = _m; if (m == 0) m = 1; b = ((lll)1 << 64) / m; } FastMod(int _m) { init(_m); }int operator()(ll a) { ll q = ((lll)a * b) >> 64; a -= q * m; if (a >= m) a -= m; return a; } } Mod(mod);
template<const int mod, typename T> 
struct z {
    T val;
    z(T _val = 0) : val(Mod(Mod(_val) + mod)) {}
    const T operator * () const { return val; }
    T& operator * () { return val; }
    inline friend bool operator < (const z & a, const z & b) { return *a < *b; }
    inline friend bool operator < (const z & a, const T & b) { return *a < b; }
    inline friend bool operator <= (const z & a, const z & b) { return *a <= *b; }
    inline friend bool operator <= (const z & a, const T & b) { return *a <= b; }
    inline friend bool operator > (const z & a, const z & b) { return *a > *b; }
    inline friend bool operator > (const z & a, const T & b) { return *a > b; }
    inline friend bool operator >= (const z & a, const z & b) { return *a >= *b; }
    inline friend bool operator >= (const z & a, const T & b) { return *a >= b; }
    inline friend bool operator == (const z & a, const z & b) { return *a == *b; }
    inline friend bool operator == (const z & a, const T & b) { return *a == b; }
    inline friend bool operator != (const z & a, const z & b) { return *a != *b; }
    inline friend bool operator != (const z & a, const T & b) { return *a != b; }
    z& operator = (const z & b) { val = *b; return *this; }
    z& operator = (const T & b) { val = b; return *this; }
    void operator ++ () { ++ val; } void operator -- () { -- val; } 
	z operator - () { return z(val == 0 ? val : mod - val); } 

    friend z operator + (const z & a, const z & b) { z ret(*a + *b); if (*ret >= mod) *ret -= mod; return ret; }
    z operator += (const z & b) { val += *b; if (val >= mod) val -= mod; return *this; }
    friend z operator - (const z & a, const z & b) { z ret(*a - *b + mod); if (*ret >= mod) *ret -= mod; return ret; }
    z operator -= (const z & b) { val -= *b; if (val < 0) val += mod; return *this; }
    friend z operator * (const z & a, const z & b) { z ret(Mod(1ll * *a * *b)); return ret; }
    z operator *= (const z & b) { val = Mod(1ll * val * *b); return *this; }
    friend z operator / (const z & a, const z & b) { z ret(*a / *b); return ret; }
    z operator /= (const z & b) { val = val / *b; return *this; }
    friend z operator % (const z & a, const z & b) { z ret(*a % *b); return ret; }
    z operator %= (const z & b) { val = val % *b; return *this; }
    friend z operator + (const z & a, const T & b) { z ret(*a + b); if (*ret >= mod) *ret -= mod; return ret; }
    z operator += (const T & b) { val += b; if (val >= mod) val -= mod; return *this; }
    friend z operator - (const z & a, const T & b) { z ret(*a - b + mod); if (*ret >= mod) *ret -= mod; return ret; }
    z operator -= (const T & b) { val -= b; if (val < 0) val += mod; return *this; }
    friend z operator * (const z & a, const T & b) { z ret(Mod(1ll * *a * b)); return ret; }
    z operator *= (const T & b) { val = Mod(1ll * val * b); return *this; }
    friend z operator / (const z & a, const T & b) { z ret(*a / b); return ret; }
    z operator /= (const T & b) { val = val / b; return *this; }
    friend z operator % (const z & a, const T & b) { z ret(*a % b); return ret; }
    z operator %= (const T & b) { val = val % b; return *this; }

    friend istream& operator >> (istream& in, z& b) { in >> *b; return in; }
    friend ostream& operator << (ostream& out, z b) { out << *b; return out; }
}; using mint = z<mod, int>;
template <typename T1, typename T2> T1 qp(T1 a, T2 b) { T1 ret = 1; for (; b > 0; a = a * a, b >>= 1) if (b & 1) ret = ret * a; return ret; }
mint fac[N], inv[N], ifac[N];
void init_fac(int bnd, int mask = 0b111) {
    if ((mask >> 2) & 1) { fac[0] = fac[1] = 1; rep(i, 2, bnd) fac[i] = fac[i - 1] * i; }
    if ((mask >> 1) & 1) { inv[0] = inv[1] = 1; rep(i, 2, bnd) inv[i] = (mod - mod / i) * inv[mod % i]; }
    if ((mask >> 0) & 1) { ifac[0] = ifac[1] = 1; if (!((mask >> 1) & 1)) init_fac(bnd, 0b010); rep(i,1,bnd) ifac[i] = ifac[i - 1] * inv[i]; }
}
template <typename T1, typename T2> mint C(T1 n, T2 m) { if (n < m) return 0; return fac[n] * ifac[m] * ifac[n - m]; }
template <typename T1, typename T2> mint invC(T1 n, T2 m) { if (n < m) return 0; return ifac[n] * fac[m] * fac[n - m]; }


const int n = 3, m = 4;
struct mat {
	vector<vector<mint>> a;
	mat() { a = vector<vector<mint>>(n, vector<mint>(m)); }
	int rank() {
		int rk = 0;
		rep(col, 0, m - 1) {
			int t = rk;
			while (t < n && a[t][col] == 0)
				++t;
			if (t == n)
				continue;
			swap(a[t], a[rk]);
			mint inv = qp(a[rk][col], mod - 2);
			rep(i, col, m - 1)
				a[rk][i] *= inv;
			rep(i, 0, n - 1) {
				if (i == rk) continue;
				mint mul = a[i][col];
				rep(j, col, m - 1)
					a[i][j] -= a[rk][j] * mul;
			}
			++rk;
		}
		return rk;
	}
};

using ll = long long;
ll result[1000], cnt[1000];

void dfs(int pos, vector<int>& cur) {
	if (pos == n * m) {
		mat M;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < m; j++)
				M.a[i][j] = cur[i * m + j];
		int k = M.rank();
		++ cnt[k];
		int countNonZero = 0;
		for (int val : cur)
			if (val != 0)
				countNonZero++;
		result[k] += countNonZero;
		return;
	}
	for (int v = 0; v < mod; v++) {
		cur[pos] = v;
		dfs(pos + 1, cur);
	}
}

int main() {
	vector<int> arr(n * m, 0);
	dfs(0, arr);
	rep(k,0,min(n,m)) {
		cout << k << ' ';
		ll g = __gcd(result[k], cnt[k]);
		cout << result[k] / g << "/" << cnt[k] / g << endl;
	}
	return 0;
}

\(k = 1\) 的部分分：

容易知道答案即

\[\frac{\sum _{i=1}^m \binom{m}{i} \sum _{j=1}^n\binom{n}{j} i j (p-1)^{i-1} (p-1)^j }{\sum _{i=1}^m \binom{m}{i} \sum _{j=1}^n\binom{n}{j} (p-1)^{i-1} (p-1)^j } = \frac{nm (1-1/p)^2}{(1-1/p^m)(1-1/p^n)} \]

公式由 @H_Kaguya 提供。她说如果她没做出 \(k = 2\) 就不用放 credit，所以题解里没有提到她 /qd