BZOJ4818 [SDOI2017] 序列计数 【矩阵快速幂】

题目分析：

一个很显然的同类项合并。注意到p的大小最大为100，考虑把模p意义下相同的求出来最后所有的减去没有质数的做矩阵快速幂即可。

代码：

#include<bits/stdc++.h>
using namespace std;

const int maxn = 20000005;

const int mod = 20170408;

int n,m,q;

int base[105][2];

int flag[maxn],prime[maxn/10],num;

int mat[105][105],Res[105][105],po[105][105];

void init(){
    flag[1] = 1;
    for(int i=2;i<=m;i++){
    if(!flag[i]) prime[++num] = i;
    for(int j=1;j<=num&&i*prime[j]<=m;j++){
        flag[i*prime[j]] = 1;
        if(i%prime[j] == 0) break;
    }
    }
    for(int i=1;i<=m;i++){
    if(flag[i]) base[i%q][1]++;
    base[i%q][0]++;
    }
}

void BuildMatrix(int now){
    memset(mat,0,sizeof(mat));
    for(int i=0;i<q;i++){
    for(int j=0;j<q;j++){
        int nw = j-i;if(nw < 0) nw += q;
        mat[j][nw] += base[i][now];
    }
    }
}

void fast_pow(int now){
    if(now == 1){
    for(int i=0;i<q;i++) for(int j=0;j<q;j++) Res[i][j]=mat[i][j];
    }else{
    fast_pow(now/2);
    memset(po,0,sizeof(po));
    for(int k=0;k<q;k++)
        for(int i=0;i<q;i++)
        for(int j=0;j<q;j++){
            po[i][j] += (1ll*Res[i][k]*Res[k][j])%mod;
            po[i][j] %= mod;
        }
    for(int i=0;i<q;i++)for(int j=0;j<q;j++)Res[i][j]=po[i][j];
    if(now & 1){
        memset(po,0,sizeof(po));
        for(int k=0;k<q;k++)
        for(int i=0;i<q;i++)
            for(int j=0;j<q;j++){
            po[i][j] += (1ll*Res[i][k]*mat[k][j])%mod;
            po[i][j] %= mod;
            }
        for(int i=0;i<q;i++)for(int j=0;j<q;j++)Res[i][j]=po[i][j];
    }
    }
}

void work(){
    BuildMatrix(0);
    fast_pow(n);
    int ans = Res[0][0];
    BuildMatrix(1);
    fast_pow(n);
    ans -= Res[0][0]; 
    if(ans < 0) ans += mod;
    printf("%d",ans);
}

int main(){
    scanf("%d%d%d",&n,&m,&q);
    init();
    work();
    return 0;
}