-数论-杜教筛-莫比乌斯反演- [LOJ 6491]「XXOI 2018」简单的最大公约数

题目描述

Link

Solution

先推一波柿子：

$$\begin{align}\sum_{i_1=1}^m \sum_{i_2=1}^m \sum_{i_3=1}^m ...\sum_{i_n=1}^m \gcd(i_1,i_2,i_3,...,i_n) &= \sum_{d=1}^m \sum_{i_1=1}^m \sum_{i_2=1}^m \sum_{i_3=1}^m ...\sum_{i_n=1}^m [\gcd(i_1,i_2,i_3,...,i_n) = d]\\\\ &= \sum_{d=1}^m \sum_{i_1=1}^{\lfloor \frac{m}{d}\rfloor} \sum_{i_2=1}^{\lfloor \frac{m}{d}\rfloor} \sum_{i_3=1}^{\lfloor \frac{m}{d}\rfloor} ...\sum_{i_n=1}^{\lfloor \frac{m}{d}\rfloor} [\gcd(i_1,i_2,i_3,...,i_n) = 1] \\\\&= \sum_{d=1}^m d\sum_{i_1=1}^{\lfloor \frac{m}{d}\rfloor} \sum_{i_2=1}^{\lfloor \frac{m}{d}\rfloor} \sum_{i_3=1}^{\lfloor \frac{m}{d}\rfloor} ...\sum_{i_n=1}^{\lfloor \frac{m}{d}\rfloor} \sum_{t|i_1,t|i_2,t|i_3,...,t|i_n} \mu(t) \\\\ &= \sum_{d=1}^m \lfloor \frac{m}{T}\rfloor^n \sum_{d|T} d \mu(\frac{T}{d}) \\\\ &= \sum_{d=1}^m \lfloor \frac{m}{T}\rfloor^n \varphi(T) \end{align}​$$

然而这道题的数据范围是 $10 \le n,m \le 10^{11}​$，所以单纯用整除分块是不够的，所以还需要用杜教筛。

Code

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <tr1/unordered_map>
#define MAXN 10000010
#define ull unsigned long long
using namespace std;
tr1::unordered_map <ull, ull> h;
ull phi[MAXN], p[MAXN], tot, sum[MAXN], n, m, ans;
bool chk[MAXN];
ull qpow(ull a, ull b) {
    ull res = 1LL;
    while (b) {
        if (b & 1) res = res * a;
        b >>= 1;
        a = a * a;
    }
    return res;
}
void seive(ull n) {
    phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!chk[i]) {
            p[++tot] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= tot && i * p[j] <= n; j++) {
            chk[i * p[j]] = 1;
            if (i % p[j]) {
                phi[i * p[j]] = phi[i] * phi[p[j]];
            }
            else {
                phi[i * p[j]] = phi[i] * p[j];
                break;
            }
        }
    }
    for (int i = 1; i <= n; i++) {
        sum[i] = sum[i - 1] + phi[i];
    }
}
ull djseive(ull n) {
    if (n <= 10000000) return sum[n];
    if (h[n]) return h[n];
    ull ans = (n & 1) ? (n + 1) / 2 * n : n / 2 * (n + 1);
    for (ull l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans -= (r - l + 1LL) * djseive(n / l);
    }
    return h[n] = ans;
}
int main() {
    scanf("%llu%llu", &n, &m);
    seive(10000000);
    for (ull l = 1, r; l <= m; l = r + 1) {
        r = m / (m / l);
        ans += (djseive(r) - djseive(l - 1)) * qpow(m / l, n); 
    }
    printf("%llu", ans);
}