「SDOI2017」数字表格

知识点：莫比乌斯反演

学到了 \(\prod\) 的交换。
另外莫比乌斯反演的代码真好写（
只用预处理一吨东西就做完了（

简述

原题面：Luogu，时限 5s 的 Loj。

设 \(f_{i}\) 表示斐波那契数列的第 \(i\) 项。
\(T\) 组数据，每次给定参数 \(n,m\)，求：

\[\prod_{i=1}^{n}\prod_{j=1}^{m}f_{\gcd(i,j)} \pmod {10^9 + 7} \]

\(1\le T\le 10^3\)，\(1\le n,m\le 10^6\)。
2S~3S，128MB。

这里的时限是 Luogu 的时限，需要略微卡常才能过。
不想卡常可以移步 Loj。

分析

以下钦定 \(n\ge m\)。
大力化式子，先套路地枚举 \(\gcd(i,j)\)，用初中知识把两个 \(\prod\) 化到指数位置，原式等于：

\[\large\begin{aligned} &\prod_{d = 1}^{n}\prod_{i=1}^{n}\prod_{j=1}^{m}f_{d}[\gcd(i,j)=d]\\ =&\prod_{d = 1}^{n}f_{d}^{\left(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=d]\right)} \end{aligned}\]

分母套路一波，有：

\[\begin{aligned} &\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j) = d]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[\gcd(i,j) = 1]\\ =& \sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum_{k\mid \gcd(i,j)}\mu (k)\\ =& \sum_{k=1}\mu(k)\left\lfloor\dfrac{n}{kd}\right\rfloor\left\lfloor\dfrac{m}{kd}\right\rfloor \end{aligned}\]

代回原式，原式等于：

\[\large \begin{aligned} &\prod_{d = 1}^{n}f_{d}^{\left(\sum\limits_{k=1}\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor\right)}\\ =& \prod_{d = 1}^{n}\left(f_{d}^{{\sum\limits_{k=1}\mu(k)}}\right)^{\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor} \end{aligned}\]

考虑再暴力拆一波，原式等于：

\[\large \begin{aligned} &\prod_{d = 1}^{n}\left(f_{d}^{{\sum\limits_{k=1}\mu(k)}}\right)^{\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor}\\ =& \prod_{d = 1}^{n}\left(\prod_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f_{d}^{\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor}\right) \end{aligned}\]

做不动了，但发现变量仅有 \(k,d,kd\)，考虑更换枚举对象改为枚举 \(t = kd\) 与 \(d\)，则原式等于：

\[\large\prod_{t=1}^{n}\left(\prod_{d | t}^{n}f_{d}^{{\mu(\frac{t}{d})}}\right)^{\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor} \]

枚举对象变成了约数形式。从后面的式子推前面的式子是比较显然的，可以认为这种枚举 \(t=kd\) 的形式是一种逆向操作。

设：

\[\large g(t)=\prod_{d | t}^{n}f_{d}^{{\mu(\frac{t}{d})}} \]

\(g(t)\) 可以用类似埃氏筛的方法 \(O(n\log n)\) 地预处理出来。再把 \(g\) 代回原式，原式等于：

\[\large\prod_{t=1}^{n}g(t)^{\left\lfloor\frac{n}{t}\right\rfloor\left\lfloor\frac{m}{t}\right\rfloor} \]

可以考虑预处理 \(g(t)\) 的前缀积，数论分块枚举指数求解即可。

总时间复杂度 \(O(n\log n + T\sqrt n)\)，轻微卡常，多预处理一些东西可以过。

实现

//知识点：莫比乌斯反演 
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cstdio>
#include <cstring>
#define LL long long
const LL mod = 1e9 + 7;
const int kN = 1e6;
//=============================================================
LL n, m, ans;
int p_num, p[kN + 10];
bool vis[kN + 10];
LL mu[kN + 10], f[kN + 10], g[kN + 10], prod[kN + 10];
LL invf[kN + 10], invp[kN];
//=============================================================
inline int read() {
  int f = 1, w = 0;
  char ch = getchar();
  for (; !isdigit(ch); ch = getchar())
    if (ch == '-') f = -1;
  for (; isdigit(ch); ch = getchar()) {
    w = (w << 3) + (w << 1) + (ch ^ '0');
  }
  return f * w;
}
void Chkmax(int &fir_, int sec_) {
  if (sec_ > fir_) fir_ = sec_;
}
void Chkmin(int &fir_, int sec_) {
  if (sec_ < fir_) fir_ = sec_;
}
LL QPow(LL x_, LL y_) {
  x_ %= mod;
  LL ret = 1;
  for (; y_; y_ >>= 1ll) {
    if (y_ & 1) ret = ret * x_ % mod;
    x_ = x_ * x_ % mod;
  }
  return ret;
}
void Euler() {
  vis[1] = true, mu[1] = 1;
  for (int i = 2; i <= kN; ++ i) {
    if (! vis[i]) {
      p[++ p_num] = i;
      mu[i] = -1;
    }
    for (int j = 1; j <= p_num && i * p[j] <= kN; ++ j) {
      vis[i * p[j]] = true;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = -mu[i];
    }
  }
}
void Prepare() {
  g[1] = g[2] = 1;
  f[1] = f[2] = 1;
  invf[1] = invf[2] = 1;
  for (int i = 3; i <= kN; ++ i) {
    g[i] = 1;
    f[i] = (f[i - 1] + f[i - 2]) % mod;
    invf[i] = QPow(f[i], mod - 2);
  }
  
  Euler();
  for (int d = 1; d <= kN; ++ d) {
    for (int j = 1; d * j <= kN; ++ j) {
      if (mu[j] == 1) {
        g[d * j] = g[d * j] * f[d] % mod;
      } else if (mu[j] == -1) {
        g[d * j] = g[d * j] * invf[d] % mod;
      }
    }
  }
  invp[0] = prod[0] = 1;
  for (int i = 1; i <= kN; ++ i) {
    prod[i] = prod[i - 1] * g[i] % mod;
    invp[i] = QPow(prod[i], mod - 2);
  }
}
//=============================================================
int main() {
  Prepare();
  int T = read();
  while (T -- ){
    n = read(), m = read(), ans = 1;
    if (n < m) std::swap(n, m);
    for (LL l = 1, r = 1; l <= m; l = r + 1) {
      r = std::min(n / (n / l), m / (m / l));
      ans = (ans * QPow(prod[r] * invp[l - 1] % mod, 1ll * (n / l) * (m / l))) % mod;
    }
    printf("%lld\n", ans);
  }
  return 0;
}