luogu 3702 [SDOI2017]序列计数 矩阵乘法+容斥

现在看来这道题真的不难啊~

正着求不好求，那就反着求：答案=总-全不是质数 

这里有一个细节要特判：1不是质数，所以在算全不是质数的时候要特判1 

code: 

#include <bits/stdc++.h>           
#define N 104 
#define M 20000002 
#define mod 20170408    
#define ll long long  
#define setIO(s) freopen(s".in","r",stdin) 
using namespace std;  
int n,m,p,tot; 
int cnt[N];
bool vis[M];
int prime[200000],cnt2[N];  
struct matrix 
{   
    ll a[100][100];               
    void re() { memset(a,0,sizeof(a));}  
    void I() { re(); for(int i=0;i<p;++i) a[i][i]=1ll;  }
    ll*operator[](int x) { return a[x]; }                       
}A,B;   
matrix operator*(matrix a,matrix b) 
{    
    matrix c;c.re();           
    int i,j,k; 
    for(i=0;i<p;++i) 
    {
        for(j=0;j<p;++j)
            for(k=0;k<p;++k) 
                c[i][j]=(c[i][j]+a[i][k]*b[k][j]%mod+mod)%mod;   
    }   
    return c;      
}
matrix operator^(matrix a,ll k) 
{   
    matrix tmp;                    
    for(tmp.I();k;k>>=1,a=a*a) if(k&1) tmp=tmp*a;   
    return tmp;      
}
int main() 
{ 
    // setIO("input");           
    int i,j; 
    scanf("%d%d%d",&n,&m,&p);        
    for(i=1;i<=m;++i) cnt[i%p]++;        
    for(i=2;i<=m;++i) 
    {   
        if(!vis[i]) prime[++tot]=i;  
        for(j=1;j<=tot&&prime[j]*i<=m;++j) 
        {
            vis[i*prime[j]]=1;  
            if(i%prime[j]==0) break;       
        }         
    }    
    for(i=1;i<=m;++i) if(vis[i]) ++cnt2[i%p]; 
    ++cnt2[1%p];                            
    A.re();          
    B.re();    
    for(i=0;i<p;++i) 
    {
        for(j=0;j<p;++j) 
        {
            int pp=(j-i+p)%p;        
            A[i][j]=cnt[pp];                            
            B[i][j]=cnt2[pp];                                 
        }
    }                      
    A=A^n;       
    B=B^n;          
    // printf("%lld %lld\n",A[0][0],B[0][0]);   
    printf("%lld\n",(A[0][0]-B[0][0]+mod)%mod);        
    return 0;
}