闲话 25.3.(4.5)

闲话

大学生活？退休生活！

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

推歌：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})}

$\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}$ 使得

P A Q = [\begin{matrix} E_{r} & 0_{r \times (n - r)} \\ 0_{(m - r) \times r} & 0_{(m - r) \times (n - r)} \end{matrix}]

$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{matrix} E_{r} & 0_{r \times (n - r)} \\ 0_{(m - r) \times r} & 0_{(m - r) \times (n - r)} \end{matrix}] Q^{- 1} \\ = P^{- 1} [\begin{matrix} E_{r} \\ 0_{(m - r) \times r} \end{matrix}] [\begin{matrix} E_{r} & 0_{r \times (n - r)} \end{matrix}] Q^{- 1} \\ = B_{n \times r} C_{r \times m} \end{aligned}

$\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}

$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]

$\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}

$\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$ 维子空间的方案数为

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

${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{{[\binom{m - 1}{k}]}_{p}}{{[\binom{m}{k}]}_{p}} = \frac{p^{m} - p^{m - k}}{p^{m} - 1}

$\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] = (\frac{p^{n} - p^{n - 1}}{p^{n} - 1}) (\frac{p^{m} - p^{m - k}}{p^{m} - 1}) = \frac{(1 - 1 / p) (1 - 1 / p^{k})}{(1 - 1 / p^{n}) (1 - 1 / p^{m})}

$\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{n m (1 - 1 / p)^{2}}{(1 - 1 / p^{m}) (1 - 1 / p^{n})}

$\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

posted @ 2025-03-05 00:18 joke3579 阅读(36) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 闲话 22.11.8

· 闲话 23.1.4

· 2025.3.2 闲话

· 多项式全家桶

· 多项式小寄

Dr. Zhou

joke3579

闲话 25.3.(4.5)

闲话

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

$k = 1$ 的部分分：

最新随笔

我的标签

阅读排行榜

评论排行榜

推荐排行榜

joke3579

闲话 25.3.(4.5)

闲话

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

k=1k=1k = 1 的部分分：

最新随笔

我的标签

阅读排行榜

评论排行榜

推荐排行榜

$k = 1$ 的部分分：