# [BZOJ 3994] 约数个数和

## 题意

$\sum_{i=1}^n\sum_{j=1}^md(ij)$

$$n,m\le 5\times 10 ^ 4, q\le 5\times 10^4$$

## 题解

\begin{align} d(ij)&=\sum_{a|i}\sum_{b|j}[a\perp b] \\ \text{Ans}&=\sum_i\sum_jd(ij)\\ &=\sum_i\sum_j\sum_{a|i}\sum_{b|j}[i\perp j] \\ &=\sum_i\sum_j\sum_{a|i}\sum_{b|j}\sum_{k|a,k|b}\mu(k)\\ &=\sum_k\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\sum_a^{\lfloor \frac n {ki} \rfloor}\sum_b^{\lfloor \frac m {kj} \rfloor}\mu(k) \\ &=\sum_k\mu(k)\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\sum_a^{\lfloor \frac n {ki} \rfloor}\sum_b^{\lfloor \frac m {kj} \rfloor}1 \\ &=\sum_k\mu(k)\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\left\lfloor\frac n {ki}\right\rfloor\left\lfloor\frac m {kj}\right\rfloor \\ &=\sum_k\mu(k)\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\left\lfloor\frac {\left\lfloor\frac n k\right\rfloor} {i}\right\rfloor\left\lfloor\frac {\left\lfloor\frac m k\right\rfloor} {j}\right\rfloor\\ %&=\sum_k\sum_i^{\lfloor \frac n k \rfloor}\sum_j^{\lfloor \frac m k \rfloor}\left\lfloor\frac {\left\lfloor\frac n k\right\rfloor} {i}\right\rfloor\left\lfloor\frac {\left\lfloor\frac m k\right\rfloor} {j}\right\rfloor\mu(k) &=\sum_k\mu(k)\left(\sum_i^{\lfloor \frac n k \rfloor}\left\lfloor\frac {\left\lfloor\frac n k\right\rfloor} {i}\right\rfloor\right)\left(\sum_j^{\lfloor \frac m k \rfloor}\left\lfloor\frac {\left\lfloor\frac m k\right\rfloor} {j}\right\rfloor\right)\\ \end{align}

## 代码实现

#include <bits/stdc++.h>

const int MAXN=5e4+10;

int cnt;
int mu[MAXN];
int pr[MAXN];
bool npr[MAXN];
long long g[MAXN];

long long Calc(int);
void EulerSieve(int);

int main(){
int T;
scanf("%d",&T);
EulerSieve(5e4);
while(T--){
int n,m;
scanf("%d%d",&n,&m);
if(n>m)
std::swap(n,m);
long long ans=0;
for(int i=1,j;i<=n;i=j+1){
j=std::min(n/(n/i),m/(m/i));
ans+=(mu[j]-mu[i-1])*Calc(n/i)*Calc(m/i);
}
printf("%lld\n",ans);
}
return 0;
}

long long Calc(int x){
if(g[x]!=0)
return g[x];
else{
for(int i=1,j;i<=x;i=j+1){
j=x/(x/i);
g[x]+=(j-i+1)*(x/i);
}
return g[x];
}
}

void EulerSieve(int n){
npr[0]=npr[1]=true;
mu[1]=1;
for(int i=2;i<=n;i++){

if(!npr[i]){
pr[cnt++]=i;
mu[i]=-1;
}
for(int j=0,t;j<cnt&&(t=i*pr[j])<=n;j++){
npr[t]=true;
if(i%pr[j])
mu[t]=-mu[i];
else{
mu[t]=0;
break;
}
}
}
for(int i=1;i<=n;i++)
mu[i]+=mu[i-1];
}



posted @ 2019-01-08 09:20  rvalue  阅读(255)  评论(0编辑  收藏  举报