【bzoj4818】[Sdoi2017]序列计数 矩阵乘法

原文地址：http://www.cnblogs.com/GXZlegend/p/6825132.html

Alice想要得到一个长度为n的序列，序列中的数都是不超过m的正整数，而且这n个数的和是p的倍数。Alice还希望，这n个数中，至少有一个数是质数。Alice想知道，有多少个序列满足她的要求。

输入

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

1<=n<=10^9,1<=m<=2×10^7,1<=p<=100

输出

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

样例输入

3 5 3

样例输出

33

矩阵乘法

至少有1个质数的方案数=总方案数-不含质数的方案数

可以先在O(m)的时间内把1~m的质数筛出来。

然后考虑：从mod p=0，到mod p=0，可以由mod p=a和mod p=p-a两个阶段组成。

可以设f[i][j]表示从mod p=i到mod p=j的方案数，不难看出这是一个矩阵，而且自乘m次就能得到答案。

所以只要处理出2种f即可。

对于每个数可以更新所有的f[i]（0≤i<n），不过这样会TLE

其实这样做没有必要，因为f都是循环出现的，只需要求f[0]即可，再据此推其余的f。

代码中，我将mod p=0当作了p来处理，其实看作0也是一样的。

#include <cstdio>
#include <cstring>
#include <algorithm>
#define mod 20170408
using namespace std;
typedef long long ll;
bool notprime[20000010];
int pri[5000010] , top;
struct matrix
{
	int n , m;
	ll num[110][110];
	matrix()
	{
		n = m = 0;
		memset(num , 0 , sizeof(num));
	}
	matrix operator*(const matrix a)const
	{
		matrix ans;
		ans.n = n , ans.m = a.m;
		int i , j , k;
		for(i = 1 ; i <= ans.n ; i ++ )
			for(j = 1 ; j <= ans.m ; j ++ )
				for(k = 1 ; k <= m ; k ++ )
					ans.num[i][j] = (ans.num[i][j] + num[i][k] * a.num[k][j]) % mod;
		return ans;
	}
}A , B;
matrix pow(matrix x , int y)
{
	matrix ans;
	int i;
	ans.n = x.n , ans.m = x.m;
	for(i = 1 ; i <= ans.n ; i ++ )
			ans.num[i][i] = 1;
	while(y)
	{
		if(y & 1) ans = ans * x;
		x = x * x , y >>= 1;
	}
	return ans;
}
int main()
{
	int n , m , p , i , j;
	scanf("%d%d%d" , &n , &m , &p);
	notprime[1] = 1;
	for(i = 2 ; i <= m ; i ++ )
	{
		if(!notprime[i]) pri[++top] = i;
		for(j = 1 ; j <= top && i * pri[j] <= m ; j ++ )
		{
			notprime[i * pri[j]] = 1;
			if(i % pri[j] == 0) break;
		}
	}
	A.n = A.m = B.n = B.m = p;
	for(i = 1 ; i <= m ; i ++ )
	{
		A.num[p][(i - 1) % p + 1] ++ ;
		if(notprime[i]) B.num[p][(i - 1) % p + 1] ++ ;
	}
	for(i = p - 1 ; i >= 1 ; i -- )
		for(j = 1 ; j <= p ; j ++ )
			A.num[i][j] = A.num[i + 1][j % p + 1] , B.num[i][j] = B.num[i + 1][j % p + 1];
	printf("%lld\n" , (pow(A , n).num[p][p] - pow(B , n).num[p][p] + mod) % mod);
	return 0;
}