SP20174 DIVCNT3 - Counting Divisors (cube)
Rigel
·
·
题解
前置知识：Min_25 筛。
默认读者已经掌握 Min_25 筛。
思路
设 f(n) = \sigma(n^{\textcolor{C74541}{3}})。显然 f 为积性函数。
考虑 Min_25 筛。
-
对于质数 p，f(p) = \textcolor{C74541}{4}，符合 f(p) 是可快速求值的完全积性函数的要求。
-
对于质数 p，f(p^c) = \textcolor{C74541}{3}c + 1，符合 f(p^c) 可快速求值的要求。
设 g(p) = 1，则 f(p) = \textcolor{C74541}{4} g(p)。
然后就可以用 Min_25 筛做了。
首先预处理出前 \lfloor \sqrt{n} \rfloor 个质数与它们对应的 G_{\text{prime}} 值。
由 g(p) = 1，有 G_{\text{prime}}(p_i) = G_{\text{prime}}(p_{i-1}) + 1。
然后用一次整除分块找到所有有效点 n'。
显然 G_0(n') = n' - 1。
对于所有 n'，从 G_0(n') 递推到 G_{\ell(n')}(n')，即得 G_\text{prime}(n')。
$$\begin{aligned}
G_k(n) &= G_{k-1}(n) - g(p_k)\left( G_{k-1}(n/p_k) - G_{\text{prime}}(p_{k-1}) \right)\\
&= G_{k-1}(n) - \left( G_{k-1}(n/p_k) - G_{\text{prime}}(p_{k-1}) \right).
\end{aligned}$$
最后计算 $F_1(n)$。
$F_k(n)$ 的递推式为：
$$\begin{aligned}
F_k(n) &= F_{\text{prime}}(n) - F_{\text{prime}}(p_{k-1}) + \sum\limits_{\substack{i \ge k\\p_i^2 \le n}}\sum\limits_{\substack{c \ge 1\\p_i^{c+1} \le n}}\left( f(p_i^{c+1}) + f(p_i^c)F_{i+1}(n/p_i^c) \right)\\
&= \textcolor{C74541}{4} \left( G_{\text{prime}}(n) - G_{\text{prime}}(p_{k-1}) \right) + \sum\limits_{\substack{i \ge k\\p_i^2 \le n}}\sum\limits_{\substack{c \ge 1\\p_i^{c+1} \le n}}\left( \textcolor{C74541}{3}(c+1) + 1 + (\textcolor{C74541}{3}c + 1)F_{i+1}(n/p_i^c) \right).
\end{aligned}$$
递归计算即可。
$F_1(n) + 1$ 即为所求的 $S_3(n)$。
时间复杂度为 $\Theta(n^{1 - \epsilon})$。
## Code
```cpp line-numbers
#include<bits/stdc++.h>
#define int long long
const int B=3.2e5+10;
using namespace std;
int T;
int n,b,ans;
int cnt,p[B];
bool vis[B];
int s[B];
int tot,w[B<<1];
int id1[B],id2[B];
int g[B<<1];
inline int read(){
int ret=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0'&&ch<='9')ret=ret*10+(ch&15),ch=getchar();
return ret*f;
}
int _sqrt(int x){
int ret=floor(sqrtl(x));
while((ret+1)*(ret+1)<=x)ret++;
return ret;
}
int get_id(int _n){
if(_n<=b)return id1[_n];
return id2[n/_n];
}
void _init(){
cnt=tot=0;
for(int i=1;i<=b;i++)p[i]=vis[i]=s[i]=id1[i]=id2[i]=0;
for(int i=1;i<=b*2;i++)w[i]=g[i]=0;
for(int i=2;i<=b;i++){
if(!vis[i]){
p[++cnt]=i;
s[cnt]=s[cnt-1]+1;
}
for(int j=1;j<=cnt&&i*p[j]<=b;j++){
vis[i*p[j]]=1;
if(!(i%p[j]))break;
}
}
for(int l=1,r;l<=n;l=r+1){
r=n/(n/l);
int _n=n/l;
w[++tot]=_n;
if(_n<=b)id1[_n]=tot;
else id2[n/_n]=tot;
g[tot]=_n-1;
}
for(int k=1;k<=cnt;k++){
for(int i=1;i<=tot&&p[k]*p[k]<=w[i];i++){
int id=get_id(w[i]/p[k]);
g[i]-=g[id]-s[k-1];
}
}
}
int F(int k,int _n){
if(_n<p[k])return 0;
int id=get_id(_n);
int ret=(g[id]-s[k-1])*4;
for(int i=k;i<=cnt&&p[i]*p[i]<=_n;i++){
for(int c=1,_p=p[i];_p*p[i]<=_n;c++,_p*=p[i]){
ret+=((c+1)*3+1)+(c*3+1)*F(i+1,_n/_p);
}
}
return ret;
}
void _solve(){
n=read();
b=_sqrt(n);
_init();
ans=F(1,n)+1;
printf("%lld\n",ans);
}
signed main(){
T=read();
while(T--)_solve();
return 0;
}
```
## 后记
双倍经验：**[DIVCNTK](https://www.luogu.com.cn/article/63nfpshh)**。
事实上 **[DIVCNT2](https://www.luogu.com.cn/article/87fdh0gf)** 也可以用 Min_25 筛做，但是跑得比较慢。