BZOJ4816【SDOI2017】数字表格 <莫比乌斯反演>

Problem

【SDOI2017】数字表格

$\mathrm{Time\;Limit:\;50\;Sec}$
$\mathrm{Memory\;Limit:\;128\;MB}$

Description

$\mathrm{Doris}$刚刚学习了$\mathrm{fibonacci}$数列。用$f[i]$表示数列的第$i$项，那么

$$f[n]= \begin{cases} 0&n=0\\ 1&n=1\\ f[n-1]+f[n-2]&n\ge2 \end{cases}$$

$\mathrm{Doris}$用老师的超级计算机生成了一个$n\times m$的表格，第$i$行第$j$列的格子中的数是$f[\gcd(i,j)]$，其中$\gcd(i,j)$表示$i,j$的最大公约数。
$\mathrm{Doris}$的表格中共有$n\times m$个数，她想知道这些数的乘积是多少。答案对$10^9+7$取模。

Input

有多组测试数据。
第一个一个数$T$，表示数据组数。
接下来$T$行，每行两个数$n,m$。

Output

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

Sample Input

3
2 3
4 5
6 7

Sample Output

1
6
960

HINT

$T\le1000,\;1\le n,m\le10^6$

Source

鸣谢infinityedge上传

标签：莫比乌斯反演

Solution

转换题目求和的角度为枚举$\gcd(i,j)$，求$f[\gcd(i,j)]$对答案的贡献。
那么有

$$\begin{aligned} Ans&=\prod_{d=1}^{n}f(d)^{\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=d]}\\ &=\prod_{d=1}^{n}f(d)^{\sum_{t=1}^{\lfloor\frac{\min(n,m)}{d}\rfloor}\mu(t)\lfloor\frac{n}{dt}\rfloor\lfloor\frac{m}{dt}\rfloor}\\ &=\prod_{T=1}^{n}\prod_{d|T}f(d)^{\mu(\frac{T}{d})\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}\\ &=\prod_{T=1}^{n}(\prod_{d|T}f(d)^{\mu(\frac{T}{d})})^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}\\ \end{aligned}$$

将中间$\prod_{d|T}f(d)^{\mu(\frac{T}{d})}$单独分开，设$g(T)=\prod_{d|T}f(d)^{\mu(\frac{T}{d})}$，那么$Ans=\prod_{T=1}^{n}g(T)^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}$，预处理出$\lbrace g\rbrace$的前缀积后数论分块即可。

发现对于每个$d\in[1,n]$，$f(d)^{\mu(\frac{T}{d})}$最多只会对$\frac{n}{d}$个$g(T)$的值产生贡献，枚举$d$累加贡献的时间复杂度是调和级数。于是枚举$d$，枚举$d$在$[1,n]$内的倍数$t$，将$f(d)^{\mu(\frac{t}{d})}$乘到$g(t)$中即可处理出所有$g$。而$\mu$只能取$\pm1$，于是需要预处理出$f(d)$和$f(d)^{-1}$即$f(d)$在模意义下的逆元。处理出$g$后再处理$g$的前缀积即可。

时间复杂度$O(n\log{n}+T(\sqrt{n}+\sqrt{m})\log\mathrm{MOD})$。

Code

#include <bits/stdc++.h>
#define MAX_N 1000000
#define MOD 1000000007
using namespace std;
typedef long long lnt;
template <class T> inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
bool NotPri[MAX_N+5]; int cnt, pri[MAX_N+5]; lnt ans, gp;
lnt mu[MAX_N+5], f[MAX_N+5], g[MAX_N+5], inv[MAX_N+5];
lnt Pow(lnt x, lnt k) {
    lnt ret = 1LL;
    for (; k; k >>= 1, (x *= x) %= MOD)
        if (k&1) (ret *= x) %= MOD;
    return ret;
}
void init() {
    mu[1] = f[1] = inv[1] = 1;
    for (int i = 0; i <= MAX_N; i++) g[i] = 1;
    for (int i = 2; i <= MAX_N; i++) {
        if (!NotPri[i]) pri[cnt++] = i, mu[i] = -1;
        for (int j = 0; j < cnt; j++) {
            if (i*pri[j] > MAX_N) break;
            NotPri[i*pri[j]] = true;
            if (i%pri[j]) mu[i*pri[j]] = -mu[i];
            else {mu[i*pri[j]] = 0; break;}
        }
    }
    for (int i = 2; i <= MAX_N; i++)
        f[i] = (f[i-1]+f[i-2])%MOD, inv[i] = Pow(f[i], MOD-2);
    for (int i = 1; i <= MAX_N; i++) if (mu[i])
        for (int j = i; j <= MAX_N; j += i)
            if (mu[i] == 1) (g[j] *= f[j/i]) %= MOD;
            else (g[j] *= inv[j/i]) %= MOD;
    for (int i = 2; i <= MAX_N; i++)
        (g[i] *= g[i-1]) %= MOD;
}
int main() {
    int T; read(T), init();
    while (T--) {
        int n, m; read(n), read(m), ans = 1LL;
        for (int l = 1, r; l <= min(n, m); l = r+1)
            r = min(n/(n/l), m/(m/l)), 
            gp = g[r]*Pow(g[l-1], MOD-2)%MOD, 
            gp = Pow(gp, 1LL*(n/l)*(m/l)%(MOD-1)), 
            (ans *= gp) %= MOD;
        printf("%lld\n", ans);
    }
    return 0;
}