[Sdoi2017]序列计数

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

思路

非常容易想到 O(n*m*p)的转移 不再赘述

发现不需要枚举每一个m值 而是可以在一开始就对m分类  (mod p的剩余系) 分类之后的结果就用c数组记录 这样就可以从i向j转移即  s[(i+j)%p]+=s[i]*c[j]; O(n*p*p); 还是会T

那么我们可以构造一个 s[i][j]的(p*p)的一个矩阵 表示从i向j转移的变化 那么可以定义一个初始向量 f 表示第一个位置的选择方案

明显的是 f[i]=c[i]; 即分类之后的选择方案  那么 f*s就是第二个位置选择之后每一个剩余系的方案数 现在我们可以对s做矩阵快速幂(n-1)次 最后再用f去乘s 之后的 f[0][0]就是不考虑质数的答案ans1

然后在把质数从c中抛出去  再做一遍 表示没有质数的方案数ans2  那么ans1-ans2就是答案

代码

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
typedef long long ll;
const ll MOD=20170408;
const int maxn=200+5;
const int N=20000000+5;
int prime[N],tot;
bool visit[N];
int n,m,p;
ll s[maxn][maxn],c[maxn];
struct matrix
{
	int n,m,a[maxn][maxn];//n行 m列 
	void clear()
	{
		n=m=0; memset(a,0,sizeof(a));
	}
	matrix operator * (const matrix &p) const
	{
		matrix tmp; tmp.clear();
		for(int i=0; i<n; i++)
		{
			for(int j=0; j<p.m; j++)
			{
				for(int k=0; k<m; k++)
				{
					tmp.a[i][j]+=(1LL*a[i][k]*p.a[k][j])%MOD;
					tmp.a[i][j]%=MOD;
				}
			}
		}
		tmp.n=n; tmp.m=p.m;
		return tmp;
	}
};
void getprime()
{
	for(int i=2; i<=m; i++)
	{
		if(!visit[i])
		{
			prime[++tot]=i;
		}
		for(int j=1; j<=tot,prime[j]*i<=m; j++)
		{
			visit[prime[j]*i]=true;
			if(i%prime[j]==0) break;
		}
	}
}
matrix quick(int n)
{
	matrix e; e.clear(); e.clear();
	e.n=e.m=p; matrix ans; 
	ans.n=ans.m=p;
	for(int i=0; i<p; i++)
	for(int j=0; j<p; j++)
	{
		e.a[i][j]=s[i][j];
	}
	for(int i=0; i<p; i++) ans.a[i][i]=1;
	while(n)
	{
		if(n&1) ans=ans*e;
		e=e*e;
	/*	for(int i=0;i<p;i++)
		{
			for(int j=0;j<p;j++)
			{
				cout<<e.a[i][j]<<' ';
			}
			cout<<endl;
		}
		puts("**********************************");*/
		n>>=1;
	}
	return ans;
}
int main()
{
	freopen("test.in","r",stdin);
	freopen("test.out","w",stdout);	
	scanf("%d %d %d",&n,&m,&p);
	for(int i=0; i<p; i++)
		c[i]+=m/p;
	for(int j=1; j<=m%p; j++)
		c[j]++;
	for(int i=0; i<p; i++)
	{
		for(int j=0; j<p; j++)
		{
			s[i][j]=c[(j-i+p)%p]%MOD;
		}
	}
	matrix f1;
	f1.n=1,f1.m=p;
	for(int i=0; i<p; i++)
	{
		f1.a[0][i]=c[i];
	}
	matrix ans1=f1*quick(n-1);
	getprime();
	for(int i=1; i<=tot; i++)
	{
//		cout<<prime[i]<<endl;
		c[prime[i]%p]--;
	}
	for(int i=0; i<p; i++)
	{
		c[i]=((c[i]%MOD)+MOD)%MOD;
	}
	matrix f2;
	f2.n=1,f2.m=p;
	for(int i=0; i<p; i++)
	{
		f2.a[0][i]=c[i];
	}
	for(int i=0; i<p; i++)
	{
		for(int j=0; j<p; j++)
		{
			s[i][j]=c[(j-i+p)%p]%MOD;
//			cout<<s[i][j]<<' ';
		}
//		cout<<endl;
	}
	matrix ans2=f2*quick(n-1);
	cout<<(ans1.a[0][0]-ans2.a[0][0]+MOD)%MOD<<endl;
	return 0;
}