BZOJ4816 数字表格

4816: [Sdoi2017]数字表格

Time Limit: 50 Sec  Memory Limit: 128 MB

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

---

题解

这道题很好地延续了SDOI的优良传统，考了一道莫比乌斯反演以供娱乐。

由于我们一眼发现了这是一道莫比乌斯反演水题，正如做所有的莫比乌斯反演一样，我们把要求的式子先写出来并推导

\begin{eqnarray*}
ans & = & \prod_{i}^{n}\prod _{j}^{m}f( \gcd(i,j) ) \\
& = & \prod_{k}^{n}f(k) ^ {\sum_{i}^{\frac{n}{k}}\sum_{j}^{\frac{m}{k}}[\gcd(i,j)=1]} \\
& = & \prod_{k}^{n}f(k) ^ {\sum_{i}^{\frac{n}{k}}\sum_{j}^{\frac{m}{k}}\sum_{x\mid{i}\&x\mid{j}}~~\mu(x)} \\
& = & \prod_{k}^{n}\prod_{x}^{\frac{n}{k}}(f(k) ^ {\mu(x)})^{\frac{n}{kx}\frac{m}{kx}}
\end{eqnarray*}

为了能够把\(f(k) ^ {\mu(x)}\)提出来，显然，我们可以设\(T=kx\)，\(g(T)=\prod_{k\mid{T}}f(k) ^ {\mu(\frac{T}{k})}\)

化简得到\(ans = \prod_{T}^{n}g(T)^{\frac{n}{T}\frac{m}{T}}\)

求出$f$以及$f$的逆元，线性筛求$\mu$，Dirichlet卷积求出$g$，然后计算$g$的前缀积$g'$以及$g'$的逆元

查询使用大众喜闻乐见的分块，至此，我们切掉了这道水题。

代码

#include<bits/stdc++.h>
using namespace std;
template <class _T> inline void read(_T &_x) {
    int _t; bool flag = false;
    while ((_t = getchar()) != '-' && (_t < '0' || _t > '9')) ;
    if (_t == '-') _t = getchar(), flag = true; _x = _t - '0';
    while ((_t = getchar()) >= '0' && _t <= '9') _x = _x * 10 + _t - '0';
    if (flag) _x = -_x;
}
typedef long long LL;
const int maxn = 1000010;
const int mod = 1000000007;
int f[maxn], f_re[maxn], mu[maxn], g[maxn], g_re[maxn];
bool vis[maxn]; int prime[maxn / 10], pcnt;
#define trans(x) ((int)((LL)x % mod))
#define reg register
inline void update(int &a, reg int b) {
    a = trans(a * b);
    if (a < 0) a += mod;
}
inline int qpower(int a, reg LL b) {
    int ret = 1;
    while (b) {
        if (b & 1) update(ret, a);
        update(a, a), b >>= 1;
    }
    return ret;
}
inline int calc(reg int a, reg int b) {
    if (b == 0) return 1;
    return b < 0 ? f_re[a] : f[a];
}
void init() {
    f[0] = 0, f[1] = f_re[1] = g[1] = mu[1] = 1;
    reg int i, j;
    for (i = 2; i < maxn; ++i) {
        f[i] = f[i - 1] + f[i - 2];
        if (f[i] >= mod) f[i] -= mod;
        f_re[i] = qpower(f[i], mod - 2);
    }
    for (i = 2; i < maxn; ++i) {
        g[i] = 1;
        if (!vis[i]) {
            prime[++pcnt] = i;
            mu[i] = -1;
        }
        for (j = 1; j <= pcnt && prime[j] * i < maxn; ++j) {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    for (i = 1; i * i < maxn; ++i) {
        update(g[i * i], calc(i, mu[i]));
        for (j = i + 1; i * j < maxn; ++j)
            update(g[i * j], trans(calc(i, mu[j]) * calc(j, mu[i])));
    }
    g[0] = g_re[0] = 1;
    for (i = 1; i < maxn; ++i) {
        update(g[i], g[i - 1]);
        g_re[i] = qpower(g[i], mod - 2);
    }
}
inline int query(reg int a, reg int b) {
    if (a > b) {int t = a; a = b, b = t; }
    int ret = 1;
    for (reg int i = 1, j, x, y; i <= a; i = j + 1) {
        x = a / i, y = b / i;
        j = min(a / x, b / y);
        update(ret, qpower(trans(g[j] * g_re[i - 1]), (LL)x * y));
    }
    return ret;
}
int main() {
    //freopen(".in", "r", stdin);
    //freopen(".out", "w", stdout);
    init();
    int T, a, b; read(T);
    while (T--) {
        read(a), read(b);
        printf("%d\n", query(a, b));
    }
    return 0;
}