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

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

题意

$$求解\sum _{i=1}^a\sum _{j=1}^b\sum _{k=1}^{c}d(i\cdot j\cdot k)$$

思路

$前置技能：d(i\cdot j)={\sum }_{x|i}{\sum }_{y|j}[gcd(x,y)=1]$

$$\sum _{i=1}^a\sum _{j=1}^b\sum _{k=1}^{c}d(i\cdot j\cdot k)$$

$$\sum _{i=1}^a\sum _{j=1}^b\sum _{k=1}^c\sum _{x|i}\sum _{y|j}\sum _{z|k}[gcd(i,j)=1][gcd(i,k)=1][gcd(j,k)=1]$$

$改变枚举顺序：$

$$\sum _{x=1}^a\sum _{y=1}^b\sum _{z=1}^c\lfloor \frac{a}{x}\rfloor \lfloor \frac{b}{y}\rfloor \lfloor \frac{c}{z}\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\lfloor \frac{a}{i}\rfloor \lfloor \frac{b}{j}\rfloor \lfloor \frac{c}{k}\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\lfloor \frac{a}{i}\rfloor \lfloor \frac{b}{j}\rfloor \lfloor \frac{c}{k}\rfloor [gcd(i,j)=1][gcd(i,k)=1]\sum _{d|jd|k}\mu (d)$$

$枚举d：$

$$\sum _{i=1}^a\lfloor \frac{a}{i}\rfloor \sum _{d=1}^{min(b,c)}\mu (d)\sum _{j=1}^{\lfloor \frac{b}{d}\rfloor }\lfloor \frac{b}{jd}\rfloor [gcd(i,jd)=1]\sum _{k=1}^{\lfloor \frac{c}{d}\rfloor }\lfloor \frac{c}{kd}\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）

#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();
}