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;
}
posted @ 2019-03-19 15:46  ImagineC  阅读(...)  评论(...编辑  收藏