BZOJ4816: [Sdoi2017]数字表格(莫比乌斯反演)

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

解题思路：

让你求：${prod_{i=1}^{N}}{prod_{j=1}^{M}}{Fib[gcd(i,j)]}$

不妨设$N{leq}M$,Fib[i]表示斐波那契数列第i项。

接下来就是喜闻乐见的公式推导了QAQ

$Ans={prod_{i=1}^{N}}{prod_{j=1}^{M}}{Fib[gcd(i,j)]}$

$={prod_{d=1}^{N}}{prod_{d|i}^{N}}{prod_{d|j}^{M}}{Fib[d](gcd(i,j)==d)}$

$={prod_{d=1}^{N}}{prod_{i-1}^{left lfloor {frac{N}{d}} ight floor}}{prod_{i-1}^{left lfloor {frac{M}{d}} ight floor}}{Fib[d](gcd(i,j)==1)}$

$={prod_{d=1}^{N}}{Fib[d]}^{{prod_{k=1}^{left lfloor {frac{N}{d}} ight floor}}{mu(k)}{left lfloor {frac{N}{dk}} ight floor}{left lfloor {frac{M}{dk}} ight floor}}$

设T=dk，则

$Ans={prod_{T=1}^{N}}{prod_{d|T}}{{Fib[d]}^{{mu(frac{T}{d})}{left lfloor {frac{N}{T}} ight floor}{left lfloor {frac{M}{T}} ight floor}}}$

那个${prod_{d|T}}{Fib[d]}^{mu(frac{T}{d})}$对于某个T来说是一定的，那么预处理一下就好了

代码：

  1 #include<cstdio>
  2 #include<cstring>
  3 #include<algorithm>
  4 typedef long long lnt;
  5 const int N=1000010;
  6 const lnt mod=(lnt)(1e9+7);
  7 int prime[N];
  8 lnt miu[N];
  9 lnt fib[N];
 10 lnt ifi[N];
 11 lnt ansp[N];
 12 bool vis[N];
 13 int cnt;
 14 int T;
 15 lnt Inv(lnt x)
 16 {
 17     lnt ans=1;
 18     lnt y=mod-2;
 19     while(y)
 20     {
 21         if(y&1)
 22             ans=ans*x%mod;
 23         x=x*x%mod;
 24         y=y/2;
 25     }
 26     return ans;
 27 }
 28 lnt ksm(lnt x,lnt y)
 29 {
 30     lnt ans=1;
 31     while(y)
 32     {
 33         if(y&1)
 34             ans=ans*x%mod;
 35         x=x*x%mod;
 36         y=y/2;
 37     }
 38     return ans;
 39 }
 40 void gtp(void)
 41 {
 42     miu[1]=1;
 43     for(int i=2;i<N;i++)
 44     {
 45         if(!vis[i])
 46         {
 47             prime[++cnt]=i;
 48             miu[i]=-1;
 49         }
 50         for(int j=1;j<=cnt&&prime[j]*i<N;j++)
 51         {
 52             vis[i*prime[j]]=true;
 53             if(i%prime[j]==0)
 54             {
 55                 miu[i*prime[j]]=0;
 56                 break;
 57             }
 58             miu[i*prime[j]]=-miu[i];
 59         }
 60     }
 61     ansp[0]=ansp[1]=ansp[2]=1;
 62     fib[1]=1;
 63     ifi[1]=Inv(fib[1]);
 64     fib[2]=1;
 65     ifi[2]=Inv(fib[2]);
 66     for(int i=3;i<N;i++)
 67     {
 68         fib[i]=(fib[i-1]+fib[i-2])%mod;
 69         ansp[i]=1;
 70         ifi[i]=Inv(fib[i]);
 71     }
 72     for(int p=1;p<N;p++)
 73     {
 74         if(!miu[p])
 75             continue;
 76         for(int j=1;j*p<N;j++)
 77         {
 78             lnt tmp=(miu[p]==1)?fib[j]:ifi[j];
 79             ansp[j*p]=(ansp[j*p]*tmp)%mod;
 80         }
 81     }
 82     for(int i=1;i<N;i++)
 83         ansp[i]=(ansp[i]*ansp[i-1])%mod;
 84     return ;
 85 }
 86 int main()
 87 {
 88     gtp();
 89     scanf("%d",&T);
 90     while(T--)
 91     {
 92         lnt n,m;
 93         scanf("%lld%lld",&n,&m);
 94         if(n>m)
 95             std::swap(n,m);
 96         lnt ans=1;
 97         for(lnt i=1,j;i<=n;i=j+1)
 98         {
 99             j=std::min(n/(n/i),m/(m/i));
100             lnt tmp=ansp[j]*Inv(ansp[i-1])%mod;
101             ans=ans*ksm(tmp,(lnt)(n/i)*(lnt)(m/i)%(mod-1))%mod;
102         }
103         printf("%lld
",(ans%mod+mod)%mod);
104     }
105     return 0;
106 }