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

题目链接

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

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

题解

参考照搬了 wxh 的博客。

为了方便，下文用 \((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]\)。证明如下：

令 \(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) = \prod_\limits{r = 1}^t(a_r + b_r + c_r + 1)\)。

考虑计算 \(\sum_\limits{x |i}\sum_\limits{y | j}\sum_\limits{z | k}[(x, y) = 1][(y, z) = 1][(x, z) = 1]\) 的值。不难发现，满足 \(x | i, y | j, z | k\) 且 \((x, y) = 1, (y, z) = 1, (x, z) = 1\) 的 \(x, y, z\) 必然有：\(\forall r(1 \leq r \leq t)\)，\(x, y, z\) 中至少有两个数满足其唯一分解式中 \(p_r\) 的指数为 \(0\)。当 \(x, y, z\) 中至少有两个数满足其唯一分解式中 \(p_r\) 的指数为 \(0\) 时，\(x, y, z\) 三个数的唯一分解式中 \(p_r\) 的指数共有 \(a_r + b_r + c_r + 1\) 种情况。显然，对于各个 \(r\)，情况是独立的，因此根据乘法原理，满足 \(x | i, y | j, z | k\) 且 \((x, y) = 1, (y, z) = 1, (x, z) = 1\) 的 \(x,y, z\) 共有 \(\prod_\limits{r = 1}^t (a_r + b_r + c_r + 1)\) 组。其数量恰好等于 \(d(ijk)\) 的值。

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