P3327 约数的个数和 [约数函数性质,数论分块]

P3327 约数的个数和

题意

$d(x)$为约数的个数,对于每个询问,回答${\sum }_{i=1}^n{\sum }_{j=1}^{m}d(ij)$.

题解

这个题推得我头皮发麻,然后还没推出来,后来发现要做这题的先知道一个性质:

$d(ij)={\sum }_{x|i}{\sum }_{y|j}[gcd(x,y)=1]$

通过这个性质,我们把原式写成
${\sum }_{i=1}^n{\sum }_{j=1}^m{\sum }_{x|i}{\sum }_{y|j}[gcd(x,y)=1]$

我们知道${\sum }_{d|x}\mu (d)=[x=1]$,代换进去,就得到了:

${\sum }_{i=1}^n{\sum }_{j=1}^m{\sum }_{x|i}{\sum }_{y|j}{\sum }_{d|gcd(x,y)}\mu (d)$

变枚举$i,j$为枚举$x,y$:

${\sum }_{x=1}^n{\sum }_{y=1}^m\lfloor \frac{n}{x}\rfloor \lfloor \frac{m}{y}\rfloor {\sum }_{d|gcd(x,y)}\mu (d)$

再转为枚举$d$,得到:

${\sum }_{d=1}\mu (d){\sum }_{x=1}^{n/d}{\sum }_{y=1}^{m/d}\lfloor \frac{n}{xd}\rfloor \lfloor \frac{m}{yd}\rfloor$

也即

${\sum }_{d=1}\mu (d)({\sum }_{x=1}^{n/d}\lfloor \frac{n}{xd}\rfloor )({\sum }_{y=1}^{m/d}\lfloor \frac{m}{yd}\rfloor )$

记$f(x)={\sum }_{i=1}^x\lfloor \frac{x}{i}\rfloor$,则原式:

${\sum }_{d=1}\mu (d)f(\lfloor \frac{n}{d}\rfloor )f(\lfloor \frac{m}{d}\rfloor )$

若$f(x)$可以$O(1)$查询的话,上面的式子就可以$O(\sqrt{n})$数论分块求出.

显然,$f(x)$可以用$O(n\sqrt{n})$的时间复杂度预处理出来,方法也是数论分块.

代码

// luogu-judger-enable-o2
#include <iostream>
#include <algorithm>
#include <cstring>
#define pr(x) std::cout << #x << ':' << x << std::endl
#define rep(i,a,b) for(int i = a;i <= b;++i)
typedef long long LL;
const int N = 50010;
int n,m,T;
int prime[N+10],mu[N+10],pcnt,zhi[N+10],low[N+10];
void sieve() {
    mu[1] = zhi[1] = 1;
    for(int i = 2;i <= N;++i) {
        if(!zhi[i]) {
            prime[pcnt++] = i;
            mu[i] = -1;
        }
        for(int j = 0;j < pcnt && i * prime[j] <= N;++j) {
            zhi[i*prime[j]] = 1;
            if(i % prime[j] == 0) {
                mu[i*prime[j]] = 0;
                break;
            }
            else{
                mu[i*prime[j]] = -mu[i];
            }
        }
    }
}
LL F[N+10];
int main() {
    std::ios::sync_with_stdio(false);
    std::cin >> T;
    sieve();
    for(int i = 1;i <= N;++i) {
        mu[i] += mu[i-1];
    }
    for(int i = 1;i <= N;++i) {
        for(int x = 1,last;x <= i;x = last+1) {
            last = i/(i/x);
            F[i] += (last-x+1)*(i/x);
        }
    }
    while(T--) {
        std::cin >> n >> m;
        LL ans = 0;
        int lim = n > m?m:n;
        for(int x = 1,nx1,nx2,nxt;x <= lim;x = nxt+1) {
            nx1 = n/(n/x);
            nx2 = m/(m/x);
            nxt = nx1>nx2?nx2:nx1;
            ans += (mu[nxt]-mu[x-1])*F[n/x]*F[m/x];
        }
        std::cout << ans << std::endl;
    }
    return 0;
}