【noip2016十连测round3】T3 涂色游戏 【矩阵快速幂优化dp】

涂色游戏

题解：
推一推公式。
我们让f[i][j]表示第i列有j种颜色的方案总数，k表示i-1列用了多少种颜色，l表示第i列用了多少种没有在i-1列出现的颜色,G(i,j)表示i个格子涂j种颜色的方案总数。
则f[i][j]=Σ(f[i−1][k]∗C(k,i−l)∗C(p−k,l)∗G(n,j))$f[i][j]=\Sigma (f[i-1][k]\ast C(k,i-l)\ast C(p-k,l)\ast G(n,j))$。其中要满足k+l>=q$k+l>=q$。
G(i,j)怎么求？
可以得到G(i,j)=ji−Σ(g[i][k]∗c[j][k])$G(i,j)=j^i-\Sigma (g[i][k]\ast c[j][k])$，其中k=1..j−1$k=\mathrm{1..}j-1$。这个式子代表先get到所有的涂法，再减去少于j种颜色的方案总数。
但是我们发现列数为1e9，怎么办？
太显然了吧！^ ^直接上矩阵快速幂优化。
在最后orz orz %%%zjr大佬。老中医保佑！
代码：

#include<cstdio>
#include<cstring>
#pragma GCC optimize (3)
typedef long long ll;
const int N=105;
const ll mod=998244353;
int n,m,p,q;
ll ans,c[N][N],g[N][N];
ll fast_pow(ll a,int x){
    ll res=1;
    while(x){
        if(x&1){
            res*=a;
            res%=mod;
        }
        x>>=1;
        a*=a;
        a%=mod;
    }
    return res;
}
struct matrix{
    ll a[N][N];
    matrix(){
        memset(a,0,sizeof(a));
    }
    matrix operator * (const matrix &b) const{
        matrix c;
        for(int i=1;i<=p;i++){
            for(int j=1;j<=p;j++){
                for(int k=1;k<=p;k++){
                    c.a[i][j]+=a[i][k]*b.a[k][j];
                    c.a[i][j]%=mod;
                }
            }
        }
        return c;
    }
}bg,x;
int main(){
    scanf("%d%d%d%d",&n,&m,&p,&q);
    if(p<q){
        puts("0");
        return 0;
    }
    for(int i=0;i<=p;i++){
        c[i][0]=c[i][i]=1;
    }
    for(int i=1;i<=p;i++){
        for(int j=1;j<=i;j++){
            c[i][j]=c[i-1][j-1]+c[i-1][j];
            c[i][j]%=mod;
        }
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=i&&j<=p;j++){
            g[i][j]=fast_pow(j,i);
            for(int k=1;k<j;k++){
                g[i][j]-=g[i][k]*c[j][k]%mod;
                g[i][j]=(g[i][j]%mod+mod)%mod;
            }
        }
    }
    for(int i=1;i<=p&&i<=n;i++){
        for(int j=1;j<=p&&j<=n;j++){
            for(int k=0;k<=i;k++){
                if(j+k<q){
                    continue;
                }
                x.a[j][i]+=c[j][i-k]*c[p-j][k]%mod*g[n][i]%mod;
                x.a[j][i]%=mod;
            }
        }
    }
    for(int i=1;i<=p;i++){
        bg.a[1][i]=c[p][i]*g[n][i]%mod;
    }
    m--;
    while(m){
        if(m&1){
            bg=bg*x;
        }
        m>>=1;
        x=x*x;
    }
    for(int i=1;i<=p;i++){
        ans+=bg.a[1][i];
        ans%=mod;
    }
    printf("%lld\n",ans);
    return 0;
}