題目鏈接
以下除法均指下取整
\[\prod_{i=1}^n \prod_{j=1}^mf(gcd(i,j))\=\prod_{x}f(x)^{\sum_i \sum_j [gcd(i,j)=x] } \=\prod_{x}f(x)^{ \sum_{x|d} \mu(\frac{d}{x})\frac{n}{d} \frac{m}{d} }\=\prod_{d}(\prod_{x|d} f(x)^{\mu(\frac{d}{x})} )^{\frac{n}{d}\frac{m}{d}}\]
然後把中間那個括號裏的東西看成\(g(d)\)預處理,即可做到\(O(nlogn)\)預處理\(O(\sqrt n *logn)\)

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<vector>
#include<algorithm>
#include<cmath>
#define P puts("lala")
#define cp cerr<<"lala"<<endl
#define ln putchar(‘\n‘)
#define pb push_back
#define fi first
#define se second
 

#define mkp make_pair
using namespace std;
inline int read()
{
    char ch=getchar();int g=1,re=0;
    while(ch<‘0‘||ch>‘9‘) {if(ch==‘-‘)g=-1;ch=getchar();}
    while(ch<=‘9‘&&ch>=‘0‘) re=(re<<1)+(re<<3)+(ch^48),ch=getchar();
    return re*g;
}
typedef long long ll;
typedef
 
 pair<int,int> pii;

const int N=1000050;
const int inf=0x3f3f3f3f;
const int mod=1e9+7;
inline ll qpow(ll a,ll n)
{
    ll ans=1;
    for(;n;n>>=1,a=a*a%mod) if(n&1) ans=ans*a%mod;
    return ans;
}

int f[N],finv[N],g[N],ginv[N],prime[N],cnt=0,mu[N];
bool isnotprime[N];

void init()
{
    int n=N-50;
    f[1]=1; g[1]=1;
    for(int i=2;i<=n;++i) f[i]=(f[i-1]+f[i-2])%mod,g[i]=1;
    for(int i=1;i<=n;++i) finv[i]=qpow(f[i],mod-2);

    isnotprime[1]=1; mu[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!isnotprime[i]) prime[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&i*prime[j]<=n;++j)
        {
            isnotprime[i*prime[j]]=1;
            if(!(i%prime[j])) break;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;++i) for(int j=i;j<=n;j+=i)
    {
        if(mu[j/i]>0) g[j]=1ll*g[j]*f[i]%mod;
        else if(!mu[j/i]) ;
        else g[j]=1ll*g[j]*finv[i]%mod;
    }
    g[0]=1; ginv[0]=1;
    for(int i=1;i<=n;++i) g[i]=1ll*g[i]*g[i-1]%mod;
    for(int i=1;i<=n;++i) ginv[i]=qpow(g[i],mod-2);
}

int main()
{
#ifndef ONLINE_JUDGE
    freopen("1.in","r",stdin);freopen("1.out","w",stdout);
#endif
    init();
    int T=read();
    for(int cas=1;cas<=T;++cas)
    {
        int n=read(),m=read();
        int M=min(n,m),ans=1;
        for(int i=1,j;i<=M;i=j+1)
        {
            j=min(n/(n/i),m/(m/i));
            ans=1ll*ans*qpow(1ll*g[j]*ginv[i-1]%mod,1ll*(n/i)*(m/i))%mod;
        }
        printf("%d\n",ans);
    }
    return 0;
}

bzoj4816[SDOI2017]數字表格