#矩阵乘法#洛谷 3702 [SDOI2017]序列计数

题目链接

---

分析

考虑容斥，用总方案减去全是合数的方案数，
可以发现 (n) 很大，(p) 很小，直接用矩阵乘法转移即可

---

代码

#include <cstdio>
#include <cctype>
#include <algorithm>
#define rr register
using namespace std;
const int mod=20170408; bool v[mod];
struct maix{int p[100][100];}A0,A1,ANS0,ANS1;
int n,m,k,prime[1300011],Cnt,c0[100],c1[100];
inline signed iut(){
	rr int ans=0; rr char c=getchar();
	while (!isdigit(c)) c=getchar();
	while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
	return ans;
}
inline signed mo(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline void mul(maix &A,maix B){
	rr maix C;
    for (rr int i=0;i<k;++i)
	for (rr int j=0;j<k;++j){
    	C.p[i][j]=0;
    	for (rr int o=0;o<k;++o)
    	    C.p[i][j]=mo(C.p[i][j],1ll*A.p[i][o]*B.p[o][j]%mod);
	}
	for (rr int i=0;i<k;++i)
	for (rr int j=0;j<k;++j)
	    A.p[i][j]=C.p[i][j];
}
signed main(){
	scanf("%d%d%d",&n,&m,&k),v[1]=1;
	for (rr int i=2;i<=m;++i){
		if (!v[i]) prime[++Cnt]=i;
		for (rr int j=1;j<=Cnt&&prime[j]<=m/i;++j){
			v[i*prime[j]]=1;
			if (i%prime[j]==0) break;
		}
	}
	for (rr int i=1;i<=m;++i) ++c0[i%k];
	for (rr int i=1;i<=m;++i) if (v[i]) ++c1[i%k];
	for (rr int i=0;i<k;++i)
	for (rr int j=0;j<k;++j)
	    A0.p[i][(i+j)%k]=c0[j],A1.p[i][(i+j)%k]=c1[j];
	ANS0.p[0][0]=ANS1.p[0][0]=1;
	for (;n;n>>=1,mul(A0,A0),mul(A1,A1))
	    if (n&1) mul(ANS0,A0),mul(ANS1,A1);
	return !printf("%d",mo(ANS0.p[0][0],mod-ANS1.p[0][0]));
}