太才了
(prod limits_{i=1}^{n}prodlimits_{j=1}^{m}f[gcd(i,j)])
(prodlimits_{k=1}^{n}f[k]^{sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)==k]})
幂我们很熟悉
就是
(g(x)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i,j)==x]=sumlimits_{i=1}^{frac{n}{x}}sumlimits_{j=1}^{frac{m}{x}}[gcd(i,j)==1])
(f(x)=sumlimits_{x|d}^ng(d)=frac{n}{x}frac{m}{x})
(g(x)=sumlimits_{x|d}^nmu(frac{d}{x})f(d))
(g(x)=sumlimits_{x|d}^nmu(frac{d}{x})frac{n}{d}frac{m}{d})
(g(x)=sumlimits_{d=1}^{n/x}mu(d)frac{n}{xd}frac{m}{xd})
带回去
(prodlimits_{k=1}^{n}f[k]^{sumlimits_{d=1}^{n/k}mu(d)frac{n}{dk}frac{m}{dk}})
令i=d*k
(prodlimits_{k=1}^{n}f[k]^{sumlimits_{k|i}^{n}mu(i/k)frac{n}{i}frac{m}{i}})
(prodlimits_{k=1}^{n}prodlimits_{k|i}^{n}f[k]^{mu(i/k)frac{n}{i}frac{m}{i}})
(prodlimits_{i=1}^{n}prodlimits_{k|i}^{n}f[k]^{mu(i/k)frac{n}{i}frac{m}{i}})
设(Au(k)=prodlimits_{k|i}^{n}f[k]^{mu(i/k)})
(ans=prodlimits_{i=1}^{n}Au(k)^{frac{n}{i}frac{m}{i}})
注意,错误
幂次是不能直接模除的,用(x^{p-1}equiv 1(mod p))
代码
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll N=1e6+7,mod=1e9+7;
ll read() {
ll x=0,f=1;char s=getchar();
for(;s>'9'||s<'0';s=getchar()) if(s=='-') f=-1;
for(;s>='0'&&s<='9';s=getchar()) x=x*10+s-'0';
return x*f;
}
ll f[N],mu[N],pri[N],Au[N],cnt;
bool vis[N];
ll q_pow(ll a,ll b) {
ll ans=1;
while(b) {
if(b&1) ans=ans*a%mod;
a=a*a%mod;
b>>=1;
}
return ans;
}
void Euler() {
mu[1]=vis[1]=1;
for(ll i=2;i<=1000000;++i) {
if(!vis[i]) {pri[++cnt]=i,mu[i]=-1;}
for(ll j=1;j<=cnt&&i*pri[j]<=1000000;++j) {
vis[i*pri[j]]=1;
if(i%pri[j]==0) {
mu[i*pri[j]]=0;
break;
} else mu[i*pri[j]]=-mu[i];
}
}
f[1]=1;
for(ll i=2;i<=1000000;++i) f[i]=(f[i-1]+f[i-2])%mod;
for(ll i=0;i<=1000000;++i) Au[i]=1;
for(ll i=1;i<=1000000;++i) {
ll inv=q_pow(f[i],mod-2);
for(ll j=i;j<=1000000;j+=i) {
if(mu[j/i]==-1) Au[j]=Au[j]*inv%mod;
else if(mu[j/i]==1) Au[j]=Au[j]*f[i]%mod;
}
}
for(ll i=1;i<=1000000;++i) Au[i]=(Au[i]*Au[i-1])%mod;
}
int main() {
Euler();
ll T=read();
while(T--) {
ll n=read(),m=read();
if(n>m) swap(n,m);
ll ans=1;
for(ll l=1,r=1;l<=n;l=r+1) {
r=min(n/(n/l),m/(m/l));
ans=ans*q_pow(Au[r]*q_pow(Au[l-1],mod-2)%mod,(n/l)*(m/l)%(mod-1))%mod;
}
cout<<ans<<"
";
}
return 0;
}