【bzoj4816】[Sdoi2017]数字表格

题目链接

Description

Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项，那么
f[0]=0
f[1]=1
f[n]=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取模。

Input

有多组测试数据。

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

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

Sample Input

2 3

4 5

6 7

Sample Output

960

题解

莫比乌斯反演。
求

∏x=1n∏y=1mfib(gcd(x,y))

\prod_{x=1}^n\prod_{y=1}^m fib(gcd(x,y))
如果我们能找到一个 f(x) f(x)使 fib(x)=∏d|xf(d) fib(x) = \prod_{d|x}f(d)。
那么原式可以转化为：

∏x=1n∏y=1m∏d|gcd(x,y)f(d)

\prod_{x=1}^n\prod_{y=1}^m \prod_{d|gcd(x,y)}f(d)

⇒∏x=1n∏y=1m∏d|x,d|yf(d)

\Rightarrow \prod_{x=1}^n\prod_{y=1}^m \prod_{d|x,d|y}f(d)

⇒∏d=1min(n,m)f(d)⌊nd⌋⌊md⌋

\Rightarrow \prod_{d=1}^{min(n,m)}f(d)^{\left \lfloor\frac nd \right \rfloor\left \lfloor\frac md \right \rfloor}
然后相同的 ⌊nd⌋⌊md⌋ \left \lfloor\frac nd \right \rfloor\left \lfloor\frac md \right \rfloor分块计算即可。
考虑 f(x) f(x)怎么求。
显然

f(x)=fib(x)∏d|x,d≠xf(d)

f(x)=\frac {fib(x)}{\prod_{d|x,d\not= x}f(d)}
那么就可以递推求解了。

#include<bits/stdc++.h>
using namespace std;inline int read(){int x = 0, f = 1; char c = getchar();while(!isdigit(c)) { if(c == '-') f = -1; c = getchar(); }while(isdigit(c)) { x = x * 10 + c - '0'; c = getchar(); }return x * f;
}typedef long long ll;
const int N = 1000000 + 10, mo = 1000000007;
ll f[N], fv[N];
int n, m;ll poww(ll a, ll b){ll cnt = 1;while(b){if(b & 1) cnt  = (cnt * a) % mo;a = (a * a) % mo;b >>= 1;}return cnt;
}void init(){f[1] = 1;for(int i = 2; i < N; i++)f[i] = (f[i-1] + f[i-2]) % mo;for(int i = 1; i < N; i++){ll inv = poww(f[i], mo - 2);for(int j = i + i; j < N; j += i)f[j] = (f[j] * inv) % mo;}f[0] = 1;for(int i = 1; i < N; i++) f[i] = (f[i] * f[i-1]) % mo;for(int i = 0; i < N; i++) fv[i] = poww(f[i], mo - 2);
}void work(){int t = read();ll ans;while(t--){int n = read(), m = read();if(n > m) swap(n, m);ans = 1;for(int i = 1, p; i <= n; i = p + 1){p = min(n/(n/i), m/(m/i));ans = ans * poww(poww((f[p] * fv[i-1] % mo), (ll)n/i), (ll)m/i) % mo;}printf("%lld\n", ans);}
}int main(){init();work();return 0;
}