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

题目链接

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

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

题解

参考照搬了 wxh 的博客。

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

先分析 LOJ2476。

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

[egin{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)，那么原式等于：

[egin{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)，原式等于：

[egin{aligned}&sum_{{ m lcm}(p, q) leq n} f(E)(p) imes f(F)(q) imes g(C)({ m 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 leftlfloorfrac{n}{d} ight floor, (x, y) = 1} f(E)(xd) imes f(F)(yd) imes g(C)(xyd)sum_{ix leq leftlfloorfrac{n}{d} ight floor} sum_{jy leq leftlfloorfrac{n}{d} ight floor} A_{ixd} B_{jyd} D_{d imes (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 = leftlfloorfrac{n}{d} ight floor)，原式等于：

[egin{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}) 可以直接得到。考虑如何求后面的部分：

[egin{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)，那么上式等于：

[egin{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) imes 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} imes 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 Oleft( left(frac{n}{i} ight) log left(frac{n}{i} ight) + left(frac{n}{i} ight)^{1.5} log logleft(frac{n}{i} ight) ight)) (= 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) 的唯一分解式：

[egin{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]) 代入原式得：

[egin{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]leftlfloorfrac{A}{x} ight floor leftlfloorfrac{B}{y} ight floor leftlfloorfrac{C}{z} ight floorend{aligned} ]

定义函数 (a_i = leftlfloorfrac{A}{i} ight floor, b_i = leftlfloorfrac{B}{i} ight floor, c_i = leftlfloorfrac{C}{i} ight floor, 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
", 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
", answer);
  }
  return 0;
}