Description
Solution
一开始的时候我是这么推的（(f(n))表示斐波那契数列的第(n)项）
[egin{aligned}
Ans&=prod_{x=1}^{min(n,m)}f(x)^{(sum_{i=1}^nsum_{j=1}^m[gcd(i,j)=x])}\
&=prod_{x=1}^{min(n,m)}f(x)^{sum_{e=1}^{min(lfloorfrac nx
floor,lfloor frac mx
floor)}mu(e)lfloorfrac n {ex}
floorlfloorfrac m {ex}
floor}\
end{aligned}
]
然后我想，根号分段套根号分段，(mathcal O(sqrt n(log+sqrt n)))解决！嗯不错，一看数据组数1000........不过还是能拿70的分。
题解的思路非常神。
假设我们能构造一个函数(g)，使得
[f(n)=prod_{d|n}g(d)
]
那么答案就变成
[egin{aligned}
Ans&=prod_{i=1}^nprod_{j=1}^mf({gcd(i,j))}\
&=prod_{i=1}^nprod_{j=1}^mprod_{d|i,d|j}g(d)\
&=prod_{d=1}^{min(n,m)}g(d)^{lfloorfrac nd
floorlfloorfrac md
floor}
end{aligned}
]
这样就可以在(mathcal O(2sqrt n))的时间内处理每一个询问了。前提是我们知道(g)及其前缀积。
考虑式子(f(n)=prod_{d|n}g(d))十分像莫比乌斯反演，能否用类似的形式反演出(g)呢？
在(sum)的意义下
[g_n=sum_{d|n}mu(d)f({frac nd})
]
反演的本质是通过加减来容斥出所需要的组合。而在乘法的意义下，不就是通过乘除来容斥出所需要的组合吗？所以有：
[g_n=prod_{d|n}f({frac nd})^{mu(d)}=prod_{d|n}f(d)^{mu(frac nd)}
]
因此我们可以在(mathcal O(n lg n))的时间内处理出(g)的取值和前缀积。那么上面的根号分段也就迎刃而解了。
总体思路是仿造莫比乌斯反演构造一个可求函数，利用该函数化简式子使得答案求和式变得简明且复杂度较低，再通过传统根号分段求解。
Code
#include <cstdio>
using namespace std;
const int N=1e6+10,MOD=1e9+7,PMOD=MOD-1;
int n,m;
int fib[N],ifib[N],g[N],ig[N];
bool vis[N];
int p[N],pcnt,mu[N];
inline int min(int x,int y){return x<y?x:y;}
inline int fmi(int x,int y){
int res=1;
for(;y;x=1LL*x*x%MOD,y>>=1)
if(y&1) res=1LL*res*x%MOD;
return res;
}
void sieve(){
mu[1]=1;
for(int i=2;i<=1e6;i++){
if(!vis[i]) p[++pcnt]=i,mu[i]=-1;
for(int j=1;j<=pcnt&&i*p[j]<=1e6;j++){
int x=i*p[j];
vis[x]=true;
if(i%p[j]==0){
mu[x]=0;
break;
}
mu[x]=-mu[i];
}
}
}
void prework(){
sieve();
fib[0]=0; fib[1]=1; ifib[1]=1;
for(int i=2;i<=1e6;i++){
fib[i]=(fib[i-2]+fib[i-1])%MOD;
ifib[i]=fmi(fib[i],MOD-2);
}
for(int i=1;i<=1e6;i++) g[i]=1;
for(int d=1;d<=1e6;d++)
for(int n=d;n<=1e6;n+=d)
g[n]=1LL*g[n]*(mu[n/d]==1?fib[d]:(mu[n/d]==-1?ifib[d]:1))%MOD;
g[0]=ig[0]=1;
for(int i=1;i<=1e6;i++){
g[i]=1LL*g[i]*g[i-1]%MOD;
ig[i]=fmi(g[i],MOD-2);
}
}
int main(){
prework();
int T,up,ans;
scanf("%d",&T);
while(T--){
scanf("%d%d",&n,&m);
up=min(n,m);
ans=1;
for(int i=1,j;i<=up;i=j+1){
j=min(n/(n/i),m/(m/i));
ans=1LL*ans*fmi(1LL*g[j]*ig[i-1]%MOD,1LL*(n/i)*(m/i)%PMOD)%MOD;
}
printf("%d
",ans<0?ans+MOD:ans);
}
return 0;
}