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

题目传送门


Solution

首先将题目转换,容易发现一个点(x,y)kgcd(x,y)1

然后就变成了求

i=1nj=1m2gcd(i,j)1

它等于
(2i=1nj=1mgcd(i,j))nm

于是我们只需求
i=1nj=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)r1kr

观察最右边的公式,似乎和
(a+b)n=r=0nCrnarbnr

这个公式形式相同。
所以
d|nμ(d)=(1+1)k=0(n>1)

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

因为F(n)=d|nf(d)=fI,μI=e,所以Fμ=fIμ=fe=f, 即f(n)=d|nμ(d)f(nd)

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

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

然后开始变形:
原式=mg=1gngi=1mgj=1e(gcd(i,j))

=mg=1gngi=1mgj=1t|gcd(i,j)μ(t)

=mg=1gmgt=1μ(t)ngtmgt

s=gt
=ms=1nsmst|sμ(t)st

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

μid=μϕI=ϕe=ϕ

接下来:

Ans=s=1mnsmsϕ(s)

大道至简。

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

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

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

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

然后对于nimi一对一组合,时间复杂度是O(n+m)。处理多组询问就没问题了。

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

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

f(i)=nid=1g(di)=>g(i)=nid=1f(di)μ(d)

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

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

首先,设f(d)表示有多少个gcd(i,j)=dg(d)表示有多少个gcd(i,j)d的倍数,g(d)=ndmd

因为g(d)=ndi=1f(di),

则根据公式得
f(d)=ndi=1g(di)μ(i)


f(d)=ndi=1μ(i)ndimdi

而答案就等于
ng=1gndi=1μ(i)ndimdi(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[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;
}

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

发布了177 篇原创文章 · 获赞 109 · 访问量 6万+
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 精致技术 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览