SPOJ DIVCNT2 - Counting Divisors (square)
题意：求
\[\sum_{i=1}^n\sigma_0(i^2)
\]
\]
好棒啊！
带着平方没法做，考虑用其他函数表示\(\sigma_0(i^2)\)，把平方消去。
\(\sigma_0(n) = (1*1)(n) = \sum_{d\mid n}1\)
我们考虑那些\(n^2\)有而\(n\)没有的因子，\(n=\prod p_i^{a_i}\)，那么这些因子里一定有\(p_i^c:c>a_i\)。
对于因子\(d\)，他的每个质因子都可以指数加上\(a_i\)成为\(n^2\)独有的因子，贡献为\(2^{\omega(d)}\)，其中\(\omega(n)\)表示不同的质因子个数。
\(2^{\omega(n)} = \sum_{d\mid n}\mu^2(d)\)
\[\sigma_0(n^2) = \sum_{d\mid n} 2^{\omega(d)} = \sum_{d\mid n} \sum_{e\mid d} \mu^2(e) = ((\mu^2 * 1) * 1) (n)
\]
\]
我们就是要求\((\mu * 1) * 1 = \mu * (1*1) = \mu * \sigma_0\)的前缀和
\[\begin{align}
ans &= \sum_{i=1}^n \sum_{d\mid i} \mu^2(d) \cdot \sigma_0(\frac{i}{d}) \\
&= \sum_{i=1}^n \mu^2(i) \sum_{j=1}^{\lfloor \frac{n}{i} \rfloor} \sigma_0(j)
\end{align}
\]
ans &= \sum_{i=1}^n \sum_{d\mid i} \mu^2(d) \cdot \sigma_0(\frac{i}{d}) \\
&= \sum_{i=1}^n \mu^2(i) \sum_{j=1}^{\lfloor \frac{n}{i} \rfloor} \sigma_0(j)
\end{align}
\]
不用杜教筛，我们也可以求。
我们只要得到\(\mu^2\)和\(\sigma_0\)的前缀和就可以整除分块了。
\(\sum_{i=1}^n \mu^2(i) = \sum_{i=1}^{\sqrt{n}}\mu(i)\lfloor \frac{n}{i^2} \rfloor\) 就是无平方因子数的个数呀
\(\sum_{i=1}^n\sigma_0(i) = \sum_{i=1}^n \lfloor \frac{n}{i} \rfloor\) 也可以整除分块
同时我们使用线性筛预处理前\(O(n^{\frac{2}{3}})\)的前缀和，剩下的部分用上面两个式子\(O(\sqrt{n})\)计算
复杂度分析和杜教筛类似，
\[T(n) =O(k + \sum_{i=1}^{\frac{n}{k}}\sqrt{\frac{n}{i}})=O(k + \frac{n}{\sqrt{k}})
\]
\]
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=1e8+5;
int U=1e8;
inline ll read(){
char c=getchar(); ll x=0,f=1;
while(c<'0' || c>'9') {if(c=='-')f=-1; c=getchar();}
while(c>='0' && c<='9') {x=x*10+c-'0'; c=getchar();}
return x*f;
}
bool notp[N]; int p[N/10], mu[N], lp[N], mu2[N]; ll si[N];
void sieve(int n) {
mu[1]=1; si[1]=1; mu2[1]=1;
for(int i=2; i<=n; i++) {
if(!notp[i]) p[++p[0]] = i, mu[i] = -1, si[i] = lp[i] = 2;
for(int j=1; j <= p[0] && i*p[j] <= n; j++) {
int t = i*p[j];
notp[t] = 1;
if(i%p[j] == 0) {
mu[t] = 0;
lp[t] = lp[i] + 1;
si[t] = si[i] / lp[i] * lp[t];
break;
}
mu[t] = -mu[i];
lp[t] = 2;
si[t] = si[i] * 2;
}
mu2[i] = mu2[i-1] + mu[i] * mu[i];
si[i] += si[i-1];
}
}
inline ll sum_mu2(ll n) {
if(n <= U) return mu2[n];
int m = sqrt(n); ll ans=0;
for(int i=1; i<=m; i++) if(mu[i]) ans += mu[i]>0 ? (n / ((ll)i*i)) : -(n / ((ll)i*i));
return ans;
}
inline ll sum_si(ll n) {
if(n <= U) return si[n];
ll ans=0, r;
for(ll i=1; i<=n; i=r+1) {
r = n/(n/i);
ans += (r-i+1) * (n/i);
}
return ans;
}
ll solve(ll n) {
ll ans=0, r, last=0, now;
for(ll i=1; i<=n; i=r+1, last=now) {
r = n/(n/i); //printf("begin %lld\n", r);
now = sum_mu2(r); //printf("[%lld, %lld]\n", i, r);
ans += (now - last) * sum_si(n/i);
}
return ans;
}
ll q[10005], mx;
int main() {
freopen("in", "r", stdin);
int T=read();
for(int i=1; i<=T; i++) q[i]=read(), mx = max(mx, q[i]);
U = pow(mx,2/3.0);
if(mx >= 1e6 && mx <= 1e8+1) U=1e8;
sieve(U);
for(int i=1; i<=T; i++) printf("%lld\n", solve(q[i]));
}