好菜啊,现在才会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<=mf(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<jf(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<jf(pri[k])。
- 至此我们就能求出 ∑ i < = m f ( i ) \sum_{i<=m}f(i) ∑i<=mf(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<jf(pri[k])+∑k=jpri[k]2<=n∑q=1pri[k]q+1<=nS(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);
}