BZOJ4818 LOJ2002 SDOI2017 序列计数 【矩阵快速幂优化DP】*

BZOJ4818 LOJ2002 SDOI2017 序列计数

---

Description

Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。
Alice还希望，这n个数中，至少有一个数是质数。
Alice想知道，有多少个序列满足她的要求。

Input

一行三个数，n,m,p。
1<=n<=10^9,1<=m<=2×10^7,1<=p<=100

Output

一行一个数，满足Alice的要求的序列数量，答案对20170408取模。

Sample Input

3 5 3

Sample Output

33

---

#include<bits/stdc++.h>
using namespace std;
#define M 20000010
#define N 110
#define Mod 20170408
#define LL long long
int n,m,p,pri[M];
LL s[N],tot=0;
bool vis[M];
struct Matrix{
    LL t[N][N];
    Matrix(){memset(t,0,sizeof(t));}
    Matrix operator *(const Matrix b)const{
        Matrix c;
        for(int i=0;i<p;i++)
            for(int k=0;k<p;k++)
                for(int j=0;j<p;j++){
                    c.t[i][j]+=t[i][k]*b.t[k][j]%Mod;
                    if(c.t[i][j]>Mod)c.t[i][j]-=Mod;
                }
        return c;
    }
};
void init(){
    vis[1]=1;   
    for(int i=2;i<M;i++){
        if(!vis[i])pri[++tot]=i;
        for(int j=1;j<=tot&&i*pri[j]<M;j++){
            vis[i*pri[j]]=1;
            if(i%pri[j]==0)break;
        }
    }
}
Matrix fast_pow(Matrix a,int b){
    Matrix ans;
    for(int i=0;i<p;i++)ans.t[i][i]=1;
    while(b){
        if(b&1)ans=ans*a;
        a=a*a;
        b>>=1;
    }
    return ans;
}
int solve(){
    Matrix res;
    for(int i=0;i<p;i++)
        for(int j=0;j<p;j++)res.t[i][(i+j)%p]=s[j];
    res=fast_pow(res,n);
    return res.t[0][0]; 
}
int main(){
    init();
    scanf("%d%d%d",&n,&m,&p);
    for(int i=1;i<=m;i++)++s[i%p];
    LL ans=solve();
    for(int i=1;i<=m;i++)if(!vis[i])--s[i%p];
    ans-=solve();
    printf("%lld",(ans+Mod)%Mod);
//  system("pause");
    return 0;
}