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

最新推荐文章于 2020-09-28 17:12:28 发布

Lynstery

最新推荐文章于 2020-09-28 17:12:28 发布

阅读量450

点赞数

分类专栏： --------------数学相关------------ 矩阵乘法

本文链接：https://blog.csdn.net/chhnz/article/details/78693590

版权

102 篇文章 0 订阅

订阅专栏

矩阵乘法

5 篇文章 0 订阅

订阅专栏

注意到 $p$ 很小，容易想到矩乘。 $f_i \to f_{(i+j)\text{%}P}$ 。
关于那个质数的限制，只需要 $1$ 到 $m$ 的方案减去 $1$ 到 $m$ 扣去所有质数的方案即可。
构造 $1$ 到 $m$ 转移矩阵 $O(P^2)$ 搞， $\sum_{k=1}^m [k\text{%}P=i]\ =\ \lfloor \frac{m-i}k \rfloor+[i>0]$ 。
扣去质数直接暴力 $O(m/\ln m*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;
}