「SPOJ」DIVCNT2-Counting Divisors(square)

Problem

Description

S_{2}(n)=\sum_{i=1}^{n}\tau(i^{2})=\sum_{i=1}^{n}\sigma_{0}(i^{2})

S_{2}(n)

Constraints

1\le n\le 10^{12}

Solution

Analysis

考虑如何计算 \tau(n^{2})
对于 n 的每个约数 d ,可以在 d^{2} 中去除 d 的一个质因子集合。
可以证明这样可以不重复不遗漏地得到 n^{2} 的所有约数。
而对于 d 的每个质因子都有去掉与保留两种选择,总方案数即为 2^{\omega(d)}
即:

\tau(n^{2})=\sum_{d | n}2^{\omega(d)}

考虑如何计算 2^{\omega(n)} ,可以发现其意义即为 n 的无平方因子的约数个数,即:

2^{\omega(n)}=\sum_{d | n}\mu^{2}(d)

\tau(n^{2})=\sum_{d | n}2^{\omega(d)}=\sum_{d | n}\sum_{e | d}\mu^{2}(e)=(\mu^{2}\ast 1\ast 1)(n)=(\mu^{2}\ast\tau)(n)

故:

S_{2}(n)=\sum_{i=1}^{n}\tau(i^{2})=\sum_{i=1}^{n}(\mu^{2}\ast\tau)(i)=\sum_{i=1}^{n}\sum_{d | i}\mu^{2}(d)\tau\left(\frac{i}{d}\right)=\sum_{d=1}^{n}\mu^{2}(d)\sum_{id \le n}\tau(i)

故需要求出 \mu^{2}(n) \tau(n) 的前缀和。
易知 \mu^{2}(n) 的前缀和即为 [1, n] 中没有大于 1 的完全平方数因子的数的个数。
设:

f(n)=\max\left\{d^{2} | d\in\mathbb{Z^{+}},d^{2} | n\right\}

则:

\sum_{i=1}^{n}\mu^{2}(i)=\sum_{i=1}^{n}\left[f(i)=1\right]=\sum_{i=1}^{n}\sum_{d | f(i)}\mu(d)

发现 \mu(d)\ne 0 当且仅当 d 无大于 1 的完全平方数因子。
又由 f(n) 的定义知 f(i) 必为完全平方数,则对于所有有贡献的 \mu(d) ,均有 d^{2} | f(i)
即:

\begin{aligned} \sum_{i=1}^{n}\mu^{2}(i)&=\sum_{i=1}^{n}\sum_{d^{2} | f(i)}\mu(d)\\ &=\sum_{i=1}^{n}\sum_{d^{2} | i}\mu(d)\\ &=\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{n}\left[d^{2} | i\right]\\ &=\sum_{d=1}^{n}\mu(d)(n / d^{2})\\ &=\sum_{d=1}^{\sqrt{n}}\mu(d)(n / d^{2}) \end{aligned}

又易知

\sum_{i=1}^{n}\tau(i)=\sum_{i=1}^{n}\sum_{d=1}^{n}\left[d | i\right]=\sum_{d=1}^{n}\sum_{i=1}^{n}\left[d | i\right]=\sum_{d=1}^{n}(n / d)

考虑杜教筛的过程,类似地,可以使用线性筛预处理前 n^{\frac{2}{3}} 个前缀和,其余的部分使用以上公式可以 O\left(\sqrt{n}\right) 计算。

复杂度 O\left(n^{\frac{2}{3}}\right)

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#include<cmath>
#include<cstdio>

using i64=long long;

template<typename _Tp>
inline void chkMax(_Tp&x,const _Tp&y)
{(x<y)&&(x=y);}

inline int gi(){
int s=0;char c;
do c=getchar();while(c<48||c>57);
do s=s*10+c-48,c=getchar();while(c>47&&c<58);
return s;
}

inline i64 GetLL(){
i64 s=0;char c;
do c=getchar();while(c<48||c>57);
do s=s*10+c-48,c=getchar();while(c>47&&c<58);
return s;
}

constexpr int maxq(10000);
const int maxn(100000000);
const int p(10000000);
/*maxn/ln(maxn)=5428681.0249554...*/

bool notp[maxn+1];
int prime[_p];
int cnt;

short mu[maxn+1];
int Mu2[maxn+1];

i64 tau[maxn+1];
short Minnu[maxn+1];

void Sieve(const int&_L){
mu[1]=1;
tau[1]=1;
Mu2[1]=1;
for(int i=2;i<=_L;++i){
if(!notp[i])
prime[++cnt]=i,mu[i]=-1,tau[i]=Minnu[i]=2;

i64 v;
for(int j=1;j<=cnt;++j){
v=(i64)i*prime[j];
if(v>_L)
break;

notp[v]=true;
if(i%prime[j])
mu[v]=-mu[i],tau[v]=tau[i]<<1,Minnu[v]=2;
else{
mu[v]=0;
Minnu[v]=Minnu[i]+1;
tau[v]=tau[i]/Minnu[i]*Minnu[v];
break;
}
}

Mu2[i]=Mu2[i-1]+mu[i]*mu[i];
tau[i]+=tau[i-1];
}
}

int _lim;

i64 calcMu2(const i64&n){
if(n<=_lim)
return Mu2[n];
i64 ans=0;
for(int d=1,sqrn=sqrt(n);d<=sqrn;++d)
if(mu[d])
~mu[d]?(ans+=n/((i64)d*d)):(ans-=(n/((i64)d*d)));
return ans;
}

inline i64 calcTau(const i64&n){
if(n<=_lim)
return tau[n];
i64 ans=0;
for(i64 d=1,nxd;d<=n;d=nxd+1){
nxd=n/(n/d);
ans+=(nxd-d+1)*(n/d);
}return ans;
}

inline i64 calc(const i64&n){
i64 ans=0;
for(i64 d=1,nxd,las=0,cur;d<=n;d=nxd+1,las=cur){
nxd=n/(n/d);
cur=calcMu2(nxd);
ans+=(cur-las)*calcTau(n/d);
}return ans;
}

i64 Q[maxq+1];

int main(){
int T=gi();
i64 Maxn=0;
for(int i=1;i<=T;++i)
Q[i]=GetLL(),chkMax(Maxn,Q[i]);

_lim=pow(Maxn,2./3.);
if(1e6<=_lim&&_lim<=1e8+1)
_lim=1e8;

Sieve(_lim);//precalc(_lim);

for(int i=1;i<=T;++i)
printf("%lld\n",calc(Q[i]));

return 0;
}