1. 程式人生 > >BZOJ4816: [Sdoi2017]數字表格(莫比烏斯反演)

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}} \right \rfloor}}{\prod_{i-1}^{\left \lfloor {\frac{M}{d}} \right \rfloor}}{Fib[d](gcd(i,j)==1)}$

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

設T=dk,則

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

那個${\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\n",(ans%mod+mod)%mod);
104     }
105     return 0;
106 }