[BZOJ 4816][SDOI2017]数字表格

danihao123 posted @ 2018年3月05日 07:31 in 题解 with tags 数论 SDOI bzoj 莫比乌斯反演 狄利克雷卷积 , 621 阅读
转载请注明出处:http://danihao123.is-programmer.com/

当年在考场上一点反演都不会啊啊啊啊啊啊

这题虽然牵扯到了对积的反演,但其实也不是很难。先让我们大力推导一波(逃

考虑枚举柿子中的\(f[(i, j)]\)中的最大公约数(设为\(k\)),然后式子变成了(这里不妨偷个懒,钦定\(n\leq m\)):

\[\prod_{k = 1}^n f[k]^{\sum_{i = 1}^n \sum_{j = 1}^m\:[(i, j) = k]}\]

然后发现上面那个指数就是Problem b的式子,显然可化为\(\sum_{d = 1}^n \mu(d) \lfloor\frac{n}{kd}\rfloor \lfloor\frac{m}{kd}\rfloor\)。然后似乎化简没有余地了,但是根据直觉,我们可以钦定\(T = kd\),然后枚举\(T\)。

然后柿子变成了:

\[\prod_{T = 1}^n (\prod_{k\mid T} f[k]^{\mu(\tfrac{T}{k})})^{\lfloor\tfrac{n}{T}\rfloor \lfloor\tfrac{m}{T}\rfloor}\]

柿子中的\(\prod_{k\mid T} f[k]^{\mu(\tfrac{T}{k})}\)是一个类似狄利克雷卷积的东西(取完对数就是了),可以枚举约数预处理。并且这个东西很不错的一点就是上面的指数是莫比乌斯函数,可以取到的值只有那三种,不需要快速幂。

然后剩下的部分就和通常的反演题一般无二了……

代码(卡线过的蒟蒻……求轻喷……但我在LOJ上跑的还挺快的):

/**************************************************************
    Problem: 4816
    User: danihao123
    Language: C++
    Result: Accepted
    Time:50840 ms
    Memory:33048 kb
****************************************************************/
 
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <utility>
#include <cctype>
#include <cassert>
const int N = 1000000;
const int maxn = N + 5;
typedef long long ll;
const ll ha = 1000000007LL;
bool p_flag = false;
ll pow_mod(ll a, ll b) {
  if(!b) return 1LL;
  ll ans = pow_mod(a, b / 2LL);
  ans = (ans * ans) % ha;
  if(1LL & b) ans = (ans * a) % ha;
#ifdef LOCAL
  if(p_flag) printf("%lld^%lld : %lld\n", a, b, ans);
  if(b == ha - 2LL) p_flag = false;
#endif
  return ans;
}
ll inv(ll x) {
  // if(x == 1LL) p_flag = true;
  return pow_mod(x, ha - 2LL);
}
int miu[maxn];
void sieve() {
  static int prm[maxn]; int cnt = 0;
  static bool vis[maxn];
  miu[1] = 1;
  for(int i = 2; i <= N; i ++) {
    if(!vis[i]) {
      miu[i] = -1;
      prm[cnt ++] = i;
    }
    for(int j = 0; j < cnt && prm[j] * i <= N; j ++) {
      int v = prm[j] * i;
      vis[v] = true;
      if(i % prm[j] == 0) {
        miu[v] = 0;
        break;
      } else {
        miu[v] = miu[i] * -1;
      }
    }
  }
}
ll f[maxn], inv_d[maxn], F[maxn];
void process() {
  sieve();
  f[1] = 1LL; inv_d[1] = 1LL;
  for(int i = 2; i <= N; i ++) {
    f[i] = (f[i - 1] + f[i - 2]) % ha;
    inv_d[i] = inv(f[i]);
    assert(f[i] != 0LL); assert(inv_d[i] != 0LL);
  }
  for(int i = 1; i <= N; i ++) F[i] = 1LL;
  for(int i = 1; i <= N; i ++) {
    for(int j = i; j <= N; j += i) {
      int res = j / i;
      if(miu[res] == 1) {
        F[j] = (F[j] * f[i]) % ha;
#ifdef LOCAL
        if(j <= 3) printf("f[%d](%lld) -> F[%d]\n", i, f[i], j);
#endif
        continue;
      }
      if(miu[res] == -1) {
        F[j] = (F[j] * inv_d[i]) % ha;
#ifdef LOCAL
        if(j <= 3) printf("inv_f[%d](%lld) -> F[%d]\n", i, inv_d[i], j);
#endif
      }
    }
  }
  F[0] = 1LL;
  bool flag = true;
  for(int i = 1; i <= N; i ++) {
    F[i] = (F[i] * F[i - 1]) % ha;
    if(flag && F[i] == 0LL) {
      printf("SF[%d] has been 0.\n", i);
      flag = false;
    }
  }
}
 
ll calc(ll n, ll m) {
  if(n > m) std::swap(n, m);
  ll ans = 1LL;
  for(ll i = 1; i <= n;) {
    ll nx = std::min(n / (n / i), m / (m / i));
#ifdef LOCAL
    // printf("nx of %lld : %lld\n", i, nx);
#endif
    ll mulv = pow_mod((F[nx] * inv(F[i - 1])) % ha, (n / i) * (m / i));
    ans = (ans * mulv) % ha;
    i = nx + 1LL;
  }
  return ans;
}
 
int main() {
  process();
  int T; scanf("%d", &T);
  while(T --) {
    int n, m; scanf("%d%d", &n, &m);
    printf("%lld\n", calc(n, m));
  }
  return 0;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter