[BZOJ4816][SDOI2017]数字表格（莫比乌斯反演）

最新推荐文章于 2019-11-06 20:09:08 发布

xyz32768

最新推荐文章于 2019-11-06 20:09:08 发布

阅读量317

点赞数

CC 4.0 BY-SA版权

分类专栏： BZOJ、UOJ、LOJ 文章标签： gcd 莫比乌斯反演

本文链接：https://blog.youkuaiyun.com/xyz32768/article/details/78289077

BZOJ、UOJ、LOJ 专栏收录该内容

341 篇文章

订阅专栏

本文介绍了一种利用莫比乌斯反演解决特定数学问题的方法。通过转换思考角度，将问题转化为枚举最大公约数，进而推导出求解公式。文中详细解释了从原始表达式到最终形式的转换过程，并提供了完整的代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

容易看出，此题用莫比乌斯反演求解。
首先，要把 $\prod_{i=1}^n\prod_{j=1}^mf[\gcd(i,j)]$ 换一个方向去思考，即不枚举 $i,j$ ，而是枚举 $d=\gcd(i,j)$ 。可以得到，在这个表格里， $f[d]$ 出现的次数就是 $\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]$ 。
所以：
$Ans=\prod_df[d]^{\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]}$
$=\prod_df[d]^{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]}$
而 $\sum_{i=1}^x\sum_{j=1}^y[\gcd(i,j)=1]$ 是一个经典问题，它等于 $\sum_d\lfloor\frac{x}{d}\rfloor\lfloor\frac{y}{d}\rfloor\mu(d)$ 。
即原式等于 $\prod_df[d]^{\sum_k\lfloor\frac{n}{dk}\rfloor\lfloor\frac{m}{dk}\rfloor\mu(k)}$ 。
继续考虑。可以发现， $\lfloor\frac{n}{dk}\rfloor$ 和 $\lfloor\frac{m}{dk}\rfloor$ 这两个式子不容易直接分块。所以这里令 $u=dk$ ，原式化为：
$\prod_u\prod_{k|u}f[\frac{u}{k}]^{\mu(k)\lfloor\frac{n}{u}\rfloor\lfloor\frac{m}{u}\rfloor}$
$=\prod_u(\prod_{k|u}f[\frac{u}{k}]^{\mu(k)})^{\lfloor\frac{n}{u}\rfloor\lfloor\frac{m}{u}\rfloor}$ 。
这样就可以将 $\lfloor\frac{n}{u}\rfloor$ 和 $\lfloor\frac{m}{u}\rfloor$ 的值分块了。令 $F[u]=\prod_{k|u}f[\frac{u}{k}]^{\mu(k)}$ ，那么就可以预处理的前缀积，计算区间的积时使用逆元计算。
代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
const int MaxN = 1e6, PYZ = 1e9 + 7, N = MaxN + 5;
int tot, f[N], g[N], pri[N], miu[N], F[N], prod[N];
bool mark[N];
int qpow(int a, int b) {
    int res = 1;
    while (b) {
        if (b & 1) res = 1ll * res * a % PYZ;
        a = 1ll * a * a % PYZ;
        b >>= 1;
    }
    return res;
}
void sieve() {
    int i, j; f[0] = 0; f[1] = miu[1] = g[1] = F[1] = prod[0] = 1;
    mark[0] = mark[1] = 1;
    for (i = 2; i <= MaxN; i++) {
        f[i] = (f[i - 1] + f[i - 2]) % PYZ;
        g[i] = qpow(f[i], PYZ - 2); F[i] = 1;
        if (!mark[i]) pri[++tot] = i, miu[i] = -1;
        for (j = 1; j <= tot; j++) {
            if (1ll * i * pri[j] > MaxN) break;
            mark[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
            else miu[i * pri[j]] = -miu[i];
        }
    }
    for (i = 1; i <= MaxN; i++) {
        if (miu[i] != 0) for (j = i; j <= MaxN; j += i)
            F[j] = 1ll * F[j] * (miu[i] == 1 ? f[j / i] : g[j / i]) % PYZ;
        prod[i] = 1ll * prod[i - 1] * F[i] % PYZ;
    }
}
int solve(int a, int b) {
    int i, n = min(a, b), ans = 1;
    for (i = 1; i <= n;) {
        int nxt = min(a / (a / i), b / (b / i)), pro;
        pro = 1ll * prod[nxt] * qpow(prod[i - 1], PYZ - 2) % PYZ;
        ans = 1ll * ans * qpow(pro, 1ll * (a / i) * (b / i) % (PYZ - 1)) % PYZ;
        i = nxt + 1;
    }
    return ans;
}
int main() {
    int a, b, T = read(); sieve();
    while (T--) a = read(), b = read(),
        printf("%d\n", solve(a, b));
    return 0;
}