Luogu P3327 [SDOI2015]约数个数和 题解
题目大意
求
\[\sum_{i=1}^n\sum_{j=1}^m d(ij)
\]
其中 \(n,m,T\le 5\cdot 10^4\)，多组测试数据，\(T\) 为测试组数。
题目解析
众所周知，考虑莫比乌斯反演的时候，需要出现 \(\gcd\)，然后枚举 \(\gcd\)，展开 \([x=1]\)。
所以我们考虑咋么出现 \(\gcd\)。
根据 人类智慧，我们可以得到约数函数 \(d\) 具有如下性质：
\[\begin{aligned}
& d(ij) \\
= & \sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]\\
= & \sum_{x|i}\sum_{y|j}\sum_{k|x,k|y}\mu(k)\\
= & \sum_{k|i,k|j}\mu(k)\sum_{k|x,x|i}\sum_{k|y,y|j}1\\
= & \sum_{k|i,k|j}\mu(k)d\left(\dfrac{i}{k}\right)d\left(\dfrac{j}{k}\right)
\end{aligned}
\]
这样的好处就是，这里出现了 \(\dfrac{i}{k}\)，所以最后我们枚举 \(\dfrac{i}{k}\) 的时候，约数函数里面就不会有最外面枚举的 \(k\) 了。
然后就是套路了，直接回带。
\[\begin{aligned}
& \sum_{i=1}^n\sum_{j=1}^md(i,j)\\
= & \sum_{i=1}^n\sum_{j=1}^m\sum_{k|i,k|j}\mu(k)d\left(\dfrac{i}{k}\right)d\left(\dfrac{j}{k}\right)\\
= & \sum_{k=1}^{\min\{n,m\}}\mu(k)\sum_{k|i,i\le n}\sum_{k|j,j\le m} d\left(\dfrac{i}{k}\right)d\left(\dfrac{j}{k}\right)\\
= & \sum_{k=1}^{\min\{n,m\}}\mu(k)\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor m/k\rfloor}d(i)d(j)\\
= & \sum_{k=1}^{\min\{n,m\}}\mu(k) S(\lfloor n/k \rfloor)S(\lfloor m/k \rfloor)
\end{aligned}
\]
其中
\[S(n)=\sum_{i=1}^nd(i)
\]
线性筛+整除分块，\(\Theta(n+T\sqrt n)\)。
其实可以算是多余的代码：
#include<cstdio>
#define db double
#define gc getchar
#define pc putchar
#define U unsigned
#define ll long long
#define ld long double
#define ull unsigned long long
#define Tp template<typename _T>
#define Me(a,b) memset(a,b,sizeof(a))
Tp _T mabs(_T a){ return a>0?a:-a; }
Tp _T mmax(_T a,_T b){ return a>b?a:b; }
Tp _T mmin(_T a,_T b){ return a<b?a:b; }
Tp void mswap(_T &a,_T &b){ _T tmp=a; a=b; b=tmp; return; }
Tp void print(_T x){ if(x<0) pc('-'),x=-x; if(x>9) print(x/10); pc((x%10)+48); return; }
#define EPS (1e-7)
#define INF (0x7fffffff)
#define LL_INF (0x7fffffffffffffff)
#define N 50000
#define maxn 50039
#define maxm
#define MOD
#define Type int
#ifndef ONLINE_JUDGE
//#define debug
#endif
using namespace std;
Type read(){
char c=gc(); Type s=0; int flag=0;
while((c<'0'||c>'9')&&c!='-') c=gc(); if(c=='-') c=gc(),flag=1;
while('0'<=c&&c<='9'){ s=(s<<1)+(s<<3)+(c^48); c=gc(); }
if(flag) return -s; return s;
}
int n,m,isp[maxn],pr[maxn],cnt; ll mu[maxn],d[maxn],g[maxn];
void init(){
int i,j; d[1]=1; mu[1]=1;
for(i=2;i<=N;i++){
if(!isp[i]) pr[++cnt]=i,d[i]=2,mu[i]=-1,g[i]=1;
for(j=1;j<=cnt&&i*pr[j]<=N;j++){
isp[i*pr[j]]=1; mu[i*pr[j]]=-mu[i]; d[i*pr[j]]=d[i]*2; g[i*pr[j]]=1;
if(i%pr[j]==0){ mu[i*pr[j]]=0; g[i*pr[j]]=g[i]+1; d[i*pr[j]]=d[i]/(g[i]+1)*(g[i]+2); break; }
}
}
for(i=1;i<=N;i++) mu[i]+=mu[i-1],d[i]+=d[i-1];
}
void work(){
n=read(); m=read(); int l=1,r; ll ans=0;
while(l<=n&&l<=m){
r=mmin(n/(n/l),m/(m/l));
ans+=1ll*(mu[r]-mu[l-1])*d[n/l]*d[m/l];
l=r+1;
} print(ans),pc('\n');
}
int main(){
//freopen(".in","r",stdin);
//freopen(".out","w",stdout);
init(); int T=read(); while(T--) work(); return 0;
}
