56pts求调qwq
查看原帖
56pts求调qwq
99643
陈刀仔楼主2021/7/20 21:15

调不动了救救孩子。。。加了点注释,应该勉强能看吧qwq

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int maxn=2e5+50;
int n,K,mod,sqr,pri[maxn],prik[maxn],tot,prek[maxn],g[maxn];
//g 辅助数组,带上了滚动 
int w[maxn],ind1[maxn],ind2[maxn],cnt;//离散化记录有限的取值个数 
int fy[maxn],inv[maxn];//拉格朗日插值的点值,和预处理逆元
bool notp[maxn];//质数标记 
inline int ksm(int a,int p){
	int res=1;
	a=(a%mod+mod)%mod;
	while(p){
		if(p&1)res=res*a%mod;
		p>>=1;a=a*a%mod;
	}
	return res;
}
void init(){//线性筛,预处理拉格朗日插值要用的点值 
	notp[0]=notp[1]=1;
	for(int i=2;i*i<=n;i++){
		if(!notp[i]){
			pri[++tot]=i;prik[tot]=ksm(i,K);
			prek[tot]=prek[tot-1]+prik[tot];prek[tot]%=mod;
		}
		for(int j=1;j<=tot&&pri[j]*i<=sqr;j++){
			notp[pri[j]*i]=1;
			if(i%pri[j]==0)break;
		}
	}
	fy[1]=1;
	for(int i=2;i<=K+2;i++){
		fy[i]=fy[i-1]+ksm(i,K);
		fy[i]%=mod;
	}
	for(int i=0;i<=2*K+4;i++)inv[i]=ksm(i-K-2,mod-2);
}
inline int getsum(int x){ //O(K^2)拉格朗日插值 
	int ans=0;
	for(int i=1;i<=K+2;i++){
		int tmp=fy[i];
		for(int j=1;j<=K+2;j++){
			if(j==i)continue;
			tmp=tmp*(x-j)%mod*inv[i-j+K+2]%mod;
		}
		ans=(ans+tmp)%mod;
		if(ans<0)ans+=mod;
	}
	return ans;
}
void getg(){
	for(int i=1,las;i<=n;i=las+1){
		las=n/(n/i);
		w[++cnt]=n/i;
		g[cnt]=(getsum(w[cnt])-1+mod)%mod;
		if(w[cnt]<=sqr)ind1[w[cnt]]=cnt;
		else ind2[las]=cnt;
	}
	for(int i=1;i<=tot;i++){
		for(int j=1;j<=cnt&&pri[i]*pri[i]<=w[j];j++){
			int x=w[j];
			int t=(x/pri[i]<=sqr?ind1[x/pri[i]]:ind2[n/(x/pri[i])]);
			g[j]-=prik[i]*(g[t]-prek[i-1]+mod)%mod;
			if(g[j]<0)g[j]+=mod;
		}
	}
}
signed main(){
//	freopen("prime.in","r",stdin);
//	freopen("prime.out","w",stdout);
	cin>>n>>K>>mod;
	sqr=sqrt(n);
	init();
	getg();
	int ans=0;
	for(int i=1;i<=sqr;i++){
		int tmp=n/i,k;
		if(tmp<=sqr)k=ind1[tmp];
		else k=ind2[n/tmp];
		ans+=i*i%mod*g[k]%mod;
		ans%=mod;
//		cout<<"** "<<tmp<<' '<<g[k]<<endl;
	}
	printf("%lld\n",ans);
	return 0;
}


2021/7/20 21:15
加载中...