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

题目链接

题解

为了方便，下文用 \((x, y)\) 表示 \({\rm gcd}(x, y)\)。

先分析 LOJ2476。

注意到对于任意一个数组 \(a\)，第 \(x\) 项的值 \(a_x\) 可以展开写成 \(\sum_\limits{i = 1}^{x} a_i[i = x]\)，进一步地，有：

\[\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}
\]

对于数组 \(a\)，不妨设 \(f(a) = a * \mu\)，那么原式等于：

\[\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}
\]

对于每个 \(d\)，我们可以先预处理出数组 \(P, Q, R, S, T, W\)。这样，如果 \(x\) 与 \(y\) 也已知，那么 \(P_xQ_yR_{xy}\) 可以直接得到。考虑如何求后面的部分：

\[\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)\)。

先忽略所有求 \(f(a), g(a), h(a)\) 等函数的复杂度。首先我们需要枚举 \(d\)。对于每一个 \(d\)，我们需要暴力枚举所有可能的 \(x, y\)，对于每一组合法的 \(x, y\)（即满足 \(xy \leq m, (x, y) = 1\)），再暴力枚举不超过 \(m\) 的 \(y\) 的倍数。经程序验证，当 \(n = 10^5\) 时，总枚举量为 \(98380871\)，在可接受的范围内。

对于一个长度为 \(n\) 的数组 \(a\)，求 \(f(a)\) 的时间复杂度为 \(O(n \log n)\)，求 \(g(a)\) 与 \(h(a)\) 的时间复杂度均为 \(O(n \log \log n)\)。首先，我们需要预处理出 \(f(E), f(F)\) 与 \(g(C)\)，时间复杂度为 \(O(n \log n)\)。对于每一个 \(d\)，我们需要求出 \(P, Q, R, S, T, W\) 与 \(f(W)\)，时间复杂度为 \(O(m \log m)\)。当 \(d\) 确定时，对于每一个 \(x\)，我们需要求出 \(f_1, g(f_1), f_2, h(f_2)\)，由于 \(x \leq \sqrt m\)，因此时间复杂度为 \(O(m \sqrt m \log \log m)\)。

最终时间复杂度为 \(O(n \log n)+ \sum_\limits{i = 1}^n O\left( \left(\frac{n}{i}\right) \log \left(\frac{n}{i}\right) + \left(\frac{n}{i}\right)^{1.5} \log \log\left(\frac{n}{i}\right)\right)\) \(= O(n \sqrt n \log \log n)\)。

接下来分析 LOJ2565。

首先有 \(d(ijk) = \sum_\limits{x |i}\sum_\limits{y | j}\sum_\limits{z | k}[(x, y) = 1][(y, z) = 1][(x, z) = 1]\)。证明如下：

将 \(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}
\]

定义函数 \(a_i = \left\lfloor\frac{A}{i}\right\rfloor, b_i = \left\lfloor\frac{B}{i}\right\rfloor, c_i = \left\lfloor\frac{C}{i}\right\rfloor, f_i = [i = 1]\)。那么原式即为：

\[\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];

          }

        }

      }

    }

  }

  printf("%llu\n", answer);

  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) {

      add(a[j], a[j * primes[i]]);

    }

  }

}

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) {

      add(a[j * primes[i]], a[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);

    int 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] = 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) {

              add(s1, u[j]);

              add(s2, v[j]);

            }

            add(answer, mul(s1, p[x], q[y], r[x * y]));

            if (x != y) {

              add(answer, mul(s2, p[y], q[x], r[x * y]));

            }

          }

        }

      }

    }

    printf("%d\n", answer);

  }

  return 0;

}