bzoj4816: [Sdoi2017]数字表格

Posted 2020-10-31 lrj998244353

Description

Doris刚刚学习了\(fibonacci\)数列。用\(f[i]\)表示数列的第\(i\)项，那么

\[\begin{align}f[0]&=0\\f[1]&=1\\f[n]&=f[n-1]+f[n-2],n>=2\end{align}\]

Doris用老师的超级计算机生成了一个\(n\times m\)的表格，第\(i\)行第\(j\)列的格子中的数是\(f[\mathrm{gcd}(i,j)]\)，其中\(\mathrm{gcd}(i,j)\)表示\(i\),\(j\)的最大公约数。Doris的表格中共有\(n\times m\)个数，她想知道这些数的乘积是多少。答案对\(10^9+7\)取模。

Input

有多组测试数据。第一个一个数\(T\)，表示数据组数。

接下来\(T\)行，每行两个数\(n,m\)

\(T \leq 1000\),\(1 \leq n\),\(m \leq 10^6\)

Output

输出\(T\)行，第\(i\)行的数是第\(i\)组数据的结果

Sample Input

3
2 3
4 5
6 7

Sample Output

1
6
960

题解

要求的东西是
\[ ans=\prod_{i=1}^n\prod_{j=1}^mf[\mathrm{gcd}(i,j)] \]
假设\(n<m\)，根据常见套路，枚举约数
\[ ans=\prod_{d=1}^nf(d)^{\sum_{i=1}^n\sum_{j=1}^m[\mathrm{gcd}(i,j)==d]} \]
把指数拿出来化简
\[ \begin{align} &\sum_{i=1}^n\sum_{j=1}^m[\mathrm{gcd}(i,j)==d]\=&\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}[\mathrm{gcd}(i,j)==1]\=&\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}\sum_{k|\mathrm{gcd}(i,j)}\mu(k)\=&\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(k)\left \lfloor \frac{n}{kd} \right \rfloor\left \lfloor \frac{m}{kd} \right \rfloor \end{align} \]
令\(T=kd\)，代入上式
\[ \sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(k)\left \lfloor \frac{n}{kd} \right \rfloor\left \lfloor \frac{m}{kd} \right \rfloor=\sum_{T=1}^n\mu(\frac{T}{d})\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor\\]
把\(T\)放到外面枚举
\[ ans=\prod_{T=1}^n\prod_{d|T}f(d)^{\mu(\frac{T}{d})\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor} \]
令函数\(g(x)=\prod_{d|x}f(d)^{\mu(\frac{x}{d})}\)
\[ ans=\prod_{T=1}^ng(T)^{\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor} \]
只要我们求出了\(g(x)\)的前缀积以及对应的逆元就可以整除分块了。

只要预处理出莫比乌斯函数和\(f\)以及\(f\)的逆元，\(g(x)\)就可以通过枚举每个数的筛它的倍数求出，复杂度是调和级数。

问题就是如何求所有\(f\)以及对应的逆元\(invf\)。

这里我们需要一个辅助的东西\(mulf\)，表示\(f\)的前缀积，即\(mulf(x)=\prod_{i=1}^x f(i)\).

假设我们知道了\(mulf(i+1)\)的逆元\(invm(i+1)\)，那么有
\[ \begin{align} invm(i+1) &\equiv invm(i)\times invf(i+1) (\mathrm{mod}\ P)\invm(i)&\equiv invm(i+1)\times f(i+1)(\mathrm{mod}\ P) \end{align} \]
而\(invf(i)\equiv invm(i)*mulf(i-1)(\mathrm{mod} P)\),这样我们就可以求出\(f\)以及\(invf\)了。

用同样的方法，我们也可以求出\(g(x)\)的前缀积以及对应的逆元。

代码

#include<bits/stdc++.h>
#define MAXN 1000010
#define P 1000000007
#define LL long long
namespace IO{
    char buf[1<<15],*fs,*ft;
    inline char gc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
    inline int qr(){
        int x=0,rev=0,ch=gc();
        while(ch<'0'||ch>'9'){if(ch=='-')rev=1;ch=gc();}
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=gc();}
        return rev?-x:x;}
}using namespace IO;
using namespace std;
LL f[MAXN],invm[MAXN],invf[MAXN],g[MAXN],invg[MAXN],mulf[MAXN],mulg[MAXN],ans;
int p[MAXN],cnt,u[MAXN],T,N,M;
bool np[MAXN]; 
inline LL qp(LL x,LL y){
    LL ret=1;
    while(y){
        if(y&1)ret=ret*x%P;
        x=x*x%P;y>>=1; 
    }
    return ret;
}
int main(){
    #ifndef ONLINE_JUDGE
    freopen("product.in","r",stdin);
    freopen("product.out","w",stdout);
    #endif
    f[0]=0;
    f[1]=u[1]=np[1]=g[1]=mulf[1]=mulg[1]=invm[0]=1;
    for(int i=2;i<=1000000;i++){
        f[i]=(f[i-1]+f[i-2])%P;
        mulf[i]=mulf[i-1]*f[i]%P;
        g[i]=1;
        if(!np[i])p[++cnt]=i,u[i]=-1;
        for(int j=1;j<=cnt;j++){
            int t=p[j]*i;
            if(t>1000000)break;
            np[t]=1;
            if(i%p[j]==0){u[t]=0;break;}
            u[t]=-u[i];
        }
    }
    invm[1000000]=qp(mulf[1000000],P-2);
    invf[1000000]=qp(f[1000000],P-2);
    for(int i=999999;i>=1;i--){
        invm[i]=invm[i+1]*f[i+1]%P;
        invf[i]=mulf[i-1]*invm[i]%P;
    } 
    for(int i=2;i<=1000000;i++){
        for(int j=1,k;(k=j*i)<=1000000;j++){
            if(u[j]==1)g[k]=g[k]*f[i]%P;
            else if(u[j]==-1)g[k]=g[k]*invf[i]%P;
        }
        mulg[i]=mulg[i-1]*g[i]%P;
    }
    invg[1000000]=qp(mulg[1000000],P-2);
    for(int i=999999;i>=0;i--)invg[i]=invg[i+1]*g[i+1]%P;
    T=qr();
    while(T--){
        N=qr();M=qr();ans=1;
        if(N>M)swap(N,M);
        for(int i=1,j;i<=N;i=j+1){
            j=min(N/(N/i),M/(M/i));
            ans=ans*qp(mulg[j]*invg[i-1]%P,1LL*(N/i)*(M/i)%(P-1))%P;
        }
        printf("%lld\n",ans);
    }
    return 0;
}