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