1. 程式人生 > 實用技巧 >luogu P5325 【模板】Min_25篩

luogu P5325 【模板】Min_25篩

程式碼:

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>

using namespace std;

typedef long long LL;
const int M=1000000007,N=1000009;
int a[N],cnt,inv6,ind1[N],ind2[N],tot;
LL sp1[N],sp2[N],Sqr,n,p[N],g2[N],g1[N],w[N];

void prework(int fuck)
{
	for (LL i=2;i<=fuck;i++)
	{
		if(!a[i])
			a[i]=i,p[++cnt]=i,sp1[cnt]=(sp1[cnt-1]+i)%M,sp2[cnt]=(sp2[cnt-1]+i*i%M)%M;
		for (int j=1;j<=cnt;j++)
		{
			if(p[j]>a[i]||p[j]>fuck/i)
				break;
			a[p[j]*i]=p[j];
		}
	}
}

int ksm(int a,int b)
{
	int res=1;
	while(b)
	{
		if(b&1)
			res=1ll*res*a%M;
		b>>=1,a=1ll*a*a%M;
	}
	return res;
}

void init()
{
	scanf("%lld",&n);
	inv6=ksm(6,M-2);
	Sqr=(LL)sqrt(n);
	prework(Sqr);
}

LL S(LL x,int y)
{
	if(p[y]>=x) return 0;
	int k=x<=Sqr?ind1[x]:ind2[n/x];
	LL res=(g2[k]-g1[k]-sp2[y]+sp1[y])%M;
	for (int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++)
	{
		LL pe=p[i];
		for (int e=1;pe<=x;e++,pe*=p[i])
		{
			LL xx=pe%M;
			res=(res+xx*(xx-1)%M*(S(x/pe,i)+(e!=1))%M)%M;
		}
	}
	return res;
}

void work()
{
	for (LL i=1;i<=n;i++)
	{
		LL j=n/(n/i);
		w[++tot]=n/i;
		g1[tot]=w[tot]%M;
		g2[tot]=(g1[tot]*(g1[tot]+1)%M*(g1[tot]+g1[tot]+1)%M*inv6%M-1)%M;
		g1[tot]=(g1[tot]+1)*g1[tot]/2%M-1;
		if(n/i<=Sqr) ind1[n/i]=tot;
		else ind2[n/(n/i)]=tot;
		i=j;
	}
	for (int i=1;i<=cnt;i++)
		for (int j=1;j<=tot&&p[i]*p[i]<=w[j];j++)
		{
			int k=w[j]/p[i]<=Sqr?ind1[w[j]/p[i]]:ind2[n/(w[j]/p[i])];
			g1[j]=(g1[j]-p[i]*(g1[k]-sp1[i-1])%M)%M;
			g2[j]=(g2[j]-p[i]*p[i]%M*(g2[k]-sp2[i-1])%M)%M;
		}
	printf("%lld\n",(S(n,0)+1+M)%M);
}

int main()
{
	init();
	work();
	return 0;
}