[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} ight floor} sum_{j=1}^{left lfloor frac{m}{k} ight floor} [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} ight floor} sum_{j=1}^{left lfloor frac{m}{k} ight floor} sum_{d|gcd(i,j)}mu(d) } $$

　　对每个$d$统计$mu(d)$被算了几次，则

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

　　这个式子还是不好算，我们自然希望能把$d$拖出来，于是令$p=kd$

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

　　这个式子就舒服多了，发现对于每个$p$，$prod_{k|p} f[k]^{mu(frac{p}{k})}$的值是固定的，与$n$和$m$无关，于是我们先用筛法预处理出每个$p$对应的这个式子的值，前缀积一下，利用$ left lfloor frac{n}{p} ight floor left lfloor frac{m}{p} ight floor $只有 $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
",ans);
    }
}