莫比乌斯反演题目结(上)

首先是喜闻乐见的几个入门基础题,连题面基本都是一样的,流程是预处理mu函数,得到输入数据后整除分块(如果时间复杂度需要的话),思路主要是套用下图第二个公式:

(截图来自电科bilibili上的算法讲堂大家可以去看哦)

洛谷2257

题意:全部gcd(1~n,1~m)为质数的个数。

刚开始不用管他要求的数有什么特性,就列一般式子:

(截图来自pengym的博客题解)

但这个题用第二条推出的式子还得变一变,就直观感受去变!d的范围?n是质数?把F(d)提出去?

下图中T即为上图中的d,t即为n。

式子写到这里就可以写代码了,先把mu全加上,这样写:

1   rep(j, 1, cnt) {
2         for (int i = 1; i * primes[j] <= n; i++) {
3             g[primes[j] * i] += mu[i];
4         }
5     }

然后整除分块是某一块乘上的(n/T)*(m/T)是不变的,所以这一块一起算,代码是这样的:

1   rep(i, 1, n)    sum[i] = sum[i-1] + g[I];//这个是预处理
1         ll ans = 0ll;
2         for (int l = 1, r; l <= n; l = r + 1) {
3             r = min(n / (n / l), m / (m / l));
4             ans += (sum[r] - sum[l - 1]) * (n / l) * (m / l);
5         }

总代码(部分):

 1 void pre(int n) {
 2     mu[1] = 1;
 3     rep(i, 2, n) {
 4         if (!vis[i]) {
 5             mu[i] = -1;
 6             primes[++cnt] = i;
 7         }
 8         for (int j = 1; j <= cnt && primes[j] * i <= n; j++) {
 9             vis[primes[j] * i] = true;
10             if (i % primes[j] == 0)    break;
11             mu[primes[j] * i] = -mu[i];    
12         }
13     }
14 
15     rep(j, 1, cnt) {
16         for (int i = 1; i * primes[j] <= n; i++) {
17             g[primes[j] * i] += mu[i];
18         }
19     }
20     rep(i, 1, n)    sum[i] = sum[i-1] + g[i];
21 }
22 
23 int main() {
24     pre(maxn);
25     for (T = ri; T; T--) {
26         n = ri, m = ri;
27         if (n > m)    swap(n, m);
28 
29         ll ans = 0ll;
30         for (int l = 1, r; l <= n; l = r + 1) {
31             r = min(n / (n / l), m / (m / l));
32             ans += (sum[r] - sum[l - 1]) * (n / l) * (m / l);
33         }
34         printf("%lld
", ans);
35     }
36     return 0;
37 }

HDU1695

(复习一下式子)

题意:gcd(1~b,1~d)== k的个数。

跟上一题很像,甚至这道题根本不用整除分块,也不用进一步推式子……就是求f(k),直接拿第二个式子扫一遍加和就行了。另外去个重,本题(x,y)和(y,x)算一种。

 1 void pre(int n) {
 2     mu[1] = 1;
 3     rep(i, 2, n) {
 4         if (!vis[i]) {
 5             primes[++cnt] = i;
 6             mu[i] = -1;
 7         }
 8         for (int j = 1; j <= cnt && primes[j] * i <= n; j++) {
 9             vis[primes[j] * i] = true;
10             if (i % primes[j] == 0)    break;
11             mu[primes[j] * i] = -mu[i];
12         }
13     }
14 }
15 
16 ll cal(int b, int d) {
17     ll ret = 0ll;
18     for (int i = 1; i * k <= b; i++) {
19         int t = i * k;
20         ret += (ll)mu[i] * (b / t) * (d / t);
21     }
22     return ret;
23 }
24 
25 int main() {
26     pre(maxn - 5);
27     for (int T = ri, i = 1; i <= T; i++) {
28         a = ri, b = ri, c = ri, d = ri, k = ri;
29         if (b > d)    swap(b, d);
30         printf("Case %d: %lld
", i, k ? cal(b, d) - cal(b, b) / 2 : 0);
31     }
32     return 0;
33 }

BZOJ1101

n次询问gcd(1~a, 1~b) == d。同一道题。

 1 #pragma comment(linker, "/STACK:1024000000,1024000000")
 2 #include <cstdio>
 3 #include <cstring>
 4 #include <cstdlib>
 5 #include <cmath>
 6 #include <ctime>
 7 #include <cctype>
 8 #include <climits>
 9 #include <iostream>
10 #include <iomanip>
11 #include <algorithm>
12 #include <string>
13 #include <sstream>
14 #include <stack>
15 #include <queue>
16 #include <set>
17 #include <map>
18 #include <vector>
19 #include <list>
20 #include <fstream>
21 #include <bitset>
22 #define init(a, b) memset(a, b, sizeof(a))
23 #define rep(i, a, b) for (int i = a; i <= b; i++)
24 #define irep(i, a, b) for (int i = a; i >= b; i--)
25 using namespace std;
26 
27 typedef double db;
28 typedef long long ll;
29 typedef unsigned long long ull;
30 typedef pair<int, int> P;
31 const int inf = 0x3f3f3f3f;
32 const ll INF = 1e18;
33 
34 template <typename T> void read(T &x) {
35     x = 0;
36     int s = 1, c = getchar();
37     for (; !isdigit(c); c = getchar())
38         if (c == '-')    s = -1;
39     for (; isdigit(c); c = getchar())
40         x = x * 10 + c - 48;
41     x *= s;
42 }
43 
44 template <typename T> void write(T x) {
45     if (x < 0)    x = -x, putchar('-');
46     if (x > 9)    write(x / 10);
47     putchar(x % 10 + '0');
48 }
49 
50 template <typename T> void writeln(T x) {
51     write(x);
52     puts("");
53 }
54 
55 const int maxn = 5e4 + 5;
56 int n, a, b, d;
57 int mu[maxn], primes[maxn], tot, sum[maxn];
58 bool vis[maxn];
59 
60 void pre(int n) {
61     mu[1] = 1;
62     rep(i, 2, n) {
63         if (!vis[i]) {
64             primes[++tot] = i;
65             mu[i] = -1;
66         }
67         for (int j = 1; j <= tot && primes[j] * i <= n; j++) {
68             vis[primes[j] * i] = true;
69             if (i % primes[j] == 0)    break;
70             mu[primes[j] * i] = - mu[i];
71         }
72     }
73     rep(i, 1, n)    sum[i] = sum[i - 1] + mu[i];
74 }
75 
76 ll solve(int n, int m) {
77     if (n > m)    swap(n, m);
78     ll ret = 0ll;
79     for (int l = 1, r; l <= n; l = r + 1) {
80         r = min(n / (n / l), m / (m / l));
81         ret += (ll)(sum[r] - sum[l - 1]) * (n / l) * (m / l);
82     }
83     return ret;
84 }
85 
86 int main() {
87     pre(maxn - 5);
88     for (read(n); n; n--) {
89         read(a), read(b), read(d);
90         writeln(solve(a / d, b / d));
91     }
92     return 0;
93 }
bzoj1101

BZOJ2301

跟HDU1695的区别是变成了gcd(a~b,c~d)== k的个数。

还是那么做,但是ans容斥一下:

1 printf("%lld
", cal(b, d) - cal(b, c) - cal(a, d) + cal(a, c));

唔,实际上应该是cal(b,c-1)、cal(a-1,d)……主要是上一行代码直接操作掉了:

1 int main() {
2     pre(maxn - 5);
3     for (T = ri; T; T--) {
4         a = ri, b = ri, c = ri, d = ri, k = ri;
5         a = (a - 1) / k, b = b / k, c = (c - 1) / k, d = d / k;
6         printf("%lld
", cal(b, d) - cal(b, c) - cal(a, d) + cal(a, c));
7     }
8     return 0;
9 }

除以k有降低复杂度的功效,gcd和范围同时除以k使得问题变成求f(1)!也就是所有的mu[d]*F(d)加和,不用预处理k的倍数了。上一题也可以这么写的。

先对照一下公式,旧书不厌百回读:

 1 void pre(int n) {
 2     mu[1] = 1;
 3     rep(i, 2, n) {
 4         if (!vis[i]) {
 5             primes[++cnt] = i;
 6             mu[i] = -1;
 7         }
 8         for (int j = 1; j <= cnt && primes[j] * i <= n; j++) {
 9             vis[primes[j] * i] = true;
10             if (i % primes[j] == 0)    break;
11             mu[primes[j] * i] = -mu[i];
12         }
13     }
14     rep(i, 1, n)    sum[i] = sum[i - 1] + mu[i];
15 }
16 
17 ll cal(int n, int m) {
18     if (n > m)    swap(n, m);
19     ll ans = 0;
20 
21     for (int l = 1, r; l <= n; l = r + 1) {
22         r = min(n / (n / l), m / (m / l));
23         ans += (sum[r] - sum[l - 1]) * (n / l) * (m / l);
24     }
25     return ans;
26 }

POJ3904

题意:首先,这次范围为给定数列,不是从1~n、从1~m所有的数这种大范围了;然后抽取四个数使他们gcd为1。也就是共有多少四元组的gcd是1,还是计数~

通过这个题可以得知套路公式中的F(n)一般是根据莫比乌斯而跟f(d)有关的一个式子,另外又会等于一个与题目相关的、直观感受的、可以直接求出的式子。

比如本题还是使用套路公式,但是F(d)不再是(n/d)*(m/d),而是数列中d的倍数的这些数中,随便抽4个,有多少种取法:

1 ll ans = 0ll;
2 rep(i, 1, maxn - 5) {
3     ans += mu[i] * C(c[i], 4);//gcd是1的莫比乌斯公式
4 }  

总代码:

 1 void mobi(int n) {
 2     mu[1] = 1;
 3     rep(i, 2, n) {
 4         if (!vis[i]) {
 5             primes[++cnt] = i;
 6             mu[i] = -1;
 7         }
 8         for (int j = 1; j <= cnt && primes[j] * i <= n; j++) {
 9             vis[primes[j] * i] = true;
10             if (i % primes[j] == 0)    break;
11             mu[primes[j] * i] = -mu[i];
12         } 
13     }
14 }
15 
16 void get_cnt() {
17     init(c, 0);
18     rep(i, 1, n) {
19         int x = a[i];
20         for (int j = 1; j * j <= x; j++) {
21             if (x % j == 0) {
22                 c[j]++;
23                 if (j * j != x)    c[x / j]++;
24             }
25         }
26     }
27 }
28 
29 ll C(int n, int k) {
30     if (n < k)    return 0;
31 
32     ll ret = 1ll;
33     rep(i, 1, k) {
34         ret = ret * (n - i + 1) / i;
35     }
36     return ret;
37 }
38 
39 int main() {
40     ios_base::sync_with_stdio(false);
41     cin.tie(0);
42 
43     mobi(maxn - 5);
44     while (cin >> n) {
45         rep(i, 1, n)    cin >> a[i];
46         get_cnt();
47 
48         ll ans = 0ll;
49         rep(i, 1, maxn - 5) {
50             ans += mu[i] * C(c[i], 4);
51         }
52         cout << ans << endl;
53     }
54     return 0;
55 }

BZOJ2005

做完以上几题以后此题可独立做出。还是个gcd(1~n,1~m),但这次不是个数加和而是值加和。参考代码:

 1 #pragma comment(linker, "/STACK:1024000000,1024000000")
 2 #include <cstdio>
 3 #include <cstring>
 4 #include <cstdlib>
 5 #include <cmath>
 6 #include <ctime>
 7 #include <cctype>
 8 #include <climits>
 9 #include <iostream>
10 #include <iomanip>
11 #include <algorithm>
12 #include <string>
13 #include <sstream>
14 #include <stack>
15 #include <queue>
16 #include <set>
17 #include <map>
18 #include <vector>
19 #include <list>
20 #include <fstream>
21 #include <bitset>
22 #define init(a, b) memset(a, b, sizeof(a))
23 #define rep(i, a, b) for (int i = a; i <= b; i++)
24 #define irep(i, a, b) for (int i = a; i >= b; i--)
25 using namespace std;
26 
27 typedef double db;
28 typedef long long ll;
29 typedef unsigned long long ull;
30 typedef pair<int, int> P;
31 const int inf = 0x3f3f3f3f;
32 const ll INF = 1e18;
33 
34 template <typename T> void read(T &x) {
35     x = 0;
36     int s = 1, c = getchar();
37     for (; !isdigit(c); c = getchar())
38         if (c == '-')    s = -1;
39     for (; isdigit(c); c = getchar())
40         x = x * 10 + c - 48;
41     x *= s;
42 }
43 
44 template <typename T> void write(T x) {
45     if (x < 0)    x = -x, putchar('-');
46     if (x > 9)    write(x / 10);
47     putchar(x % 10 + '0');
48 }
49 
50 template <typename T> void writeln(T x) {
51     write(x);
52     puts("");
53 }
54 
55 const int maxn = 1e5 + 5;
56 int n, m;
57 int mu[maxn], primes[maxn], tot;
58 bool vis[maxn];
59 ll ans;
60 
61 void pre(int n) {
62     mu[1] = 1;
63     rep(i, 2, n) {
64         if (!vis[i]) {
65             primes[++tot] = i;
66             mu[i] = -1;
67         }
68         for (int j = 1; j <= tot && i * primes[j] <= n; j++) {
69             vis[i * primes[j]] = true;
70             if (i % primes[j] == 0)    break;
71             mu[i * primes[j]] = -mu[i];
72         }
73     }
74 }
75 
76 ll cal(int n, int m) {
77     ll ret = 0ll;
78     rep(i, 1, n) {
79         ret += (ll)mu[i] * (n / i) * (m / i);
80     }
81     return ret;
82 }
83 
84 int main() {
85     read(n), read(m);
86     if (n > m)    swap(n, m);
87     pre(maxn - 5);
88     for (int l = 1, r; l <= n; l = r + 1) {
89         r = min(n / (n / l), m / (m / l));
90         ans += (ll)(l + r) * (r - l + 1) / 2 * cal(n / l, m / l);
91     }
92     writeln(2 * ans - (ll)n * m);
93     return 0;
94 }

BZOJ2154

这个我就推不出来了Orz。这式子还蛮复杂的,推荐这个博客,证明得蛮清楚的。

核心:

#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cctype>
#include <climits>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <string>
#include <sstream>
#include <stack>
#include <queue>
#include <set>
#include <map>
#include <vector>
#include <list>
#include <fstream>
#include <bitset>
#define init(a, b) memset(a, b, sizeof(a))
#define rep(i, a, b) for (int i = a; i <= b; i++)
#define irep(i, a, b) for (int i = a; i >= b; i--)
using namespace std;

typedef double db;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> P;
const int inf = 0x3f3f3f3f;
const ll INF = 1e18;

template <typename T> void read(T &x) {
    x = 0;
    int s = 1, c = getchar();
    for (; !isdigit(c); c = getchar())
        if (c == '-')    s = -1;
    for (; isdigit(c); c = getchar())
        x = x * 10 + c - 48;
    x *= s;
}

template <typename T> void write(T x) {
    if (x < 0)    x = -x, putchar('-');
    if (x > 9)    write(x / 10);
    putchar(x % 10 + '0');
}

template <typename T> void writeln(T x) {
    write(x);
    puts("");
}

const int maxn = 1e7 + 5;
const int mod = 20101009;
int n, m;
ll ans;
int mu[maxn], primes[maxn], tot;
ll sum[maxn];
bool vis[maxn];

void pre(int n) {
    mu[1] = 1;
    rep(i, 2, n) {
        if (!vis[i]) {
            primes[++tot] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= tot && primes[j] <= n / i; j++) {
            vis[primes[j] * i] = true;
            if (i % primes[j] == 0)    break;
            mu[primes[j] * i] = -mu[i];
        }
    }
    rep(i, 1, n) {
        if (mu[i] < 0)    mu[i] = mod - 1;
        sum[i] = sum[i - 1] + (ll)i * i % mod * mu[i] % mod;
        sum[i] %= mod;
    }
}

ll S(ll x, ll y) {
    return (x * (x + 1) / 2 % mod) * (y * (y + 1) / 2 % mod) % mod;
}

ll F(int x, int y) {
    ll ret = 0ll;
    for (int l = 1, r; l <= x; l = r + 1) {
        r = min(x / (x / l), y / (y / l));
        ret += (sum[r] - sum[l - 1] + mod) % mod * S(x / l, y / l);
        ret %= mod;
    }
    return ret;
}

int main() {
    read(n), read(m);
    if (n > m)    swap(n, m);
    pre(m);
    for (int l = 1, r; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans += (ll)(r - l + 1) * (l + r) / 2 % mod * F(n / l, m / l);
        ans %= mod;
    }
    writeln(ans);
    return 0;
}

莫比乌斯反演题目结(下)

原文地址:https://www.cnblogs.com/AlphaWA/p/10478075.html