HDU 6706 - huntian oy

发表于 更新于 分类于 阅读次数：

题目链接：HDU 6706

思路

\[ {\rm gcd}(i,j) = 1, i > j \rightarrow {\rm gcd}(i^a - j^a, i^b - j^b) = i^{ {\rm gcd} (a,b)} - j^{ {\rm gcd} (a,b)} \]

\[ \begin{aligned} f(n, a, b) = & \sum_{i = 1}^{n} \sum_{j = 1}^{i} {\rm gcd} (i^a - j^a, i^b - j^b)[ {\rm gcd} (i,j)=1] \\ = & \sum_{i = 1}^{n} \sum_{j=1}^{i} (i-j)[{\rm gcd}(i, j)=1] \\ = & \sum_{i = 1}^{n} \sum_{j = 1}^{i} i [{\rm gcd}(i, j) = 1] - \sum_{i = 1}^{n} \sum_{j = 1}^{i} j[{\rm gcd}(i, j) = 1]  \end{aligned} \]

被减数 \(\sum_{i = 1}^{n} \sum_{j=1}^{i} = \sum_{i = 1}^{n} i \sum_{j = 1}^{i} [{\rm gcd}(i,j) = 1] = \sum_{i = 1}^{n} i \varphi(i)\)。

减数考虑\({\rm gcd}(i,j) = {\rm gcd}(i, i - j)\)，故如果存在一个有贡献的\(j'\)，一定存在一个有贡献的\(i - j'\)，他们的贡献和是\(i\)。对于每一个这样的\((i)\)都有\(\varphi(i)\)组\((i,j)\)，去重之后可以表示为\(\sum_{i = 1}^{n} (\frac{i \varphi(i)}{2} + [i=1])\)。

所以我们可以推得：

\[f(n, a, b) = \frac{\sum_{i = 1}^{n} i \varphi(i) - 1}{2}\]

\(f(x) = x \varphi(x)\)是积性函数，考虑使用筛法。设\(g(x) = x\)，令\(f\)和\(g\)做Direclet卷积：

\[ f * g (x) = \sum_{d|x} \frac{x}{d} \varphi(\frac{x}{d}) d = \sum_{d|n} x \varphi(\frac{x}{d}) = x \sum_{d|n} \varphi(\frac{x}{d}) = x^2  \]

可求前缀和，记\(\sum_{x=1}^{n} f(x) = F(n)\)，也就是我们要求的值，那么

\[ \begin{aligned} \frac{n(n+1)(2n + 1)}{6} = \sum_{x = 1}^{n} f * g(x) =&  \sum_{x = 1}^{n} \sum_{d|x} x \varphi(\frac{x}{d}) \\ = & \sum_{d=1}^{n}\sum_{x=1}^{n} x \varphi(\frac{x}{d}) [d | x] \\ = & \sum_{d=1}^{n}  \sum_{x = 1}^{\lfloor \frac{n}{d} \rfloor} xd\varphi(x) \\ = & \sum_{d=1}^{n} d \sum_{x = 1}^{\lfloor \frac{n}{d} \rfloor} x \varphi(x) \\ = & F(n) + \sum_{d = 2}^{n} d F(\lfloor \frac{n}{d} \rfloor) \end{aligned} \]

所以我们可以得到

\[F(n) = \frac{n(n+1)(2n + 1)}{6} -\sum_{d = 2}^{n} d F(\lfloor \frac{n}{d} \rfloor)\]

杜教筛即可。

实测结果预处理前\(n \leq 3.5 \times 10^6\)的\(F(n)\)是最快的，在HDU的老爷机上跑了998 MS。

预处理线性筛，对于积性函数\(f(x) = x \varphi(x)\)：

  • \(f(p) = p(p - 1)\)
  • \(f(p^k) = p^k \varphi(p^k)\)，\(f(p^{k+1}) = p^{k+1} \varphi(p^{k+1}) = p^{k+1} \varphi(p^k) p = p^{k+2} \varphi(p^k)\)，所以\(f(p^{k+1}) = p^2 f(p^k)\)

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int MAX_LEN = int(3E6) + int(5E5) + 5;
const int P = int(1E9) + 7;
LL getPow(LL x, LL y) {
    LL ret = 1;
    while (y) {
        if (y & 1) ret = ret * x % P;
        x = x * x % P;
        y >>= 1;
    }
    return ret;
}
LL inv_6, inv_2;
LL n, T;
bool vis[MAX_LEN];
LL tot, p[MAX_LEN], low[MAX_LEN];
LL f[MAX_LEN];
void linearSieve() {
    vis[0] = vis[1] = 1; low[1] = f[1] = 1;
    for (LL i = 2; i < MAX_LEN; ++i) {
        if (!vis[i]) {
            low[i] = p[++tot] = i;
            f[i] = 1LL * i * (i - 1) % P;
        }
        for (LL j = 1; j <= tot && 1LL * i * p[j] < MAX_LEN; ++j) {
            vis[i * p[j]] = 1;
            if (i % p[j] == 0) {
                low[i * p[j]] = low[i] * p[j];
                if (low[i] == i) f[i * p[j]] = ((1LL * p[j] * p[j]) % P * f[i]) % P;
                else f[i * p[j]] = (f[i / low[i]] * f[p[j] * low[i]]) % P;
                break;
            } else {
                f[i * p[j]] = (f[i] * f[p[j]]) % P;
                low[i * p[j]] = p[j];
            }
        }
    }
    for (LL i = 1; i < MAX_LEN; ++i) f[i] = (f[i] + f[i - 1]) % P;
}
unordered_map <LL, LL> hashS;
LL sum(LL x) {
    if (x < MAX_LEN) return f[x];
    if (hashS.find(x) != hashS.end()) return hashS[x];
    LL ret = (1LL * x * (x + 1) % P) * (1LL * (2 * x + 1) * inv_6 % P) % P;
    LL l, r;
    for (l = 2; l <= x; l = r + 1) {
        r = x / (x / l);
        ret -= (((l + r) * (r - l + 1) % P * inv_2) % P * sum(x / l));
        ret = ((ret % P) + P) % P;
    }
    return hashS[x] = ret;
}
int main() {
    scanf("%lld", &T);
    linearSieve();
    inv_6 = getPow(6, P - 2);
    inv_2 = getPow(2, P - 2);
    int ta, tb;
    for (int cs = 1; cs <= T; ++cs) {
        scanf("%lld%lld%lld", &n, &ta, &tb);
        LL ans = ((sum(n) - 1) % P + P) % P * inv_2 % P;
        printf("%lld\n", ans);
    }
    return 0;
}