【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

3

2 3

4 5

6 7

Sample Output

1

6

960

题解

莫比乌斯反演。
求

$$\prod_{x=1}^n\prod_{y=1}^m fib(gcd(x,y))$$
如果我们能找到一个 $f(x)$使 $fib(x) = \prod_{d|x}f(d)$。
那么原式可以转化为： $$\prod_{x=1}^n\prod_{y=1}^m \prod_{d|gcd(x,y)}f(d)$$
$$\Rightarrow \prod_{x=1}^n\prod_{y=1}^m \prod_{d|x,d|y}f(d)$$
$$\Rightarrow \prod_{d=1}^{min(n,m)}f(d)^{\left \lfloor\frac nd \right \rfloor\left \lfloor\frac md \right \rfloor}$$
然后相同的 $\left \lfloor\frac nd \right \rfloor\left \lfloor\frac md \right \rfloor$分块计算即可。
考虑 $f(x)$怎么求。
显然 $$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;
}