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

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

1、枚举gcd，转化为ans=∏nd=1f(d)h(d)，其中h(d)=∑ni=1∑mj=1[gcd(i,j)==d]；

2、h(d)是莫比乌斯反演的模板题，直接上结论：h(d)=∑⌊nd⌋k=1μ(k)⌊ndk⌋⌊mdk⌋，所以转化为ans=∏nd=1f(d)h(d)=∑⌊nd⌋k=1μ(k)⌊ndk⌋⌊mdk⌋；

做到这里通过枚举约数的方法就能够得到答案，这个方法的时间复杂度是O((n+m)34)，但是乘上一个T就有可能超时。这样至少就是60分。

3、设T=dk，并且改变一下运算顺序（因为全部都是乘法），转换为ans=∏nT=1(∏d|Tf(d)μ(Td))⌊nT⌋⌊mT⌋；

4、记g(x)=∏d|Tf(d)μ(Td)，现在如果我们已经知道了g(x)的前缀积，我们就可以在O(N−−√logN)的复杂度内求出每次询问，其中根号是分块的近似复杂度，log是快速幂的复杂度；

5、那么如何求g(x)呢？可以通过枚举d，对于每个满足d|T的g(T)都乘上f(d)μ(Td)，只需提前预处理μ(x)、f(x)和f(x)的逆元即可，根据调和级数的相关证明可以得到这部分的时间复杂度是O(lnNlogN)；

6、那么如何求f(x)的逆元呢？显然不能O(P)求出所有数的逆元。可以选择扩欧来做，时间复杂度为O(NlnN)，不过网上有O(N)的做法：设Si=∏i1ai，然后求出Sn的逆元invSn，然后通过递推式invSi−1=ai∗invSi求出invS数组，再通过invai−1=invSi∗Si−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
,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
= quick_multi(Sf
,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
= quick_multi(Sg
,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;
}