Luogu P3327 [SDOI2015]约数个数和

又是恶心的莫比乌斯反演，蒟蒻我又是一脸懵逼的被CXR dalao狂虐。

题目要求(ans=sum_{i=1}^n sum_{j=1}^m d(ij))，其中(d(ij))表示数(x)的约数个数

这道题的一大难点就在于(d(ij))这个函数，它有一个重要的性质：

[d(ij)=sum_{x|i}sum_{y|i}[gcd(i,j)=1] ]

大致的证明思路就是对于(i,j)的所有约数，为了避免重复计算，我们只取互质的一对。

知道了这个就是反演的套路了（如果不知道为什么要设去看Luogu P2257 YY的GCD's Sol）：

[f(d)=sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=d] ]

[F(n)=sum_{n|d} f(d) ]

还是由莫比乌斯反演定理得到：

[f(n)=sum_{n|d}mu(lfloorfrac{d}{n} floor)F(d) ]

然后对于原式(ans=sum_{i=1}^n sum_{j=1}^m d(ij))有：

[ans=sum_{i=1}^n sum_{j=1}^msum_{x|i}sum_{y|i}[gcd(i,j)=1] ]

然后根据莫比乌斯函数的性质：(sum_{d|n} mu(d)=[n=1])，我们可以把这个性质带进去：

[ans=sum_{i=1}^n sum_{j=1}^msum_{x|i}sum_{y|i}sum_{d|gcd(x,y)}mu(d) ]

我去，5个Sigma，我们更换一下枚举项，由枚举(d|gcd(x,y))变为直接枚举(d)

[ans=sum_{i=1}^n sum_{j=1}^msum_{x|i}sum_{y|i}sum_{d=1}^{min(n,m)}mu(d)cdot[d|gcd(x,y)] ]

把(mu(d))提到前面来：

[ans=sum_{d=1}^{min(n,m)}mu(d)sum_{i=1}^n sum_{j=1}^msum_{x|i}sum_{y|i}cdot[d|gcd(x,y)] ]

重点来了，我们原来是在枚举(i,j)和他们的约数。这样对计算很不利。

我们现在直接枚举它们的约数再乘上这些约数的倍数的个数。

为什么？因为每一个约数都只会对它的倍数产生贡献，所以有：

[ans=sum_{d=1}^{min(n,m)}mu(d)sum_{x=1}^n sum_{y=1}^m[d|gcd(x,y)]lfloorfrac{n}{x} floorlfloorfrac{m}{y} floor ]

这时候干扰计算的只剩下([d|gcd(x,y)])这个条件了，我们考虑如何消去。

注意到当([d|gcd(x,y)])时意味着(d|x)且(d|y)，所以我们直接枚举(x,y)关于(d)的倍数关系，则有：

[ans=sum_{d=1}^{min(n,m)}mu(d)sum_{x=1}^{lfloorfrac{n}{d} floor} sum_{y=1}^{lfloorfrac{m}{d} floor}lfloorfrac{n}{dx} floorlfloorfrac{m}{dy} floor ]

由于(lfloorfrac{n}{dx} floor)和(lfloorfrac{m}{dy} floor)互不干扰，我们拆开考虑：

[ans=sum_{d=1}^{min(n,m)}mu(d)(sum_{x=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{dx} floor) (sum_{y=1}^{lfloorfrac{m}{d} floor}lfloorfrac{m}{dy} floor) ]

这时候我们惊喜的发现(sum_{x=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{dx} floor)就是一裸的除法分块式子，我们可以直接预处理出(g(x)=sum_{i=1}^xlfloorfrac{x}{i} floor)，这样显然是(O(nsqrt n))的。

然后再回头看一眼式子，把(sum_{x=1}^{lfloorfrac{n}{d} floor}lfloorfrac{n}{dx} floor)用(g(lfloorfrac{n}{d} floor))代进去（另一边同理）得：

[ans=sum_{d=1}^{min(n,m)}mu(d)cdot g(lfloorfrac{n}{d} floor)cdot g(lfloorfrac{m}{d} floor) ]

再看上去就很舒服了，我们发现这个式子可以再来一遍除法分块，就可以做到单次询问(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('
')); RI ptop=0;
            while (x) pt[++ptop]=x%10,x/=10; while (ptop) pc(pt[ptop--]+48); pc('
');
        }
        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;
};