BZOJ4816 [Sdoi2017]数字表格 【莫比乌斯反演】

题目

Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项，那么
f[0]=0
f[1]=1
f
=f[n-1]+f[n-2],n>=2
Doris用老师的超级计算机生成了一个n×m的表格，第i行第j列的格子中的数是f[gcd(i,j)]，其中gcd(i,j)表示i,
j的最大公约数。Doris的表格中共有n×m个数，她想知道这些数的乘积是多少。答案对10^9+7取模。

输入格式

有多组测试数据。

第一个一个数T，表示数据组数。
接下来T行，每行两个数n,m
T<=1000,1<=n,m<=10^6

输出格式

输出T行，第i行的数是第i组数据的结果

输入样例

3

2 3

4 5

6 7

输出样例

1

6

960

题解

一道满是套路的莫比乌斯反演题

我们要求：
\[\prod\limits_{i=1}^{N}\prod\limits_{j=1}^{M} f[gcd(i,j)]\]

根据莫比乌斯反演的套路，我们将gcd改为枚举d
\[\prod\limits_{d=1}^{N} f[d]^{\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{M} 1 [gcd(i,j) == d]}\]

然后把指数中的\(d\)提掉
\[\prod\limits_{d=1}^{N} f[d]^{\sum\limits_{i=1}^{\lfloor \frac{N}{d} \rfloor}\sum\limits_{j=1}^{\lfloor \frac{M}{d} \rfloor} 1 [gcd(i,j) == 1]}\]

然后指数部分就是经典的莫比乌斯反演
\[\prod\limits_{d=1}^{N} f[d]^{\sum\limits_{i=1}^{\lfloor \frac{N}{d} \rfloor} \mu(i) * \lfloor \frac{N}{i * d} \rfloor * \lfloor \frac{M}{i * d} \rfloor}\]

如果只有一组询问，这样接近\(O(n)\)可以过，但是多组询问，我们考虑继续优化
我们有一个枚举\(i * d\)的套路
我们记\(T = i * d\)
有
\[\prod\limits_{T=1}^{N} \prod\limits_{d|T} f[d]^{\mu(\frac{T}{d}) * \lfloor \frac{N}{T} \rfloor * \lfloor \frac{M}{T} \rfloor}\]

划分一下：
\[\prod\limits_{T=1}^{N} ( \prod\limits_{d|T} f[d]^{\mu(\frac{T}{d})} ) ^ {\lfloor \frac{N}{T} \rfloor * \lfloor \frac{M}{T} \rfloor}\]

容易发现，内层的东西就是关于\(T\)的因子的积，根据经验，这玩意可以通过枚举\(O(nlogn)\)预处理
外层是只与\(T\)有关的整除，可以\(O(\sqrt{N})\)分块

那我们就用\(O(nlogn + T\sqrt{N})\)的复杂度做完了

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
#define res register
using namespace std;
const int maxn = 1000005,maxm = 100005,INF = 1000000000,P = 1000000007,md = 1000000006;
int fv[maxn],f[maxn],g[maxn],gv[maxn];
int mu[maxn],p[maxn],fac[maxn],pi;
int isn[maxn];
int qpow(int a,LL b){
b = (b % md + md) % md;
int ans = 1;
for (; b; b >>= 1,a = (LL)a * a % P)
if (b & 1) ans = (LL)ans * a % P;
return ans;
}
void init(){
mu[1] = 1;
for (res int i = 2; i < maxn; i++){
if (!isn[i]) p[++pi] = i,mu[i] = -1;
for (res int j = 1; j <= pi && i * p[j] < maxn; j++){
isn[i * p[j]] = true;
if (i % p[j] == 0){
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
f[0] = 0; f[1] = 1;
fv[0] = 1; fv[1] = 1;
for (res int i = 2; i < maxn; i++){
f[i] = (f[i - 1] + f[i - 2]) % P;
fv[i] = qpow(f[i],P - 2);
}
for (res int i = 1; i < maxn; i++) g[i] = 1;
for (res int i = 1; i < maxn; i++){
for (int j = i; j < maxn; j += i){
if (mu[j / i] == 1) g[j] = (LL)g[j] * f[i] % P;
if (mu[j / i] == -1) g[j] = (LL)g[j] * fv[i] % P;
}
}
g[0] = gv[0] = 1;
gv[1] = qpow(g[1],P - 2);
for (res int i = 2; i < maxn; i++){
g[i] = (LL)g[i] * g[i - 1] % P;
gv[i] = qpow(g[i],P - 2);
}
}
int main(){
init();
int n,m,T;
LL ans;
cin >> T;
while (T--){
cin >> n >> m;
ans = 1;
if (n > m) swap(n,m);
for (res int i = 1,nxt; i <= n; i = nxt + 1){
nxt = min(n / (n / i),m / (m / i));
ans = (LL)ans * qpow((LL)g[nxt] * gv[i - 1] % P,(LL)(n / i) * (m / i)) % P;
}
cout << ans << endl;
}
return 0;
}