解题：SDOI 2017 数字表格

题面

反演题，推式子么=。=

$prodlimits_{d=1}^{min(n,m)}prodlimits_{i=1}^nprodlimits_{j=1}^m[gcd(i,j)==d]fib[d]$

把$fib[d]$前提，前面的连乘就跑到指数上去了

$prodlimits_{d=1}^{min(n,m)}fib[d]^{sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)==d]}$

开始反演那坨指数，等等这玩意不是做过么=。=

$sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)==d]$

$sumlimits_{i=1}^{leftlfloorfrac{n}{d} ight floor}sumlimits_{j=1}^{leftlfloorfrac{m}{d} ight floor}[gcd(i,j)==1]$

$sumlimits_{i=1}^{min(leftlfloorfrac{n}{d} ight floor,leftlfloorfrac{m}{d} ight floor)}μ(i)leftlfloorfrac{n}{id} ight floorleftlfloorfrac{m}{id} ight floor$

于是把$id$捉出来，在原来的整个式子里枚举$id$(不是那个$id$，都懂)

$prodlimits_{k=1}^{min(n,m)}(prodlimits_{d|k}fib[d]^{μ(frac{k}{d})})^{leftlfloorfrac{n}{k} ight floorleftlfloorfrac{m}{k} ight floor}$

停，可以做了

对于$prodlimits_{d|k}fib[d]^{μ(frac{k}{d})}$，预处理，大力把每个数乘到倍数上去，复杂度$O(nlog n)$

对于$leftlfloorfrac{n}{k} ight floorleftlfloorfrac{m}{k} ight floor$这个指数，可以数论分块，这样再加个快速幂每次回答复杂度就是$O(sqrt nlog mod)$了，可能有点卡常？我倒是一次过了

#include<cstdio>
#include<cstring> 
#include<algorithm>
using namespace std;
const int N=1000005,M=1000000,mod=1e9+7;
int fib[N],ifb[N],pfb[N],ipf[N];
int pri[N],npr[N],mul[N];
int T,n,m,x,y,mn,cnt,ans;
int qpow(int x,int k)
{
    if(k==1) return x;
    int tmp=qpow(x,k/2);
    return k%2?1ll*tmp*tmp%mod*x%mod:1ll*tmp*tmp%mod;
}
void exGCD(int a,int b,int &x,int &y)
{
    if(!b) {x=1,y=0; return ;} 
    exGCD(b,a%b,y,x); y-=a/b*x;
}
int Inv(int b)
{
    exGCD(b,mod,x,y);
    return (x+mod)%mod;
}
void prework()
{
    register int i,j;
    npr[1]=true,mul[1]=1,fib[1]=1,ifb[1]=1,pfb[1]=1;
    for(i=2;i<=M;i++)
    {
        fib[i]=(fib[i-1]+fib[i-2])%mod;
        ifb[i]=Inv(fib[i]),pfb[i]=1;
        if(!npr[i])     
            pri[++cnt]=i,mul[i]=-1;
        for(j=1;j<=cnt&&i*pri[j]<=M;j++)
        {
            npr[i*pri[j]]=true;
            if(i%pri[j]==0) break;
            else mul[i*pri[j]]=-mul[i];
        }
    }
    for(i=1;i<=M;i++)
        for(j=i;j<=M;j+=i)
            if(mul[j/i]) pfb[j]=1ll*pfb[j]*((~mul[j/i])?fib[i]:ifb[i])%mod;
    pfb[0]=ipf[0]=1;
    for(i=1;i<=M;i++)
    {
        ipf[i]=Inv(pfb[i]);
        pfb[i]=1ll*pfb[i-1]*pfb[i]%mod;
        ipf[i]=1ll*ipf[i-1]*ipf[i]%mod;
    }
}
int main()
{
    register int i,j;
    scanf("%d",&T),prework();
    while(T--)
    {
        scanf("%d%d",&n,&m);
        mn=min(n,m),ans=1;
        for(i=1;i<=mn;i=j+1)
        {
            j=min(n/(n/i),m/(m/i));
            ans=1ll*ans*qpow(1ll*pfb[j]*ipf[i-1]%mod,1ll*(n/i)*(m/i)%(mod-1))%mod;
        }
        printf("%d
",ans);
    }
    return 0;
}