【SDOI2017】序列计数

题目描述

题目解析

呵呵，实际这是NOIP题你信吗？

如果不考虑包含质数的情况，因为p很小，我们不妨构建p*p的矩阵，然后对于(i,j)填有多少个数x使得(i+x)%p=j，最后矩阵快速幂一下就行了。

如果考虑质数我们便可以求出所有方案数再减去不包含质数的方案数。

代码

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<cstdlib>
#include<algorithm>
#include<queue>
#include<stack>
using namespace std;

#define MAXN
#define MAXM 20000000
#define INF 0x3f3f3f3f
typedef long long int LL;

template<class T>
void Read(T &x){
    x=;char c=getchar();bool flag=;
    while(c<'0'||'9'<c){if(c=='-')flag=;c=getchar();}
    while('0'<=c&&c<='9'){x=x*+c-'0';c=getchar();}
    if(flag)x=-x;
}

const int MOD = ;
struct Matrix{
    int n;
    LL data[][];
    Matrix(){}
    Matrix(int _n){n=_n;memset(data,,sizeof(data));}
    Matrix operator + (const Matrix &x)const{
        Matrix ret=Matrix(n);
        for(int i=;i<n;++i)
            for(int j=;j<n;++j)
                ret.data[i][j]=data[i][j]+x.data[i][j];
        return ret;
    }
    Matrix operator - (const Matrix &x)const{
        Matrix ret=Matrix(n);
        for(int i=;i<n;++i)
            for(int j=;j<n;++j)
                ret.data[i][j]=data[i][j]-x.data[i][j];
        return ret;
    }
    Matrix operator * (const Matrix &x)const{
        Matrix ret=Matrix(n);
        for(int i=;i<n;++i)
            for(int j=;j<n;++j)
                for(int k=;k<n;++k)
                    ret.data[i][j]+=data[i][k]*x.data[k][j];
        return ret;
    }
    Matrix operator % (const int mod)const{
        Matrix ret=Matrix(n);
        for(int i=;i<n;++i)
            for(int j=;j<n;++j)
                ret.data[i][j]=data[i][j]%mod;
        return ret;
    }
};

Matrix quick_power(Matrix a,int p,int mod){
    Matrix ret=Matrix(a.n);
    for(int i=;i<ret.n;++i)ret.data[i][i]=;

    while(p){
        if(p&)ret=ret*a%mod;
        a=a*a%mod;
        p>>=;
    }
    return ret;
}

bool isprime[MAXM+];
int prime[MAXM+],cnt;
void init(int ed){
    memset(isprime,,sizeof(isprime));
    cnt=;

    isprime[]=;
    for(int i=;i<=ed;++i){
        if(isprime[i])prime[++cnt]=i;
        for(int j=;j<=cnt&&prime[j]*i<=ed;++j){
            isprime[prime[j]*i]=;
            if(i%prime[j]==)break;
        }
    }
}
int change[];

int main(){
    //freopen("count.in","r",stdin);
    //freopen("count.out","w",stdout);

    int n,m,p;
    Read(n),Read(m),Read(p);
    init(m);

    Matrix origin=Matrix(p);
    memset(change,,sizeof(change));
    for(int i=;i<=m;++i)++change[(i%p)];
    for(int i=;i<p;++i)
        for(int j=;j<p;++j)origin.data[i][j]=change[((j-i)%p+p)%p]%MOD;
    Matrix ans1=quick_power(origin,n,MOD);

    origin=Matrix(p);
    memset(change,,sizeof(change));
    for(int i=;i<=m;++i)if(!isprime[i])++change[(i%p)];
    for(int i=;i<p;++i)
        for(int j=;j<p;++j)origin.data[i][j]=change[((j-i)%p+p)%p]%MOD;
    Matrix ans2=quick_power(origin,n,MOD);

    printf("%I64d\n",(ans1.data[][]-ans2.data[][]+MOD)%MOD);
}
/*
1 2 1
*/