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 100010

using namespace std;

typedef long long LL;
int n, m, cnt;
bool vis
;
int prime
;
LL phi
, sum
;

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;
}

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