BZOJ 2005: [Noi2010]能量采集(莫比乌斯反演)

题目传送门

---

Solution

首先将题目转换，容易发现一个点(x,y)的k为gcd(x,y)−1。

然后就变成了求

∑i=1n∑j=1m2∗gcd(i,j)−1

它等于

(2∗∑i=1n∑j=1mgcd(i,j))−n∗m

于是我们只需求

∑i=1n∑j=1mgcd(i,j)

这是一个莫比乌斯反演的经典问题。

在解决这个问题之前，先来复习一下莫比乌斯反演的有关知识：

设F(n)和f(n)是两个定义在正整数域上的数论函数，且F(n)=∑d|nf(d),那么有

f(n)=∑d|nμ(d)F(nd)

f(n)=∑n|dμ(dn)F(d)

如何证明？下面引用LYW大神的证明：

倒数第二步如何推到最后一步呢？我们知道狄利克雷卷积的元函数e(n)=[n=1]。而μ∗I=e，就是∑d|nμ(d)=e(n)。

这是莫比乌斯函数的一个重要性质。
如何证明？

n=1明显成立。
其他情况把n分解质因数得到：
n=px11px22px33...pxkk
因为d是n的约数，那么d的所有质因数也只有p1到pk。
因为只要有一个质因数的指数大于1，μ(d)就是0，所以系数只可能是0或者1。
如果有r个质因子的指数是1，那么μ(d)=(−1)r
在k个质因子中选出r个，所以d有Crk个。
那么总的公式就是：

∑d|nμ(d)=∑r=0kCrk∗(−1)r

我们再乘上一个1，得到

∑d|nμ(d)=∑r=0kCrk∗(−1)r∗1k−r

观察最右边的公式，似乎和

(a+b)n=∑r=0nCrn∗ar∗bn−r

这个公式形式相同。
所以

∑d|nμ(d)=(−1+1)k=0(n>1)

得到这个后，其实还有另一种证明前面的公式的方法，这里借用了狄利克雷卷积的一些变换(狄利克雷卷积的定义等省略)。

因为F(n)=∑d|nf(d)=f∗I,μ∗I=e,所以F∗μ=f∗I∗μ=f∗e=f, 即f(n)=∑d|nμ(d)f(nd)。

到这里，我们简要复习了一下莫比乌斯反演。

回到题目(假设m<n)，先枚举gcd，然后就是
∑mg=1g∑⌊ng⌋i=1∑⌊mg⌋j=1e(gcd(i,j))
这里的做法类似于统计gcd为g的数对的数量。

然后开始变形：
原式=∑mg=1g∑⌊ng⌋i=1∑⌊mg⌋j=1e(gcd(i,j))

=∑mg=1g∑⌊ng⌋i=1∑⌊mg⌋j=1∑t|gcd(i,j)μ(t)

=∑mg=1g∑⌊mg⌋t=1μ(t)⌊ngt⌋⌊mgt⌋

令 s=gt，
=∑ms=1⌊ns⌋⌊ms⌋∑t|sμ(t)st

到这里，引入第二条式子ϕ(n)=∑d|nμ(d)nd
这个证明就用∑d|nϕ(d)=n(想想就知道),即ϕ∗I=id, 然后就

μ∗id=μ∗ϕ∗I=ϕ∗e=ϕ

接下来：

Ans=∑s=1m⌊ns⌋⌊ms⌋ϕ(s)

大道至简。

在此之前，时间复杂度为O(mlogm)
后来就是O(m)

然而多组数据时，我们要加上分块优化。就是将下底函数进行分组。

当我们计算∑ni=1⌊ni⌋时，可以在根号n的时间内搞定。

因为随着i长，i<=n−√，有n−√个取值，当i变大，⌊ni⌋也只有n−√个取值，所以最多是2n−√个取值。

然后对于⌊ni⌋与⌊mi⌋一对一组合，时间复杂度是O(n−√+m−−√)。处理多组询问就没问题了。

对于同样的⌊ni⌋与⌊mi⌋，我们预处理出ϕ的前缀和，然后就可以一起算了。这就是算出每一块的开头结尾，同一块的一起算的的分块优化方法。

说了那么多，不如多熟记几条公式。以下必记：

f(i)=∑⌊ni⌋d=1g(d∗i)=>g(i)=∑⌊ni⌋d=1f(d∗i)μ(d)

f(x)=∑x|d,d<=ng(d)=>g(x)=∑x|d,d<=nf(d)μ(dx)

我们可以直接使用上面的公式直接跳到中间的某步运算过程：

首先，设f(d)表示有多少个gcd(i,j)=d，g(d)表示有多少个gcd(i,j)是d的倍数，g(d)=⌊nd⌋⌊md⌋

因为g(d)=∑⌊nd⌋i=1f(d∗i),

则根据公式得
f(d)=∑⌊nd⌋i=1g(d∗i)μ(i)

即
f(d)=∑⌊nd⌋i=1μ(i)⌊ndi⌋⌊mdi⌋

而答案就等于
∑ng=1g∑⌊nd⌋i=1μ(i)⌊ndi⌋⌊mdi⌋(n<m)

这里就到了“令s=gt”之前了，然后就可以推下去，记住几条常见的公式推得会更快。

好了，就到这里了(数学公式敲到我快吐血)。

---

Code

#include <iostream>#include <cstdio>#include <cstdlib>#include <cstring>#include <cmath>#include <algorithm>#define N 100010using namespace std;typedef long long LL;int n, m, cnt;bool vis[N];int prime[N];LL phi[N], sum[N];void Get_Phi(){    vis[1] = true;    phi[1] = 1ll;    for(int i = 2; i <= n; i++){      if(!vis[i]){        prime[++cnt] = i;        phi[i] = i - 1;      }      for(int j = 1; j <= cnt && i * prime[j] <= n; j++){        vis[i * prime[j]] = true;        if(i % prime[j] == 0){          phi[i * prime[j]] = phi[i] * prime[j];          break;        }        else  phi[i * prime[j]] = phi[i] * (prime[j] - 1);      }    }    sum[0] = 0ll;    for(int i = 1; i <= n; i++)  sum[i] = sum[i-1] + phi[i];}LL Get_Ans(){    LL res = 0ll;    int last;    for(int i = 1; i <= n; i = last+1){      last = min(n/(n/i), m/(m/i));      res += (LL)(n/i) * (LL)(m/i) * (sum[last] - sum[i-1]);    }    return res * 2 - (LL)n * (LL)m; }int main(){    freopen("bzoj2005.in", "r", stdin);    freopen("bzoj2005.out", "w", stdout);    scanf("%d%d", &n, &m);    if(n > m)  swap(n, m);    Get_Phi();    printf("%lld\n", Get_Ans());    return 0;}

---

多希望我的人生是大梦一场，梦醒时分你依旧在我身旁。