「bzoj2671」Calc - 数论+莫比乌斯反演

题目大意

    给出$N$，统计满足下面条件的数对$(a,b)$的个数：

1. $1\le a\lt b\le N$
2. $a+b$ 整除 $a\cdot b$

---

题目分析

做这道题时跟着CQzhangyu学，然后发现他给出了假的式子和假的反演，代码却是对的。

前置技能

• 使用$(a,b)$表示$\gcd(a,b)$
• 若$(a,b)=1$，则$(a+b,ab)=1$
  证明（来自Hineven）：
  反证法。
  设$(a+b,ab)=k\gt1$，因为$(a,b)$互质，故$k\mid a$或$k\mid b$。
  设$k\mid a$，则： $$a=a'k \Rightarrow a+b = a'k+b \Rightarrow k|(a'k+b) \Rightarrow k|b \Rightarrow k|\gcd(a,b)$$
• $$\left\lfloor\frac{\lfloor\frac ab\rfloor}{c}\right\rfloor=\left\lfloor\frac a{bc}\right\rfloor$$

正式分析

设$d=\gcd(a,b)$
则：

$$a+b\mid ab\Rightarrow a'd+b'd\mid a'b'd^2\Rightarrow a'+b'\mid a'b'd\Rightarrow a'+b'\mid d$$

问题就变成了问有多少对$a,b$满足$(a,b)=1$且$b\times d\le n$满足$a+b\mid d$，又因为$d\gt b$所以b最多只能是$m=\sqrt n$。

$$\sum_{b=1}^m\sum_{a=1}^{b-1}[\gcd(a,b)=1]\left\lfloor \frac{n}{(a+b)b}\right\rfloor$$

最后的$\lfloor \frac{n}{(a+b)b}\rfloor$是因为，$b\times d\le n$，故$d\le\lfloor\frac nb\rfloor$，又因为$a+b\mid d$，则合法的数对有$\left\lfloor\frac{\lfloor\frac nb\rfloor}{a+b}\right\rfloor=\left\lfloor\frac{n}{(a+b)b}\right\rfloor$个。
反演一波：

$$\begin{aligned} ans&=\sum_{b=1}^m\sum_{a=1}^{b-1}[\gcd(a,b)=1]\left\lfloor \frac{n}{(a+b)b}\right\rfloor \\ &=\sum_{b=1}^m\sum_{a=1}^{b-1}\sum_{d\mid a,d\mid b}\mu(d)\left\lfloor \frac{n}{(a+b)b}\right\rfloor \\ &=\sum_d\mu(d)\sum_{d\mid b}^m\sum_{d\mid a}^{b-1}\left\lfloor \frac{n}{(a+b)b}\right\rfloor \\ &=\sum_d\mu(d)\sum_{b=1}^{\lfloor \frac md\rfloor}\sum_{a=1}^{b-1}\left\lfloor \frac{n}{d^2b(a+b)}\right\rfloor \end{aligned}$$

不难发现，最后这一部分：

$$\sum_{a=1}^{b-1}\left\lfloor \frac{n}{d^2b(a+b)}\right\rfloor=\sum_{a=1}^{b-1}\left\lfloor\frac{\left\lfloor\frac{n}{d^2b}\right\rfloor}{a+b}\right\rfloor$$

可以枚举$a+b$整除分块。

时间复杂度由简单积分可知为$O(n^{\frac34}-n^{\frac12})$。

---

代码

#include<bits/stdc++.h>

using namespace std;

typedef long long LL; 

inline int Get_Int() {
	int num=0,bj=1;
	char x=getchar();
	while(!isdigit(x)) {if(x=='-')bj=-1;x=getchar();}
	while(isdigit(x)) {num=num*10+x-'0';x=getchar();}
	return num*bj;
}

const int maxn=1000005;

int vst[maxn],p[maxn],u[maxn],cnt=0;

void Mobius_Table(int n) {
	u[1]=1;
	for(int i=2; i<=n; i++) {
		if(!vst[i]) {p[++cnt]=i;u[i]=-1;}
		for(int j=1; j<=cnt&&i*p[j]<=n; j++) {
			vst[i*p[j]]=1;
			if(i%p[j]==0) {u[i*p[j]]=0;break;}
			u[i*p[j]]=-u[i];
		}
	}
}

LL Cal(int n,int m) {
	LL ans=0;
	for(int b=1; b<=m; b++) {
		int up=n/b;
		for(int i=b+1,last; i<(b<<1)&&i<=up; i=last+1) {
			last=min((b<<1)-1,up/(up/i));
			ans+=1ll*(last-i+1)*(up/i);
		}
	}
	return ans;
}

int n,m;

int main() {
	n=Get_Int();
	m=sqrt(n);
	Mobius_Table(m);
	LL ans=0;
	for(int d=1; d<=m; d++)ans+=u[d]*Cal(n/d/d,m/d);
	printf("%lld\n",ans);
	return 0;
}