洛谷P3704 [SDOI2017] 数字表格题解

题目链接：P3704 [SDOI2017] 数字表格

题意：

$q$ 组询问，给定 $n,m$ ，求
$\mathrm{Ans} = \prod_{i=1}^n\prod_{j=1}^{m}F_{\gcd(i,j)}$
其中 $F$ 的幂次可以快速计算（在本题中是斐波那契数列，$F_0=0,F_1=1$ ）

$1\le q \le 10^3,1 \le n,m \le 10^6$

莫比乌斯反演 + 数论分块。

既然 $F$ 可以快速计算，那就是叫我们统计 $F(k)$ 是怎么贡献的了。

$\prod_{i=1}^n\prod_{j=1}^{m}F_{\gcd(i,j)} = \prod_{k=1}^nF_k^{\sum_{i=1}^n \sum_{j=1}^m[\operatorname{gcd}(i, j)=k]}$

把上面那一堆提出来

$\begin{aligned} \sum_{i=1}^n \sum_{j=1}^m[\operatorname{gcd}(i, j)=k] &= \sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[\operatorname{gcd}(i, j)=1] \\[6pt]&= \sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum_{d \mid \operatorname{gcd}(i, j)} \mu(d) \\[6pt]&= \sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum_{d=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \mu(d)\cdot [d \mid \operatorname{gcd}(i, j)] \\[6pt]&= \sum_{d=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \mu(d)\sum_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}[d \mid \operatorname{gcd}(i, j)] \end{aligned}$

因为能产生贡献的 $d$ 一定为 $\gcd(i,j)$ 的因数，反过来能产生贡献的 $i,j$ 一定是 $d$ 的倍数。

所以我们可以转为枚举 $d$ 的倍数

$\begin{aligned} &= \sum_{d=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \mu(d)\sum_{i=1}^{\left\lfloor\frac{n}{kd}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{kd}\right\rfloor}[d \mid \operatorname{gcd}(di, dj)] \\[6pt] &= \sum_{d=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \mu(d)\sum_{i=1}^{\left\lfloor\frac{n}{kd}\right\rfloor} \sum_{j=1}^{\left\lfloor\frac{m}{kd}\right\rfloor}1 \\[6pt] &= \sum_{d=1}^{\left\lfloor\frac{n}{k}\right\rfloor} \mu(d)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor \end{aligned}$

这里需要用到洛谷P2257 YY的GCD 中的技巧，把 $kd$ 看作一个整体枚举

在本题中我们直接把上面那一坨东西扔回去再枚举 $T=kd$

$\begin{aligned} \mathrm{Ans}& =\prod_{k=1}^N F_k\left(\sum_{d=1}^{\left\lfloor\frac{N}{k}\right\rfloor} \mu(d)\left\lfloor\frac{N}{k d}\right\rfloor\left\lfloor\frac{M}{k d}\right\rfloor\right) \\ & =\prod_{T=1}^N\left(\prod_{k \mid T} F_k{ }^{\mu\left(\frac{T}{k}\right)}\right)^{\left\lfloor\frac{N}{T}\right\rfloor\left\lfloor\frac{M}{T}\right\rfloor} \end{aligned}$

令 $f(n) = \prod_{d \mid n} F_d^{\mu\left(\frac{n}{d}\right)}$ ，则

$\mathrm{Ans} = \prod_{T=1}^N f(T)^{\left\lfloor\frac{N}{T}\right\rfloor\left\lfloor\frac{M}{T}\right\rfloor}$

其中 $f(T)$ 可以暴力预处理出来，每个数直接往它的倍数贡献，复杂度是 $\mathcal{O}(n \log n)$

注意数论分块的 $f(T)$ 的乘积要预处理，不要直接用快速幂算，否则会多一个 log 。

总时间复杂度 $\mathcal{O}(n \log nM + q\sqrt{N}\log M)$ ，其中 $M$ 是模数。

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
const int mod = 1e9 + 7;
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(1e6 + 15))

char ck[N];
int pcnt, prime[N], mu[N], f[N], fr[N];
int qpow(int a, int b)
{
    int r = 1;
    while(b) {
        if(b & 1) r = r * a % mod;
        a = a * a % mod; b >>= 1;
    }
    return r;
}
int inv(int x) { return qpow(x, mod - 2); }
void init()
{
    ck[1] = true; mu[1] = 1;
    for(int i = 2; i <= N - 5; i++)
    {
        if(!ck[i]) { prime[++pcnt] = i; mu[i] = -1; }
        for(int j = 1; j <= pcnt && i * prime[j] <= N - 5; j++)
        {
            int pos = prime[j] * i; ck[pos] = true;
            if(i % prime[j]) { mu[pos] = -mu[i]; } else break;
        }
    }
    for(int i = 1; i <= N - 5; i++) { f[i] = 1, fr[i] = 1; }
    int A = 1, B = 0;
    for(int i = 1; i <= N - 5; i++)
    {
        B = (A + B) % mod; A = (B - A + mod) % mod;
        int g[3] = {inv(B), 1, B}; 
        for(int j = i, k = 1; j <= N - 5; j += i, ++k)
        {
            f[j] = f[j] * g[1 + mu[k]] % mod;
            fr[j] = fr[j] * g[1 - mu[k]] % mod;
        }
    }
    f[0] = fr[0] = 1;
    for(int i = 1; i <= N - 5; i++)
    {
        f[i] = f[i - 1] * f[i] % mod;
        fr[i] = fr[i - 1] * fr[i] % mod;
    }
}
void solve()
{
    int n, m; cin >> n >> m; if(n > m) swap(n, m);
    int res = 1;
    for(int l = 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        res = res * qpow(f[r] * fr[l - 1] % mod, (n / l) * (m / l)) % mod;
    }
    cout << res << '\n';
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    init(); int qwq; cin >> qwq; while(qwq--) solve();
    return 0;
}