正題

題目鏈接:https://www.luogu.com.cn/problem/P5325

---

題目大意

定義一個積性函數滿足$f(p^k)=p^k(p^k?1)$
求${\sum }_{i=1}^{n}f(i)$

---

解題思路

首先我們可以把$f(p^k)$是質數的情況拆成一個$2$階的多項式$f(x)=x^2?x$。

然后就是$\text{Min25}$篩了。我們先需要求出質數有關的答案，考慮把這個多項式的多個階分開來求。

有關質數的求法，我們定義$dp$數組
$$g_k(n,j)=\sum _{i=1}^n[i\in Priorminp(i)>p_j]i^k$$（上面$Pri$表示質數集，$minp(x)$表示$x$的最小質因子，$p_j$表示$\le \sqrt{n}$的第$j$個質因子，$k$表示多項式的某個階，后面為了方便$g_k$寫作$g$，后文中的$g(n,j)$中的$n$都與輸入的$n$無關，只是一個變量）

具體的說就是如果$i$是質數或者最小質因子超過$p_j$，那當某個$p_j$最接近但小于$\sqrt{n}$時，$g(n,j)$就是我們需要的答案了。

因為一個合數$x$滿足$minp(x)\le \sqrt{x}$，所以$p_j$的數量級很小，可以利用這個性質。

遞推這個數組時對于每次$j$往上加相當于多加上了一些限制條件，也就是$g(n,j)$相對于$g(n,j?1)$我們需要減去一些多余的答案。

這些多余的答案就是最小質因數是$p_j$的數，這些數就是$p_j^k?(g(\frac{n}{p_j},j?1)?g(p_{j?1},j?1))$。
$g(\frac{n}{p_j},j?1)?g(p_{j?1},j?1))$表示枚舉一個乘上$p_j$之后不會超過$n$且最小質因子超過$p_j$的數，這些數乘上$p_j$后就是不重不漏的表示最小質因子是$p_j$的數，因為$y=x^k$是一個完全積性函數，所以可以直接乘上一個$p_j^k$。

現在我們就有遞推式
$$g(n,j)=g(n,j?1)+p_j^k(g(\frac{n}{p_j},j?1)?g(p_{j?1},j?1))$$
快速處理$g(n,0)$然后遞推，可以注意到$g(p_{j?1},j?1)$就是前$j?1$個質數的$k$次方和，因為$p_j$數量級只有$\sqrt{n}$所以可以直接預處理，因為后面還要用就定義$s{p}_i={\sum }_{i=1}^{n}f(p_i)$

還有發現無論如何每個$g(x,j)$的取值都只與$?\frac{n}{x}?$有關，所以我們對于$g$數組的每一個$j$都只有$\sqrt{n}$個連續的范圍內的同一取值。所以我們可以直接整除分塊來大大縮減數量級的空間和時間。

給出一個$x$如何快速得到$g(x,j)$的儲存地址？一個比較巧妙的方法是如果$x\le \sqrt{n}$那么$ind{1}_x$表示$x$的儲存位置，否則用$ind{2}_{\frac{n}{x}}$表示$x$的儲存地址。

然后$j$那個維度我們只需要用到$j=cnt$（$cnt$表示$\sqrt{n}$以內的質數數量）值，所以我們直接滾動來遞推就好了。

現在我們知道了所有的$g(i,cnt)$，這道題的話具體的值就是$g_2(i,cnt)?g_1(i,cnt)$就可以表示前綴質數$f$的函數和了。下面簡寫$g(i)=g_2(i,cnt)?g_1(i,cnt)$。

然后就是要求答案了，同樣使用$dp$的技巧，設
$$S(n,j)=\sum _{i=2}^n[minp(i)>p_j]f(i)$$
然后對于$S(n,j)$的答案我們可以分為質數和非質數兩部分，質數部分顯然就是$g(n)?s{p}_n$（不記得$sp$的去看前面定義）。
然后合數部分我們考慮遞歸下去處理就是${\sum }_{k>j}^{p_k^e\le n}f(p_k^e)(S(\frac{n}{p_k^e},k)+[e=1])$。$p_k^e$是枚舉質因子就不多說了，那個$[e=1]$是因為$p_k$已經作為質數被計算過了，但是后面的$p_k^e$還沒有。

現在就有$S(n,j)$的遞推式了
$$S(n,j)=g(n)?s{p}_n+\sum _{k>j}^{p_k^e\le n}f(p_k^e)(S(\frac{n}{p_k^e},k)+[e=1])$$

遞歸下去做就好了，說是不知道因為啥定理所以是不用記憶化的。

一個細節就是$f(1)$最好不要統計進去，最好加上就好了

據說時間復雜度是$O(\frac{n^{\frac{3}{4}}}{\log n})$的，反正這題我沒開$O2$的話$1e10$只跑了七百多毫秒。

---

code

#include<cstdio> #include<cstring> #include<algorithm> #include<cmath> #define ll long long using namespace std; const ll N=1e6+10,P=1e9+7; ll n,T,cnt,tot,pri[N],sp1[N],sp2[N],ind1[N],ind2[N],g1[N],g2[N],w[N]; bool v[N]; void init(ll n){for(ll i=2;i<=n;i++){if(!v[i]){pri[++cnt]=i;sp1[cnt]=(sp1[cnt-1]+i)%P;sp2[cnt]=(sp2[cnt-1]+i*i)%P;}for(ll j=1;j<=cnt&&i*pri[j]<=n;j++){v[i*pri[j]]=1;if(i%pri[j]==0)break;}}return; } ll S(ll x,ll y){if(pri[y]>=x)return 0;ll pos=(x>T)?ind2[n/x]:ind1[x];ll ans=g2[pos]-g1[pos]-(sp2[y]-sp1[y]);ans=(ans%P+P)%P;for(ll i=y+1;i<=cnt&&pri[i]*pri[i]<=x;i++){for(ll j=1,p=pri[i];p<=x;j++,p=p*pri[i]){ll w=p%P;(ans+=w*(w-1)%P*(S(x/p,i)+(j!=1))%P)%=P;}}return ans; } signed main() {scanf("%lld",&n);T=sqrt(n);init(T);ll inv2=(P+1)/2,inv6=(P+1)/6;for(ll l=1,r;l<=n;l=r+1){ll x=n/l;r=n/(n/l);w[++tot]=n/l;x%=P;g1[tot]=x*(x+1)%P*inv2%P-1;g2[tot]=x*(x+1)%P*(2*x+1)%P*inv6%P-1;if(n/l<=T)ind1[n/l]=tot;else ind2[n/(n/l)]=tot;}for(ll i=1;i<=cnt;i++)for(ll j=1;j<=tot&&pri[i]*pri[i]<=w[j];j++){ll k=(w[j]/pri[i]);k=(k>T)?ind2[n/k]:ind1[k];(g1[j]+=P-pri[i]*(g1[k]-sp1[i-1])%P)%=P;(g2[j]+=P-pri[i]*pri[i]%P*(g2[k]-sp2[i-1])%P)%=P;}printf("%lld\n",(S(n,0)+1)%P);return 0; }

總結

以上是生活随笔為你收集整理的P5325-[模板]Min_25筛的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。