Description
Solution
一般求填色方案的题目都要用到动态规划。
看到数据范围有一个超级大的数，那么很显然就是矩阵乘法，往那个地方去想。
首先要填颜色：
所以设f[i][j]表示i个格子填上j个颜色的方案数，那么
f[i][j]=f[i−1][j−1]∗(p−(j−1))+f[i−1][j]∗j
f[i][j]太麻烦了，有用的只有f[n][j]，所以我们设g[i]=f[n][i]
我们考虑从第i-1列的j种颜色推到第i列的k中颜色，因为重复的情况不是很好的去处理，所以我们可以枚举用到的颜色总数为x，那么构成的方案数为：
g[k]Ckp∗∑x=min(j,k,q)min(j+k,p)Cj+k−xj∗Cx−jp−j
这个东西画一个维恩图看一看就好了，因为要考虑j和k集合是那些重复，然后还要考虑，不重复的地方是那些数搞出来的。但是g[k]为什么要除以C(p,k)呢？因为g[k]是p个数的，现在要确定这些数之后的情况（每个颜色至少有一个，其实就是第二类斯特林数）
然后设
dp[i][j]表示第i列填j种颜色的情况，那么
dp[i][j]=∑dp[i−1][j]∗zhuan[j][k]
zhuan[j][k]=g[k]Ckp∗∑min(j+k,p)min(j,k,q)Cj+k−xj∗Ck−xp−j
所以，这个东西只用打一个矩阵乘法就好了。
Code
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fod(i,a,b) for(i=a;i>=b;i--)
using namespace std;
const int mo=998244353;
typedef long long ll;
const int maxn=107;
ll i,j,k,l,t,n,m,ans,p,q;
ll f[maxn][maxn],g[maxn];
ll fact[maxn],ni[maxn],dp[maxn][maxn],x;
struct node{
ll ju[maxn][maxn];
node friend operator *(node x,node y){
int i,j,k;node z;
memset(z.ju,0,sizeof(z.ju));
fo(k,1,100)fo(i,1,100)fo(j,1,100){
z.ju[i][j]=(z.ju[i][j]+x.ju[i][k]*y.ju[k][j]%mo)%mo;
}
return z;
}
}a,b;
node qsm(node x,ll y){
node z;int i;memset(z.ju,0,sizeof(z.ju));
fo(i,1,100)z.ju[i][i]=1;
for(;y;y/=2,x=x*x)if(y&1)z=z*x;
return z;
}
ll qsm1(ll x,ll y){
ll z=1;
for(;y;y/=2,x=x*x%mo)if(y&1)z=z*x%mo;
return z;
}
ll c(ll x,ll y){
if(x<y)return 0;
return fact[x]*ni[y]%mo*ni[x-y]%mo;
}
ll chu(ll x,ll y){
ll o=qsm1(y,mo-2);
return x*o%mo;
}
int main(){
freopen("color.in","r",stdin);
freopen("color.out","w",stdout);
fact[0]=ni[0]=1;
fo(i,1,100)fact[i]=fact[i-1]*i%mo;
ans=c(3,1);
ni[100]=qsm1(fact[100],mo-2);
fod(i,99,1)ni[i]=ni[i+1]*(i+1)%mo;
scanf("%lld%lld%lld%lld",&n,&m,&p,&q);
f[0][0]=1;
fo(i,1,n)fo(j,1,min(i,p))f[i][j]=(f[i-1][j-1]*(p-(j-1))%mo+f[i-1][j]*j%mo)%mo;
fo(i,1,p)b.ju[1][i]=f[n][i],g[i]=chu(f[n][i],c(p,i));
fo(i,1,p)fo(j,1,p){
fo(x,max(j,max(i,q)),min(i+j,p)){
a.ju[i][j]=(a.ju[i][j]+c(i,i+j-x)*c(p-i,x-i)%mo)%mo;
}
a.ju[i][j]=a.ju[i][j]*g[j]%mo;
}
b=b*qsm(a,m-1);
fo(i,1,p)ans=(ans+b.ju[1][i])%mo;
printf("%lld\n",ans);
}
