☆【线性规划】【容斥原理】【NOI2010】能量采集

Description
栋栋有一块长方形的地，他在地上种了一种能量植物，这种植物可以采集太阳光的能量。在这些植物采集能量后，栋栋再使用一个能量汇集机器把这些植物采集到的能量汇集到一起。 栋栋的植物种得非常整齐，一共有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

【数据规模和约定】
对于10%的数据：1 ≤ n, m ≤ 10；
对于50%的数据：1 ≤ n, m ≤ 100；
对于80%的数据：1 ≤ n, m ≤ 1000；
对于90%的数据：1 ≤ n, m ≤ 10,000；
对于100%的数据：1 ≤ n, m ≤ 100,000。

首先有一个公式：(0, 0)到(x, y)的线段（不包括(0, 0)）上整点的个数为gcd(x, y)。
那么对于这道题，朴素算法自然是暴力累加，把所有2gcd(i, j) - 1(1<=i<=n, 1<=j<=m)累加起来即为答案。

显然这个算法只能过8组。

于是寻找快速算出所有gcd(i, j)的方法。

可以用容斥原理进行优化。
                                            │ n │  │ m │
令f[i]为所有质因数为i的数的组数，于是f[i] = │---│·│---│ - sigma(j·i <= max(n, m)) f[i·j]。
                                            └ i ┘  └ i ┘
在此基础上还可以优化。
                            │ n │  │ m │
可以发现，对于一段连续的i， │---│和│---│的取值都是一样的,所以可以放在一起考虑（这点还没想清楚。）
                            └ i ┘  └ i ┘
最后的结果，直接把所有的i·f[i]累加再乘以2再减去n·m即可。
Accode：

#include <cstdio>
typedef long long int64;
int64 f[100010], n, m, ans, tmp;
int main()
{
    freopen("energy.in", "r", stdin);
    freopen("energy.out", "w", stdout);
    scanf("%I64d%I64d", &n, &m);
    int64 k = n > m ? n : m;
    for (int i = k; i; --i)
    {
        if ((n / i) * (m / i) == tmp)
            f[i] = f[i + 1];
        else
        {
            tmp = f[i] = (n / i) * (m / i);
            for (int j = 2; j <= k / i; ++j)
                f[i] -= f[i * j];
        }
        ans += f[i] * i;
    }
    printf("%I64d\n", (ans <<= 1) -= n * m);
    return 0;
}

再贴个朴素的代码：

#include <cstdio>
inline int gcd(int a, int b)
{
    while (b)
    {
        int tmp = a;
        a = b;
        b = tmp % b;
    }
    return a;
}

int main()
{
    freopen("energy.in", "r", stdin);
    freopen("energy.out", "w", stdout);
    int n = 0, m = 0, ans = 0;
    scanf("%d%d", &n, &m);
    for (int i = 1; i < n + 1; ++i)
    for (int j = 1; j < m + 1; ++j)
        --ans += gcd(i, j) << 1;
    printf("%d\n", ans);
    return 0;
}

第二次做：

#include <cstdio>
long long f[100010], ans, tmp, n, m, K;
int main()
{
    freopen("energy.in", "r", stdin);
    freopen("energy.out", "w", stdout);
    scanf("%I64d%I64d", &n, &m); K = n > m ? n : m;
    for (long long i = K; i; --i)
    {
        f[i] = (n / i) * (m / i); //
        if (f[i] == tmp) f[i] = f[i + 1];
        else
        {
            tmp = f[i];
            for (long long j = i << 1; j <= K; j += i)
                f[i] -= f[j];
        }
        ans += f[i] * i; //
    }
    printf("%I64d\n", (ans << 1) - n * m);
        //这里若n, m为int类型的话会超范围。 
    return 0;
}