矩阵乘法
我们有一个 DP 式,它的转移系数相对固定,不受 DP 值的变化而变化,可以递推
且它通常有一或两维状态,可以分为很多阶段,但每个阶段中状态数不多
现在我们要递推很多次,是线性复杂度接受不了的
便把 DP 的式子写为矩阵的形式,一般在 \(O(w^3\log n)\) 复杂度内计算(\(w\) 为矩阵的大小)
它不仅应用于优化递推,还可以做到求某个区间的 递推/DP 值,此时再用类似线段树的数据结构维护每个区间的矩阵乘积
例题
预处理转移矩阵
这个题,看 \(m,k\) 小且不变,考虑加入状态中
以 \(m=3\) 为例,剩下的把状态减小即可
设 \(f(i,a,b,c)\) 表示当前攻击了 \(i\) 次,此时随从中生命值为 \(1\) 的有 \(a\) 个,为 \(2\) 的有 \(b\) 个,为 \(3\) 的有 \(c\) 个
此时通过枚举 \((a,b,c)\) 发现合法的三元组不超过 170 个
但是 \(n\) 很大,想到矩阵乘法
初始和最终求出的为行向量,初始局面对应三元组概率为 \(1\)
我们列出 \(170\times 170\) 的矩阵,每次枚举两个三元组编号为 \(i,j\),判断它们能否一步到达,分类讨论是攻击的谁,矩阵中 \(a_{i,j}\) 乘以相应概率
但要求的是期望,只需把矩阵加一行一列,行向量加一格代表期望,当不攻击随从,即 \(i=j\) 时贡献期望,为 \(1\times P[\text{攻击boss}]\)
但是多组询问中,每次都矩阵快速幂,\(O(170^3T\log n)\) 显然过不去
用到一个技巧:
发现转移矩阵不变,且注意到转移矩阵乘行/列向量复杂度为 \(O(w^2)\)
考虑把所有转移矩阵的 \(2^0,2^1\dots,2^{\log n}\) 预处理,矩阵乘法有结合律,把行向量与需要乘的矩阵相乘
复杂度 \(O((w^3+Tw^2)\log n),w=170\),卡个常就过了
代码:
#include<bits/stdc++.h>
#define reg register
using namespace std;
typedef long long ll;
const ll mod = 998244353;
struct matrix
{
ll a[170][170];
inline matrix() {memset(a, 0, sizeof(a));}
}f, I, kong, b, pw[65], ans;
struct group
{
ll p, q, r;
}lin[170];
inline ll read()
{
reg char ch = getchar(); reg ll x = 0;
while(ch < '0' || ch > '9') ch = getchar();
while(ch >= '0' && ch <= '9') x = (x << 1) + (x << 3) + ch - '0', ch = getchar();
return x;
}
inline void print(ll x)
{
if(x / 10) print(x / 10);
putchar(x % 10 + '0');
}
ll t, n, m, k, cnt, targ, pr, invf[110];
inline ll qmi(ll a, ll b)
{
reg ll res = 1;
while(b)
{
if(b & 1) res = res * a % mod;
b >>= 1, a = a * a % mod;
}
return res;
}
inline ll inv(ll x) {return qmi(x, mod - 2);}
inline ll add(ll a, ll b) {return (a + b >= mod) ? (a + b - mod) : (a + b);}
inline matrix mul(matrix x, matrix y)
{
matrix z;
for(reg ll i = 1; i <= cnt; ++i)
for(reg ll u = 1; u <= cnt; ++u)
if(x.a[i][u])
for(reg ll j = 1; j <= cnt; ++j)
z.a[i][j] = add(z.a[i][j], x.a[i][u] * y.a[u][j] % mod);
return z;
}
inline matrix hmul(matrix x, matrix y)
{
matrix res;
for(reg ll i = 1; i <= cnt; ++i)
for(reg ll j = 1; j <= cnt; ++j)
res.a[1][i] = add(res.a[1][i], x.a[1][j] * y.a[j][i] % mod);
return res;
}
inline void work1()
{
for(reg ll i = 1; i < cnt; ++i)
for(reg ll j = 1; j < cnt; ++j)
{
if(lin[i].p == lin[j].p + 1)
f.a[i][j] = invf[lin[i].p + 1] * lin[i].p % mod;
else if(i == j) f.a[i][i] = f.a[i][cnt] = invf[lin[i].p + 1] % mod;
}
}
inline void work2()
{
for(reg ll i = 1; i < cnt; ++i)
for(reg ll j = 1; j < cnt; ++j)
{
if(lin[i].p == lin[j].p + 1 && lin[i].q == lin[j].q)
f.a[i][j] = invf[lin[i].p + lin[i].q + 1] * lin[i].p % mod;
else if(lin[i].p + lin[i].q < k && lin[i].p == lin[j].p - 1 && lin[i].q == lin[j].q)
f.a[i][j] = invf[lin[i].p + lin[i].q + 1] * lin[i].q % mod;
else if(lin[i].p + lin[i].q == k && lin[i].p == lin[j].p - 1 && lin[i].q == lin[j].q + 1)
f.a[i][j] = invf[lin[i].p + lin[i].q + 1] * lin[i].q % mod;
else if(i == j) f.a[i][i] = f.a[i][cnt] = invf[lin[i].p + lin[i].q + 1] % mod;
}
}
inline void work3()
{
for(reg ll i = 1; i < cnt; ++i)
for(reg ll j = 1; j < cnt; ++j)
{ // j 是下一步的可能状态,i 是上一步,a[i][j] 为 i->j 的概率
if(lin[i].p == lin[j].p + 1 && lin[i].q == lin[j].q && lin[i].r == lin[j].r) // 打生命值为 1
f.a[i][j] = invf[lin[i].p + lin[i].q + lin[i].r + 1] * lin[i].p % mod;
else if(lin[i].p + lin[i].q + lin[i].r < k && lin[i].p == lin[j].p - 1 && lin[i].q == lin[j].q + 1 && lin[i].r == lin[j].r - 1)
f.a[i][j] = invf[lin[i].p + lin[i].q + lin[i].r + 1] * lin[i].q % mod; // 打生命值为 2
else if(lin[i].p + lin[i].q + lin[i].r == k && lin[i].p == lin[j].p - 1 && lin[i].q == lin[j].q + 1 && lin[i].r == lin[j].r)
f.a[i][j] = invf[lin[i].p + lin[i].q + lin[i].r + 1] * lin[i].q % mod;
else if(lin[i].p + lin[i].q + lin[i].r < k && lin[i].p == lin[j].p && lin[i].q == lin[j].q - 1 && lin[i].r == lin[j].r)
f.a[i][j] = invf[lin[i].p + lin[i].q + lin[i].r + 1] * lin[i].r % mod; // 打生命值为 3
else if(lin[i].p + lin[i].q + lin[i].r == k && lin[i].p == lin[j].p && lin[i].q == lin[j].q - 1 && lin[i].r == lin[j].r + 1)
f.a[i][j] = invf[lin[i].p + lin[i].q + lin[i].r + 1] * lin[i].r % mod;
else if(i == j) // 打 boss
f.a[i][i] = f.a[i][cnt] = invf[lin[i].p + lin[i].q + lin[i].r + 1];
}
}
int main()
{
t = read(), m = read(), k = read();
for(reg ll x = 0; x <= 8; ++x)
for(reg ll y = 0; y <= 8; ++y)
for(reg ll z = 0; z <= 8; ++z)
if(x + y + z <= k && (m > 2 || !z) && (m > 1 || (!y && !z)))
{
++cnt, lin[cnt] = (group){x, y, z};
if(m == 1 && x == 1) targ = cnt;
else if(m == 2 && x == 0 && y == 1) targ = cnt;
else if(m == 3 && x == 0 && y == 0 && z == 1) targ = cnt;
}
// printf("%lld", cnt); 最多共 165 种合法状态
++cnt;
for(reg int i = 1; i <= 100; ++i) invf[i] = qmi(i, mod - 2);
for(reg ll i = 1; i <= cnt; ++i) I.a[i][i] = 1;
b.a[1][targ] = f.a[cnt][cnt] = 1;
if(m == 1) work1();
else if(m == 2) work2();
else work3();
pw[0] = f;
for(reg ll i = 1; i <= 61; ++i) pw[i] = mul(pw[i - 1], pw[i - 1]); // 预处理矩阵的 2^i 次方
while(t--)
{
n = read(), ans = b;
// printf("%lld \n", f.a[cnt][2]);
for(reg ll i = 0; i <= 61; ++i)
if((n >> i) & 1) ans = hmul(ans, pw[i]); // 行向量乘矩阵, O(166^2)(不会 T)
print(ans.a[1][cnt] % mod), putchar('\n'); // 最后一格表示期望
// printf("%lld %lld\n", ans.a[][], ans.a[][]);
}
return 0;
}
标签:cnt,lin,矩阵,&&,乘法,ll,mod
From: https://www.cnblogs.com/KellyWLJ/p/18015718