SDOI2017 数字表格
题意：
题解：
答案的式子大致是这样的：
[prod_{i = 1} ^ n prod_{j = 1} ^ m f_{gcd(i, j)}
]
然后大力反演一波（这里假设(n leq m)）：
[prod_{d = 1}^nprod_{i = 1} ^ {lfloor frac{n}{d}
floor} prod_{j = 1} ^ {lfloor frac{m}{d}
floor} f_d [gcd(i, j) == 1]
]
[prod_{d = 1}^n f_d ^ {sum_{i = 1}^ {lfloor frac{n}{d}
floor} sum_{j = 1}^ {lfloor frac{m}{d}
floor} [gcd(i, j) == 1]}
]
发现指数这部分很熟悉：
[prod_{d = 1}^n f_d ^ {sum_{i = 1} ^ {lfloor frac{n}{d}
floor} mu(i) lfloor frac{n}{id}
floor lfloor frac{m}{id}
floor}
]
这样预处理出(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}
floor lfloor frac{m}{id}
floor}
]
记(g(n) = sum_{d | n} f_{frac{n}{d}}^{mu(d)})
最后的答案式就是这样：
[prod_{x = 1}^n g(x) ^{lfloor frac{n}{id}
floor lfloor frac{m}{id}
floor}
]
发现(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
", ans);
}
return 0;
}