min_25筛基础

min_25筛

作用及使用条件

可以得到积性函数的单点前缀和。时间复杂度为：

[O(frac{n^{3/4}}{log n}+n^{1-epsilon}) ]

由2018年某篇集训队论文证明。具体而言就是当(n)趋于无穷时，时间复杂度趋于(O(n))。(n)较小时时间复杂度为前者。

使用条件：

• 我们要找一些(f'(p) , f'(p))为可以快速求前缀和的完全积性函数。并且其在质数处的取值为原函数在质数处的取值相等或是一部分。并不要求其与原函数或其一部分相等，也就是说，构造一个，只要在质数处的点值与原来相同即可。
• (f(p^c))可以快速计算。

Step 1 求出质数处的函数值之和

我们设一个(g(n,j))表示小于等于(n)且最小质因数大于(p_j)的数的(f'(i))之和。

显然根据定义，(g(n,j))是无论如何包含不了1处的点值的。

对于(p_j^2 > n)的，显然有(S(n,j)=S(n,j-1))。

对于(p_j^2 leq n)的，然后我们有如下的式子：

[large g(n,j)=g(n,j-1)-f'(p_j)(g(frac{n}{p_j},j-1)-sum_{i=1}^{j-1}f'(p_j)) ]

(g(n,j-1)-g(n,j)=)小于等于(n)的最小质因数等于(p_j)的(f'(x))的值之和。

在(g(frac{n}{p_j},j-1))中，不符合条件的只有那些最小质因数(leq p_{j-1})的数，但是根据(S)的定义只有质数，所以减去那些质数。

对于边界情况，我们要求的答案就是(g(n,0))，这个是一个可以快速求前缀和的完全积性函数。

我们做了这么多,只为了求一样东西：(Sp(n)=g(n,|P|),n in )。

由于我们都只要所有的(g(i,|P|),P)为小于(n)的质数集大小（根据前面的分类，只要筛(O(sqrt{n}))的就行了）。所以我们使用滚动数组完成DP。

Step 2 组织答案

我们设一个(S(n,j))表示小于等于(n)且最小质因数大于(p_j)的数的(f(i))之和。

于是我们有

[large G(n,j)=sum_{j}(Sp_j(n)-sum_{i=1}^{j}f_j(i))+sum_{k > j}^{p_k leq sqrt{n}}sum_{e=1}^{p_k^eleq n}f(p_k^e)(S(frac{n}{p_k^e},k)+[e eq 1]) ]

为什么有一个([e eq 1])呢,因为我们在质数(>sqrt{n})时的点值是只有通过预处理得出。如果在这里要算的话，那么第一个(sum_k)的枚举次数直飚(pi(n))！而且由于没有合数最小质因子大于(sqrt{n})，所以后面的(S(frac{n}{p_k^e},k)=0)，只有质数本身有个点值。完全做不了。

但是有了预处理过后，我们可以只枚举到(O(sqrt{n}))

这就是为什么我们想尽办法算出了质数处的点值。

模板：P5325 【模板】Min_25筛

#include <bits/stdc++.h>
using namespace std;
const int S=1000005,mod=1000000007,inv3=333333336;
long long p[S],sp1[S],sp2[S],nn,n,w[S],g1[S],g2[S];
bool pp[S];
int ind1[S],ind2[S],tot=0;
void getp(int n)
{
	for (int i=2;i<=n;++i)
	{
		if (!pp[i])
		{
			p[++p[0]]=i;
			sp1[p[0]]=(sp1[p[0]-1]+i)%mod;
			sp2[p[0]]=(sp2[p[0]-1]+1ll*i*i)%mod;
		}
		for (int j=1;j<=p[0] && i*p[j]<=n;++j)
		{
			pp[i*p[j]]=true;
			if (i%p[j]==0) break;
		}
	}
}
long long func(long long x,int y)
{
	if (p[y]>=x) return 0;
	int k=x<=nn?ind1[x]:ind2[n/x];
	long long ret=(g2[k]-g1[k]+mod-(sp2[y]-sp1[y])+mod)%mod;
	for (int i=y+1;i<=p[0] && p[i]*p[i]<=x;++i)
	{
		long long pe=p[i],xx;
		for (int e=1;pe<=x;++e,pe*=p[i])
		{
			xx=pe%mod;
			ret=(ret+xx*(xx-1)%mod*(func(x/pe,i)+(e!=1)))%mod;
		}
	}
	return ret;
}
int main()
{
	scanf("%lld",&n);
	nn=sqrt(n);
	getp(nn);
	for (long long i=1;i<=n;)
	{
		long long r=n/(n/i);
		w[++tot]=n/i;
		long long o=w[tot]%mod;
		g1[tot]=o*(o+1)/2;
		g2[tot]=o*(o+1)/2%mod*(o<<1|1)%mod*inv3%mod;
		--g1[tot];--g2[tot];
		if (w[tot]<=nn) ind1[w[tot]]=tot;
		else ind2[r]=tot;
		i=r+1;
	}
	for (int i=1;i<=p[0];++i)
	{
		for (int j=1;j<=tot && p[i]*p[i]<=w[j];++j)
		{
			long long k=w[j]/p[i]<=nn?ind1[w[j]/p[i]]:ind2[n/(w[j]/p[i])];
			g1[j]-=p[i]*(g1[k]-sp1[i-1]+mod);
			g2[j]-=p[i]*p[i]%mod*(g2[k]-sp2[i-1]+mod);
			g1[j]%=mod;g2[j]%=mod;
			if (g1[j]<0) g1[j]+=mod;
			if (g2[j]<0) g2[j]+=mod;
		}
	}
	printf("%lld
",(func(n,0)+1)%mod);
	return 0;
}