「BZOJ 3994」「SDOI 2015」约数个数和「莫比乌斯反演」

题意

设(d(x))为(x)的约数个数，求(sum_{i=1}^{n}sum_{j=1}^{m}d(ij))。

题解

首先证个公式：

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

可以这么考虑：利用唯一分解定理把(i,j)分解，即：

$i=prod_{k = 1}^{m} p_k^{c_k},j=prod_{k=1}^m p_k^{d_k} $

那等式左边显然为(prod(c_k+d_k+1))，

然后考虑等式右边在干什么事情：约数的最大公约数为(1)。

我们把(x，y)分解以后(p_k)的指数分别记为(a_k,b_k)

这就是说对于每个(p_k)，(a_k,b_k)最多有一个大于(0)。

那每一位分成三种情况，(a_k=b_k=0)；(a_kin[1,c_k])且(b_k=0)；(a_k=0)且(b_kin[1,d_k])

然后根据乘法原理乘起来，得到的结果和左边一样，就证出来了

然后带进去化简式子：（假设(nleq m)）

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

[=sum_{x=1}^n sum_{y=1}^m [gcd(x, y)=1]lfloor frac{n}{x} floorlfloor frac{m}{y} floor ]

[=sum_{x=1}^n sum_{y=1}^m (sum_{d|x,d|y} mu(d) )lfloor frac{n}{x} floorlfloor frac{m}{y} floor ]

[=sum_{d=1}^n mu(d) sum_{x=1}^{lfloor frac{n}{d} floor} sum_{y=1}^{lfloor frac{m}{d} floor} lfloor frac{n}{xd} floor lfloor frac{m}{xd} floor ]

这里利用整除的小性质：(lfloorfrac{lfloor frac{a}{b} floor}{c} floor=lfloorfrac{lfloorfrac{a}{c} floor}{b} floor)

[=sum_{d=1}^n mu(d) sum_{x=1}^{lfloor frac{n}{d} floor} sum_{y=1}^{lfloor frac{m}{d} floor} lfloor frac{lfloorfrac{n}{d} floor}{x} floor lfloor frac{lfloorfrac{m}{d} floor}{y} floor ]

记(f(n)=sum_{i=1}^nlfloor frac{n}{i} floor)，则答案为(sum_{d=1}^n mu(d)f(lfloor frac{n}{d} floor) f(lfloor frac{m}{d} floor))

这里如果(f)通过数论分块预处理，虽然可以过但是显然不够优美，可以换个角度考虑

这个式子的意义是，对于每个(i)统计它(leq n)的所以倍数。那么每个数会被统计(d(i))次（(d)为约数个数），因此：

(f(n)=sum_{i=1}^n d(i))

就可以线性筛预处理(d)的前缀和，还要预处理(mu)的前缀和，然后每次数论分块。

#include <algorithm>
#include <cstdio>
using namespace std;

typedef long long ll;

const int N = 5e4 + 10;

ll sum[N], mu[N], d[N];
int p[N], tot, e[N];
bool tag[N];

void init(int n) {
	mu[1] = d[1] = 1;
	for(int i = 2; i <= n; i ++) {
		if(!tag[i]) {
			p[tot ++] = i; mu[i] = -1; d[i] = 2; e[i] = 1;
		}
		for(int j = 0; j < tot && i * p[j] <= n; j ++) {
			tag[i * p[j]] = 1;
			if(i % p[j] == 0) {
				mu[i * p[j]] = 0;
                e[i * p[j]] = e[i] + 1;
                d[i * p[j]] = d[i] / (e[i] + 1) * (e[i] + 2);
				break ;
			}
			mu[i * p[j]] = - mu[i];
            e[i * p[j]] = 1;
            d[i * p[j]] = d[i] * 2;
		}
        mu[i] += mu[i - 1];
        d[i] += d[i - 1];
	}
}

ll calc(int n, int m) {
	ll ans = 0;
	for(int i = 1, j; i <= n; i = j + 1) {
		j = min(n / (n / i), m / (m / i));
		ans += (mu[j] - mu[i - 1]) * d[n / i] * d[m / i];
	}
	return ans;
}

int main() {
	init(50000);
	int test; scanf("%d", &test);
	for(int n, m; test --; ) {
		scanf("%d%d", &n, &m);
		if(n > m) swap(n, m);
		printf("%lld
", calc(n, m));
	}
	return 0;
}