[矩阵乘法] LOJ#2002. 「SDOI2017」序列计数

注意到 p 很小，容易想到矩乘。fi→f(i+j)%P。
关于那个质数的限制，只需要 1 到 m的方案 减去 1 到 m 扣去所有质数的方案即可。
构造 1 到 m 转移矩阵 O(P2) 搞，∑mk=1[k%P=i] = ⌊m−ik⌋+[i>0]。
扣去质数直接暴力 O(m/lnm∗P)。

#include<cstdio>#include<cstring>#include<algorithm>using namespace std;typedef long long LL;const int maxn=20000005,MOD=20170408;int ans,n,m,P,p[maxn],vis[maxn];void getP(){    for(int i=2;i<=m;i++){        if(!vis[i]) p[++p[0]]=i;        for(int j=1;j<=p[0]&&(LL)i*p[j]<=m;j++){            vis[i*p[j]]=true;            if(i%p[j]==0) break;        }    }}struct Matrix{    int n,m,a[105][105];     Matrix(){ memset(a,0,sizeof(a)); n=m=0; }    Matrix operator * (const Matrix &b){        Matrix c; c.n=n; c.m=b.m;        for(int i=0;i<c.n;i++)         for(int j=0;j<c.m;j++)           for(int k=0;k<m;k++) (c.a[i][j]+=(LL)a[i][k]*b.a[k][j]%MOD)%=MOD;        return c;    }} T,Ans;Matrix Pow(Matrix a,int b){    Matrix res; res.n=res.m=P; for(int i=0;i<P;i++) res.a[i][i]=1;    for(;b;b>>=1,a=a*a) if(b&1) res=res*a;    return res;}int main(){    freopen("loj2002.in","r",stdin);    freopen("loj2002.out","w",stdout);    scanf("%d%d%d",&n,&m,&P);    getP();    T=Matrix(); T.n=T.m=P;    for(int i=0;i<=P-1&&i<=m;i++)     for(int j=0;j<=P-1;j++) (T.a[j][(j+i)%P]+=(m-i)/P+(i>0))%=MOD;    Ans=Matrix(); Ans.n=1; Ans.m=P; Ans.a[0][0]=1;    Ans=Ans*Pow(T,n); ans=Ans.a[0][0];    for(int i=1;i<=p[0];i++)     for(int j=0;j<=P-1;j++) T.a[j][(j+p[i])%P]--;    Ans=Matrix(); Ans.n=1; Ans.m=P; Ans.a[0][0]=1;    Ans=Ans*Pow(T,n); ans=(ans-Ans.a[0][0])%MOD;    printf("%d\n",(ans+MOD)%MOD);    return 0;}