[BZOJ4816]数字表格

题目链接：BZOJ4816

题目大意

求∏1≤i≤n,1≤j≤mF(gcd(i,j))，答案对1e9+7取模，其中F为fibonacci数列，F(0)=0,F(1)=1。

分析

1. 枚举gcd，转化为求∏nd=1F(d)h(d)，其中h(d)=∑ni=1∑mj=1[gcd(i,j)==d]。

2. 而求h(d)作为莫比乌斯反演的经典问题，我们知道结论，于是变成了∏nd=1F(d)(∑⌊nd⌋k=1μ(k)⌊ndk⌋⌊mdk⌋)。

3. 设T=dk，同时外移T，就变成了∏nT=1(∏d∣TF(d)μ(Td))⌊nT⌋⌊mT⌋。

4. 设g(T)=∏d∣TF(d)μ(Td)，那么预处理出g(x)的前缀积以及前缀积的逆元即可N−−√logN求出每次询问。

5. 那么如何求g(x)呢？首先初始化g(x)=1，我们枚举倍数d，对所有d∣T，g(T)=g(T)∗F(d)μ(Td)；因此，我们在求出g之前，需要先筛出μ，求出F，以及F的逆元。

6. 然而因为模数过大，不能O(P)预处理逆元，于是我们需要另一种方法求数列的逆元：对于数列a1,a2,a3...an，我们设Si为前i项积；那么我们可以O(logP)求出Sn的逆元invs[n]=(a1a2a3...an)−1，那么an的逆元inva[n]=Sn−1⋅invs[n]=a−1n，而且invs[n−1]=invs[n]⋅an，这样就可以递推出数列元素的逆元以及元素前缀积的逆元了。

7. 总结：预处理筛出μ，求出F以及其前缀积，然后求逆元；枚举因数求出g；然后求前缀积，求逆元；最后就可以O(TN−−√logN)过掉了。虽然步骤多，但是代码不长。

上代码

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long LL;
const int N = 1e6 + 10;
const int MOD = 1e9 + 7;

int n, m;
inline int read() {
register char ch = getchar();
register int ans = 0, neg = 1;
for (; !isdigit(ch); ch = getchar())
if (ch == '-') neg = -1;
for (; isdigit(ch); ch = getchar())
ans = ans * 10 + ch - '0';
return ans * neg;
}

inline int pls(int a, int b) {
LL ans = (LL)a + b;
return ans >= (LL)MOD ? ans - MOD : ans;
}
inline int mul(int a, int b) {
return (LL)a * b % MOD;
}
inline int power(int a, int b) {
int ans = 1;
while (b) {
if (b & 1) ans = mul(ans, a);
b >>= 1, a = mul(a, a);
}
return ans;
}

bool notPrime
;
int prime
, mu
;
int fi
, pf
, invfi
, invpf
;
int gi
, pg
, invpg
;
void preCalc(int a) {
int cnt = 0, mut = 1;
mu[1] = 1, fi[0] = 0, fi[1] = pf[1] = 1;
for (int i = 2; i <= a; ++i) {
if (!notPrime[i]) prime[++cnt] = i, mu[i] = -1;
for (int j = 1, k; j <= cnt && (k = i * prime[j]) <= a; ++j) {
notPrime[k] = true;
if (i % prime[j] == 0) {
mu[k] = 0; break;
}
mu[k] = -mu[i];
}
pf[i] = mul(pf[i - 1], fi[i] = pls(fi[i - 1], fi[i - 2]));
}
invpf[a] = power(pf[a], MOD - 2);
for (int i = a; i >= 1; --i) {
gi[i] = 1;
invpf[i - 1] = mul(invpf[i], fi[i]);
invfi[i] = mul(pf[i - 1], invpf[i]);
}
for (int i = 2; i <= a; ++i)
for (int j = 1, k; (k = i * j) <= a; ++j)
if (mu[j] == 1) gi[k] = mul(gi[k], fi[i]);
else if (mu[j] == -1) gi[k] = mul(gi[k], invfi[i]);
pg[1] = 1;
for (int i = 2; i <= a; ++i)
pg[i] = mul(pg[i - 1], gi[i]);
invpg[a] = power(pg[a], MOD - 2);
for (int i = a; i >= 1; --i)
invpg[i - 1] = mul(invpg[i], gi[i]);
}

int main() {
int T = read();
preCalc(1e6);
while (T--) {
n = read(), m = read();
if (n > m) swap(n, m);
int ans = 1;
for (int i = 1, nxt; i <= n; i = ++nxt) {
nxt = min(n / (n / i), m / (m / i));
ans = mul(ans, power(mul(pg[nxt], invpg[i - 1]), (LL)(n / i) * (m / i) % (MOD - 1)));
}
printf("%d\n", ans);
}
return 0;
}

以上