- 给定(n)，求有多少对(a,b)满足(1le a<ble n)且((a+b)|ab)。
- (nle2^{31}-1)
初始推式子
说实话，最近翻了翻数论小蓝本，于是看到这种式子首先就想到令(d=gcd(a,b))得到(a_0=frac ad,b_0=frac bd)互质。
然后就发现原式相当于是：
[(a_0+b_0)d|a_0b_0d^2Leftrightarrow (a_0+b_0)|a_0b_0d
]
由于(gcd(a_0+b_0,a_0)=gcd(b_0,a_0)=1)（同理(gcd(a_0+b_0,b_0)=1)），所以我们可以直接消去整除符号右边的(a_0,b_0)，得到一个非常漂亮的式子：
[(a_0+b_0)|d
]
考虑我们现在的几个限制：
[egin{cases}
b=b_0dle nRightarrow dlelfloorfrac n{b_0}
floor\
(a_0+b_0)|dRightarrow a_0+b_0le d
end{cases}
]
综合两个条件得到(a_0+b_0lelfloorfrac n{b_0} floor)。
所以说，我们发现，(a_0,b_0,d)都应该是不超过(sqrt n)的。
莫比乌斯反演
根据上面的结论，答案的计算式就应该是：
[sum_{a_0=1}^{sqrt n}sum_{b_0=a_0+1}^{sqrt n}lfloorfrac{lfloorfrac{n}{b_0}
floor}{a_0+b_0}
floor[gcd(a_0,b_0)=1]
]
经典莫比乌斯反演，枚举(gcd)得到：
[sum_{d=1}^{sqrt n}mu(d)sum_{a_0=1}^{lfloorfrac{sqrt n}d
floor}sum_{b_0=a_0+1}^{lfloorfrac{sqrt n}d
floor}lfloorfrac{lfloorfrac n{b_0d}
floor}{(a_0+b_0)d}
floor
]
其中大分母中的(d)可以移到分子的分母中：
[sum_{d=1}^{sqrt n}mu(d)sum_{a_0=1}^{lfloorfrac{sqrt n}d
floor}sum_{b_0=a_0+1}^{lfloorfrac{sqrt n}d
floor}lfloorfrac{lfloorfrac n{b_0d^2}
floor}{a_0+b_0}
floor
]
我们发现(a_0)只出现在分母(a_0+b_0)这一项中，因此我们不如不枚举(a_0)，而是去枚举分母，得到：
[sum_{d=1}^{sqrt n}mu(d)sum_{b_0=2}^{lfloorfrac{sqrt n}d
floor}sum_{t=b_0+1}^{2b_0-1}lfloorfrac{lfloorfrac n{b_0d^2}
floor}{t}
floor
]
发现此时的(sum_{t=b_0+1}^{2b_0-1}lfloorfrac{lfloorfrac n{b_0d^2} floor}{t} floor)已经成为一个标准的除法分块形式。
那么我们只要以共计(O(sqrt nlogsqrt n))的复杂度枚举(d,b_0)，然后除法分块(O(sqrt{sqrt n}))枚举(t)即可解决此题。
代码：(O(n^{frac34}logsqrt n))
#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define SN 47000
#define LL long long
using namespace std;
int n,sn;struct LinearSieve
{
int Pt,P[SN+5],mu[SN+5];I int operator [] (CI x) {return mu[x];}
I LinearSieve()//线性筛预处理μ
{
mu[1]=1;for(RI i=2,j;i<=SN;++i) for(!P[i]&&(mu[P[++Pt]=i]=-1),
j=1;j<=Pt&&i*P[j]<=SN;++j) if(P[i*P[j]]=1,i%P[j]) mu[i*P[j]]=-mu[i];else break;
}
}Mu;
I LL Calc(CI L,CI R,CI n)//除法分块
{
RI l,r;LL t=0;for(l=L;l<=min(n,R);l=r+1) r=min(R,n/(n/l)),t+=1LL*(n/l)*(r-l+1);return t;
}
int main()
{
RI i,j;LL t=0;for(scanf("%d",&n),sn=sqrt(n),i=1;i<=sn;++i)//枚举d
for(j=2;i*j<=sn;++j) t+=Mu[i]*Calc(j+1,2*j-1,n/i/i/j);return printf("%lld
",t),0;//枚举b0，然后除法分块
}