洛谷 P3327 [SDOI2015]约数个数和 (莫比乌斯反演)

题目描述

设d(x)为x的约数个数，给定N、M，求 ∑Ni=1∑Mj=1d(ij)

---

输入输出格式

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

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

---

输入输出样例

输入样例#1：

2
7 4
5 6

输出样例#1：

110
121

---

说明

1<=N, M<=50000

1<=T<=50000

---

题目分析

数论题对我这种智障来说就是要命的。这题我想了一会还是没有想出来。

某Cla口中的水题对蒟蒻我来说都是难题，所谓的好题就是我根本不会做的题。。。

本题必须用到的结论：

d(ij)=∑x|i∑y|j[(x,y)=1]

证明很简单：

将一个数唯一分解为pk11pk22pk33...pknn，则其约数个数为(k1+1)(k2+1)(k3+1)...(kn+1)。因为对于任意i≠j，(pi,pj)=1。

然后考虑一个质数p对d(ij)的贡献。假设i的质因数分解中有k个p，j的质因数分解中有q个p，那么对d(ij)产生的贡献就是k+q+1。接下来这条关键的式子又不容易想到(反正我想不到)：

k+q+1=∑x=0k∑y=0q[(px,py)=1]

按照某Cla的说法，这就相当于一个(k+1)∗(q+1)的矩形，只取了左下角的”L”字型的一块。 (然而我还不知道这有什么用)

然后根据前面的规律，我们设i=pk11pk22pk33...pknn，j=pq11pq22pq33...pqnn，逐一考虑每个质数pi，根据乘法原理，有

d(ij)=∑x1=0k1∑y1=0q1[(px11,py11)=1]∑x2=0k2∑y2=0q2[(px22,py22)=1]...∑xn=0kn∑yn=0qn[(pxnn,pynn)=1]

我们令x=px11px22px33...pxnn，y=py11py22py33...pynn，就得到了一开始的结论

d(ij)=∑x|i∑y|j[(x,y)=1]

其实很好理解，x就代表了pxii的总出现状态，分开的式子相当于分组计数再乘法合并，合起来相当于直接枚举总状态计数，本质就是在统计数对的数量(显然所有质因数不同时出现等价于原数互质)。

---

现在就变成求

∑i=1n∑j=1m∑x|i∑y|j[(x,y)=1]

根据数论中切换枚举次序的套路以及莫比乌斯反演对布尔条件框的替换，我们得到

原式

=∑x=1n∑y=1m[(x,y)=1]⌊nx⌋⌊my⌋

=∑x=1n∑y=1mμ(gcd(x,y))⌊nx⌋⌊my⌋

将gcd用枚举去掉

=∑x=1n∑y=1m∑k|x,k|yμ(k)⌊nx⌋⌊my⌋

最后将约数k的枚举提到外面，内部枚举k的互补约数

=∑k=1min(n,m)μ(k)∑x=1⌊nk⌋⌊nkx⌋∑y=1⌊mk⌋⌊mky⌋

然后我们设g(n)=∑ni=1⌊ni⌋，那么μ(k)右边的部分变成g(⌊nk⌋)g(⌊mk⌋)。

这样对于连续一段k，g函数的参数相同，函数值也相同，于是预处理μ的前缀和，将第一个枚举改成下底函数分块。

g函数也用下底函数分块来去处理，算出所有g的时间是O(nn√)。

还有一种做法，就是g函数其实是约数个数的前缀和(g(n)=∑ni=1∑d|i1=∑nd=1⌊nd⌋)。做线性筛的同时能顺带求出约数个数，最后做一次前缀和就行了。(某Cla想到的)

---

代码

#include <cstdio>#include <cstdlib>#include <cstring>#include <algorithm>#include <cmath>#include <iostream>#define N 50010using namespace std;typedef long long LL;int T, n, m;LL g[N], Ans;int prime[N], miu[N], cnt;bool vis[N];void Get_miu(){    for(int i = 1; i < N; i++){      int last;      for(int j = 1; j <= i; j = last+1){        last = i/(i/j);        g[i] += (LL)(i/j)*(LL)(last-j+1);      }    }    vis[1] = true;    miu[1] = 1;    for(int i = 2; i < N; i++){      if(!vis[i]){        prime[++cnt] = i;        miu[i] = -1;      }      for(int j = 1; j <= cnt && i * prime[j] < N; j++){        vis[i * prime[j]] = true;        if(i % prime[j] == 0){          miu[i * prime[j]] = 0;          break;        }        else  miu[i * prime[j]] = -miu[i];      }    }    for(int i = 1; i < N; i++)  miu[i] += miu[i-1];}int main(){    Get_miu();    scanf("%d", &T);    while(T --){      scanf("%d%d", &n, &m);      if(n > m)  swap(n, m);      Ans = 0;      int last;      for(int i = 1; i <= n; i = last+1){        last = min(n/(n/i), m/(m/i));        Ans += (LL)(miu[last] - miu[i-1]) * g[n/i] * g[m/i];      }      printf("%lld\n", Ans);    }    return 0;} 

---

给所有的现充施以破坏的铁锤吧!