bzoj 2005 能量采集 - 容斥原理

栋栋有一块长方形的地，他在地上种了一种能量植物，这种植物可以采集太阳光的能量。在这些植物采集能量后，

栋栋再使用一个能量汇集机器把这些植物采集到的能量汇集到一起。 栋栋的植物种得非常整齐，一共有n列，每列

有m棵，植物的横竖间距都一样，因此对于每一棵植物，栋栋可以用一个坐标(x, y)来表示，其中x的范围是1至n，

表示是在第x列，y的范围是1至m，表示是在第x列的第y棵。 由于能量汇集机器较大，不便移动，栋栋将它放在了

一个角上，坐标正好是(0, 0)。 能量汇集机器在汇集的过程中有一定的能量损失。如果一棵植物与能量汇集机器

连接而成的线段上有k棵植物，则能量的损失为2k + 1。例如，当能量汇集机器收集坐标为(2, 4)的植物时，由于

连接线段上存在一棵植物(1, 2)，会产生3的能量损失。注意，如果一棵植物与能量汇集机器连接的线段上没有植

物，则能量损失为1。现在要计算总的能量损失。 下面给出了一个能量采集的例子，其中n = 5，m = 4，一共有20

棵植物，在每棵植物上标明了能量汇集机器收集它的能量时产生的能量损失。 在这个例子中，总共产生了36的能

量损失。

Input

仅包含一行，为两个整数n和m。

Output

仅包含一个整数，表示总共产生的能量损失。

Sample Input

【样例输入1】
5 4
【样例输入2】
3 4

Sample Output

【样例输出1】
36
【样例输出2】
20
对于100%的数据：1 ≤ n, m ≤ 100,000。

---

　　题目大意 求

　　因为n和m单个比较小，而且gcd(i, j)不超过min(n, m)，因此可以想到去枚举gcd(i, j)的取值，然后统计个数，再乘一乘即可。

　　至于这个统计个数有一个很好的方法就是容斥原理，计算i和j都是d的倍数时的对数，这个很好算，直接搬结论吧，是，然后再减去gcd(i, j)的值为2d, 3d, ...的时候。

　　这样显然倒着算能够更简单，这样的时间复杂度是O(nlog₂n)。

Code

 1 /**
 2  * bzoj
 3  * Problem#2005
 4  * Accepted
 5  * Time:24ms
 6  * Memory:2068k
 7  */
 8 #include <iostream>
 9 #include <cstdio>
10 #include <ctime>
11 #include <cmath>
12 #include <cctype>
13 #include <cstring>
14 #include <cstdlib>
15 #include <fstream>
16 #include <sstream>
17 #include <algorithm>
18 #include <map>
19 #include <set>
20 #include <stack>
21 #include <queue>
22 #include <vector>
23 #include <stack>
24 #ifndef WIN32
25 #define Auto "%lld"
26 #else
27 #define Auto "%I64d"
28 #endif
29 using namespace std;
30 typedef bool boolean;
31 const signed int inf = (signed)((1u << 31) - 1);
32 const signed long long llf = (signed long long)((1ull << 61) - 1);
33 const double eps = 1e-6;
34 const int binary_limit = 128;
35 #define smin(a, b) a = min(a, b)
36 #define smax(a, b) a = max(a, b)
37 #define max3(a, b, c) max(a, max(b, c))
38 #define min3(a, b, c) min(a, min(b, c))
39 template<typename T>
40 inline boolean readInteger(T& u){
41     char x;
42     int aFlag = 1;
43     while(!isdigit((x = getchar())) && x != '-' && x != -1);
44     if(x == -1) {
45         ungetc(x, stdin);    
46         return false;
47     }
48     if(x == '-'){
49         x = getchar();
50         aFlag = -1;
51     }
52     for(u = x - '0'; isdigit((x = getchar())); u = (u << 1) + (u << 3) + x - '0');
53     ungetc(x, stdin);
54     u *= aFlag;
55     return true;
56 }
57 
58 int n, m;
59 
60 inline void init() {
61     readInteger(n);
62     readInteger(m);
63     if(n > m)    swap(n, m);
64 }
65 
66 long long res = 0;
67 long long f[100005];
68 inline void solve() {
69     for(int d = n; d; d--) {
70         f[d] = (n / d) * 1LL * (m / d);
71         for(int i = (d << 1); i <= n; i += d)
72             f[d] -= f[i];
73         res += d * 1LL * f[d];
74     }
75     res = (res << 1) - m * 1LL * n;
76     printf(Auto, res);
77 }
78 
79 int main() {
80     init();
81     solve();
82     return 0;
83 }

 　　但是另外用点黑科技，可以让它变成O(n)。

　　记得有一个关于欧拉函数的结论。于是有:

　　

　　现在等同于求有多少个数对(i, j)使得i能够整除d，j能够整除d，这是一个刚刚就解决了的问题，于是我们得到了下面这个优美的式子:

　　对于求phi函数的值这个事情就交给已经闲置了一会儿的线性筛去做吧。

Code

/**
 * bzoj
 * Problem#2005
 * Accepted
 * Time:4ms
 * Memory:2172k
 */
#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>
#include <cctype>
#include <cstring>
#include <cstdlib>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <map>
#include <set>
#include <stack>
#include <queue>
#include <vector>
#include <list>
#ifndef WIN32
#define Auto "%lld"
#else
#define Auto "%I64d"
#endif
using namespace std;
typedef bool boolean;
const signed int inf = (signed)((1u << 31) - 1);
const signed long long llf = (signed long long)((1ull << 61) - 1);
const double eps = 1e-6;
const int binary_limit = 128;
#define smin(a, b) a = min(a, b)
#define smax(a, b) a = max(a, b)
#define max3(a, b, c) max(a, max(b, c))
#define min3(a, b, c) min(a, min(b, c))
template<typename T>
inline boolean readInteger(T& u){
    char x;
    int aFlag = 1;
    while(!isdigit((x = getchar())) && x != '-' && x != -1);
    if(x == -1) {
        ungetc(x, stdin);    
        return false;
    }
    if(x == '-'){
        x = getchar();
        aFlag = -1;
    }
    for(u = x - '0'; isdigit((x = getchar())); u = (u << 1) + (u << 3) + x - '0');
    ungetc(x, stdin);
    u *= aFlag;
    return true;
}

const int limit = 1e5;

int n, m;
int num = 0;
int prime[limit + 1];
int phi[limit + 1];
boolean vis[limit + 1];
long long res = 0;

inline void init() {
    readInteger(n);
    readInteger(m);
    if(n > m)    swap(n, m);
}

inline void Euler() {
    phi[1] = 1;
    memset(vis, false, sizeof(boolean) * (n + 1));
    for(int i = 2; i <= n; i++) {
        if(!vis[i])    prime[num++] = i, phi[i] = i - 1;
        for(int j = 0; j < num && i * 1LL * prime[j] <= n; j++) {
            int c = i * prime[j];
            vis[c] = true;
            if((i % prime[j]) == 0) {
                phi[c] = prime[j] * phi[i];
                break;
            } else {
                phi[c] = phi[prime[j]] * phi[i];
            }
        }
    }
}

inline void solve() {
    for(int i = 1; i <= n; i++)
        res += (n / i) * 1LL * (m / i) * phi[i];
    printf(Auto"\n", (res * 2) - (long long)n * m);
}

int main() {
    init();
    Euler();
    solve();
    return 0;
}