Luogu 3704 [SDOI2017]数字表格
列一下式子:
$\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} \right \rfloor}\sum_{j = 1}^{\left \lfloor \frac{m}{d} \right \rfloor}[gcd(i, j) == 1]}$
右上角的那个东西:
$\sum_{i = 1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j = 1}^{\left \lfloor \frac{m}{d} \right \rfloor}[gcd(i, j) == 1]$
太熟悉了。
$\sum_{d = 1}^{min(n, m)}\sum_{t = 1}^{min(\left \lfloor \frac{n}{d} \right \rfloor,\left \lfloor \frac{m}{d} \right \rfloor)}\mu (t) * \left \lfloor \frac{n}{td} \right \rfloor * \left \lfloor \frac{m}{td} \right \rfloor$
代回去之后按照套路枚举$T = dt$:
$\prod_{T = 1}^{min(n, m)}\sum_{d | T}fib^{\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor \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} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor}$。
发现外面已经可以整除分块了。
然而这个$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("\n"); */ 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("\n"); */ 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("\n"); */ 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\n", ans); } return 0; }