Min_25筛学习小计

好菜啊，现在才会Min_25

简单介绍

• 由Min_25在他的博客中介绍的做法，又称Min_25筛。
• 用于求积性函数 f ( n ) f(n) f(n)的前缀和，其中要求 f ( p ) f(p) f(p)可以表示成多项式，并且 f ( p k ) f(p^k) f(pk)可以快速算出。
• 能够在 O ( n 0.75 l o g ) O(\frac{n^{0.75}}{log}) O(logn0.75​)的时间内算出。

基本思路

• 使用 D P DP DP先求出所有素数的 f ( p ) f(p) f(p)的前缀和，然后再通过枚举最小的质因子计算出合数的前缀和。

求素数的前缀和

• 下面的所有除法都为整除。
• 我们之后要求对所有的 m = n i m=\frac{n}{i} m=in​求 ∑ i ⊆ p r i m e , i < = m f ( i ) \sum_{i\subseteq prime,i<=m}f(i) ∑i⊆prime,i<=m​f(i)
• 为了方便计算我们将 f ( i ) f(i) f(i)拆成几个完全积性函数的和。使得这些积性函数在素数 p p p的位置的和与 f ( i ) f(i) f(i)一致。
• 接下来考虑一个完全积性函数 f ( i ) f(i) f(i)在素数位置的和怎么求。
• 那么用容斥，考虑用所有的和减去合数的和。
• g ( n , j ) = ∑ i = 1 n [ i ⊆ p r i m e ∣ m i n p ( i ) > p r i [ j ] ] f ( i ) g(n,j)=\sum_{i=1}^n[i\subseteq prime|minp(i)>pri[j]]f(i) g(n,j)=∑i=1n​[i⊆prime∣minp(i)>pri[j]]f(i)

• 从 g ( n , j ) g(n,j) g(n,j)求 g ( n , j + 1 ) g(n,j+1) g(n,j+1)的过程相当于是一个埃氏筛法（即找出 n \sqrt n n

  ​下的质数，将它的倍数筛掉，相当于现在筛到的数的最小因子为当前质数）。

• 那么如果 p r i [ j ] 2 > n pri[j]^2>n pri[j]2>n，则 p r i [ j ] pri[j] pri[j]不能继续筛了，所以：

  g ( n , j ) = g ( n , j − 1 ) g(n,j)=g(n,j-1) g(n,j)=g(n,j−1)。

• 否则 p r i [ j ] pri[j] pri[j]可以筛，即为：

  g ( n , j ) = g ( n , j − 1 ) − p r i [ j ] ( g ( n / p r i [ j ] , j − 1 ) − ∑ k < j f ( p r i [ k ] ) ) g(n,j)=g(n,j-1)-pri[j](g(n/pri[j],j-1)-\sum_{k<j}f(pri[k])) g(n,j)=g(n,j−1)−pri[j](g(n/pri[j],j−1)−∑k<j​f(pri[k]))

• 因为 g g g里面有小于 p r i [ j ] pri[j] pri[j]的质数，但是它们乘上 p r i [ j ] pri[j] pri[j]的最小因子并不是 p r i [ j ] pri[j] pri[j]，不满足当前筛的数的范围——最小质因子为 p r i [ j ] pri[j] pri[j]，所以减去。
• 最后求得 g ( n , ∣ P ∣ ) g(n,|P|) g(n,∣P∣)即为素数的和了。

亿点点细节

• 一般使用 p k p^k pk作为完全积性函数，初值 g ( n , 0 ) g(n,0) g(n,0)就是一个自然数幂求和。
• 这个东西用滚动数组求。
• 由于要求所有 g ( n i , ∣ P ∣ ) g(\frac{n}{i},|P|) g(in​,∣P∣)，查找的时候可以设一个 i d 1 , i d 2 id1,id2 id1,id2，如果 i < = s q r t ( n ) i<=sqrt(n) i<=sqrt(n)就直接存在该位置，否则存在 n / i n/i n/i的位置，可以证明一一对应。然后离散化去储存和转移。
• 还需要在线筛时预先求出 ∑ k < j f ( p r i [ k ] ) \sum_{k<j}f(pri[k]) ∑k<j​f(pri[k])。
• 至此我们就能求出 ∑ i < = m f ( i ) \sum_{i<=m}f(i) ∑i<=m​f(i)，这里的 f ( i ) f(i) f(i)是题意中的函数。
• 所有计算都不计算1的贡献，因为筛的时候比较方便

求所有数的前缀和

• 设 S ( n , j ) = ∑ i = 1 n [ m i n p ( i ) > = p r i [ j ] ] f ( i ) S(n,j)=\sum_{i=1}^n[minp(i)>=pri[j]]f(i) S(n,j)=∑i=1n​[minp(i)>=pri[j]]f(i)
• S ( n , j ) = g ( n , ∣ P ∣ ) − ∑ k < j f ( p r i [ k ] ) + ∑ k = j p r i [ k ] 2 < = n ∑ q = 1 p r i [ k ] q + 1 < = n S ( n / p r i [ k ] q , k + 1 ) + f ( p r i [ k ] q + 1 ) S(n,j)=g(n,|P|)-\sum_{k<j}f(pri[k])+\sum_{k=j}^{pri[k]^2<=n}\sum_{q=1}^{pri[k]^{q+1}<=n}S(n/pri[k]^q,k+1)+f(pri[k]^{q+1}) S(n,j)=g(n,∣P∣)−∑k<j​f(pri[k])+∑k=jpri[k]2<=n​∑q=1pri[k]q+1<=n​S(n/pri[k]q,k+1)+f(pri[k]q+1)
• 前面的质数的和，后面合数的和也是枚举最小的质因子的次幂去筛它，注意最后一项算不到要补上。
• 出口即 n < p r i [ j ] , n = 1 n<pri[j],n=1 n<pri[j],n=1返回 0 0 0。
• 因为整除的数是有限且不重复的，所以 S ( n , j ) S(n,j) S(n,j)并不需要像杜教筛一样记忆化。
• 听说时间是 O ( n 0.75 l o g ) O(\frac{n^{0.75}}{log}) O(logn0.75​)的，最快能过 1 0 12 10^{12} 1012。

• 相较杜教筛普适性更高了，但是速度也更慢了。

  luoguP3235

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define maxn 500005
#define ll long long 
#define mo 1000000007
using namespace std;

ll n,sqr,i,j,k,m,inv6;
ll id1[maxn],id2[maxn],w[maxn];
ll bz[maxn],tot,pri[maxn],sg1[maxn],sg2[maxn];
ll s[maxn],g1[maxn],g2[maxn];

ll ksm(ll x,ll y){
	ll s=1;
	for(;y;y/=2,x=x*x%mo) if (y&1)
		s=s*x%mo;
	return s;
}

void getpri(){
	for(i=2;i<=sqr;i++){
		if (!bz[i]) {
			pri[++tot]=i;
			sg1[tot]=(sg1[tot-1]+i)%mo;
			sg2[tot]=(sg2[tot-1]+1ll*i*i)%mo;
		}
		for(j=1;j<=tot&&i*pri[j]<=sqr;j++){
			bz[i*pri[j]]=1;
			if (i%pri[j]==0) break;
		}
	}
}

void prepare(){
	getpri();
	for(i=1;i<=n;i=j+1){
		j=n/(n/i),w[++m]=n/i;
		if (w[m]<=sqr) id1[w[m]]=m; else
			id2[n/w[m]]=m;
		ll tmp=w[m]%mo; 
		g1[m]=(tmp+1)*tmp/2%mo-1;
		g2[m]=tmp*(tmp+1)%mo*(tmp*2+1)%mo*inv6%mo-1;
	}
	for(j=1;j<=tot;j++){
		for(i=1;i<=m&&pri[j]*pri[j]<=w[i];i++){
			int k=(w[i]/pri[j]<=sqr)?id1[w[i]/pri[j]]:id2[n/(w[i]/pri[j])];	
			(g1[i]-=pri[j]*(g1[k]-sg1[j-1]))%=mo;
			(g2[i]-=pri[j]*pri[j]%mo*(g2[k]-sg2[j-1]))%=mo;
		}
	}
}

ll S(ll x,int j){
	if (x==1||pri[j]>x) return 0;
	int k=(x<=sqr)?id1[x]:id2[n/x];
	ll sum=g2[k]-g1[k]-sg2[j-1]+sg1[j-1];
	for(k=j;k<=tot&&pri[k]*pri[k]<=x;k++){
		ll p=pri[k],q=p%mo;
		while (p*pri[k]<=x){
			sum+=S(x/p,k+1)*(q*(q-1)%mo)%mo;
			p=p*pri[k],q=q*pri[k]%mo,sum+=q*(q-1)%mo;
		}
	}
	return sum%mo;
}

int main(){
	scanf("%lld",&n),sqr=sqrt(n),inv6=ksm(6,mo-2);
	prepare();
	printf("%lld",S(n,1)+1);
}