列一下式子:
$prod_{i = 1}^{n}prod_{j = 1}^{m}fib_{gcd(i, j)}$
很套路的变成这样:
$prod_{d = 1}^{min(n, m)}fib_{d}^{sum_{i = 1}^{left lfloor frac{n}{d} ight floor}sum_{j = 1}^{left lfloor frac{m}{d} ight floor}[gcd(i, j) == 1]}$
右上角的那个东西:
$sum_{i = 1}^{left lfloor frac{n}{d} ight floor}sum_{j = 1}^{left lfloor frac{m}{d} ight floor}[gcd(i, j) == 1]$
太熟悉了。
$sum_{d = 1}^{min(n, m)}sum_{t = 1}^{min(left lfloor frac{n}{d} ight floor,left lfloor frac{m}{d} ight floor)}mu (t) * left lfloor frac{n}{td} ight floor * left lfloor frac{m}{td} ight floor$
代回去之后按照套路枚举$T = dt$:
$prod_{T = 1}^{min(n, m)}sum_{d | T}fib^{left lfloor frac{n}{T} ight floorleft lfloor frac{m}{T} ight floor mu (frac{T}{d})}_{d}$
这样子的话我们记$h(i) = sum_{d | T}fib^{mu (frac{T}{d})}_{d}$
原式就变为:
$prod_{T = 1}^{min(n, m)}h(T)^{left lfloor frac{n}{T} ight floorleft lfloor frac{m}{T} ight floor}$。
发现外面已经可以整除分块了。
然而这个$h(i)$怎么办,这玩意...不能线性筛的呀。
喂喂,不能线性筛就暴力算吧,暴力...似乎并不慢啊,其实是一个$O(nlogn)$。
注意到$mu (i) == -1$的时候其实是乘上一个逆元。
时间复杂度$O(MaxNlogMaxN + T sqrt{n} )$。
复杂度写的并不严格。
Code:
#include <cstdio> #include <cstring> using namespace std; typedef long long ll; const int N = 1e6 + 5; const ll P = 1e9 + 7; int testCase, pCnt = 0, pri[N]; ll mu[N], fib[N], h[N]; bool np[N]; template <typename T> inline void read(T &X) { X = 0; char ch = 0; T op = 1; for(; ch > '9' || ch < '0'; ch = getchar()) if(ch == '-') op = -1; for(; ch >= '0' && ch <= '9'; ch = getchar()) X = (X << 3) + (X << 1) + ch - 48; X *= op; } inline ll pow(ll x, ll y) { ll res = 1LL; for(; y > 0; y >>= 1) { if(y & 1) res = res * x % P; x = x * x % P; } return res; } ll inv(ll x) { return pow(x, P - 2); } void sieve() { fib[1] = 1LL; for(int i = 2; i < N; i++) fib[i] = (fib[i - 1] + fib[i - 2]) % P; /* for(int i = 1; i <= 20; i++) printf("%lld ", fib[i]); printf(" "); */ mu[1] = 1; for(int i = 2; i < N; i++) { if(!np[i]) pri[++pCnt] = i, mu[i] = -1; for(int j = 1; j <= pCnt && i * pri[j] < N; j++) { np[i * pri[j]] = 1; if(i % pri[j] == 0) { mu[i * pri[j]] = 0; break; } mu[i * pri[j]] = -mu[i]; } } /* for(int i = 1; i <= 20; i++) printf("%d ", mu[i]); printf(" "); */ for(int i = 0; i < N; i++) h[i] = 1LL; for(int i = 1; i < N; i++) { if(!mu[i]) continue; for(int j = i; j < N; j += i) h[j] = h[j] * ((mu[i] == 1) ? fib[j / i] : inv(fib[j / i])) % P; } /* for(int i = 1; i <= 20; i++) printf("%lld ", h[i]); printf(" "); */ for(int i = 1; i < N; i++) h[i] = 1LL * h[i] * h[i - 1] % P; } inline int min(int x, int y) { return x > y ? y : x; } int main() { sieve(); for(read(testCase); testCase--; ) { int n, m; read(n), read(m); int rep = min(n, m); ll ans = 1LL; for(int l = 1, r; l <= rep; l = r + 1) { r = min((n / (n / l)), (m / (m / l))); ans = ans * pow(h[r] * inv(h[l - 1]) % P, 1LL * (n / l) * (m / l) % (P - 1)) % P; } printf("%lld ", ans); } return 0; }