LOJ2325「清华集训 2017」小Y和恐怖的奴隶主（期望概率+矩阵快速幂）

LOJ2325「清华集训 2017」小Y和恐怖的奴隶主

原题地址：https://loj.ac/problem/2325

题意：

"A fight? Count me in!" 要打架了，算我一个。

"Everyone, get in here!" 所有人，都过来！

小Y是一个喜欢玩游戏的OIer。一天，她正在玩一款游戏，要打一个Boss。

虽然这个Boss有$10^{100}$​​点生命值，但它只带了一个随从——一个只有 $m$ 点生命值的“恐怖的奴隶主”。

这个“恐怖的奴隶主”有一个特殊的技能：每当它被扣减生命值但没有死亡（死亡即生命值 $\leq 0$），且Boss的随从数量小于上限 $k$，便会召唤一个新的具有 $m$点生命值的“恐怖的奴隶主”。

现在小Y可以进行 $n$次攻击，每次攻击时，会从Boss以及Boss的所有随从中的等概率随机选择一个，并扣减$1$ 点生命值，她想知道进行 $n$次攻击后扣减Boss的生命值点数的期望。为了避免精度误差，你的答案需要对$998244353$取模。

数据范围
$1≤T≤1000,1≤n≤10^{18} ,1≤m≤3,1≤k≤8$

题解：

最朴素的DP，利用概率算期望，记录第i轮，分别有j，k，l个血量为1，2，3的奴隶主的概率。
$DP[i+1][j'][k'][l']+=\frac{k/j/l}{j+k+l+1}DP[i][j][k][l]$，
每个状态都有$\frac{1}{j+k+l+1}$攻击boss，有$\frac{1}{j+k+l+1}*dp[i][j][k][l]$的贡献。

n太大了，考虑矩阵快速幂。
如何把后面三维压成一维呢？
发现所有合法状态最多$C_{10}^2+C_9^2+...+C_2^2=165$个，
可以直接把所有状态之间的转移关系列出来，然后跑矩阵快速幂。
但是，复杂度为$O(tot^3logn)$，要T。

复杂度的瓶颈在于T组数据，每次矩阵*矩阵都是$tot^3$,
但其实我们只需要答案那个$1 \times tot$的数组。
那么，预处理出2的幂次方的矩阵，由于矩阵乘法的结合率，
$Ans=Ini*Max^n=Ini*Max^{(101...0)_2}=Ini*Max^{2^{s}}*Max^{2^{s-1}}*...Max^{2^{0}}$
一项一项，每次都是Ans向量乘矩阵，复杂度$O(tot^3logn+T*tot^2logn)$

学习：

• 对于和的期望不好求，可以算每个状态的概率*它的贡献。（不过逆推可以较容易地求得和的期望）
• 表示清晰的多维DP数组难以矩阵快速幂，但是可以看看所有状态有多少，把状态间的关系列出来，就可以转移。
• 对于在贡献在中间统计的情况，例如这个$\frac{1}{j+k+l+1}*dp[i][j][k][l]$的处理，矩阵多开一位当计数器即可，既继承上一个，又加上新增的。
• 当然这道题可以逆推，即$dp[i][j][k][l]$表示从第i轮到n轮，期望的攻击boss次数。$DP[i-1][j'][k'][l']+=\frac{k/j/l}{j+k+l+1}(DP[i][j][k][l]+w)$，攻击boss时w为1，否则为0。
  顺推：$now[1][i]=\sum pre[1][j]*M[j][i]$
  倒推：$now[i][1]=\sum M[i][j]*pre[j][1]$

优化：

• 预处理
• 利用好矩阵的运算律。
• 少模一点。可以设置一个很大的lim，超了就减回来，最后再模。

代码：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#define LL unsigned long long
using namespace std;
const int mod=998244353;
const LL lim=16940360401038606353llu;
const int N=10;
const int MXN=170;
int T,m,K,tot=0,id[N][N][N];
LL n,inv[N],ans[170],tmp[MXN];
struct Matrix
{
    LL M[MXN][MXN];
    void init() {memset(M,0,sizeof(M));}
    Matrix operator*(const Matrix &A)
    {
        Matrix ret; ret.init();
        for(int i=1;i<=tot+1;i++)
        for(int j=1;j<=tot+1;j++)
        {
            for(int k=1;k<=tot+1;k++)
            {
                ret.M[i][j]+=M[i][k]*A.M[k][j];
                if(ret.M[i][j]>=lim) ret.M[i][j]-=lim;
            }
            ret.M[i][j]%=mod;
        }
        return ret;
    }
}f[65];
inline LL modpow(LL A,LL B)
{
    LL ret=1; LL base=A;
    for(;B;B>>=1) 
    {
        if(B&1) ret=ret*base%mod; 
        base=base*base%mod;
    }
    return ret;
}
void muti(Matrix &A)
{
    memset(tmp,0,sizeof(tmp));
    for(int i=1;i<=tot+1;i++)  //tmp[1][i]=ans[1][j]*M[j][i]
    {
        for(int j=1;j<=tot+1;j++)
        {
            tmp[i]+=ans[j]*A.M[j][i];
            if(tmp[i]>=lim) tmp[i]-=lim;
        }
        tmp[i]%=mod;
    }
    memcpy(ans,tmp,sizeof(tmp));
}
int main()
{
    scanf("%d%d%d",&T,&m,&K);
    for(int i=0;i<=K;i++)
    for(int j=0;j<=((m>1)?K-i:0);j++)
    for(int k=0;k<=((m>2)?K-i-j:0);k++)
    id[i][j][k]=++tot;
    for(int i=0;i<=K+1;i++) inv[i]=modpow(i,mod-2);
    for(int i=0;i<=63;i++) f[i].init();
    for(int i=0;i<=K;i++)
    for(int j=0;j<=((m>1)?K-i:0);j++)
    for(int k=0;k<=((m>2)?K-i-j:0);k++)
    {
        int cur=id[i][j][k],nk=(i+j+k)<K; LL iv=inv[i+j+k+1];
        if(m==1) f[0].M[cur][id[i-1][j][k]]=iv*i%mod;
        if(m==2)
        {
            if(i) f[0].M[cur][id[i-1][j][k]]=iv*i%mod;
            if(j) f[0].M[cur][id[i+1][j-1+nk][k]]=iv*j%mod;
        }
        if(m==3)
        {
            if(i) f[0].M[cur][id[i-1][j][k]]=iv*i%mod;
            if(j) f[0].M[cur][id[i+1][j-1][k+nk]]=iv*j%mod;
            if(k) f[0].M[cur][id[i][j+1][k-1+nk]]=iv*k%mod;
        }
        f[0].M[cur][cur]=f[0].M[cur][tot+1]=iv;
    }
    f[0].M[tot+1][tot+1]=1;
    for(int i=1;i<=63;i++) f[i]=f[i-1]*f[i-1];
    while(T--) 
    {
        scanf("%lld",&n);
        memset(ans,0,sizeof(ans));
        if(m==1) ans[id[1][0][0]]=1;
        else if(m==2) ans[id[0][1][0]]=1;
        else if(m==3) ans[id[0][0][1]]=1;
        for(int i=0;n;n>>=1,i++) if(n&1) muti(f[i]);
        printf("%lld\n",ans[tot+1]);
    }
    return 0;
}