d(x)d(x)xx的约数个数,给定NMN、M,求i=1Nj=1Md(ij)\sum_{i=1}^N \sum_{j=1}^M d(ij)

链接🔗

「SDOI2015」约数和个数

题解

约数个数公式 (不会证鸭)

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

原式等于

i=xnj=ymixjy[gcd(i,j)==1]\sum_{i=x}^n \sum_{j=y}^m \sum_{i|x} \sum_{j|y} [gcd(i,j) == 1]

将枚举因数的位置调换一下

i=1nj=1mnimj[gcd(i,j)=1]\sum_{i=1}^n \sum_{j=1}^m \left \lfloor \frac{n}{i} \right \rfloor \left \lfloor \frac{m}{j} \right \rfloor [gcd(i,j) = 1]

然后我们设

f(x)=i=1nj=1mnimj[gcd(i,j)=x]f(x) = \sum_{i=1}^n \sum_{j=1}^m \left \lfloor \frac{n}{i} \right \rfloor \left \lfloor \frac{m}{j} \right \rfloor [gcd(i,j) = x]

再设

g(x)=xdf(d)g(x) = \sum_{x|d} f(d)

等式展开

g(x)=i=1nj=1mnimj[xgcd(i,j)]g(x) = \sum_{i=1}^n \sum_{j=1}^m \left \lfloor \frac{n}{i} \right \rfloor \left \lfloor \frac{m}{j} \right \rfloor [x|gcd(i,j)]

再将xx提出来

g(x)=i=1nxj=1mxnixmjxg(x) = \sum_{i=1}^{\left \lfloor \frac{n}{x} \right \rfloor} \sum_{j=1}^{\left \lfloor \frac{m}{x} \right \rfloor} \left \lfloor \frac{n}{ix} \right \rfloor \left \lfloor \frac{m}{jx} \right \rfloor

拆开来

g(x)=i=1nxnixj=1mxmjxg(x) = \sum_{i=1}^{\left \lfloor \frac{n}{x} \right \rfloor} \left \lfloor \frac{n}{ix} \right \rfloor \sum_{j=1}^{\left \lfloor \frac{m}{x} \right \rfloor}\left \lfloor \frac{m}{jx} \right \rfloor

我们设 h(x)=i=1nxih(x) = \sum_{i=1}^n \left \lfloor \frac{x}{i} \right \rfloor

那么g(x)=h(nx)h(mx)g(x) = h(\left \lfloor \frac{n}{x} \right \rfloor) h(\left \lfloor \frac{m}{x} \right \rfloor)

我们要求的东西是 f(1)f(1) ,然后通过莫比乌斯反演,可以得出我们要求得东西是这个玩意

f(x)=xdμ(dxg(d))f(x) = \sum_{x|d} \mu( \frac{d}{x} g(d))

f(1)=1dnμ(d1)h(nd)h(md)f(1) = \sum_{1|d}^n \mu( \frac{d}{1}) h(\left \lfloor \frac{n}{d} \right \rfloor) h(\left \lfloor \frac{m}{d} \right \rfloor)

μ\mu 是可以预先筛出来的,然后我们可以预处理楚h(i)h(i) 的值,然后就可以数论分块,将每次查询降为 O(n)O(\sqrt{n})

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int maxd = 5e4+10;
int m,n,ans,k,t;
int mu[maxd], prime[maxd],summu[maxd],cnt;
long long f[maxd];
bool vis[maxd];
void init() 
{
    mu[1] = 1;
    for(int i=2;i<maxd;i++)
    {
        if(!vis[i]) prime[++cnt] = i, mu[i] = -1;
        for(int j=1;j<=cnt && i*prime[j] <= maxd;j++)
        {
            vis[i*prime[j]] = 1;
            if(i%prime[j] == 0) break;
            mu[i*prime[j]] = -mu[i]; 
        }
    }
    for(int i=1;i<maxd;i++) summu[i] = summu[i-1] + mu[i];
    for(int i=1;i<=5e4;i++)   //预处理h(i) 
        for(int l=1,r;l<=i;l = r+1) r = i/(i/l),f[i]+=1LL * (r-l+1) * (i/l);
}
long long sum(int m,int n)  //数论分块
{
    long long t = min(n, m), ans = 0;
    for (int i = 1, j = 1; i <= t; i = j + 1) 
    {
        j = min(n / ( n / i), m / (m / i));
        ans += (summu[j] - summu[i - 1]) * f[n / i] * f[m / j];
    }
    return ans;
}
int main()
{
    // freopen("a.in","r",stdin);
    // freopen("k.out","w",stdout);
    init();scanf("%d",&t);
    for(int i=1;i<=t;i++)
    {
        scanf("%d %d",&n,&m);
        printf("%lld\n",sum(n,m));
    }
    return 0;
}