SPOJ DIVCNT2 [我也不知道是什么分类了反正是数论]

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} \]

不用杜教筛，我们也可以求。

我们只要得到\(\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]));
}