【BZOJ 4816】【SDOI 2017】数字表格

考虑到$gcd(i,j)$的这个形式是一个常见的莫比乌斯反演，尝试构造函数。

1、枚举gcd，转化为$ans = \prod_{d=1}^{n}f(d)^{h(d)}$，其中$h(d)=\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)==d]$；
2、$h(d)$是莫比乌斯反演的模板题，直接上结论：$h(d)=\sum _{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(k)\left \lfloor \frac{n}{dk} \right \rfloor\left \lfloor \frac{m}{dk} \right \rfloor$，所以转化为$ans = \prod_{d=1}^{n}f(d)^{h(d)=\sum _{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(k)\left \lfloor \frac{n}{dk} \right \rfloor\left \lfloor \frac{m}{dk} \right \rfloor}$；
做到这里通过枚举约数的方法就能够得到答案，这个方法的时间复杂度是$O((n+m)^{\frac{3}{4}})$，但是乘上一个$T$就有可能超时。这样至少就是60分。
3、设$T=dk$，并且改变一下运算顺序（因为全部都是乘法），转换为$ans=\prod _{T=1}^{n}(\prod _{d|T}f(d)^{\mu(\frac{T}{d})})^{\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor}$；
4、记$g(x)=\prod _{d|T}f(d)^{\mu(\frac{T}{d})}$，现在如果我们已经知道了$g(x)$的前缀积，我们就可以在$O(\sqrt {N}logN)$的复杂度内求出每次询问，其中根号是分块的近似复杂度，$log$是快速幂的复杂度；
5、那么如何求$g(x)$呢？可以通过枚举d，对于每个满足$d|T$的$g(T)$都乘上$f(d)^{\mu(\frac{T}{d})}$，只需提前预处理$\mu(x)$、$f(x)$和$f(x)$的逆元即可，根据调和级数的相关证明可以得到这部分的时间复杂度是$O(lnNlogN)$；
6、那么如何求$f(x)$的逆元呢？显然不能$O(P)$求出所有数的逆元。可以选择扩欧来做，时间复杂度为$O(NlnN)$，不过网上有$O(N)$的做法：设$S_{i}=\prod_{1}^{i}a_i$，然后求出$S_n$的逆元$invS_n$，然后通过递推式$invS_{i-1}=a_{i}*invS_{i}$求出$invS$数组，再通过$inva_{i-1}=invS_{i}*S_{i-1}$就可以求出所有的逆元；
7、至此所有的问题都已经解决了，注意最后快速幂求答案的时候指数模一下mod-2不然会有神奇的错误。。。

#include<cmath>
#include<cstdio>
#include<vector>
#include<queue>
#include<cstring>
#include<iomanip>
#include<stdlib.h>
#include<iostream>
#include<algorithm>
#define ll long long
#define inf 1000000000
#define mo 1000000007
#define N 1000000
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
int mbs[N+5],p[N],vis[N+5];
int f[N+5],Sf[N+5],invf[N+5],invS[N+5];
int g[N+5],Sg[N+5],invg[N+5];
int T,n,m,res,i,last;
int quick_multi(int x,int a)
{
    int res = 1,t = x;
    while (a) {if (a&1) res=(1ll*res*t)%mo; t=(1ll*t*t)%mo; a>>=1;}
    return res;
}
void mobius()
{
    int i,j; mbs[1] = 1;
    fo(i,2,N)
        {
            if (!vis[i]) {p[++p[0]] = i; mbs[i] = -1;}
            for (j = 1;j <= p[0] && i * p[j] <= N; j++)
                {
                    vis[i*p[j]] = 1;
                    if (i % p[j] == 0) {mbs[i*p[j]] = 0; break;}
                    mbs[i*p[j]] = - mbs[i];
                }
        }
}
void fibonacci()
{
    int i,j; f[0] = 0; f[1] = 1; Sf[1] = 1;
    fo(i,2,N)
        {
            f[i] = f[i-1] + f[i-2];
            f[i] -= f[i]>mo?mo:0;
            Sf[i] = 1ll*Sf[i-1]*f[i]%mo; 
        }
    invS[N] = quick_multi(Sf[N],mo-2); invf[1] = 1;
    fd(i,N,2) invS[i-1] = 1ll*invS[i]*f[i] % mo;
    fd(i,N,2) invf[i] = 1ll*Sf[i-1]*invS[i] % mo;
    invf[0] = 1;
}
void pre_g()
{
    int i,d;
    fo(i,1,N) g[i] = 1;
    fo(d,1,N)
        {
            for (i = 1;i * d <= N; i++)
                if (mbs[i] == 1) g[i*d] = 1ll*g[i*d]*f[d]%mo; else
                if (mbs[i] == -1) g[i*d] = 1ll*g[i*d]*invf[d]%mo;
        }
    Sg[1] = 1;
    fo(i,2,N) Sg[i] = 1ll*Sg[i-1]*g[i]%mo;
    invg[N] = quick_multi(Sg[N],mo-2);
    fd(i,N,2) invg[i-1] = 1ll*invg[i]*g[i]%mo;
    invg[0] = 1;
}
int main()
{
    mobius();
    fibonacci();
    pre_g();
    /*n = 100;
    fo(i,1,n) cout<<invg[i]<<" "; cout<<endl;*/
    scanf("%d",&T);
    while (T--)
        {
            scanf("%d%d",&n,&m);
            if (n > m) swap(n,m);
            res = 1;
            for (i = 1;i <= n; i=last+1)
                {
                    last = min(n/(n/i),m/(m/i));
                    res = 1ll*res*quick_multi(1ll*Sg[last]*invg[i-1]%mo,1ll*(n/i)*(m/i)%(mo-1))%mo;
                }
            printf("%d\n",res);
        }
    return 0;
}