一年前看的反演神题……（然后现在才A）

考虑把原来柿子里的那个lcm展开，然后就得到了：

\[\sum_{j = 1}^n \gcd(i, j)^{c - d} i^d j^d x_j\equiv b_i\pmod{p}\]

为方便起见，令\(F(x) = x^{c-d}\)，\(G(x) = x^d\)，然后再变换一下柿子就变成了：

\[\sum_{j = 1}^n F(\gcd(i, j)) G(j)x_j\equiv \frac{b_i}{G(i)}\pmod{p}\]

然后gcd上面挂一个指数非常的烦哎……那用反演能消掉它吗？

这里我们做大多数反演题时使用的技巧都不管用了……我们这里只能采取真·莫比乌斯反演了……考虑令\(f = F\ast\mu\)（预处理莫比乌斯函数之后这个显然可以\(O(n\ln n)\)求），则有\(F = f\ast 1\)，然后可以进一步化简原式

\[
\begin{aligned}
&\quad\sum_{j = 1}^n F(\gcd(i, j)) G(j)x_j\\
&=\sum_{j = 1}^n \sum_{d | i, d | j} f(d) G(j)x_j\\
&=\sum_{d | i} f(d)\sum_{d | j} G(j)x_j
\end{aligned}
\]

然后这个有点麻烦力……暂且令\(h(d) = G(d)x_d\)，\(H(d) = \sum_{d|j} h(j)\)，然后根据另一种莫比乌斯反演柿子（本质上也是容斥）可以得到\(h(d) = \sum_{d | j} H(j)\mu(\frac{j}{d})\)，然后题目中的柿子就变成了：

\[\sum_{d | i} f(d)H(d)\equiv \frac{b_i}{G(i)}\pmod{p}\]

然后……右边是个数论函数，左边是\(f(d)H(d)\)再卷上一个1……然后这事明显的莫比乌斯反演了，可以用右边卷上μ来得到\(f(d)H(d)\)，这下我们就可以求出\(H(d)\)了，那\(h(d)\)我们也知道了，\(x_d\)我们自然也就全都知道了……

至于无解的情况，原因在于可能会有个别f值在膜p意义下为0。然后我们在求\(H(d)\)的时候，如果发现\(f(d)H(d) = 0\)，那么\(H(d)\)其实可以取任意值。但如果\(f(d)H(d)\neq 0\)，那就矛盾了，无解。

代码：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cctype>
#include <algorithm>
#include <utility>
const int maxn = 100005;
using ll = long long;
const ll ha = 998244353LL;
ll ppp(ll a, ll b) {
  ll ans = 1LL, res = a % ha;
  while(b) {
    if(1LL & b) ans = (ans * res) % ha;
    res = (res * res) % ha;
    b >>= 1;
  }
  return ans;
}
ll inv(ll x) {
  return ppp(x, ha - 2LL);
}
ll pow_mod(ll a, ll b) {
  if(b < 0LL) {
    return inv(ppp(a, -b));
  } else {
    return ppp(a, b);
  }
}

ll mu[maxn];
void sieve() {
  static int prm[maxn];
  static bool vis[maxn];
  const int n = 100000;
  int cnt = 0;
  mu[1] = 1LL; vis[1] = true;
  for(int i = 2; i <= n; i ++) {
    if(!vis[i]) {
      prm[cnt ++] = i;
      mu[i] = ha - 1LL;
    }
    for(int j = 0; j < cnt && i * prm[j] <= n; j ++) {
      int v = i * prm[j];
      vis[v] = true;
      if(i % prm[j] == 0) {
        mu[v] = 0LL; break;
      } else {
        mu[v] = (mu[i] * (ha - 1LL)) % ha;
      }
    }
  }
}

int n, c, d;
ll F[maxn], f[maxn];
void process_f() {
  for(int i = 1; i <= n; i ++) {
    F[i] = pow_mod(i, c - d);
  }
  for(int i = 1; i <= n; i ++) {
    for(int j = i; j <= n; j += i) {
      f[j] += (F[i] * mu[j / i]) % ha;
      f[j] %= ha;
    }
  }
}

ll b[maxn];
ll E[maxn], e[maxn];
void process_e() {
  std::fill(E + 1, E + 1 + n, 0LL);
  std::fill(e + 1, e + 1 + n, 0LL);
  for(int i = 1; i <= n; i ++) {
    E[i] = (b[i] * pow_mod(i, -d)) % ha;
  }
  for(int i = 1; i <= n; i ++) {
    for(int j = i; j <= n; j += i) {
      e[j] += (E[i] * mu[j / i]) % ha;
      e[j] %= ha;
    }
  }
}

bool fl;
ll h[maxn], x[maxn];
void process_x() {
  std::fill(h + 1, h + 1 + n, 0LL);
  std::fill(x + 1, x + 1 + n, 0LL);
  fl = false;
  for(int i = 1; i <= n; i ++) {
    if(f[i]) {
      h[i] = (e[i] * inv(f[i])) % ha;
    } else {
      if(e[i]) {
        fl = true; return;
      } else {
        h[i] = 0LL;
      }
    }
  }
  for(int i = 1; i <= n; i ++) {
    for(int j = i; j <= n; j += i) {
      x[i] += (h[j] * mu[j / i]) % ha;
      x[i] %= ha;
    }
  }
  for(int i = 1; i <= n; i ++) {
    x[i] = (x[i] * pow_mod(i, -d)) % ha;
  }
}

int main() {
  int q; scanf("%d%d%d%d", &n, &c, &d, &q);
  sieve(); process_f();
  while(q --) {
    for(int i = 1; i <= n; i ++) {
      scanf("%lld", &b[i]);
    }
    process_e();
    process_x();
    if(fl) {
      puts("-1");
    } else {
      for(int i = 1; i <= n; i ++) {
        if(i > 1) putchar(' ');
        printf("%lld", x[i]);
      }
      putchar('\n');
    }
  }
  return 0;
}

• 相关文章
• 全局相关文章