题面

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;

}

[SDOI2017]序列计数 （矩阵加速，小容斥）的相关教程结束。