min_25 筛学习小记

min_25筛

由 dalao min_25 发明的筛子,据说时间复杂度是极其优秀的 $ O(\frac {n^{\frac 3 4}} {\log n}) $,常数还小。

1. 质数 $ k $ 次方前缀和(基础)

求 $ \sum_{p \leq n}p^k $

我们考虑一个 $ \rm DP $ 的思路:设 $ g(n,j) $ 为:

\[\sum_{i=1}^n[(\sum_{t=1}^j[p_t|i])=0] i^k \]

其实就是不大于 $ n $ 的,且不含有 $ p_1 $ ~ $ p_j $ 中任意一个质因子的数的 $ k $ 次方之和。

我们考虑从 $ g(n,j-1) $ 转移上来。

发现 $ g(n,j) $ 相对于 $ g(n,j-1) $ 减少的就是含有 $ p_k $ 质因子的数的 $ k $ 次方之和,于是我们得到了转移方程:

\[g(n,k) = g(n,k-1)-p_j^kg(\lfloor \frac n {p_j} \rfloor,j) \]

但是,这个方程似乎有点儿问题?

我们把 $ p_j $ 以内的质数算重了,所以要加上。

最终的方程就是:

\[g(n,k) = g(n,k-1)-p_j^k(g(\lfloor \frac n {p_j} \rfloor,j)-g(p_{j-1},j-1) \]

不过如果我们要维护 $ O(n) $ 个 $ g $ 复杂度开销是很大的,所以考虑 玄学优化 去掉不重要的东西。

首先我们注意到,在 $ p_j > \sqrt n $ 是,对答案的贡献只有 $ p_j^k $,所以我们可以先只考虑 $ p_j \leq \sqrt n $ 的情况。

$ g(p_{j-1},j-1) $ 显然可以在 $ O(\sqrt n) $ 的复杂度内线性筛预处理出来,就只需要考虑 $ g(\lfloor \frac n {p_j} \rfloor,j) $

然后根据整除分块的结论,我们只需要处理 $ O(\sqrt n) $ 个 $ g $ 就行了。。。

然后复杂度就降下来了,我也不会证,不过看上去似乎是 $ O(\sqrt n\log\log n) $ 的(?)

但是 $ n $ 太大需要离散化。。。

所以我们依照杜教筛的离散化,把大于 $ \sqrt n $ 的映射到 $ 1 $ ~ $ \sqrt n $。

给一下参考代码(这里只算了质数的 $ 1 $ 次方和 $ 2 $ 次方之和):

const int St=1e6+5,mod=1e9+7;
inline int Add(const int&a,const int&b){
    return a+b>=mod?a+b-mod:a+b;
}
inline int Get(const ll&k){
	return k<=sqr?id1[k]:id2[n/k];
}
inline void seive(const int&n){
    int i,j,x;zhi[1]=1;
    for(i=2;i<=n;++i){
        if(!zhi[i]){
            pri[++top]=i;
            sum1[top]=Add(sum1[top-1],i);
            sum2[top]=Add(sum2[top-1],1ll*i*i%mod);
        }
        for(j=1;j<=top&&(x=i*pri[j])<=n;++j){
            zhi[x]=1;
            if(!(i%pri[j]))break;
        }
    }
}
ll min_25(const ll&n){
    int i,j;
    seive(sqr=sqrt(n));
    for(ll L=1,R;L<=n;L=R+1){
        R=n/(n/L);
        w[++tot]=n/L;
        g1[tot]=w[tot]%mod;
        g2[tot]=g1[tot]*(g1[tot]+1)/2%mod*(g1[tot]*2+1)%mod*inv3%mod-1;
        g1[tot]=g1[tot]*(g1[tot]+1)/2%mod-1;
        if(n/L<=sqr)id1[n/L]=tot;
        else id2[n/(n/L)]=tot;
    }
    for(i=1;i<=top;++i){
        for(j=1;j<=tot&&pri[i]*pri[i]<=w[j];++j){
            ll k=Get(w[j]/pri[i]);
            g1[j]-=pri[i]*(g1[k]-sum1[i-1]+mod)%mod;
            g2[j]-=pri[i]*pri[i]%mod*(g2[k]-sum2[i-1]+mod)%mod;
            if(g1[j]<=0)g1[j]+=mod;
            if(g2[j]<=0)g2[j]+=mod;
        }
    }
}

在实现的时候将 $ g $ 的第二维滚掉了qwq

2.min_25筛

没错前面的都是铺垫(雾)

让我们看看我们要求的问题:

\[\sum_{i=1}^nf(i) \]

其中 $ f $ 是一个低阶多项式,其在质数乘方处的取值能够快速算出。

将答案分成两部分考虑:

\[\sum_{p \leq n}f(p)+\sum_{p^e \leq n,p \leq \sqrt n}f(p^e)(\sum_{i=1 \And minp > p}^{\lfloor \frac n {p^e} \rfloor}f(i)) \]

其中 $ minp $ 是 $ i $ 最小的质因子。

按照类似上面 $ \rm DP $ 的套路,设 $ S(n,x) $ 是不大于 $ n $ 的,且不含有 $ p_1 $ ~ $ p_x $ 中任何一个因子的 $ f $ 之和。

显然答案是 $ S(n,0) $。

再将答案分成两部分,一部分是大于 $ p_x $ 的质数的 $ f $ 之和,另一部分就是剩下的。

第一部分,也就是 $ g(n) - g(p_x) $。

对照着上面的式子,我们可以得到以下的转移方程

\[S(n,x) = g(n) - g(p_x) +\sum_{p_k^e \leq n,k > x}f(p_k^e)(S(\lfloor \frac n {p_k^e} \rfloor,k)+[e!=1]) \]

然后还有一些细节,就是关于 $ 1 $ 的问题。。。

由于是照着 $ \rm wucstdio $ dalao 的题解学的,所以我也没有算 $ 1 $,而是在结束时加上。

例题

板子

\[p^k(p^k-1) = p^{2k} - p^k \]

在质数处就是 $ p^2 - p $

只需要筛出质数的 $ 1 $ 次方前缀和和 $ 2 $ 次方前缀和即可。

不要问我为什么上面的代码筛的也是 $ 1 $ 次和 $ 2 $ 次前缀和,我懒总行了吧

#include<cstdio>
#include<cmath>
typedef long long ll;
const int St=1e6+5,mod=1e9+7,inv3=333333336;
int sqr,top,tot,zhi[St],id1[St],id2[St];
ll n,w[St],g1[St],g2[St],pri[St],sum1[St],sum2[St];
inline int Add(const int&a,const int&b){
    return a+b>=mod?a+b-mod:a+b;
}
inline int Get(const ll&k){
	return k<=sqr?id1[k]:id2[n/k];
}
inline void seive(const int&n){
    int i,j,x;zhi[1]=1;
    for(i=2;i<=n;++i){
        if(!zhi[i]){
            pri[++top]=i;
            sum1[top]=Add(sum1[top-1],i);
            sum2[top]=Add(sum2[top-1],1ll*i*i%mod);
        }
        for(j=1;j<=top&&(x=i*pri[j])<=n;++j){
            zhi[x]=1;
            if(!(i%pri[j]))break;
        }
    }
}
ll S(const ll&n,const int&k){
    if(pri[k]>=n)return 0;
    ll x=Get(n),ans=Add((g2[x]-g1[x]+sum1[k]-sum2[k])%mod,mod);
    for(int i=k+1;i<=top&&pri[i]*pri[i]<=n;++i){
        ll p=pri[i];
        for(int e=1;p<=n;++e,p*=pri[i]){
            ll id=p%mod;
            ans=Add(ans,id*(id-1)%mod*((e!=1)+S(n/p,i))%mod);
        }
    }
    return ans;
}
ll min_25(const ll&n){
    int i,j;
    seive(sqr=sqrt(n));
    for(ll L=1,R;L<=n;L=R+1){
        R=n/(n/L);
        w[++tot]=n/L;
        g1[tot]=w[tot]%mod;
        g2[tot]=g1[tot]*(g1[tot]+1)/2%mod*(g1[tot]*2+1)%mod*inv3%mod-1;
        g1[tot]=g1[tot]*(g1[tot]+1)/2%mod-1;
        if(n/L<=sqr)id1[n/L]=tot;
        else id2[n/(n/L)]=tot;
    }
    for(i=1;i<=top;++i){
        for(j=1;j<=tot&&pri[i]*pri[i]<=w[j];++j){
            ll k=Get(w[j]/pri[i]);
            g1[j]-=pri[i]*(g1[k]-sum1[i-1]+mod)%mod;
            g2[j]-=pri[i]*pri[i]%mod*(g2[k]-sum2[i-1]+mod)%mod;
            if(g1[j]<=0)g1[j]+=mod;
            if(g2[j]<=0)g2[j]+=mod;
        }
    }
    return Add(S(n,0),1);
}
signed main(){
    scanf("%lld",&n);
    printf("%d",min_25(n));
}

然后是一道很经典的三倍经验题DIVCNTK,做出这道题后 DIVCNT2 和 DIVCNT3 就是三倍经验。
设积性函数 $ f(p^e) = d((pe)k) = d(p^{ek}) $

在质数处的取值就是 $ k+1 $

#include<cstdio>
#include<cmath>
typedef unsigned long long ull;
const int M=2e6+5;
ull n,m,lim,sqr,tot,top,w[M],g[M],pri[M],zhi[M],id1[M],id2[M];
inline int Get(const ull&k){
	return k<=sqr?id1[k]:id2[n/k];
}
inline void sieve(const int&n){
	int i,j,x;zhi[1]=1;
	for(i=2;i<=n;++i){
		if(!zhi[i])pri[++top]=i;
		for(j=1;j<=top&&(x=i*pri[j])<=n;++j){
			zhi[x]=1;
			if(!(i%pri[j]))break;
		}
	}
}
ull S(const ull&n,const int&k){
	if(pri[k]>=n)return 0;
	ull x=Get(n),ans=g[x]-k*(m+1);
	for(ull i=k+1;i<=lim&&pri[i]*pri[i]<=n;++i){
		for(ull e=1,p=pri[i];p<=n;++e,p*=pri[i]){
			ans+=(e*m+1)*(S(n/p,i)+(e!=1));
		}
	}
	return ans;
}
ull min_25(const ull&n){
	int i,j;
	sqr=sqrt(n);tot=0;lim=1;
	while(pri[lim]*lim[pri]<=n)++lim;--lim;
	for(ull L=1,R;L<=n;L=R+1){
		R=n/(n/L);
		w[++tot]=n/L;
		g[tot]=w[tot]-1;
		if(n/L<=sqr)id1[n/L]=tot;
		else id2[n/(n/L)]=tot;
	}
	for(i=1;i<=lim;++i){
		for(j=1;j<=tot&&pri[i]*pri[i]<=w[j];++j){
			ull k=Get(w[j]/pri[i]);
			g[j]-=g[k]-i+1;
		}
	}
	for(i=1;i<=tot;++i)g[i]*=(m+1);
	return S(n,0)+1;
}
signed main(){
	int i,T;
	sieve(2e6);
	scanf("%d",&T);
	for(i=1;i<=T;++i){
		scanf("%llu%llu",&n,&m);
		printf("%llu\n",min_25(n));
	}
}
上一篇:用python创建文件夹并移动文件


下一篇:Pandas数据探索分析,分享两个神器