循环矩乘——Luogu3746/BZOJ4870 [SHOI2017]组合数问题

题面：BZOJ4870 Luogu3746
第一次接触循环矩乘。。。
首先我们可以考虑DP，$f[i][j]$表示在i个物品中选取$modk$下余j的方案数。
状态转移很好想，$f[i][j]=f[i-1][j]+f[i-1][(j-1+k)mod k]$
然后发现这个DP可以用矩乘优化，矩阵大概长这个样子：

$$\left[ \begin{matrix} 1 & 1 & 0 & ··· & 0 & 0 \\ 0 & 1 & 1 & ··· & 0 & 0 \\ 0 & 0 & 1 & ··· & 0 & 0 \\ ··· & ··· & ··· & ··· & ··· & ··· \\ 0 & 0 & 0 & ··· & 1 & 1 \\ 1 & 0 & 0 & ··· & 0 & 1 \\ \end{matrix} \right]$$
如果直接矩乘的话时间复杂度是 $O(k^3logn)$，可以通过了
但是可以发现每一行都是上一行向右移一位，这是一个循环矩阵
循环矩阵有一个性质，循环矩阵自乘还是循环矩阵
所以其实我们只需记录第一行的信息就可以了，因为下面的信息可以通过上面信息平移得到
所以我们就可以把矩阵乘法这个地方优化到 $O(k^2)$，总的时间复杂度是 $O(k^2logn)$，很优秀了

#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <iostream>
#include <ctime>
#include <map>
#include <queue>
#include <cstdlib>
#include <string>
#include <climits>
#include <set>
#include <vector>
#define int long long
using namespace std;
inline int read(){
    int k=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){k=k*10+ch-'0';ch=getchar();}
    return k*f;
}
struct juzhen{
    int a[510];
    inline void clear(){memset(a,0,sizeof a);}
}f;
int t,n,k,r,MOD;
inline juzhen cheng(juzhen a,juzhen b){
    juzhen c;c.clear();
    for(int i=0;i<k;i++)
        for(int j=0;j<k;j++)c.a[(i+j)%k]=(c.a[(i+j)%k]+a.a[i]*b.a[j]%MOD)%MOD;
    return c;
}
inline juzhen mi(juzhen a,int b){
    juzhen x,y;x=a;y=a;
    while(b){
        if(b&1)x=cheng(x,y);
        y=cheng(y,y);b>>=1;
    }
    return x;
}
signed main()
{
    n=read();MOD=read();k=read();r=read();
    f.clear();
    f.a[0]++;f.a[1%k]++;
    f=mi(f,n*k-1);
    printf("%lld\n",f.a[r]);
    return 0;
}