\(\text{Description}\)

传送门

\(\text{Solution}\)

首先有：

\[\sigma(x\times 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]\times \left \lfloor\frac{n}{x}\right \rfloor\times \left \lfloor\frac{m}{y}\right \rfloor \]

设，

\[g(k)=\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]\times \left \lfloor\frac{n}{i}\right \rfloor\times \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)\times f(i) \]

我们需要快速求出 \(f(i)\)（根据经验，应该是可以整除分块的）。

不过这次似乎无法像 \(\sum_{i=1}^n\sum_{j=1}^m [\gcd(i,j)=1],\sum_{i=1}^n\sum_{j=1}^m [\gcd(i,j)=1]\times i\times j\) 这样的柿子猛拆二重循环，因为有 \(\left \lfloor\frac{n}{i}\right \rfloor\times \left \lfloor\frac{m}{j}\right \rfloor\) 碍事的。

尝试直接化柿子：

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

显然我们无法简化 \(i,j\)，显然不能枚举 \(\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\times k}\right \rfloor\times \left \lfloor\frac{m}{j\times k}\right \rfloor \]

预处理 \(\sum_{i=1}^n \left \lfloor\frac{n}{i}\right \rfloor\) 就行了。

时间复杂度 \(\mathcal O((T+n)\times \sqrt n)\)。

\(\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;
} 

标签：约数,lfloor,right,frac,SDOI,sum,rfloor,times,2015
来源： https://www.cnblogs.com/AWhiteWall/p/14269865.html