AmazingCounters.com

[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);
    }
}
posted on 2017-04-14 10:06  ditoly  阅读(1157)  评论(3编辑  收藏