[BZOJ]4816: [Sdoi2017]數字表格

Time Limit: 50 Sec  Memory Limit: 128 MB

Description

　　Doris剛剛學習了fibonacci數列。用f[i]表示數列的第i項，那么

　　f[0]=0

　　f[1]=1

　　f[n]=f[n-1]+f[n-2],n>=2

　　Doris用老師的超級計算機生成了一個n×m的表格，第i行第j列的格子中的數是f[gcd(i,j)]，其中gcd(i,j)表示i,j的最大公約數。Doris的表格中共有n×m個數，她想知道這些數的乘積是多少。答案對10^9+7取模。

Input

　　有多組測試數據。

　　第一個一個數T，表示數據組數。

　　接下來T行，每行兩個數n,m

　　T<=1000,1<=n,m<=10^6

Output

　　輸出T行，第i行的數是第i組數據的結果

Sample Input

　　3
　　2 3
　　4 5
　　6 7

Sample Output

　　1
　　6
　　960

Solution

　　反正推推公式唄。

　　題目所求即

　　$$ Ans(n,m)= \prod_{i=1}^{n} \prod_{j=1}^{m} f[gcd(i,j)] $$

　　枚舉$k=gcd(i,j)$，則

　　$$ Ans(n,m)= \prod_{k=1}^{min(n,m)} f[k]^{ \sum_{i=1}^{\left \lfloor \frac{n}{k} \right \rfloor} \sum_{j=1}^{\left \lfloor \frac{m}{k} \right \rfloor} [gcd(i,j)=1] } $$

　　其中$[gcd(i,j)=1]$表示若$gcd(i,j)=1$，該式值為1，否則為0。

　　利用莫比烏斯函數的性質，得到

　　$$ Ans(n,m)= \prod_{k=1}^{min(n,m)} f[k]^{ \sum_{i=1}^{\left \lfloor \frac{n}{k} \right \rfloor} \sum_{j=1}^{\left \lfloor \frac{m}{k} \right \rfloor} \sum_{d|gcd(i,j)}\mu(d) } $$

　　對每個$d$統計$\mu(d)$被算了幾次，則

　　$$ Ans(n,m)= \prod_{k=1}^{min(n,m)} f[k]^{ \sum_{d=1}^{\left\lfloor\frac{min(n,m)}{k}\right\rfloor} \left \lfloor \frac{n}{kd} \right \rfloor \left \lfloor \frac{m}{kd} \right \rfloor \mu(d) } $$

　　這個式子還是不好算，我們自然希望能把$d$拖出來，於是令$p=kd$

　　$$ Ans(n,m)= \prod_{p=1}^{min(n,m)} \prod_{k|p} f[k]^{\left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor \mu(\frac{p}{k})} = \prod_{p=1}^{min(n,m)} (\prod_{k|p} f[k]^{\mu(\frac{p}{k})})^{\left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor} $$

　　這個式子就舒服多了，發現對於每個$p$，$\prod_{k|p} f[k]^{\mu(\frac{p}{k})}$的值是固定的，與$n$和$m$無關，於是我們先用篩法預處理出每個$p$對應的這個式子的值，前綴積一下，利用$ \left \lfloor \frac{n}{p} \right \rfloor \left \lfloor \frac{m}{p} \right \rfloor $只有 $O(\sqrt{n}+\sqrt{m})$種取值，我們可以在$O((\sqrt{n}+\sqrt{m})logMOD)$時間內計算並回答每次詢問（其中$O(logMOD)$為快速冪），那么總復雜度為$O((max(n,m)+T(\sqrt{n}+\sqrt{m}))logMOD)$。

Code

#include<cstdio>
#include<algorithm>
using namespace std;
inline int read()
{
    int x;char c;
    while((c=getchar())<'0'||c>'9');
    for(x=c-'0';(c=getchar())>='0'&&c<='9';)x=x*10+c-'0';
    return x;
}
#define MN 1000000
#define MOD 1000000007
int f[MN+5],fp[MN+5][3],u[MN+5],mu[MN+5],p[MN+5],pn;
int ff[MN+5],fs[MN+5],fr[MN+5];
int pw(int x,int y)
{
    int r=1;if(y<0)y+=MOD-1;
    for(;y;y>>=1,x=1LL*x*x%MOD)if(y&1)r=1LL*r*x%MOD;
    return r;
}
int main()
{
    int T=read(),n,m,i,j,ans;
    for(f[1]=mu[1]=1,i=2;i<=MN;++i)
    {
        if((f[i]=f[i-1]+f[i-2])>=MOD)f[i]-=MOD;
        if(!u[i])p[++pn]=i,mu[i]=-1;
        for(j=1;i*p[j]<=MN&&(u[i*p[j]]=1);++j)
            if(i%p[j])mu[i*p[j]]=-mu[i];else break;
    }
    for(i=1;i<=MN;++i)fp[i][0]=pw(f[i],-1),fp[i][1]=1,fp[i][2]=f[i];
    for(i=1;i<=MN;++i)ff[i]=1;
    for(i=1;i<=MN;++i)for(j=i;j<=MN;j+=i)ff[j]=1LL*ff[j]*fp[i][mu[j/i]+1]%MOD;
    for(fs[0]=i=1;i<=MN;++i)fs[i]=1LL*fs[i-1]*ff[i]%MOD;
    for(fr[i=MN]=pw(fs[MN],-1);i--;)fr[i]=1LL*fr[i+1]*ff[i+1]%MOD;
    while(T--)
    {
        n=read();m=read();ans=1;
        for(i=1;i<=n&&i<=m;i=j+1)
        {
            j=min(n/(n/i),m/(m/i));
            ans=1LL*ans*pw(1LL*fs[j]*fr[i-1]%MOD,1LL*(n/i)*(m/i)%(MOD-1))%MOD;
        }
        printf("%d\n",ans);
    }
}