WA 85，并且从实际情况来看似乎是两个sub都有错的。。。

#include<bits/stdc++.h>
using namespace std;
#define int long long
int n,m,mod,k,c[2001][2001];
struct jz{int a[201][201];}b,t,nw;
int ksm(int x,int y){
	int ans=1;while(y){
		if(y&1)ans=ans*x%mod;
		x=x*x%mod,y/=2;
	}return ans;
}
jz calc(jz x,jz y){
	memset(t.a,0,sizeof(t.a));
	for(int i=0;i<=m+1;i++)for(int k=0;k<=m+1;k++)if(x.a[i][k])for(int j=0;j<=m+1;j++)(t.a[i][j]+=x.a[i][k]*y.a[k][j])%=mod;
	return t;
}
jz ksm(jz x,int y){
	jz ans=x;y--;while(y){
		if(y&1)ans=calc(ans,x);
		x=calc(x,x);y/=2;
	}return ans;
}
int cal(int x,int y){
	if(y==1)return x;
	if(y&1)return (cal(x,y-1)+1)*x%mod;
	return cal(x,y/2)*(1+ksm(x,y/2))%mod;
}
main(){
	cin>>n>>m>>mod;k=m*m-m*3+3;k%=mod;
	if(m<=200){
		for(int i=0;i<=m;i++)for(int j=0;j<=i;j++){
			c[i][j]=j?(c[i-1][j-1]+c[i-1][j])%mod:1;
			nw.a[i][j]=k*c[i][j]%mod;
		}for(int i=0;i<=m;i++)b.a[i][0]=1;
		nw.a[m+1][m+1]=nw.a[m+1][m]=1;
		b=calc(ksm(nw,n-1),b);cout<<(b.a[m+1][0]+b.a[m][0])*ksm(2,m)%mod*m%mod*(m-1)%mod;
	}else{
		int num=0;
		for(int i=1;i<=mod;i++)num=num+ksm(i,m)%mod*ksm(k,i-1)%mod;
		num%=mod;num=num*(cal(ksm(k,mod),n/mod-1)+1)%mod;
		for(int i=n-n%mod+1;i<=n;i++)(num+=ksm(i,m)*ksm(k,i-1)%mod)%=mod;
		cout<<num*ksm(2,m)%mod*m%mod*(m-1)%mod;
	}
	return 0;
}```