[SDOI2017]序列计数 （矩阵加速，小容斥）

题面

Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。

Alice还希望，这n个数中，至少有一个数是质数。

Alice想知道，有多少个序列满足她的要求。

输入格式

一行三个数，n,m,p。

对100\%的数据，1<= n <=1e9,1<= m <= 2e7,1<= p <= 100

输出格式

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

题解

我很少写这样的矩阵快速幂的题解

由于p的范围很小，所以p^3*log是可以过的

我们先求不加质数限制的方案数，再减去所有数都为合数的方案。

这道题，其实我们要知道[1,m]中的这些数在p的剩余系中的分布就可以计算，

因为对于不加质数限制的方案，就算值域在[1 + p,m + p]它的答案也是正确的。

设f[i][j]表示对于i个数，“ 其总和 % p = j ”的方案数，

我们发现，对于f，余数是要相加的，而方案数是要相乘的，

根据这个来推一下矩阵，i = 1的情况要直接算：

然后答案矩阵 ，这里的n就是题目中的n，

再来推一下所有数都是合数的方案

其实只要用欧拉筛筛一遍合数就行，然后：

设表示对于i个合数，“ 其总和 % p = j ”的方案数，

我们发现，对于f，余数是要相加的，而方案数是要相乘的，

根据这个来推一下矩阵，i = 1的情况要直接算：

然后答案矩阵 ，这里的n就是题目中的n，

最终的答案就是。

别忘了要全程取模

CODE

#include<cstdio>
#include<iostream>
#include<cstring>
#include<vector>
#include<queue>
#include<stack>
#include<cmath>
#include<algorithm>
#define LL long long
#define MAXN 200005
#pragma GCC optimize(2)
#pragma G++ optimize(3)
#define rg register
#define DB double
using namespace std;
inline int read() {
	int f = 1,x = 0;char s = getchar();
	while(s < '0' || s > '9') {if(s == '-') f = -1;s = getchar();}
	while(s >= '0' && s <= '9') {x = x * 10 + s - '0';s = getchar();}
	return x * f;
}
int mod = 20170408; 
int n,m,q,i,j,s,o,k,t;
struct matrix{
	int n,m;
	int a[105][105];
	matrix(){memset(a,0,sizeof(a));n = m = 1;}
}m1;
matrix operator * (matrix A,matrix B) {
	matrix C;
	C.n = A.n;C.m = B.m;
	for(rg int i = 1;i <= C.n;i ++) {
		for(rg int k = 1;k <= A.m;k ++) {
			for(rg int j = 1;j <= C.m;j ++) {
				C.a[i][j] = (C.a[i][j] + A.a[i][k] *1ll* B.a[k][j] % mod) % mod;
			}
		}
	}
	return C;
}
matrix qkpow(matrix a,LL b) {
	if(b == 0) return m1;
	if(b == 1) return a;
	matrix as = qkpow(a,b>>1);
	return as*as*qkpow(a,b&1);
}
int pri[5000005],cn;
bool f[20000005];
void sieve(int n) {
	f[1] = 1;cn = 0;
	for(rg int i = 2;i <= n;i ++) {
		if(!f[i]) pri[++cn] = i;
		for(rg int j = 1;j <= cn && i * pri[j] <= n;j ++) {
			f[i * pri[j]] = 1;
			if(i % pri[j] == 0) break;
		}
	}
}
int cnt[105],cnt2[105];
int main() {
	n = read();m = read();k = read();
	for(rg int i = 1;i <= k;i ++) m1.a[i][i] = 1;
	m1.n = m1.m = k;
	sieve(m);
	for(rg int i = 1;i <= m;i ++) {
		if(f[i]) cnt2[i % k] ++;
		cnt[i % k] ++;
	}
	matrix A,B,C,D;
	A.n = 1,A.m = B.n = B.m = k;
	for(rg int i = 1;i <= k;i ++) A.a[1][i] = cnt[i-1];
	for(rg int i = 1;i <= k;i ++) {
		for(rg int j = 1;j <= k;j ++) {
			B.a[j][i] = cnt[(k + i - j) % k];
		}
	}
	C = A * qkpow(B,n-1);
	A.n = 1,A.m = B.n = B.m = k;
	for(rg int i = 1;i <= k;i ++) A.a[1][i] = cnt2[i-1];
	for(rg int i = 1;i <= k;i ++) {
		for(rg int j = 1;j <= k;j ++) {
			B.a[j][i] = cnt2[(k + i - j) % k];
		}
	}
	D = A * qkpow(B,n-1);
	printf("%d\n",(C.a[1][1] + mod - D.a[1][1]) % mod);
	return 0;
}