Min_25筛学习小计

好菜啊，现在才会Min_25

简单介绍

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

基本思路

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

求素数的前缀和

• 下面的所有除法都为整除。
• 我们之后要求对所有的$m=frac{n}{i}$求$sum_{isubseteq prime,i<=m}f(i)$
• 为了方便计算我们将$f(i)$拆成几个完全积性函数的和。使得这些积性函数在素数$p$的位置的和与$f(i)$一致。
• 接下来考虑一个完全积性函数$f(i)$在素数位置的和怎么求。
• 那么用容斥，考虑用所有的和减去合数的和。
• $g(n,j)=sum_{i=1}^n[isubseteq prime|minp(i)>pri[j]]f(i)$
• 从$g(n,j)$求$g(n,j+1)$的过程相当于是一个埃氏筛法（即找出$sqrt n$下的质数，将它的倍数筛掉，相当于现在筛到的数的最小因子为当前质数）。
• 那么如果$pri[j]^2>n$，则$pri[j]$不能继续筛了，所以：
  $g(n,j)=g(n,j-1)$。
• 否则$pri[j]$可以筛，即为：
  $g(n,j)=g(n,j-1)-pri[j](g(n/pri[j],j-1)-sum_{k<j}f(pri[k]))$
• 因为$g$里面有小于$pri[j]$的质数，但是它们乘上$pri[j]$的最小因子并不是$pri[j]$，不满足当前筛的数的范围——最小质因子为$pri[j]$，所以减去。
• 最后求得$g(n,|P|)$即为素数的和了。

亿点点细节

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

求所有数的前缀和

• 设$S(n,j)=sum_{i=1}^n[minp(i)>=pri[j]]f(i)$
• $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})$
• 前面的质数的和，后面合数的和也是枚举最小的质因子的次幂去筛它，注意最后一项算不到要补上。
• 出口即$n<pri[j],n=1$返回$0$。
• 因为整除的数是有限且不重复的，所以$S(n,j)$并不需要像杜教筛一样记忆化。
• 听说时间是$O(frac{n^{0.75}}{log})$的，最快能过$10^{12}$。
• 相较杜教筛普适性更高了，但是速度也更慢了。
  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);
}