SDOI2017 数字表格

SDOI2017 数字表格

题意：

题目传送门

题解：

答案的式子大致是这样的：

\[\prod_{i = 1} ^ n \prod_{j = 1} ^ m f_{gcd(i, j)} \]

然后大力反演一波（这里假设\(n \leq m\)）：

\[\prod_{d = 1}^n\prod_{i = 1} ^ {\lfloor \frac{n}{d} \rfloor} \prod_{j = 1} ^ {\lfloor \frac{m}{d} \rfloor} f_d [gcd(i, j) == 1] \]

\[\prod_{d = 1}^n f_d ^ {\sum_{i = 1}^ {\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^ {\lfloor \frac{m}{d} \rfloor} [gcd(i, j) == 1]} \]

发现指数这部分很熟悉：

\[\prod_{d = 1}^n f_d ^ {\sum_{i = 1} ^ {\lfloor \frac{n}{d} \rfloor} \mu(i) \lfloor \frac{n}{id} \rfloor \lfloor \frac{m}{id} \rfloor} \]

这样预处理出\(f\)的前缀积之后，利用数论分块，我们就得到了一个\(O(n^{\frac{3}{4}} + \sqrt n log(n))\)的优秀做法了。理论上来说一千组数据应该是可以跑的，但实际上似乎常数过大而导致只有60分，似乎极限点卡一卡可以卡到90分？于是我们寻找更优秀的做法。
我们记\(x = id\)，于是原来的式子可以变成这样：

\[\prod_{x = 1}^n \prod_{d | x} f_{\frac{x}{d}} ^ {\mu(d) \lfloor \frac{n}{id} \rfloor \lfloor \frac{m}{id} \rfloor} \]

记\(g(n) = \sum_{d | n} f_{\frac{n}{d}}^{\mu(d)}\)
最后的答案式就是这样：

\[\prod_{x = 1}^n g(x) ^{\lfloor \frac{n}{id} \rfloor \lfloor \frac{m}{id} \rfloor} \]

发现\(n\)只有\(1e6\)，于是我们可以在\(O(nlog(n))\)的时间内预处理出\(g\)，然后最后的复杂度就变陈过了\(O(nlog(n) + T * \sqrt n)\)。

Code

#pragma GCC optimize(3,"inline","Ofast")
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 50;
const int Md = 1e9 + 7;
typedef long long ll;

inline int Add(const int &x, const int &y) { return (x + y >= Md) ? (x + y - Md) : (x + y);}
inline int Sub(const int &x, const int &y) { return (x - y < 0) ? (x - y + Md) : (x - y);}
inline int Mul(const int &x, const int &y) { return (ll)x * y % Md;}
inline int Min(const int &x, const int &y) { return x < y ? x : y;}
int Powe(int x, int y) {
  int ans = 1;
  while(y) {
    if(y & 1) ans = Mul(ans, x);
    x = Mul(x, x);
    y >>= 1;
  }
  return ans;
}

int n, m, cnt;
int pri[N / 10], mu[N], ntp[N], f[N], sumu[N], invf[N], g[N], dg[N], invdg[N];

void Init() {
  mu[1] = 1;
  for(int i = 2; i < N; i++) {
    if(!ntp[i]) {
      mu[i] = -1;
      pri[++cnt] = i;
    }
    for(int j = 1; j <= cnt && pri[j] * i < N; j++) {
      ntp[i * pri[j]] = 1;
      if(i % pri[j] == 0) {
	    mu[i * pri[j]] = 0;
	    break;
      }
      mu[i * pri[j]] = -mu[i];
    }
  }
  f[0] = 0; f[1] = 1, invf[1] = 1;
  for(int i = 1; i < N; i++) sumu[i] = Add(sumu[i - 1], mu[i]), g[i] = 1; 
  for(int i = 2; i < N; i++) f[i] = Add(f[i - 1], f[i - 2]), invf[i] = Powe(f[i], Md - 2);
  for(int i = 1; i < N; i++) {
    for(int j = 1; j * i < N; j++) {
      int v;
      if(mu[i] == 0) v = 1;
      else if(mu[i] == -1) v = invf[j];
      else v = f[j];
      g[i * j] = Mul(g[i * j], v);
    }
  }
  dg[0] = 1, invdg[0] = 1;
  for(int i = 1; i < N; i++) dg[i] = Mul(dg[i - 1], g[i]), invdg[i] = Powe(dg[i], Md - 2);
}

int Solve() {
  if(n > m) swap(n, m);
  int nw, ans = 1;
  for(int i = 1; i <= n; i = nw + 1) {
    nw = Min(n / (n / i), m / (m / i));
    int D = Mul(dg[nw], invdg[i - 1]);
    ans = Mul(ans, Powe(D, (ll)(n / nw) * (m / nw) % (Md - 1)));
  }
  return ans;
}

int main() {
  int T;
  scanf("%d", &T);
  Init();
  while(T--) {
    scanf("%d%d", &n, &m);
    int ans = Solve();
    printf("%d\n", ans);
  }
  return 0;
}