【bzoj3994】 SDOI2015—约数个数和

http://www.lydsy.com/JudgeOnline/problem.php?id=3994 (题目链接)

题意

　　多组询问，给出${n,m}$，求${\sum_{i=1}^n\sum_{j=1}^m d(i×j)}$，${d(i×j)}$为${ij}$的约数个数。

Solution

　　看到这个式子感觉无从下手，这个${d(ij)}$比较丑，有一个比较经典的公式：$${d(nm)=\sum_{i|n}\sum_{j|m} [gcd(i,j)=1]}$$

　　这个是怎么得来的呢。每一个${nm}$的约数都可以表示成这样的形式：${i×\frac{m}{j}}$。其中${i}$是${n}$的约数，${j}$是${m}$的约数。

　　那么如果我们直接枚举${i,j}$来统计的话，可能会有重复，所以就要加上一个条件：${[gcd(i,j)=1]}$，也就是${i,j}$互质。那他们不互质为什么就不行呢，假设${i,j}$都含有一个约数${p}$，那么${\frac{m}{j}}$就表示从${m}$中拿掉${j}$用剩下的去组成约数，也就是从约数中拿掉了${p}$；而${i}$就是从${n}$中拿${i}$去组成约数，也就向约数中加入了${p}$。那么我们拿掉一个${p}$又加入一个${p}$，不是吃多了吗→_→，这样自然会算重复啦。

　　解决了这个问题，我们就很好做了。开始推式子。

\begin{aligned}  \sum_{i=1}^{n}\sum_{j=1}^{m}d(ij) =& \sum_{i=1}^n\sum_{j=1}^m\sum_{u|i}\sum_{v|j}[gcd(u,v)=1]  \\  =&\sum_{u=1}^n\sum_{v=1}^m[gcd(u,v)=1]\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{j}\rfloor  \\  =&\sum_{t=1}^nμ(t)\sum_{i=1}^{\lfloor{n/t}\rfloor}\sum_{j=1}^{\lfloor{m/t}\rfloor}\lfloor\frac{n}{ti}\rfloor\lfloor\frac{m}{tj}\rfloor  \\  =&\sum_{t=1}^nμ(t)\sum_{i=1}^{\lfloor{n/t}\rfloor}\lfloor\frac{\lfloor{n/t}\rfloor}{i}\rfloor\sum_{j=1}^{\lfloor{m/t}\rfloor}\lfloor\frac{\lfloor{m/t}\rfloor}{j}\rfloor   \end{aligned}

　　我们令${f(n)=\sum_{i=1}^n\lfloor{n/i}\rfloor}$，那么${f}$是可以分段${O(n\sqrt{n})}$的预处理出来的。$${原式=\sum_{t=1}^{n}μ(t)f(\lfloor{n/t}\rfloor)f(\lfloor{m/t}\rfloor)}$$

　　这样我们就可以预处理出${μ}$的前缀和，然后分段求解答案就可以了。

细节

　　LL，预处理不要太多会TLE。

代码

// bzoj3994
#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<vector>
#include<cmath>
#define LL long long
#define inf 2147483640
#define Pi acos(-1.0)
#define free(a) freopen(a".in","r",stdin),freopen(a".out","w",stdout);
using namespace std;

const int maxn=50010;
LL n,m,s[maxn],mu[maxn],f[maxn];
int p[maxn],vis[maxn];

int main() {
	int T;scanf("%d",&T);
	s[1]=mu[1]=1;
	for (int i=2;i<maxn;i++) {
		if (!vis[i]) p[++p[0]]=i,mu[i]=-1;
		for (int j=1;j<=p[0] && p[j]*i<maxn;j++) {
			vis[i*p[j]]=1;
			if (i%p[j]==0) {mu[i*p[j]]=0;break;}
			else mu[i*p[j]]=-mu[i];
		}
		s[i]=s[i-1]+mu[i];
	}
	for (int i=1;i<maxn;i++)
		for (int j=1,k;j<=i;j=k+1) {
			k=i/(i/j);
			f[i]+=(LL)(k-j+1)*(i/j);
		}
	while (T--) {
		scanf("%lld%lld",&n,&m);
		if (n>m) swap(n,m);
		LL ans=0;
		for (int i=1,j;i<=n;i=j+1) {
			j=min(n/(n/i),m/(m/i));
			ans+=f[n/i]*f[m/i]*(s[j]-s[i-1]);
		}
		printf("%lld\n",ans);
	}
	return 0;
}