求调
查看原帖
求调
315132
GCY_XZT楼主2023/5/5 14:37
#include<bits/stdc++.h>
#define int long long

using namespace std;

const int mod=20170408;
bool isprime[30001000];
int prime[1000010],cnt;
int cha[2001],sss[2001];
int ch1[1001],ch2[1001];

struct Matrix{
	int a[210][210];
	int n,m;
	void init0(){
		memset(a,0,sizeof(a));
	}
	void init1(){
		memset(a,0,sizeof(a));
		for(int i=0;i<n;i++){
			a[i][i]=1;
		}
	}
	Matrix operator *(const Matrix x)const{
		Matrix ans;
		ans.init0();
		for(int i=0;i<n;i++){
			for(int j=0;j<x.m;j++){
				for(int k=0;k<m;k++){
					ans.a[i][j]=(ans.a[i][j]+a[i][k]*x.a[k][j])%mod;
				}
			}
		}
		ans.n=n;
		ans.m=x.m;
		return ans;
	}
	void print(){
		for(int i=0;i<n;i++){
			for(int j=0;j<m;j++){
				cout<<a[i][j]<<" ";
			}
			cout<<endl;
		}
	}
}mat1,mat2,chu1,chu2;

Matrix qpow(Matrix a,int b){
	Matrix sum;
	sum.n=a.n;
	sum.m=a.n;
	sum.init1();
	while(b){
		if(b&1)sum=sum*a;
		a=a*a;
		b>>=1;
	}
	return sum;
} 

//void getprime(int n){
	
//}

signed main(){
	int n,m,p;
	cin>>n>>m>>p;
	for(int i=1;i<=m;i++){
		isprime[i]=1;
	}
	isprime[1]=0;
	for(int i=2;i<=m;i++){
		if(isprime[i]){
			prime[++cnt]=i;
		}
		for(int j=1;j<=cnt&&i*prime[j]<=m;j++){
			isprime[i*prime[j]]=0;
			if(i%prime[j]==0)break;
		}
	}
	mat1.n=mat1.m=mat2.n=mat2.m=p;
	mat1.init0();
	mat2.init0();
	for(int j=1;j<=m;j++){
		ch1[j%p]=(ch1[j%p]+1);
	}
	for(int j=1;j<=m;j++){
		if(!isprime[j])ch2[j%p]=(ch2[j%p]+1);
	}
	for(int i=0;i<=p-1;i++){
		for(int j=0;j*p+i<=m;j++){
			if(j*p+i==0)continue;
			sss[i]=(sss[i]+1);
		}
	}
	for(int i=0;i<mat1.n;i++){
		for(int j=0;j<mat1.n;j++){
			if(i<j){
				int x=j-i;
				mat1.a[i][j]=sss[x]%mod;
			}
			else if(i==j){
				int x=0;
				mat1.a[i][j]=sss[0]%mod;
			}
			else{
				int x=p-(i-j);
				mat1.a[i][j]=sss[x]%mod;
			}
		}
	}
	for(int i=0;i<=p-1;i++){
		for(int j=0;j*p+i<=m;j++){
			if(j*p+i==0)continue;
			if(!isprime[j*p+i])cha[i]=(cha[i]+1);
		}
	}
	for(int i=0;i<mat2.n;i++){
		for(int j=0;j<mat2.n;j++){
			if(i<j){
				int x=j-i;
				mat2.a[i][j]=cha[x]%mod;
			}
			else if(i==j){
				int x=0;
				mat2.a[i][j]=cha[0]%mod;
			}
			else{
				int x=p-(i-j);
				mat2.a[i][j]=cha[x]%mod;
			}
		}
	}
	//mat1.print();
//	mat2.print();
	mat1=qpow(mat1,n-1);
	mat2=qpow(mat2,n-1);
	chu1.n=chu2.n=1;
	chu1.m=chu2.m=p;
	chu1.init0();
	chu2.init0();
	for(int j=0;j<chu1.m;j++){
		chu1.a[0][j]=ch1[j]%mod;
	}
	for(int j=0;j<chu2.m;j++){
		chu2.a[0][j]=ch2[j]%mod;
	}	
	//chu1.print();
//	cout<<endl<<endl;
	//chu2.print();
	chu1=chu1*mat1;
	chu2=chu2*mat2;
	cout<<(((chu1.a[0][0]-chu2.a[0][0])%mod)+mod)%mod;
	return 0;
}
2023/5/5 14:37
加载中...