# LOJ2476. 「2018 集训队互测 Day 3」蒜头的奖杯 & LOJ2565. 「SDOI2018」旧试题（莫比乌斯反演）

## 题目链接

LOJ2476：https://loj.ac/problem/2476

LOJ2565：https://loj.ac/problem/2565

## 题解

\begin{aligned}a_x &= \sum_{i = 1}^{x} a_i [i = x] \\ &= \sum_{i | x}a_i[\frac{x}{i} = 1] \\ &= \sum_{i |x} a_i \sum_{p | \frac{x}{i}} \mu (p) \\ &= \sum_{p |x} \sum_{i | p} a_i \mu(\frac{p}{i}) \\ &= \sum_{p |x} (a * \mu)(p) \end{aligned}

\begin{aligned}\sum_{k = 1}^n C_k \sum_{i = 1}^n \sum_{j = 1}^n A_iB_j D_{(i, j)} \sum_{p | (i, k)} f(E)(p)\sum_{q | (j, k)}\end{aligned} f(F)(q)

$g(a)(x) = \sum_\limits{x | i}a_i$，原式等于：

\begin{aligned}&\sum_{{\rm lcm}(p, q) \leq n} f(E)(p) \times f(F)(q) \times g(C)({\rm lcm}(p, q)) \sum_{ip \leq n}A_{ip} \sum_{jq \leq n}B_{jq} D_{(ip, jq)} \\ =& \sum_{d = 1}^n \sum_{xy \leq \left\lfloor\frac{n}{d}\right\rfloor, (x, y) = 1} f(E)(xd) \times f(F)(yd) \times g(C)(xyd)\sum_{ix \leq \left\lfloor\frac{n}{d}\right\rfloor} \sum_{jy \leq \left\lfloor\frac{n}{d}\right\rfloor} A_{ixd} B_{jyd} D_{d \times (ix, jy)}\end{aligned}

$d$ 已知时，令 $P_i = f(E)(id), Q_i = f(F)(id), R_i = g(C)(id), S_i = A_{id}, T_i = B_{id}, W_i = D_{id}$。设 $m = \left\lfloor\frac{n}{d}\right\rfloor$，原式等于：

\begin{aligned} \sum_{d = 1}^n \sum_{xy \leq m, (x, y) = 1} P_x Q_y R_{xy} \sum_{ix \leq m} \sum_{jy \leq m} S_{ix} T_{jy} W_{(ix, jy)} \end{aligned}

\begin{aligned} &\sum_{ix \leq m} \sum_{jy \leq m} S_{ix} T_{jy} W_{(ix, jy)} \\ = & \sum_{ix \leq m} S_{ix} \sum_{jy \leq m} T_{jy} \sum_{r | (ix, jy)} f(W)(r) \\ = & \sum_{jy \leq m} T_{jy} \sum_{r | jy} f(W)(r) \sum_{r | ix} S_{ix}\end{aligned}

$u = ix, v= jy$，那么上式等于：

\begin{aligned}\sum_{y | v} T_{v} \sum_{r | v} f(W)(r) \sum_{ r | u} [x | u]S_{u}\end{aligned}

$h(a)(x) = \sum_\limits{i | x} a_i$，那么上式可以从右往左依次计算，具体地，令函数 $f_1(u) = [x | u]S_{u}$，那么我们可以先计算 $g(f_1)$ 来得到每一个 $r$ 对应的 $\sum_\limits{r | u}[x | u]S_{u}$，再令函数 $f_2(r) = f(W)(r) \times g(f_1)(r)$，那么我们可以计算 $h(f_2)$ 来得到每一个 $v$ 对应的 $\sum_\limits{r | v} f(W)(r) \sum_\limits{r | u} [x | u]S_{u}$。这样，上式即为枚举所有 $y$ 的倍数 $v$，求 $T_{v} \times h(f_2)(v)$

$s = ijk$。分别写出 $i, j, k, s$ 的唯一分解式：

\begin{aligned}i &= p_1^{a_1}p_2^{a_2} \cdots p_t^{a_t} \\j &= p_1^{b_1}p_2^{b_2} \cdots p_t^{b_t} \\ k &= p_1^{c_1}p_2^{c_2} \cdots p_t^{c_t} \\ s &= p_1^{a_1 + b_1 + c_1}p_2^{a_2 + b_2 + c_2} \cdots p_t^{a_t + b_t + c_t}\end{aligned}

$d(ijk) = \sum_\limits{x |i}\sum_\limits{y | j}\sum_\limits{z | k}[(x, y) = 1][(y, z) = 1][(x, z) = 1]$ 代入原式得：

\begin{aligned} \sum_{i = 1}^A \sum_{j = 1}^B \sum_{k = 1}^C d(ijk) &= \sum_{i = 1}^A \sum_{j = 1}^B \sum_{k = 1}^C \sum_{x | i} \sum_{y | j}\sum_{z|k} [(x, y) = 1][(y, z) = 1][(x, z) = 1] \\ &= \sum_{x = 1}^A \sum_{y = 1}^B \sum_{z = 1}^C [(x, y) = 1][(y, z) = 1][(x, z) = 1]\left\lfloor\frac{A}{x}\right\rfloor \left\lfloor\frac{B}{y}\right\rfloor \left\lfloor\frac{C}{z}\right\rfloor\end{aligned}

$\sum_{x = 1}^A \sum_{y = 1}^B \sum_{z = 1}^C a_x b_y c_z f_{(x, y)}f_{(y, z)}f_{(x, z)}$

## 代码

LOJ2476 代码如下：

#include<bits/stdc++.h>

using namespace std;

const int N = 1e5 + 10;

unsigned long long a[N], b[N], c[N], d[N], e[N], f[N], p[N], q[N], r[N], s[N], t[N], w[N], u[N], v[N];
int n, primes[N], all;
bool is_prime[N];

void F(unsigned long long* a, int n) {
for (int i = 1; i <= n; ++i) {
for (int j = i + i; j <= n; j += i) {
a[j] -= a[i];
}
}
}

void G(unsigned long long* a, int n) {
for (int i = 1; i <= all && primes[i] <= n; ++i) {
for (int j = n / primes[i]; j; --j) {
a[j] += a[j * primes[i]];
}
}
}

void H(unsigned long long* a, int n) {
for (int i = 1; i <= all && primes[i] <= n; ++i) {
for (int j = 1; j * primes[i] <= n; ++j) {
a[j * primes[i]] += a[j];
}
}
}

void sieve(int n) {
fill(is_prime + 1, is_prime + 1 + n, true);
for (int i = 2; i <= n; ++i) {
if (is_prime[i]) {
primes[++all] = i;
}
for (int j = 1; j <= all; ++j) {
if (primes[j] * i > n) {
break;
}
is_prime[primes[j] * i] = false;
if (i % primes[j] == 0) {
break;
}
}
}
}

int main() {
scanf("%d", &n);
for (int i = 1; i <= n; ++i) {
scanf("%llu", &a[i]);
}
for (int i = 1; i <= n; ++i) {
scanf("%llu", &b[i]);
}
for (int i = 1; i <= n; ++i) {
scanf("%llu", &c[i]);
}
for (int i = 1; i <= n; ++i) {
scanf("%llu", &d[i]);
}
for (int i = 1; i <= n; ++i) {
scanf("%llu", &e[i]);
}
for (int i = 1; i <= n; ++i) {
scanf("%llu", &f[i]);
}
sieve(n);
F(e, n);
F(f, n);
G(c, n);
unsigned long long answer = 0;
for (int i = 1; i <= n; ++i) {
int m = n / i;
for (int j = 1; j <= m; ++j) {
p[j] = e[j * i];
q[j] = f[j * i];
r[j] = c[j * i];
s[j] = a[j * i];
t[j] = b[j * i];
w[j] = d[j * i];
}
F(w, m);
for (int x = 1; x * x <= m; ++x) {
fill(u + 1, u + 1 + m, 0);
fill(v + 1, v + 1 + m, 0);
for (int j = x; j <= m; j += x) {
u[j] = s[j], v[j] = t[j];
}
G(u, m);
G(v, m);
for (int j = 1; j <= m; ++j) {
u[j] *= w[j];
v[j] *= w[j];
}
H(u, m);
H(v, m);
for (int j = 1; j <= m; ++j) {
u[j] *= t[j];
v[j] *= s[j];
}
for (int y = x; x * y <= m; ++y) {
if (__gcd(x, y) == 1) {
unsigned long long s1 = 0, s2 = 0;
for (int j = y; j <= m; j += y) {
s1 += u[j];
s2 += v[j];
}
answer += s1 * p[x] * q[y] * r[x * y];
if (x != y) {
answer += s2 * p[y] * q[x] * r[x * y];
}
}
}
}
}
return 0;
}


LOJ2565 代码如下：

#include<bits/stdc++.h>

using namespace std;

const int N = 1e5 + 10, mod = 1e9 + 7;

void add(int& x, int y) {
x += y;
if (x >= mod) {
x -= mod;
}
}

void sub(int& x, int y) {
x -= y;
if (x < 0) {
x += mod;
}
}

template<typename T>
int mul(T x) {
return x;
}

template<typename T, typename... R>
int mul(T x, R... y) {
return (long long) x * mul(y...) % mod;
}

int tt, A, B, C, a[N], b[N], c[N], d[N], e[N], f[N], p[N], q[N], r[N], s[N], t[N], w[N], u[N], v[N], primes[N], all;
bool is_prime[N];

void F(int* a, int n) {
for (int i = 1; i <= n; ++i) {
for (int j = i + i; j <= n; j += i) {
sub(a[j], a[i]);
}
}
}

void G(int* a, int n) {
for (int i = 1; i <= all && primes[i] <= n; ++i) {
for (int j = n / primes[i]; j; --j) {
}
}
}

void H(int* a, int n) {
for (int i = 1; i <= all && primes[i] <= n; ++i) {
for (int j = 1; j * primes[i] <= n; ++j) {
}
}
}

void sieve(int n) {
fill(is_prime + 1, is_prime + 1 + n, true);
all = 0;
for (int i = 2; i <= n; ++i) {
if (is_prime[i]) {
primes[++all] = i;
}
for (int j = 1; j <= all; ++j) {
if (primes[j] * i > n) {
break;
}
is_prime[primes[j] * i] = false;
if (i % primes[j] == 0) {
break;
}
}
}
}

int main() {
scanf("%d", &tt);
while (tt--) {
scanf("%d%d%d", &A, &B, &C);
int n = max(A, max(B, C));
sieve(n);
for (int i = 1; i <= n; ++i) {
a[i] = A / i;
b[i] = B / i;
c[i] = C / i;
d[i] = e[i] = f[i] = (i == 1);
}
F(e, n);
F(f, n);
G(c, n);
for (int i = 1; i <= n; ++i) {
int m = n / i;
for (int j = 1; j <= m; ++j) {
p[j] = e[j * i];
q[j] = f[j * i];
r[j] = c[j * i];
s[j] = a[j * i];
t[j] = b[j * i];
w[j] = d[j * i];
}
F(w, m);
for (int x = 1; x * x <= m; ++x) {
fill(u + 1, u + 1 + m, 0);
fill(v + 1, v + 1 + m, 0);
for (int j = x; j <= m; j += x) {
u[j] = s[j], v[j] = t[j];
}
G(u, m);
G(v, m);
for (int j = 1; j <= m; ++j) {
u[j] = mul(u[j], w[j]);
v[j] = mul(v[j], w[j]);
}
H(u, m);
H(v, m);
for (int j = 1; j <= m; ++j) {
u[j] = mul(u[j], t[j]);
v[j] = mul(v[j], s[j]);
}
for (int y = x; x * y <= m; ++y) {
if (__gcd(x, y) == 1) {
int s1 = 0, s2 = 0;
for (int j = y; j <= m; j += y) {
}