Codeforces-235E-Number Challenge 莫比乌斯反演 + 记忆化gcd

传送门：https://codeforces.ml/problemset/problem/235/E

题意

$$求解\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^cd(i\cdot j\cdot k)$$

思路

前置技能：$d(i\cdot j)=\sum_{xi}\sum_{yj}[gcd(x,y)=1]$

$$\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^cd(i\cdot j\cdot k)$$

$$\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c\sum_{xi}\sum_{yj}\sum_{zk}[gcd(i,j)=1][gcd(i,k)=1][gcd(j,k)=1]$$

改变枚举顺序：

$$\sum_{x=1}^a\sum_{y=1}^b\sum_{z=1}^c\left \lfloor \frac{a}{x} \right \rfloor \left \lfloor \frac{b}{y} \right \rfloor \left \lfloor \frac{c}{z} \right \rfloor [gcd(x,y)=1][gcd(x,z)=1][gcd(y,z)=1]$$

x,y,z看得不习惯，还是换成i,j,k吧。

$$\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c\left \lfloor \frac{a}{i} \right \rfloor \left \lfloor \frac{b}{j} \right \rfloor \left \lfloor \frac{c}{k} \right \rfloor [gcd(i,j)=1][gcd(i,k)=1][gcd(j,k)=1]$$

$$\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c\left \lfloor \frac{a}{i} \right \rfloor \left \lfloor \frac{b}{j} \right \rfloor \left \lfloor \frac{c}{k} \right \rfloor[gcd(i,j)=1][gcd(i,k)=1]\sum_{dj\;dk}\mu(d)$$

枚举d：

$$\sum_{i=1}^a\left \lfloor \frac{a}{i} \right \rfloor\sum_{d=1}^{min(b,c)}\mu(d)\sum_{j=1}^{\left \lfloor \frac{b}{d} \right \rfloor}\left \lfloor \frac{b}{jd} \right \rfloor[gcd(i,jd)=1]\sum_{k=1}^{\left \lfloor \frac{c}{d} \right \rfloor} \left \lfloor \frac{c}{kd} \right \rfloor[gcd(i,kd)=1]$$

枚举i和d，先使gcd(i,d)=1，在枚举j，使gcd(i,j)=1，即可使gcd(i,jd)=1 k部分同理。

在计算gcd的过程中，用f[][]记忆化gcd，会有很大的优化。

Code（622MS）

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
68
69
70
71
72
73
74
75
76
77
78
79
80
81

##include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef long double ld;
typedef pair<int, int> pdd;

##define INF 0x3f3f3f3f
##define lowbit(x) x & (-x)
##define mem(a, b) memset(a , b , sizeof(a))
##define FOR(i, x, n) for(int i = x;i <= n; i++)

// const ll mod = 998244353;
// const ll mod = 1e9 + 7;
// const double eps = 1e-6;
// const double PI = acos(-1);
// const double R = 0.57721566490153286060651209;

const ll mod = 1073741824;

const int N = 2005;

int mu[N]; // 莫比乌斯函数
bool is_prime[N];
int prime[N];
int cnt;

void Init() {
    mu[1] = 1; is_prime[0] = is_prime[1] = true;
    for(int i = 2;i < N; i++) {
        if (!is_prime[i]) {
            mu[i] = -1;
            prime[++cnt] = i;
        }
        for (int j = 1; j <= cnt && i * prime[j] < N; j++) {
            is_prime[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }

}
ll f[N][N];

ll gcd(ll a, ll b) {
    if(f[a][b]) return f[a][b];
    return f[a][b] = (b ? gcd(b, a % b) : a);
}

ll sum(int i, int d, int n) {
    if(!f[i][d]) gcd(i, d);
    if(f[i][d] != 1) return 0;
    ll res = 0;
    for(int j = 1;j <= n / d; j++) {
        if(!f[i][j]) gcd(i, j);
        if(f[i][j] == 1) res += n / (j * d);
    }
    return res;
}

void solve() {
    Init();
    int a, b, c;
    cin >> a >> b >> c;
    ll ans = 0;
    for(int i = 1;i <= a; i++) {
        for(int d = 1;d <= min(b, c); d++) {
            ans = (ans + (a / i) * mu[d] * sum(i, d, b) * sum(i, d, c)) % mod;
        }
    }
    cout << ans << endl;
}

signed main() {
    solve();
}