Skip to main content

Min_25 筛

参考资料​

例题​

洛谷 P5325 【模板】Min_25 筛​

洛谷 P5325 【模板】Min_25 筛

定义积性函数 f(x)f(x)，且 f(pk)=pk(pk−1)f(p^k)=p^k(p^k-1)（pp 是一个质数），求：

∑i=1nf(i)\sum_{i=1}^n f(i)

对 109+710^9+7 取模。

Code (1)
#include <bits/stdc++.h>
using namespace std;

using ll=long long;
const int N=1000005;
const int mod=1e9+7;
const int inv2=500000004;
const int inv6=166666668;
ll n,w[N<<1];
int sq,m,cnt,pri[N],sp1[N],sp2[N];
int g1[N<<1],g2[N<<1],id1[N],id2[N];
bool vis[N];
void sieve()
{
	vis[0]=vis[1]=1;
	cnt=0;
	for(int i=2;i<N;i++)
	{
		if(!vis[i])
		{
			pri[++cnt]=i;
			sp1[cnt]=(sp1[cnt-1]+i)%mod;
			sp2[cnt]=(sp2[cnt-1]+1ll*i*i)%mod;
		}
		for(int j=1;j<=cnt;j++)
		{
			if(i*pri[j]>=N)break;
			vis[i*pri[j]]=1;
			if(i%pri[j]==0)break;
		}
	}
}
int S(ll x,int y)
{
	if(x<=1||pri[y]>x)return 0;
	int k=(x<=sq)?id1[x]:id2[n/x];
	int res=(g2[k]-g1[k]+mod)%mod;
	res=(res-(sp2[y-1]-sp1[y-1]+mod)%mod+mod)%mod;
	for(int i=y;i<=cnt&&1ll*pri[i]*pri[i]<=x;i++)
	{
		ll pe=pri[i];
		for(int e=1;pe*pri[i]<=x;e++,pe*=pri[i])
		{
			int v1=pe%mod*((pe-1)%mod)%mod;
			int v2=pe*pri[i]%mod*((pe*pri[i]-1)%mod)%mod;
			res=(res+1ll*v1*S(x/pe,i+1)+v2)%mod;
		}
	}
	return res;
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	sieve();
	cin>>n;
	sq=sqrt(n);
	for(ll l=1,r;l<=n;l=r+1)
	{
		r=n/(n/l);
		w[++m]=n/l;
		ll v=w[m]%mod;
		g1[m]=(v*(v+1)%mod*inv2%mod-1+mod)%mod;
		g2[m]=(v*(v+1)%mod*(2*v+1)%mod*inv6%mod-1+mod)%mod;
		if(w[m]<=sq)id1[w[m]]=m;
		else id2[n/w[m]]=m;
	}
	for(int j=1;j<=cnt;j++)
	{
		ll p2=1ll*pri[j]*pri[j];
		for(int i=1;i<=m&&w[i]>=p2;i++)
		{
			ll v=w[i]/pri[j];
			int k=(v<=sq)?id1[v]:id2[n/v];
			g1[i]=(g1[i]-1ll*pri[j]*(g1[k]-sp1[j-1]+mod)%mod+mod)%mod;
			g2[i]=(g2[i]-1ll*pri[j]*pri[j]%mod*(g2[k]-sp2[j-1]+mod)%mod+mod)%mod;
		}
	}
	cout<<(S(n,1)+1)%mod<<'\n';
	return 0;
}

洛谷 P7884 【模板】Meissel-Lehmer​

洛谷 P7884 【模板】Meissel-Lehmer

给定整数 nn，求出 π(n)\pi(n) 的值。

π(n)\pi(n) 表示 1∼n1 \sim n 的整数中质数的个数。

Code (1)
#include <bits/stdc++.h>
using namespace std;

using ll=long long;
const double eps=1e-6;
const int N=85000010;
const int M=30000005;
int c[N],lim,sq,tot;
ll g[M],w[M];
double inv[M];
bool vis[N];
void add(int u){vis[u]=1;while(u<=lim){c[u]++;u+=u&-u;}}
ll query(int u){ll res=u;while(u){res-=c[u];u-=u&-u;}return res;}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	ll n;
	cin>>n;
	sq=sqrt(n)+1;
	lim=max((int)pow(n/log2(n),0.66),sq);
	lim=max(lim,10000);
	lim=min((ll)lim,n);
	for(int i=1;i<=sq;i++)w[i]=i-1,inv[i]=1.0/i;
	for(tot=1;1ll*lim*tot<n;tot++)g[tot]=n*inv[tot]+eps-1;
	tot--;
	add(1);
	for(int i=2;1ll*i*i<=n;i++)
	{
		if(vis[i])continue;
		ll x0=w[i-1],t=n/i,r=1ll*i*i;
		int tl=min((ll)tot,n/r);
		int tl2=min((ll)tl,ll(n/(1.0*sq*i)));
		int tl3=min(tl2,tot/i);
		for(int j=1;j<=tl3;j++)g[j]-=g[j*i]-x0;
		for(int j=tl3+1;j<=tl2;j++)g[j]-=query(t*inv[j]+eps)-x0;
		for(int j=tl2+1;j<=tl;j++)g[j]-=w[int(t*inv[j]+eps)]-x0;
		for(int j=sq;j>=r;j--)w[j]-=w[int(j*inv[i]+eps)]-x0;
		if(1ll*i*i<=lim)for(int j=i*i;j<=lim;j+=i)if(!vis[j])add(j);
	}
	if(!tot)g[1]=query(n);
	cout<<g[1]<<'\n';
	return 0;
}