Luogu P3704 [SDOI2017]数字表格

题意

(T) 组询问，每组询问给定 (n,m)，求 (prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}) 对 (10^9+7) 取模的值，(F) 是 Fibonacci 数列。

( exttt{Data Range:}1leq Tleq 10^3,1leq n,mleq 10^6)

题解

( exttt{Y})( exttt{LWang}) 是神！

考虑套路地枚举一发 (gcd)，有

[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{d=1}^{n}prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}[gcd(i,j)=d] ]

再考虑套路地把 (gcd) 提出来

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

然后我们可以将这个式子写成 (F_d) 的幂的形式

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

注意到上面其实就是一个套路反演

[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{d=1}^{n}F_d^{sumlimits_{k=1}^{leftlfloorfrac{n}{d} ight floor}mu(k)leftlfloorfrac{n}{dk} ight floorleftlfloorfrac{m}{dk} ight floor} ]

然后给你写回来

[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{d=1}^{n}prodlimits_{k=1}^{leftlfloorfrac{n}{d} ight floor}F_d^{mu(k)leftlfloorfrac{n}{dk} ight floorleftlfloorfrac{m}{dk} ight floor} ]

现在做代换 (T=dk)

[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{T=1}^{n}prodlimits_{dmid T}F_d^{muleft(frac{T}{d} ight)leftlfloorfrac{n}{T} ight floorleftlfloorfrac{m}{T} ight floor} ]

见证奇迹的时刻到了

[prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}F_{gcd(i,j)}=prodlimits_{T=1}^{n}left(prodlimits_{dmid T}F_d^{muleft(frac{T}{d} ight)} ight)^{leftlfloorfrac{n}{T} ight floorleftlfloorfrac{m}{T} ight floor} ]

括号里面的东西可以预处理出来，括号外面整除分块即可，由于可能要求很多次逆元，所以建议写个线性求逆。

代码

#include<bits/stdc++.h>
using namespace std;
typedef int ll;
typedef long long int li;
const ll MAXN=1e6+51,MOD=1e9+7; 
ll test,ptot,n,m,x;
ll np[MAXN],prime[MAXN],mu[MAXN],fib[MAXN],f[MAXN],prod[MAXN];
ll pr[MAXN],invf[MAXN],invp[MAXN];
inline ll read()
{
    register ll num=0,neg=1;
    register char ch=getchar();
    while(!isdigit(ch)&&ch!='-')
    {
        ch=getchar();
    }
    if(ch=='-')
    {
        neg=-1;
        ch=getchar();
    }
    while(isdigit(ch))
    {
        num=(num<<3)+(num<<1)+(ch-'0');
        ch=getchar();
    }
    return num*neg;
}
inline ll qpow(ll base,ll exponent)
{
    ll res=1;
    while(exponent)
    {
        if(exponent&1)
        {
            res=(li)res*base%MOD;
        }
        base=(li)base*base%MOD,exponent>>=1;
    }
    return res;
}
inline void sieve(ll limit)
{
    np[1]=mu[1]=f[1]=fib[1]=1;
    for(register int i=2;i<=limit;i++)
    {
        fib[i]=(fib[i-1]+fib[i-2])%MOD,f[i]=1;
        if(!np[i])
        {
            prime[++ptot]=i,mu[i]=-1;
        }
        for(register int j=1;i*prime[j]<=limit;j++)
        {
            np[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
}
inline ll calc(ll n,ll m)
{
    ll res=1,cur;
    for(register int l=1,r;l<=min(n,m);l=r+1)
    {
        r=min(n/(n/l),m/(m/l)),cur=1,cur=(li)prod[r]%MOD*invp[l-1]%MOD;
        res=(li)res*qpow(cur,(li)(n/l)*(m/l)%(MOD-1))%MOD;
    }
    return res;
}
int main()
{
    test=read(),sieve(1e6+10),prod[0]=pr[0]=1;
    for(register int i=1;i<=1e6;i++)
    {
        pr[i]=(li)pr[i-1]*fib[i]%MOD;
    }
    x=qpow(pr[1000000],MOD-2),invf[1]=1;
    for(register int i=1e6-1;i;i--)
    {
        invf[i+1]=(li)pr[i]*x%MOD,x=(li)x*fib[i+1]%MOD;
    }
    for(register int i=1;i<=1e6;i++)
    {
        for(register int j=1;i*j<=1e6;j++)
        {
            f[i*j]=(li)f[i*j]*(mu[j]==1?fib[i]:mu[j]==-1?invf[i]:1)%MOD;
        }
    }
    for(register int i=1;i<=1e6;i++)
    {
        prod[i]=(li)prod[i-1]*f[i]%MOD,pr[i]=(li)pr[i-1]*prod[i]%MOD;
    }
    x=qpow(pr[1000000],MOD-2),invp[0]=invp[1]=1;
    for(register int i=1e6-1;i;i--)
    {
        invp[i+1]=(li)pr[i]*x%MOD,x=(li)x*prod[i+1]%MOD;
    }
    for(register int i=0;i<test;i++)
    {
        n=read(),m=read(),printf("%d
",calc(n,m));
    }
}