莫比乌斯反演

 前言：

 　　理论部分没写完+很不严谨……民那桑将就着看吧……以后有空补充……

　　数论部分详细参考任何一本数论教材。

　　同理，组合数学部分详细参考任何一本组合数学教材。

---

---

mobius函数 μ(x) 定义为

　　          = (-1)^r  【n = p1×p2×p3... ×pr，其中pi为不同的素数】

　 μ(x)　　 = 1 　　【n = 1】

　　　　　  = 0　　 【其他，比方说 n = 4 = 2²】

 

---

---

 

积性函数：

　　对于任意互素的正整数a,b，均有f(ab) = f(a) f(b)， 则f为积性函数。

如：φ(x) 欧拉函数，μ(x) mobius函数。

 

完全积性：

　　对于任意正整数a,b，均有f(ab) = f(a) f(b)， 则f为完全积性函数。

如：Legendre符号——

（此处为本人不严谨地说明…… ）

---

---

重要结论：

1. Σ_d|nφ(d)=n

2. Σ_d|nμ(d) = 1 (n == 1)

　[Σ_d|nμ(d) = 0 (n > 1)]

---

---

 

mobius反演公式：

F(n) = Σ_d|nf(d) ↔ f(n) = Σ_d|nμ(n/d)F(d) 【注：这里的F,f为数论函数……不是随便什么函数都可以的啊…… 不过一般在竞赛中用到的貌似就是欧拉函数+莫比乌斯函数……】

---

---

碎碎念：

其实mobius函数-1,1什么的就是容斥里面+x1 -x2 ...的那个系数，当然我的意思是如果和容斥挂钩的话——因为我本人一开始草草一看就是这么理解的……所以偏好把mobius看成容斥——封装得很好的容斥。

当然，其实完全不必把mobius看成是容斥的一部分——这样感觉好像它只是容斥的一个特例的样子。

但是，就算是作为数论函数来看，mobius函数也是很重要的，意义远不止是容斥原理。

比方说积性啊什么的。

但是有些acm题如果用mobius函数的话就比较不用想容斥里面+-什么的了。比方说求1≤a≤m && 1≤b≤n &&gcd(a,b) == 1的a,b的组合数的这种类型的题目。

 

---

---

据说mobius在acm竞赛中常用的写代码技巧是：线性筛选+分块

详细参考——"贾志鹏线性筛"。

先奉上普通的方法（就是根据定义写的）：

int mobius(int n) {
  int miu = 0;
  if(n == 1) return 1;
  for(int i = 2; i < sqrt(n); i++) {
    if(n % i == 0) {
      m *= -1;
      int flag = 0;
      while(n % i == 0) {
        n /= i;
        flag = 1;
      }
      if(flag)return 0;
    }
  }
  if(n > 1) return m *= -1;
  return m;
}

 

然后是用筛选法做的（自己根据定义写的……暂时没发现bug，而且还用此过了HDU 1695）

const int maxn = 100011;
int miu[maxn];
int prime[maxn];
void mobius() {
    memset(prime,0,sizeof prime);
    miu[1] = 1;
    prime[1] = 0;
    for(int i = 2; i<maxn; i++) {
        if(prime[i] == 0) {
            long long j = (long long )i*(long long )i;
            for(; j<maxn; j+=i) {
                if(prime[j] == 0)
                    prime[j] = i;
            }
            miu[i] = -1;
            continue;
        }

        int n = i;
        miu[i] = 1;
        while(prime[n]!=0) {
            int tp = prime[n];
            int flag = 0;
            while(n % tp == 0) {
                flag++;
                if(flag==2) {
                    miu[i] = 0;
                    break;
                }
                n /= tp;
            }
            miu[i] *= -1;
            if(flag == 2 ) break;
        }
        if(n>1)miu[i] *= -1;
    }
}

 

 

例子：

---

---

HDU 1695

题意:已知a,b,c,d,k，求有多少x,y满足a≤x≤b, c≤y≤d,使得gcd(x,y) = k。

　　然后本题a = c = 1.

　　公式计算部分参考"贾"论文…… [电脑上打数学公式简直作死……]

1. 可以用容斥+欧拉函数

 

2. 可以用莫比乌斯（mobius）反演

#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
const int maxn = 100011;
int miu[maxn];
int prime[maxn];
void mobius() {
  memset(prime,0,sizeof prime);
  miu[1] = 1;
  prime[1] = 0;
  for(int i = 2; i<maxn; i++) {
    if(prime[i] == 0) {
      long long j = (long long )i*(long long )i;
      for(; j<maxn; j+=i) {
        if(prime[j] == 0)
          prime[j] = i;
      }
      miu[i] = -1;//cout<<"miu["<<i<<"]: "<<miu[i]<<endl;
      continue;
    }

    int n = i;
    miu[i] = 1;
    //if(i == 30)cout<<"i:"<<i<<" prime["<<n<<"]: "<<prime[n]<<endl;
    while(prime[n]!=0) {
//            if(i == 30)
//                    cout<<"n:"<<n<<endl;
      int tp = prime[n];
      int flag = 0;
      while(n % tp == 0) {
        flag++;
        if(flag==2) {
          miu[i] = 0;
          break;
        }
        n /= tp;
      }
//            if(i == 30)
//                    cout<<"n:"<<n<<"flag:"<<flag<<endl;
      miu[i] *= -1;
      if(flag == 2 ) break;
    }
    if(n>1)miu[i] *= -1;
    //cout<<"miu["<<i<<"]: "<<miu[i]<<endl;
  }
}
int main() {
#ifdef LOCALL
  //freopen("in","r",stdin);
  //freopen("out","w",stdout);
#endif
  int a,b,c,d,k;
  mobius();
  //cout<<miu[100000];
  int t;
//    while(cin>>t)
//    cout<<"miu["<<t<<"]:"<<miu[t]<<endl;
  cin>>t;
  for(int kase = 1; kase < t+1; kase ++) {
    cin>>a>>b>>c>>d>>k;
    if(k == 0) {
      cout<<"Case "<<kase<<": "<<"0
";
      continue;
    }
    long long ans = 0,repeat = 0;
    b /= k;
    d /= k;
    int a1 = b<d?b:d;
    for(int i = 1; i<=a1; i++) {
      ans += (long long)miu[i]*(b/i)*(d/i);
    }
    for(int i = 1; i<=a1; i++) {
      repeat += (long long)miu[i]*(a1/i)*(a1/i);
    }

    cout<<"Case "<<kase<<": "<<ans-repeat/2<<"
";
  }
  return 0;
}

---

BZOJ 2301 

题意：题意:已知a,b,c,d,k，求有多少x,y满足a≤x≤b, c≤y≤d,使得gcd(x,y) = k。

100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000

和上一题比起来就是a,c不一定等于1。。。

首先，如果我们按照上题那样做，在思路上完全没有问题，然后BZOJ这题给了50s，其实我觉得怎么想都是想给这种方法一个活路……

          公式为  answer = [1,b; 1, d] - [1,a-1; 1, d] - [1,b; 1, c-1] + [1,a-1; 1,c-1]  

　　　 【[a,b;c,d]的意思为 x∈[a,b] && y∈[c,d]的gcd满足条件的组合数】

分析一下复杂度： 1. 筛选法初始化mobius 大概O(nlog(n)log(log(n)))……额 这是Eratosthenes筛选法的复杂度…… 

　　　　　　　　 2. 最坏可能是 求x,y都属于[1 - n]的可能组合数 那就是 O(n)

　　　　　　　　 3. 可能有多个测试数据…… 然后这个数据个数也有n这么多……

　　　　　　　　话说 n 是 50000啊！

　　　　也就是这样算，去除零头什么的，大概是 O(n*n) = O(25×10⁸) ——  10⁶≈ 1s， 2500*10⁶ ≈ 2500s……

……这…… 果断活不成……  

但是我还是试了试…… 果断超时

贴个超时代码~

#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
typedef long long LL;
const int maxn = 100011;
int p,nn,plen;
int miu[maxn];
int pri[maxn];
int prime[maxn];

void mobius() {
  memset(pri,0,sizeof pri);
  memset(prime,0,sizeof prime);

  miu[1] = 1;
  pri[1] = 0;
  plen = 0;
  for(int i = 2; i<maxn; i++) {
    if(pri[i] == 0) {
      prime[plen++] = i;
      LL j = (LL )i*(LL )i;
      for(; j<maxn; j+=i) {
        if(pri[j] == 0)
          pri[j] = i;
      }
      miu[i] = -1;
      continue;
    }

    int n = i;
    miu[i] = 1;
    while(pri[n]!=0) {
      int tp = pri[n];
      int flag = 0;
      while(n % tp == 0) {
        flag++;
        if(flag==2) {
          miu[i] = 0;
          break;
        }
        n /= tp;
      }
      miu[i] *= -1;
      if(flag == 2 ) break;
    }
    if(n>1)miu[i] *= -1;
  }
}

int main() {
#ifdef LOCALL
  freopen("in","r",stdin);
  //freopen("out","w",stdout);
#endif

  mobius();
  //cout<<miu[100000];
  int t;

  scanf("%d",&t);
  int a,b,c,d,k;
  while(t--) {
    scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
    LL ans = 0, rep = 0;
    int a1 = b/k<d/k?b/k:d/k;
    for(int i = 1; i <= a1; i ++) {
      ans += (long long)miu[i]*(b/k/i)*(d/k/i);
    }
    a1 = a/k<c/k?a/k:c/k;
    for(int i = 1; i <= a1;i++){
      ans += (long long)miu[i]*((a-1)/k/i)*((c-1)/k/i);
    }

    a1 = (a-1)/k<d/k?(a-1)/k:d/k;

    for(int i = 1; i <= a1; i++) {
      rep += (long long)miu[i]*((a-1)/k/i)*(d/k/i);
    }
    a1 = (c-1)/k<b/k?(c-1)/k:b/k;
    for(int i = 1; i <= a1;i++){
      rep += (long long)miu[i]*((c-1)/k/i)*(b/k/i);
    }
    printf("%lld
",(ans - rep));
  }
  return 0;
}

然后就要用到  此报告中的分块优化了……

详细参考POI XIV Stage.1 Queries Zap 解题报告：http://wenku.baidu.com/view/fbe263d384254b35eefd34eb.html

　　分块优化的思想说起来高大上，其实是这样的：

　　1. 注意平常的mobius，在进行反演之后主要还是进行枚举：

　　　　如求x∈[1,a], y∈[1,b]的gcd(x,y) = 1的全部组合数 （假设 a>b）

for(int i = 1; i<=b; i++) {
       ans += (long long)miu[i]*(a/i)*(b/i);
}

　　我们需要枚举全部小于b的值 —— 用容斥的说法就是，需要枚举全部gcd(x,y)大于1小于d的可能组合。

　　

　　2. 我们来看一个图：

　　

　　当gcd(x,y) = d (d > 1)时，我们可以看到，当gcd(x,y) 落在 [d, new_d]这段里时，a/gcd(x,y) 都等于 a/d。

　　就是凭借这个，我们可以进行分块优化。

　　因为 i = d 到 i = new_d 的这段都对应于(a/i)*(b/i)，所以不用一个个求了，直接可以提取出"(a/i)*(b/i)"公因式，这样的话，只需要把μ(d)到μ(new_d)的值累加就好了，

因为需要μ的累加值，一开始我们也求出来[1,n]的μ值了，所以用前缀和的方法把所有的μ值的累加和先存起来……

　　嘛…… 这样优化了以后，算算复杂度：(具体复杂度推导过程，参看那篇POI解题报告+自己推算)

　　　　　　　　 O(n*(2√(n/k) + 2√(n/k)) + n * √n) ≈ O(250000* √50000) ≈ O(250000* √50000) ≈ O(55*10⁶) ≈ 50s

 　 50s也就刚好差不了…… 剧透一下：最后跑出来也就11+s而已…… 一开始还以为哪写差了，怎么这么久……（本人代码能力差……）然后随便搜了一个代码复制上去……居然也有10+s ……心理平衡了……

　　AC代码……

#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<cmath>
#define minum(a,b) (a)>(b)?(b):(a)
using namespace std;
typedef long long LL;
const int maxn = 50011;
int miu[maxn];
int pri[maxn];
int p[maxn];

void mobius() {
  memset(pri,0,sizeof pri);

  miu[1] = 1;
  pri[1] = 0;
  for(int i = 2; i<maxn; i++) {
    if(pri[i] == 0) {
      LL j = (LL )i*(LL )i;
      for(; j<maxn; j+=i) {
        if(pri[j] == 0)
          pri[j] = i;
      }
      miu[i] = -1;
      continue;
    }

    int n = i;
    miu[i] = 1;
    while(pri[n]!=0) {
      int tp = pri[n];
      int flag = 0;
      while(n % tp == 0) {
        flag++;
        if(flag==2) {
          miu[i] = 0;
          break;
        }
        n /= tp;
      }
      miu[i] *= -1;
      if(flag == 2 ) break;
    }
    if(n>1)miu[i] *= -1;
  }
}
LL sum(int a,int b) {
  int a1 = minum(a,b);
  LL ans = 0;
  int new_i;
  for(int i = 1,new_i = 0; i <= a1; i = new_i + 1) {
    new_i = minum(a/(a/i),(b/(b/i)));
    ans += (LL)(p[new_i] - p[i-1])*(a/i)*(b/i);
  }
  return ans;
}
int main() {
#ifdef LOCALL
  freopen("in","r",stdin);
  //freopen("out","w",stdout);
#endif
  mobius();
  p[0] = 0;
  for(int i = 1; i<maxn; i++) p[i] = p[i-1] + miu[i];
  int t;
  scanf("%d",&t);
  int a,b,c,d,k;
  while(t--) {
    scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
    printf("%lld
",sum(b/k,d/k) - sum((a-1)/k,d/k) - sum(b/k,(c-1)/k) + sum((a-1)/k,(c-1)/k));
  }
  return 0;
}

　　

(BTW, BZOJ 的 C++交cin | cout 会 RE …… )

---

---

还有啊……mobius反演可以用于求可重复的圆周排……

不过……可能由于我代码能力的问题…… 用mobius写就总是超时……待我写过了再来补充吧……呵呵呵呵呵……

还有啊…… 写这么多……好累啊……