【题解】「清华集训 2017」小 Y 和恐怖的奴隶主 状压/期望 DP+矩阵优化 loj2325 —Qiuly

青一的电脑，为你点赞！

卡常真恶心。

发现恐怖的奴隶主数量的上限 $k$ 和恐怖的奴隶主血量的上限 $m$ 都特别小，考虑状压当前所有恐怖的奴隶主。由于每一个恐怖的奴隶主没有不同，我们只需要记下每个不同的血量的恐怖的奴隶主的个数即可。

设 $f_{i,a,b,c}$ 表示还剩 $i$ 次攻击，当前局面下有 $a$ 个血量为 $1$ 的恐怖的奴隶主，$b$ 个血量为 $2$ 的恐怖的奴隶主和 $c$ 个血量为 $3$ 的恐怖的奴隶主时对 $\rm{Boss}$ 的期望伤害。

转移很简单，如果以 $\frac{1}{a+b+c+1}$ 的概率打中了 $\rm{Boss}$ ，那么：
$$
f_{i,a,b,c}+=\frac{1}{a+b+c+1}(f_{i-1,a,b,c}+1)
$$
然后就是打中恐怖的奴隶主的情况，将打中的恐怖的奴隶主的血量分开转移即可，注意一下如果打中的恐怖的奴隶主没有死亡就会召唤一个新的血量为 $m$ 的恐怖的奴隶主，计数的时候就需要在 $m$ 的位置多加上一个 $1$ ，当然如果 $a+b+c=k$ 了就不能召唤恐怖的奴隶主了。

$n\leq 10^{18}$ ，使用矩阵优化转移即可，将 $a,b,c$ 的状态压起来，可以通过组合数算出最多只有 $165$ 种不同的状态。

然而现在还是过不去，$T$ 很大，每一次转移矩阵的 $2^x$ 次方都要重新算，但是对于每一次算出来的却都是相同的，可以将转移矩阵的所有的 $2^x$ 次方的值放到外面预处理。

有一个卡常技巧，我们可以取一个大小和 $2\times 10^{18}$ 差不多且是 $998244353$ 的倍数的值，矩阵乘法的时候，不要直接 $\rm{mod}$ ，现将所有贡献计算起来，如果超过了这个数就减去其，最后算入答案矩阵的时候再取模。

注意矩阵的大小也不全是 $165\times 165$ ，比如说每一次的初始矩阵就是 $165\times 1$ 的，记一下矩阵的大小，又可以优化矩阵乘法的速度。

假设状态总数为 $S$ ，那么时间复杂度为：
$$
O\left(S^3\log n+TS^2\log n\right)
$$
跑不满，如果跑满了就过不去了，记矩阵大小也是让其跑不满。

Code：

#include <cstdio>
#include <string>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define freopen(x) {freopen(x".in","r",stdin);freopen(x".out","w",stdout);}

const int N=170;
const ll mod=998244353;
const ll limit=2096313141300000000;

template <typename _Tp> inline void IN(_Tp&x) {
    char ch;bool flag=0;x=0;
    while(ch=getchar(),!isdigit(ch)) if(ch=='-') flag=1;
    while(isdigit(ch)) x=x*10+ch-'0',ch=getchar();
    if(flag) x=-x;
}

struct Matrix {
    int n,m,a[N][N];
    Matrix() {}
    Matrix(int _n,int _m):n(_n),m(_m) {memset(a,0,sizeof(a));}
    friend inline Matrix operator * (const Matrix &x, const Matrix &y) {
        Matrix res=Matrix(x.n,y.m);
        for(int i=1;i<=res.n;++i)
            for(int j=1;j<=res.m;++j) {
                ll sum=0;
                for(int k=1;k<=x.m;++k) {
                    sum+=1ll*x.a[i][k]*y.a[k][j];
                    if(sum>=limit) sum-=limit;
                }
                res.a[i][j]=sum%mod;
            }
        return res;
    }
}tre[60],Begin,Ans;

ll n;
int T,m,k,tot,inv[10],ID[10][10][10];

int main() {
    // freopen("2325");
    IN(T),IN(m),IN(k);
    inv[0]=0,inv[1]=1;
    for(int i=2;i<=k+1;++i) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    for(int a=0;a<=k;++a) 
        for(int b=0;b<=k;++b) 
            for(int c=0;c<=k;++c) 
                if(a+b+c<=k) ID[a][b][c]=++tot;
    tre[0]=Matrix(tot+1,tot+1);
    Begin=Matrix(tot+1,1),Begin.a[tot+1][1]=1;
    for(int a=0;a<=k;++a) for(int b=0;b<=k;++b) for(int c=0;c<=k;++c) if(a+b+c<=k) {
        int now=ID[a][b][c],tmp=inv[a+b+c+1];
        int resa=1ll*tmp*a%mod,resb=1ll*tmp*b%mod,resc=1ll*tmp*c%mod;
        if(a) tre[0].a[now][ID[a-1][b][c]]=resa;
        if(b) {
            if(a+b+c==k) tre[0].a[now][ID[a+1][b-1][c]]=resb;
            else {
                if(m==1) tre[0].a[now][ID[a+2][b-1][c]]=resb;
                else if(m==2) tre[0].a[now][ID[a+1][b][c]]=resb;
                else tre[0].a[now][ID[a+1][b-1][c+1]]=resb;
            }
        }
        if(c) {
            if(a+b+c==k) tre[0].a[now][ID[a][b+1][c-1]]=resc;
            else {
                if(m==1) tre[0].a[now][ID[a+1][b+1][c-1]]=resc;
                else if(m==2) tre[0].a[now][ID[a][b+2][c-1]]=resc;
                else tre[0].a[now][ID[a][b+1][c]]=resc;
            }
        }
        tre[0].a[now][now]=tre[0].a[now][tot+1]=tmp;
    }
    tre[0].a[tot+1][tot+1]=1;
    for(int i=1;i<=59;++i) tre[i]=tre[i-1]*tre[i-1];
    int ans_id=ID[m==1][m==2][m==3];
    while(T--) {
        IN(n),Ans=Begin;
        for(int i=0;i<=59;++i) if((n>>i)&1) Ans=tre[i]*Ans;
        printf("%d\n",Ans.a[ans_id][1]);
    }
    return 0;
}