Luogu P3327 [SDOI2015]约数个数和
又是恶心的莫比乌斯反演,蒟蒻我又是一脸懵逼的被CXR dalao狂虐。
题目要求\(ans=\sum_{i=1}^n \sum_{j=1}^m d(ij)\),其中\(d(ij)\)表示数\(x\)的约数个数
这道题的一大难点就在于\(d(ij)\)这个函数,它有一个重要的性质:
大致的证明思路就是对于\(i,j\)的所有约数,为了避免重复计算,我们只取互质的一对。
知道了这个就是反演的套路了(如果不知道为什么要设去看Luogu P2257 YY的GCD's Sol):
还是由莫比乌斯反演定理得到:
然后对于原式\(ans=\sum_{i=1}^n \sum_{j=1}^m d(ij)\)有:
然后根据莫比乌斯函数的性质:\(\sum_{d|n} \mu(d)=[n=1]\),我们可以把这个性质带进去:
我去,5个Sigma,我们更换一下枚举项,由枚举\(d|\gcd(x,y)\)变为直接枚举\(d\)
把\(\mu(d)\)提到前面来:
重点来了,我们原来是在枚举\(i,j\)和他们的约数。这样对计算很不利。
我们现在直接枚举它们的约数再乘上这些约数的倍数的个数。
为什么?因为每一个约数都只会对它的倍数产生贡献,所以有:
这时候干扰计算的只剩下\([d|\gcd(x,y)]\)这个条件了,我们考虑如何消去。
注意到当\([d|\gcd(x,y)]\)时意味着\(d|x\)且\(d|y\),所以我们直接枚举\(x,y\)关于\(d\)的倍数关系,则有:
由于\(\lfloor\frac{n}{dx}\rfloor\)和\(\lfloor\frac{m}{dy}\rfloor\)互不干扰,我们拆开考虑:
这时候我们惊喜的发现\(\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{dx}\rfloor\)就是一裸的除法分块式子,我们可以直接预处理出\(g(x)=\sum_{i=1}^x\lfloor\frac{x}{i}\rfloor\),这样显然是\(O(n\sqrt n)\)的。
然后再回头看一眼式子,把\(\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{dx}\rfloor\)用\(g(\lfloor\frac{n}{d}\rfloor)\)代进去(另一边同理)得:
再看上去就很舒服了,我们发现这个式子可以再来一遍除法分块,就可以做到单次询问\(O(\sqrt n)\)了。
CODE
#include<cstdio>
#include<cctype>
#define RI register int
using namespace std;
const int P=50000;
int t,n,m,prime[P+5],cnt,mu[P+5],sum[P+5],lim; long long g[P+5],ans; bool vis[P+5];
class FileInputOutput
{
private:
#define tc() (A==B&&(B=(A=Fin)+fread(Fin,1,S,stdin),A==B)?EOF:*A++)
#define pc(ch) (Ftop<S?Fout[Ftop++]=ch:(fwrite(Fout,1,S,stdout),Fout[(Ftop=0)++]=ch))
#define S 1<<21
char Fin[S],Fout[S],*A,*B; int Ftop,pt[25];
public:
FileInputOutput() { A=B=Fin; Ftop=0; }
inline void read(int &x)
{
x=0; char ch; while (!isdigit(ch=tc()));
while (x=(x<<3)+(x<<1)+(ch&15),isdigit(ch=tc()));
}
inline void write(long long x)
{
if (!x) return (void)(pc(48),pc('\n')); RI ptop=0;
while (x) pt[++ptop]=x%10,x/=10; while (ptop) pc(pt[ptop--]+48); pc('\n');
}
inline void Fend(void)
{
fwrite(Fout,1,Ftop,stdout);
}
#undef tc
#undef pc
#undef S
}F;
#define Pi prime[j]
inline void Euler(void)
{
vis[1]=mu[1]=1; RI i,j; for (i=2;i<=P;++i)
{
if (!vis[i]) prime[++cnt]=i,mu[i]=-1;
for (j=1;j<=cnt&&i*Pi<=P;++j)
{
vis[i*Pi]=1; if (i%Pi) mu[i*Pi]=-mu[i]; else break;
}
}
for (i=1;i<=P;++i) sum[i]=sum[i-1]+mu[i];
for (i=1;i<=P;++i) for (RI l=1,r;l<=i;l=r+1) r=i/(i/l),g[i]+=1LL*(r-l+1)*(i/l);
}
inline int min(int a,int b)
{
return a<b?a:b;
}
int main()
{
//freopen("CODE.in","r",stdin); freopen("CODE.out","w",stdout);
for (Euler(),F.read(t);t;--t)
{
F.read(n); F.read(m); ans=0; lim=min(n,m);
for (RI l=1,r;l<=lim;l=r+1)
{
r=min(n/(n/l),m/(m/l)); ans+=1LL*(sum[r]-sum[l-1])*g[n/l]*g[m/l];
}
F.write(ans);
}
return F.Fend(),0;
};