题目地址：洛谷：【P3327】[SDOI2015]约数个数和 – 洛谷、BZOJ：Problem 3994. — [SDOI2015]约数个数和

题目描述

设d(x)为x的约数个数，给定N、M，求\sum^N_{i=1}\sum^M_{j=1} d(ij)。

输入输出格式

输入格式：
输入文件包含多组测试数据。第一行，一个整数T，表示测试数据的组数。接下来的T行，每行两个整数N、M。

输出格式：
T行，每行一个整数，表示你所求的答案。

输入输出样例

输入样例#1：

2
7 4
5 6

输出样例#1：

110
121

说明

1<=N, M<=50000
1<=T<=50000

题解

参考资料：BZOJ 3994 – [SDOI2015]约数个数和 | Sengxian’s Blog
首先我们要知道一个事，d(xy) = \sum_{i|x} \sum_{j|y} [\mathrm{gcd}(i, j) == 1]。有了它，我们就可以对要求的式子进行变形，首先是最开始的版本（[表达式]表示表达式为真时值为1，否则值为0）
\sum^n_{i=1} \sum^m_{j=1} d(ij)
代换掉d函数
\sum^n_{i=1} \sum^m_{j=1} \sum_{a|i} \sum_{b|j} [\mathrm{gcd}(a, b) == 1]
我们知道有\sum_{x|n} \mu(x) = [n == 1]，因此可以带进上面的条件
\sum^n_{i=1} \sum^m_{j=1} \sum_{a|i} \sum_{b|j} \sum_{x|\mathrm{gcd}(a, b)} \mu(x)
考虑枚举x统计每个\mu(x)的数量，上式可变换为
\sum^n_{i=1} \sum^m_{j=1} \sum_{x|i, x|j} \mu(x) d(\frac{i}{x}) d(\frac{j}{x})
把\mu(x)项放在前面
\sum^n_{i=1} \mu(x) \sum_{x|i} d(\frac{i}{x}) \sum_{x|j} d(\frac{j}{x})
后面两项换个元
\sum^n_{i=1} \mu(x) \sum_{i=1}^{\lfloor \frac{n}{x} \rfloor} d(i) \sum_{i=1}^{\lfloor \frac{m}{x} \rfloor} d(j)
我们如果能预处理出d(x), \mu(x)，再用一波除法分块，就可以把原来O(n^2)的暴力求值变成现在的O(n \sqrt{n})优秀算法。而d(x)是一个积性函数，可以跟\mu(x)一起筛出来。
我们考虑如何筛d(x)。两两互质的情况可以直接乘起来，当遇到不互质的情况，prime[j]一定是i中最小的那个质因子，我们可以求出每个数最小质因子的次数，从而获取其他质因子的选择方案数之乘积，推得i*prime[j]这个数字的d值。这是由于如果能够将x表示为p_1^{c_1}p_2^{c_2} \cdots p_k^{c_k}的形式，则d(x) = \prod_{i=1}^n (c_i + 1)。

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>  

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) 
        ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    register char c = fgc();
    while(!isdigit(c)) {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(isdigit(c)) {
        res = (res << 1) + (res << 3) + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 50005;

bool notprime[MAXN];
int prime[MAXN], ptot, mu[MAXN], d[MAXN], cnt[MAXN];

inline void sieve() {
    mu[1] = d[1] = 1; notprime[1] = true;
    for(int i = 2; i < MAXN; i++) {
        if(!notprime[i]) {
            prime[++ptot] = i; mu[i] = -1; d[i] = 2; cnt[i] = 1;
        }
        for(int j = 1; j <= ptot && i * prime[j] < MAXN; j++) {
            int t = i * prime[j];
            notprime[t] = true;
            if(i % prime[j]) {
                mu[t] = -mu[i]; d[t] = d[i] * 2; cnt[t] = 1;
            } else {
                mu[t] = 0; cnt[t] = cnt[i] + 1; d[t] = d[i] / (cnt[i] + 1) * (cnt[t] + 1);
            }
        }
    }
    for(int i = 2; i < MAXN; i++) {
        mu[i] += mu[i - 1];
        d[i] += d[i - 1];
    }
}

inline LL cal(int n, int m) {
    if(n > m) std::swap(n, m);
    LL res = 0;
    for(int i = 1; i <= n;) {
        int j = std::min(n / (n / i), m / (m / i));
        res += 1ll * (mu[j] - mu[i - 1]) * d[n / i] * d[m / i];
        i = j + 1;
    }
    return res;
}

int T, n, m;

int main() {
    sieve();
    T = readint();
    while(T--) {
        n = readint(); m = readint();
        printf("%lld\n", cal(n, m));
    }
    return 0;
}