SPOJ DIVCNT2 - Counting Divisors (square)
题意：求
[sum_{i=1}^nsigma_0(i^2)
]
好棒啊！
带着平方没法做，考虑用其他函数表示(sigma_0(i^2))，把平方消去。
(sigma_0(n) = (1*1)(n) = sum_{dmid 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_{dmid n}mu^2(d))
[sigma_0(n^2) = sum_{dmid n} 2^{omega(d)} = sum_{dmid n} sum_{emid d} mu^2(e) = ((mu^2 * 1) * 1) (n)
]
我们就是要求((mu * 1) * 1 = mu * (1*1) = mu * sigma_0)的前缀和
[egin{align}
ans &= sum_{i=1}^n sum_{dmid i} mu^2(d) cdot sigma_0(frac{i}{d}) \
&= sum_{i=1}^n mu^2(i) sum_{j=1}^{lfloor frac{n}{i}
floor} 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} floor) 就是无平方因子数的个数呀
(sum_{i=1}^nsigma_0(i) = sum_{i=1}^n lfloor frac{n}{i} floor) 也可以整除分块
同时我们使用线性筛预处理前(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
", r);
now = sum_mu2(r); //printf("[%lld, %lld]
", 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
", solve(q[i]));
}