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));
}
}