BZOJ 3309 DZY Loves Math

题目链接

题意

求

$$\sum_{i=1}^{a}\sum_{j=1}^{b}f(gcd(i, j))$$
其中 $$\begin{eqnarray} f(x) = \begin{cases} \alpha_1, x = p1^\alpha_1*p2^\alpha_2+...+pn^\alpha_n, \alpha_1 > \alpha_2,...,\alpha_n\cr 0, x = 1 \end{cases} \end{eqnarray}$$

推导

法一：艾弗森约定

记 $gcd(i, j) = k$, 则有 $gcd(\frac{i}{k}, \frac{j}{k}) = 1$,
故原式可化为

$$\sum_{k=1}^{min(a,b)}\sum_{i=1}^{\lfloor \frac{a}{k}\rfloor}\sum_{j=1}^{\lfloor \frac{b}{k}\rfloor}f(k)[gcd(i,j)=1]$$
用莫比乌斯函数将艾弗森约定展开得 $$\sum_{k=1}^{min(a,b)}f(k)\sum_{i=1}^{\lfloor \frac{a}{k}\rfloor}\sum_{j=1}^{\lfloor \frac{b}{k}\rfloor}\sum_{d|i,d|j}\mu(d)$$
即 $$\sum_{k=1}^{min(a,b)}f(k)\sum_{d=1}^{min(\lfloor\frac{a}{k}\rfloor,\lfloor\frac{b}{k}\rfloor)}\mu(d)\sum_{i=1}^{\lfloor\frac{a}{kd}\rfloor}\sum_{j=1}^{\lfloor\frac{b}{kd}\rfloor}$$
即 $$\sum_{k=1}^{min(a,b)}f(k)\sum_{d=1}^{min(\lfloor\frac{a}{k}\rfloor,\lfloor\frac{b}{k}\rfloor)}\mu(d)\lfloor\frac{a}{kd}\rfloor\lfloor\frac{b}{kd}\rfloor$$
令 $T = kd$, 得 $$\sum_{T=1}^{min(a,b)}\lfloor\frac{a}{T}\rfloor\lfloor\frac{b}{T}\rfloor\sum_{d|T}f(\frac{T}{d})\mu(d)$$
显然，接下来的路子就是求前缀和，然后分块计算，然而 关键问题 就是怎么求这个 $$\sum_{d|T}f(\frac{T}{d})\mu(d)$$ 呢?

参考PoPoQQQ可知，对于

$$g(T) = \sum_{d|T}f(\frac{T}{d})\mu(d), T=p1^\alpha_1*p2^\alpha_2*...*pn^\alpha_n$$
有 $$\begin{eqnarray} g(T)= \begin{cases} 0, \exists \alpha_i\neq\alpha_j,\cr (-1)^(n+1), otherwise \end{cases} \end{eqnarray}$$
具体证明见上附链接，要点为从 $k$ 个元素中选取奇数个元素的种数 $=$ 选取偶数个元素的种数，因为 $(1-1)^k = 0 \rightarrow C_{k}^{0}+C_{k}^{2}+... = C_{k}^{1}+C_{k}^{3}+...$

法二：莫比乌斯反演（待补）

Code

#include <bits/stdc++.h>
#define maxn 10000000
#define maxm maxn + 10
using namespace std;
typedef long long LL;
bool check[maxm];
int prime[maxm], cnt[maxm], val[maxm], g[maxm], sum[maxm];
void init() {
    int tot = 0; g[1] = 0;
    for (int i = 2; i <= maxn; ++i) {
        if (!check[i]) {
            prime[tot++] = i;
            cnt[i] = 1, val[i] = i, g[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) {
                cnt[i * prime[j]] = cnt[i] + 1, val[i * prime[j]] = val[i] * prime[j];
                int temp = (i * prime[j]) / val[i * prime[j]];
                if (temp == 1) g[i * prime[j]] = 1;
                else g[i * prime[j]] = (cnt[temp] == cnt[i * prime[j]] ? -g[temp] : 0);
                break;
            }
            cnt[i * prime[j]] = 1, val[i * prime[j]] = prime[j];
            int temp = (i * prime[j]) / val[i * prime[j]];
            if (temp == 1) g[i * prime[j]] = 1;
            else g[i * prime[j]] = (cnt[temp] == cnt[i * prime[j]] ? -g[temp] : 0);
        }
    }
    for (int i = 1; i <= maxn; ++i) sum[i] = sum[i - 1] + g[i];
}
void work() {
    LL a, b;
    scanf("%lld%lld", &a, &b);
    int le, ri, lim = min(a, b);
    LL ans = 0;
    for (int i = 1; i <= lim; i = ri + 1) {
        le = i, ri = min(a / (a / i), b / (b / i));
        ans += (a / i) * (b / i) * (sum[ri] - sum[le - 1]);
    }
    printf("%lld\n", ans);
}
int main() {
    init();
    int T;
    scanf("%d", &T);
    while (T--) work();
    return 0;
}