「SDOI2015」约数个数和
知识点: 莫比乌斯反演
题意简述
\(T\) 次询问,每次询问给定 \(n,m\)。
定义 >\(\operatorname{d}(i)\) 为 \(i\) 的约数个数,求:\[\sum_{i=1}^{n}\sum_{j=1}^m\operatorname{d}(ij) \]\(1<T,n\le 5\times 10^4\)。
分析题意
一个结论:
证明
先考虑 \(i = p^a\),\(j=p^b(p\in \mathrm{primes})\) 的情况,有:
对于等式左侧,\(p^{a+b}\) 的约数个数为 \(a+b+1\)。
对于等式右侧,在保证 \(\gcd(x,y)=1\) 成立的情况下,有贡献的数对 \((x,y)\) 只能是下列三种形式:
- \(x>0,y-0\),\(x\) 有 \(a\) 种取值方法。
- \(x=0,y>0\),\(y\) 有 \(b\) 种取值方法。
- \(x=0,y=0\)。
则等式右侧贡献次数为 \(a+b+1\) 次,等于 \(p^{a+b}\) 的约数个数。
则当 \(i = p^a\),\(j=p^b(p\in \mathrm{primes})\) 时等式成立。
又不同质因子间相互独立,上述情况可拓展到一般的情况。
对 \(\operatorname{d}(i,j)\) 进行化简,代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\),原式等价于:
调换枚举顺序,先枚举 \(d\),原式等价于:
把各项中的 \(d\) 提出来,消去整除的条件,原式等价于:
将 \(\operatorname{d}(ij)\) 代回原式,原式等价于:
调换枚举顺序,先枚举 \(d\),原式等价于:
把 \(i,j\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}, \frac{j}{d}\),消去整除的条件,原式等价于:
考虑预处理 \(S(x) = \sum\limits_{i=1}^{x}\operatorname{d}(i)\),则原式等价于:
线性筛预处理 \(\mu,\operatorname{d}\),数论分块求解即可,复杂度 \(O(n+T\sqrt{n})\)。
爆零小技巧
注意筛法中出现平方因子的处理。
代码实现
//知识点:莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const int kMaxn = 5e4 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], num[kMaxn], d[kMaxn]; //num 为最小质因子的次数
ll summu[kMaxn], sumd[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void Getmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
vis[1] = true;
mu[1] = d[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
num[i] = 1;
d[i] = 2;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
num[i * p[j]] = num[i] + 1;
d[i * p[j]] = 1ll * d[i] / num[i * p[j]] * (num[i * p[j]] + 1ll);
break;
}
mu[i * p[j]] = - mu[i];
num[i * p[j]] = 1;
d[i * p[j]] = 2ll * d[i]; //
}
}
for (int i = 1; i <= lim_; ++ i) {
summu[i] = summu[i - 1] + mu[i];
sumd[i] = sumd[i - 1] + d[i];
}
}
//=============================================================
int main() {
Euler(kMaxn - 10);
int T = read();
while (T --) {
int n = read(), m = read(), lim = std :: min(n, m);
ll ans = 0ll;
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n / (n / l), m / (m / l));
ans += 1ll * (summu[r] - summu[l - 1]) * sumd[n / l] * sumd[m / l]; //
}
printf("%lld\n", ans);
}
return 0;
}
/*
1
32 43
*/
/*
15420
*/