bzoj 4818: [Sdoi2017]序列计数（DP+矩阵快速幂）

4818: [Sdoi2017]序列计数

Time Limit: 30 Sec   Memory Limit: 128 MB
Submit: 769   Solved: 463
[ Submit][ Status][ Discuss]

Description

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

Input

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

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

Output

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

Sample Input

3 5 3

Sample Output

33

以下假设p=4，n和m未知：

dp[n][k]表示前n个数之和对p取模为k的总情况数，x[k]表示对p取模为k的数的个数（当然<=m）

那么可以列出矩阵：

中间那个矩阵的大小是p²的，所以用矩阵快速幂的话复杂度为O(p^3logn)

可是Alice还希望这n个数中有质数，这样的话还要再算一个不出现质数的dp然后相减

可以啊

#pragma comment(linker, "/STACK:102400000,102400000")
#include<stdio.h>
#include<string.h>
#define mod 20170408
#define LL long long
LL p, x[115];
int flag[20000010];
typedef struct
{
	LL i;
	LL a[105][105];
	void init()
	{
		memset(a, 0, sizeof(a));
	}
	void unit()
	{
		memset(a, 0, sizeof(a));
		for(i=1;i<=p;i++)
			a[i][i] = 1;
	}
}Matrix;
Matrix Xpow(Matrix p1, Matrix p2)
{
	Matrix pe;
	LL i, j, k;
	pe.init();
	for(i=1;i<=p;i++)
	{
		for(j=1;j<=p;j++)
		{
			for(k=1;k<=p;k++)
				pe.a[i][j] = (pe.a[i][j]+p1.a[i][k]*p2.a[k][j])%mod;
		}
	}
	return pe;
}
Matrix Powto(Matrix p, LL k)
{
	Matrix E;
	E.unit();
	while(k)
	{
		if(k%2)
			E = Xpow(E, p);
		p = Xpow(p, p);
		k /= 2;
	}
	return E;
}
int main(void)
{
	Matrix Jz;
	LL n, m, i, j, ans;
	flag[0] = flag[1] = 1;
	for(i=2;i<=20000000;i++)
	{
		if(flag[i])
			continue;
		for(j=i*i;j<=20000000;j+=i)
			flag[j] = 1;
	}
	scanf("%lld%lld%lld", &n, &m, &p);
	Jz.init();
	for(i=0;i<=p-1;i++)
	{
		for(j=i;j<=m;j+=p)
			x[i]++;
	}
	x[0]--;
	for(i=1;i<=p;i++)
	{
		for(j=1;j<=p;j++)
			Jz.a[i][j] = x[((i-j)+p)%p];
	}
	Jz = Powto(Jz, n);
	ans = Jz.a[1][1];
	Jz.init();
	memset(x, 0, sizeof(x));
	for(i=0;i<=p-1;i++)
	{
		for(j=i;j<=m;j+=p)
			x[i] += flag[j];
	}
	x[0]--;
	for(i=1;i<=p;i++)
	{
		for(j=1;j<=p;j++)
			Jz.a[i][j] = x[((i-j)+p)%p];
	}
	Jz = Powto(Jz, n);
	ans = (ans-Jz.a[1][1]+mod)%mod;
	printf("%lld\n", ans);
	return 0;
}