[SDOI 2015] 约数个数和

\(\text{Description}\)

\(\mathcal{P}\text{ortal.}\)

\(\text{Solution}\)

首先有一个前导结论

\[\sigma(x y)=\sum_{i|x}\sum_{j|y} [\gcd(i,j)=1] \]

具体证明 戳这。那么柿子就变成了

\[\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j} [\gcd(x,y)=1] \]

而莫反的常见形式是这个样子

\[\sum_{i=1}^n\sum_{j=1}^m [\gcd(i,j)=1] \]

尝试将 \(x,y\) 提出去，向常见形式靠拢。可以发现 \(\sum_{i=1}^n\sum_{j=1}^m\) 实际上就是个系数

\[\sum_{x=1}^n\sum_{y=1}^m[\gcd(x,y)=1]\cdot \left \lfloor\frac{n}{x}\right \rfloor\cdot \left \lfloor\frac{m}{y}\right \rfloor \]

于是设

\[g(k)=\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]\cdot \left \lfloor\frac{n}{i}\right \rfloor\cdot \left \lfloor\frac{m}{j}\right \rfloor \]

\[f(k)=\sum_{k|d}g(d) \]

所以

\[\text{Ans}=g(1)=\sum_{i=1}^{\min\{n,m\}}\mu(i) f(i) \]

我们需要快速求出 \(f(i)\)。不过无法像直接猛拆二重循环，因为有 \(\displaystyle \left \lfloor\frac{n}{i}\right \rfloor\cdot\left \lfloor\frac{m}{j}\right \rfloor\).

尝试直接化柿子

\[f(k)=\sum_{i=1}^n\sum_{j=1}^m[k|\gcd(i,j)]\cdot \left \lfloor\frac{n}{i}\right \rfloor\cdot \left \lfloor\frac{m}{j}\right \rfloor \]

可以预见枚举 \(\gcd\) 是没啥前途的，因为后面是整除符号。其实可以将 \(k\) 提出来，这样 \(\gcd\) 肯定是 \(k\) 的倍数

\[=\sum_{i=1}^{\left \lfloor\frac{n}{k}\right \rfloor}\sum_{j=1}^{\left \lfloor\frac{m}{k}\right \rfloor}\left \lfloor\frac{n}{i k}\right \rfloor\cdot \left \lfloor\frac{m}{j k}\right \rfloor \]

预处理 \(\displaystyle \sum_{i=1}^n \left \lfloor\frac{n}{i}\right \rfloor\) 就行了，因为 \(\displaystyle \left \lfloor\frac{n}{i k}\right \rfloor=\left \lfloor\frac{\left \lfloor\frac{n}{k}\right \rfloor}{i}\right \rfloor\).

时间复杂度 \(\mathcal O((T+n)\cdot n^{1/2})\).

\(\text{Code}\)

#include <cstdio>

#define rep(i,_l,_r) for(register signed i=(_l),_end=(_r);i<=_end;++i)
#define fep(i,_l,_r) for(register signed i=(_l),_end=(_r);i>=_end;--i)
#define erep(i,u) for(signed i=head[u],v=to[i];i;i=nxt[i],v=to[i])
#define efep(i,u) for(signed i=Head[u],v=to[i];i;i=nxt[i],v=to[i])
#define print(x,y) write(x),putchar(y)

template <class T> inline T read(const T sample) {
    T x=0; int f=1; char s;
    while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
    while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
    return x*f;
}
template <class T> inline void write(const T x) {
    if(x<0) return (void) (putchar('-'),write(-x));
    if(x>9) write(x/10);
    putchar(x%10^48);
}
template <class T> inline T Max(const T x,const T y) {if(x>y) return x; return y;}
template <class T> inline T Min(const T x,const T y) {if(x<y) return x; return y;}
template <class T> inline T fab(const T x) {return x>0?x:-x;}
template <class T> inline T gcd(const T x,const T y) {return y?gcd(y,x%y):x;}
template <class T> inline T lcm(const T x,const T y) {return x/gcd(x,y)*y;}

typedef long long ll;

const int maxn=5e4+5;

int n,m,u[maxn],p[maxn],pc;
bool is[maxn];
ll f[maxn];

void init() {
	u[1]=1;
	rep(i,2,maxn-5) {
		if(!is[i]) p[++pc]=i,u[i]=-1;
		rep(j,1,pc) {
			if(p[j]*i>maxn-5) break;
			is[p[j]*i]=1;
			if(i%p[j]==0) break;
			u[i*p[j]]=-u[i];
		}
	} 
	rep(i,2,maxn-5) u[i]+=u[i-1];
	int r; ll ret;
	rep(i,1,maxn-5) {
		ret=0;
		for(int l=1;l<=i;l=r+1) r=i/(i/l),ret+=1ll*(r-l+1)*(i/l);
		f[i]=ret;
	}
}

ll Query() {
	int lim=Min(n,m),r; ll ret=0;
	for(int l=1;l<=lim;l=r+1) {
		r=Min(n/(n/l),m/(m/l));
		ret+=1ll*(u[r]-u[l-1])*f[n/l]*f[m/l];
	}
	return ret;
}

int main() {
	init();
	for(int T=read(9);T;--T) {
		n=read(9),m=read(9);
		print(Query(),'\n');
	}
	return 0;
}