bzoj4816: [Sdoi2017]数字表格 莫比乌斯反演

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

分析

推导一波得了
设$n<m$

$$ans=\prod_d^n F[d]^{\sum_i^n\sum_j^m[gcd(i,j)==d]}$$
$$=\prod_d^n F[d]^{\sum_i^n\sum_j^m\sum_{k|i,k|j}^n\mu(k)}$$
$$=\prod_d^n F[d]^{\sum_k^n\mu(k)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor}$$
令 $D=kd$
$$=\prod_D^n \prod_{d|D}^n F[d]^{\mu(\frac{D}{d})\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor}$$
令 $g[D]=\prod_{d|D}^n F[d]^{\mu(\frac{D}{d})}$ 暴力预处理
$$ans=\prod_D^n g[D]^{\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor}$$
上一个下底函数分块可以做到 $O(\sqrt n)$
总复杂度 $O(nlogn+T\sqrt n)$

代码

#include<cstdio>
#include<algorithm>
const int N = 1e6 + 10, P = 1e9 + 7;
int f[N], mu[N], g[N], pg[N], ng[N], nf[N], pr[N], tp; bool v[N];
int Pow(int x, int k) {int r = 1; for(;k; x = 1LL * x * x % P, k >>= 1) if(k & 1) r = 1LL * r * x % P; return r;}
void Pre(int N) {
    f[0] = 0; mu[1] = f[1] = 1; 
    for(int i = 2;i <= N; ++i) f[i] = (f[i - 1] + f[i - 2]) % P;
    for(int i = 1;i <= N; ++i) nf[i] = Pow(f[i], P - 2);
    for(int i = 2;i <= N; ++i) {
        if(!v[i]) {pr[++tp] = i; mu[i] = -1;}
        for(int j = 1;j <= tp && i * pr[j] <= N; ++j) {
            v[i * pr[j]] = true;
            if(i % pr[j]) mu[i * pr[j]] = -mu[i];
            else break;
        }
    }
    for(int i = 0;i <= N; ++i) g[i] = 1;
    for(int d = 1;d <= N; ++d)
        for(int j = d, x = 1;j <= N; j += d, ++x) if(mu[x])
            g[j] = 1LL * g[j] * (~mu[x] ? f[d] : nf[d]) % P;
    pg[0] = 1; for(int i = 1;i <= N; ++i) pg[i] = 1LL * pg[i - 1] * g[i] % P; 
    ng[N] = Pow(pg[N], P - 2); for(int i = N; ~i; --i) ng[i - 1] = 1LL * ng[i] * g[i] % P;
}
int main() {
    Pre(1e6); int T, n, m;
    for(scanf("%d", &T); T--;) {
        scanf("%d%d", &n, &m); int ans = 1;
        for(int i = 1, p; i <= std::min(n, m); i = p + 1) {
            p = std::min(m / (m / i), n / (n / i));
            ans = (1LL * ans *  Pow(1LL * pg[p] * ng[i - 1] % P, 1LL * (n / i) * (m / i) % (P - 1))) % P;
        }
        printf("%d\n", ans);
    }
    return 0;
}