莫比乌斯反演小结

关于莫比乌斯反演入门，我是参考了这篇博客

http://www.cnblogs.com/chenyang920/p/4811995.html

最常用到的一个公式为

原式 :  G(n)=sigma(F(d))  (其中n|d,d<=N)  
反演公式:   F(n)=sigma(U(d/n)*G(d))  (其中n|d,d<=N)

通常用来计算与GCD计数有关的问题。

1.BZOJ 2818 Gcd

    https://vjudge.net/problem/HYSBZ-2818

     给定整数N，求1<=x,y<=N且Gcd(x,y)为素数的数对(x,y)有多少对。

    令f(i)为gcd(x,y)为i的数对的个数，则ans = ∑f(i), (i为素数)。因为本题只有一组样例，直接根据莫比乌斯反演对于每个素数i求出相应的f(i), 最终把答案累加即可。

#include<cstdio>
#include<cstring>
#include<string>
#include<cctype>
#include<iostream>
#include<set>
#include<map>
#include<cmath>
#include<sstream>
#include<vector>
#include<stack>
#include<queue>
#include<algorithm>
#define fin(a) freopen("a.txt","r",stdin)
#define fout(a) freopen("a.txt","w",stdout)
typedef long long LL;
using namespace std;
typedef pair<int, int> P;
const int dx[] = {0, 0, 0, -1, 1};
const int dy[] = {0, -1, 1, 0, 0};
const int INF = 1e8 + 10;
const int maxn = 1e7 + 10;
int prime[maxn], mu[maxn];
bool check[maxn];

void getMu() {
   mu[1] = 1;
   int tot = 0;
   for(int i = 2; i < maxn; i++) {
      if(!check[i]) { prime[tot++] = i; mu[i] = -1; }
      for(int j = 0; j < tot; j++) {
         if(i * prime[j] > maxn) break;
         check[i * prime[j]] = true;
         if(i % prime[j] == 0) {
            mu[i * prime[j]] = 0;
            break;
         }
         else {
            mu[i * prime[j]] = -mu[i];
         }
      }
   }
}

int main() {
   getMu();
   int n;
   scanf("%d", &n);
   LL ans = 0;
   for(int i = 0; prime[i] <= n; i++) {
      for(int j = prime[i]; j <= n; j += prime[i])
         ans += (LL)mu[j / prime[i]] * LL(n/j)*LL(n/j);
   }
   cout << ans << endl;
}

2. BZOJ 2301 Problem b

    https://vjudge.net/problem/HYSBZ-2301

    求有多少个数对 (x,y) ，满足 a ≤ x ≤ b ， c ≤ y ≤ d ，且 gcd(x,y) = k ， gcd(x,y) 函数为 x 和 y 的最大公约数。

     设solve(n, m)为满足1 ≤ x ≤ n , 1 ≤ y ≤ m , 且gcd(x,y) = k的数对的个数，而gcd(x,y) = k等价于gcd(x/k , y/k) = 1， 则易得最终答案为solve(b/k, d/k) - solve((a-1)/k, d/k) - solve((c-1)/k, b/k) + solve((a-1)/k, (c-1)/k)。

     因为题目有多组样例，因此若直接按照上面的方法不加优化的求解将会超时。在计算solve(n, m)时每次 ans = ∑mu[i] * (n/i) * (m/i) , (1 ≤ i ≤ min(n, m))，时间复杂度为O(n)。仔细观察发现，这部分计算中实际有一些重复的计算，在某一段i的连续区间内(n/i) * (m/i)的结果是相等的，因此我们可以给[1, min(n, m)]这段区间分块，即按照(n/i) * (m/i)的乘积进行分块，相同的分到一块中。假设某一块为[a, b]，则这一块贡献的答案为(sum[b] - sum[a-1]) * (n/i) * (m/i)，其中sum数组存储了mu数组的前缀和。

#include<cstdio>
#include<cstring>
#include<string>
#include<cctype>
#include<iostream>
#include<set>
#include<map>
#include<cmath>
#include<sstream>
#include<vector>
#include<stack>
#include<queue>
#include<algorithm>
#define fin(a) freopen("a.txt","r",stdin)
#define fout(a) freopen("a.txt","w",stdout)
typedef long long LL;
using namespace std;
typedef pair<int, int> P;
const int dx[] = {0, 0, 0, -1, 1};
const int dy[] = {0, -1, 1, 0, 0};
const int INF = 1e8 + 10;
const int maxn = 5e5 + 10;
int prime[maxn], mu[maxn], sum[maxn];
bool check[maxn];

void getMu() {
   mu[1] = 1;
   int tot = 0;
   for(int i = 2; i < maxn; i++) {
      if(!check[i]) { prime[tot++] = i; mu[i] = -1; }
      for(int j = 0; j < tot; j++) {
         if(i * prime[j] > maxn) break;
         check[i * prime[j]] = true;
         if(i % prime[j] == 0) {
            mu[i * prime[j]] = 0;
            break;
         }
         else {
            mu[i * prime[j]] = -mu[i];
         }
      }
   }
   for(int i = 1; i < maxn; i++)
      sum[i] = sum[i-1] + mu[i];
}

LL solve(int n, int m) {
   if(n > m) swap(n, m);
   LL ans = 0, pos;
   for(int i = 1; i <= n; i = pos + 1) {
      pos = min(n/(n/i), m/(m/i));
      ans += (sum[pos] - sum[i-1]) * LL(n/i) * LL(m/i);
   }
   return ans;
}

int main() {
   getMu();
   int T;
   scanf("%d", &T);
   while(T--) {
      int a, b, c, d, k;
      scanf("%d%d%d%d%d", &a, &b, &c, &d, &k);
      if(!k) { printf("0\n"); continue; }
      LL ans = solve(b/k, d/k) - solve((a-1)/k, d/k) - solve((c-1)/k, b/k) + solve((a-1)/k, (c-1)/k);
      printf("%lld\n", ans);
   }
}

3. HDU 4746 Mophues

    求满足1≤x≤n，1≤y≤m，且gcd(x，y)的质因子数不超过P的数对的个数。

     首先要预处理出每个数的质因子个数。

     如果令f(i)表示gcd(x，y)为i的数对的个数，用上一题mu数组前缀和的思想求出每个f(i)，再把所有满足i的质因子数不超过P的f(i)加起来的方法的话，很遗憾，会超时。

     怎么优化呢？再考虑分块的思想，在某个区间[a, b]内，当这个区间满足(n/i) * (m/i)相等时，我们还可以求出该区间内满足i的质因子数不超过P的mu[i]之和，这样之后，该区间只需计算一次即可。因此，我们令f(i, j)表示因子数为j时对mu[i]的贡献值，再通过求f(i, j)两维的前缀和即可得到f(i, j) ( 因子数不超过j时对mu[1]——mu[i]）的贡献值。

#include<cstdio>
#include<cstring>
#include<string>
#include<cctype>
#include<iostream>
#include<set>
#include<map>
#include<cmath>
#include<sstream>
#include<vector>
#include<stack>
#include<queue>
#include<algorithm>
#define fin(a) freopen("a.txt","r",stdin)
#define fout(a) freopen("a.txt","w",stdout)
typedef long long LL;
using namespace std;
typedef pair<int, int> P;
const int dx[] = {0, 0, 0, -1, 1};
const int dy[] = {0, -1, 1, 0, 0};
const int INF = 1e8 + 10;
const int maxn = 5e5 + 10;
int prime[maxn], mu[maxn], number[maxn];
int f[maxn][25];
bool check[maxn];

void getMu() {
   mu[1] = 1;
   int tot = 0;
   for(int i = 2; i < maxn; i++) {
      if(!check[i]) { prime[tot++] = i; mu[i] = -1; }
      for(int j = 0; j < tot; j++) {
         if(i * prime[j] > maxn) break;
         check[i * prime[j]] = true;
         if(i % prime[j] == 0) {
            mu[i * prime[j]] = 0;
            break;
         }
         else {
            mu[i * prime[j]] = -mu[i];
         }
      }
   }

}

int Divide(int x) {
   if(!check[x]) return 1;
   int len = sqrt(x+0.5);
   int cnt = 0;
   for(int i = 2; i <= len; i++) if(x % i == 0) {
       while(x && (x % i == 0)) {
          ++cnt;
          x /= i;
       }
       if(x == 1) break;
       if(!check[x]) { ++cnt; break; }
   }
   return cnt;
}

void getNumber() {
   for(int i = 2; i < maxn; i++)
      number[i] = Divide(i);
   for(int i = 1; i < maxn; i++)
      for(int j = i; j < maxn; j += i)
         f[j][number[i]] += mu[j / i];
   for(int i = 1; i < maxn; i++)
      for(int j = 1; j < 25; j++)
        f[i][j] += f[i][j-1];
   for(int i = 1; i < maxn; i++)
      for(int j = 0; j < 25; j++)
        f[i][j] += f[i-1][j];
}

int main() {
   getMu();
   getNumber();
   int T;
   scanf("%d", &T);
   while(T--) {
      int n, m, p;
      scanf("%d%d%d", &n, &m, &p);
      int mi = min(n, m);
      int pos;
      LL ans = 0;
      for(int i = 1; i <= mi; i = pos+1) {
         pos = min(n/(n/i), m/(m/i));
         ans += LL(n/i) * LL(m/i) * LL(f[pos][p] - f[i-1][p]);
      }
      printf("%lld\n", ans);
   }
}

   

Mophues