CF1864H Asterism Stream【概率 DP，矩阵优化】
给定一变量，初始为 \(1\)，每次等概率随机进行以下两种操作之一：
- 令 \(x\) 加一。
- 令 \(x\) 乘二。
求期望多少次操作之后 \(x\) 会 \(\ge n\)。
\(T\) 组数据，\(T\le 100\)，\(n\le 10^{18}\)。
对着 aw 老师的题解学的，感觉太深刻。
下文中所有表示下标的分式均为下取整。
设 \(f_i\) 表示 \(i\) 被遍历到的期望次数，则答案为 \(\sum\limits_{i = 1} ^ {n - 1} f_i\)。容易写出转移：\(f_i = \begin{cases}\frac{1}{2}(f_{i-1}+f_{i/2})&2 \mid i\\\frac{1}{2}f_{i-1}&2 \nmid i\end{cases}\)。
注意到每个状态的前驱较少，尝试用矩阵描述转移。考虑向量 \(A_i\) 里需要记录哪些信息。首先 \(f_i\) 必须有。当 \(i\) 为偶数时，还要知道 \(f_{i/2}\)。这样，如果 \(i\) 是偶数但不是 \(4\) 的倍数，那么从 \(A_{i - 1}\) 转移到 \(A_{i}\) 是容易的：\(f_i\) 的前驱是 \(f_{i - 1}\) 和 \(f_{i / 2}\)，而 \(\frac{i}{2}\) 是奇数，所以 \(f_{i / 2} = \frac 1 2 f_{i / 2 - 1} = \frac 1 2 f_{(i - 1) / 2}\)，可由 \(A_{i - 1}\) 记录的 \(f_{(i - 1) / 2}\) 直接求出。对于 \(f_{i / 2}\) 同理。但如果 \(i\) 是 \(4\) 的倍数，那就还要继续递归下去，因为 \(f_{i / 2} = \frac 1 2 (f_{i / 2 - 1} + f_{i / 4})\)。
综上，我们考虑在 \(A_i\) 中记录所有 \(f_{i / 2 ^ k}\)，其中 \(k\) 为整数且 \(0\leq k\leq \log_2 \max n\)。为了防止对于边界的讨论，根据 \(f_1 = 1\) 且 \(f_1 = \frac 1 2 f_0\)，不妨认为 \(f_0 = 2\)。这样信息就封闭了。
根据转移过程容易发现，转移矩阵只和当前数 \(i\) 所包含的 \(2\) 的幂次 \(\nu_2(i)\) 有关。一个直观的理解是把 \(i\) 写成二进制形式，那么实际上 \(\nu_2(i)\) 就是 \(i\) 最低的 \(1\) 对应的位置，这样容易证明 \(i\) 和 \(i-1\) 只有最低的 \(\nu_2(i)\) 位不同。设当前元素所包含的 \(2\) 的幂次为 \(j\) 的转移矩阵为 \(F_j\)。我们知道关于 \(2\) 的幂的一个简单性质：对于 \(n < 2 ^ k\)，\(\nu_2(n) = \nu_2(n + 2 ^ k)\)，这样就可以类似矩阵快速幂那样递推了。设 \(H_j = \prod\limits_{i=1}^{2^j} F_i\)，\(G_j = \prod\limits_{i=1}^{2^{j+1}-1} F_i\)，则有 \(H_0 = G_0=F_0\)，且 \(H_i= G_{i-1} \times F_i\)，\(G_i = H_i \times G_{i-1}\)。
最后我们需要知道 \(\sum\limits_{i = 1} ^ {n - 1} f_i\)，把这个也放进 \(A_i\) 里一起转移就行了。时间复杂度 \(\mathcal{O}((T + \log n) \log^3n)\)。
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned long long ull;
typedef pair <int, int> pi;
#define fi first
#define se second
constexpr int N = 60;
constexpr int mod = 998244353, iv = (mod + 1) >> 1;
bool Mbe;
struct mat {
ull a[N + 1][N + 1];
void init() { memset(a, 0, sizeof(a)); }
mat operator * (const mat &x) const {
mat y; y.init();
for (int i = 0; i <= N; i++) {
for (int k = 0; k <= i; k++) {
for (int j = 0; j <= k; j++) y.a[i][j] += a[i][k] * x.a[k][j];
if ((k & 15) == 15 || k == i) {
for (int j = 0; j <= i; j++) y.a[i][j] %= mod;
}
}
}
return y;
}
mat operator / (const mat &x) const {
mat y; y.init();
for (int k = 0; k <= N; k++) {
for (int j = 0; j <= k; j++) y.a[0][j] += a[0][k] * x.a[k][j];
if ((k & 15) == 15 || k == N) {
for (int j = 0; j <= N; j++) y.a[0][j] %= mod;
}
}
return y;
}
} f[N], g[N], h[N], I;
LL n;
void solve() {
cin >> n;
mat ans = I;
for (int i = 59; ~i; i--) {
if (n >> i & 1) ans = ans / h[i];
}
cout << (ans.a[0][0] + mod - 2) % mod << "\n";
}
bool Med;
int main() {
// fprintf(stderr, "%.9lfMb\n", 1.0 * (&Mbe - &Med) / 1048576);
ios :: sync_with_stdio(false);
cin.tie(0), cout.tie(0);
for (int i = 0; i < N; i++) {
f[i].a[0][0] = 1;
f[i].a[1][0] = 1;
for (int j = 0; j <= i; j++) {
int coef = iv;
for (int k = j; k < i; k++) {
f[i].a[k + 1][j + 1] = coef;
coef = 1LL * coef * iv % mod;
}
f[i].a[i + 1][j + 1] = coef;
}
for (int j = i + 1; j < N; j++) {
f[i].a[j + 1][j + 1] = 1;
}
}
h[0] = g[0] = f[0];
for (int i = 1; i < N; i++) {
h[i] = g[i - 1] * f[i];
g[i] = h[i] * g[i - 1];
}
for (int i = 1; i <= N; i++) I.a[0][i] = 2;
int t = 1;
cin >> t;
while (t--) solve();
// cerr << 1e3 * clock() / CLOCKS_PER_SEC << "ms\n";
return 0;
}
