[BZOJ2005][NOI2010]能量采集（莫比乌斯反演）

最新推荐文章于 2019-08-26 22:56:57 发布

原创最新推荐文章于 2019-08-26 22:56:57 发布 · 271 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#2010 #gcd #NOI

BZOJ、UOJ、LOJ 专栏收录该内容

341 篇文章

订阅专栏

本文探讨了如何求解特定形式的数论求和问题，给出了一种有效的算法实现，并通过代码示例展示了求解过程。针对求和∑ni=1∑mi=1(2gcd(i,j)−1)，文章详细解释了如何将其转换为更易处理的形式，并利用莫比乌斯反演等数学工具进行高效计算。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题意：给定 $n,m$ ，求 $\sum_{i=1}^n\sum_{i=1}^m(2\gcd(i,j)-1)$ 。
容易推出，原式可化为 $2\sum_{i=1}^n\sum_{i=1}^m\gcd(i,j)-nm$ 。
而关键就是求 $\sum_{i=1}^n\sum_{i=1}^m\gcd(i,j)$ 的值。考虑先枚举 $\gcd(i,j)$ 的值，就可以得到：
原式 $=\sum_{d=1}^{\min(n,m)}\sum_{i=1}^n\sum_{i=1}^m[\gcd(i,j)=d]$
$=\sum_{d=1}^{\min(n,m)}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]$ 。
而求 $\sum_{i=1}^x\sum_{i=1}^y[\gcd(i,j)=1]$ 是莫比乌斯反演中最经典的模型。这个式子等于 $\sum_{d=1}^{\min(x,y)}\mu(d)\lfloor\frac{x}{d}\rfloor\lfloor\frac{y}{d}\rfloor$ 。可以在根号的时间内求出，具体做法可查阅其他资料。
补充一点：如果有多组询问，那么：
$\sum_{i=1}^n\sum_{i=1}^m\gcd(i,j)=\sum_{u=1}^{\min(n,m)}\lfloor\frac{n}{u}\rfloor\lfloor\frac{m}{u}\rfloor\sum_{d|u}(\mu(\frac{u}{d})*d)$ 。
$\sum_{d|u}(\mu(\frac{u}{d})*d)$ 的值和前缀和可以预处理。证明由读者去思考。
代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
typedef long long ll;
const int MaxN = 1e5, N = MaxN + 5;
int n, m, tot, pri[N], miu[N], sum[N];
bool mark[N];
void sieve() {
    int i, j; mark[0] = mark[1] = 1; miu[1] = 1;
    for (i = 2; i <= MaxN; i++) {
        if (!mark[i]) pri[++tot] = i, miu[i] = -1;
        for (j = 1; j <= tot; j++) {
            if (1ll * i * pri[j] > MaxN) break;
            mark[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
            else miu[i * pri[j]] = -miu[i];
        }
    }
    for (i = 1; i <= MaxN; i++) sum[i] = sum[i - 1] + miu[i];
}
ll solve(int n, int m) {
    int i; ll ans = 0;
    for (i = 1; i <= min(n, m);) {
        int nxt = min(n / (n / i), m / (m / i));
        ans += 1ll * (sum[nxt] - sum[i - 1]) * (n / i) * (m / i);
        i = nxt + 1;
    }
    return ans;
}
int main() {
    int i; n = read(); m = read(); ll ans = 0; sieve();
    for (i = 1; i <= min(n, m); i++)
        ans += solve(n / i, m / i) * i;
    cout << (ans - 1ll * n * m) * 2 + 1ll * n * m << endl;
    return 0;
}