【BZOJ4816】【SDOI2017】数字表格 [莫比乌斯反演]

数字表格

Time Limit: 50 Sec  Memory Limit: 128 MB
[Submit][Status][Discuss]

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

Output

　　输出T行，第i行的数是第i组数据的结果

Sample Input

　　3
　　2 3
　　4 5
　　6 7

Sample Output

　　1
　　6
　　960

HINT

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

Solution

　　运用莫比乌斯反演，得到式子：

　　这样我们对于内外分块即可，复杂度为O(n^(0.75)*T)。

　　然后我们会发现在BZOJ上过不去，怎么办呢？卡常！BearChild运用了如下的卡常技巧：

　　　　1.  读入优化； 2. O(n)预处理逆元； 3. 内嵌汇编实现乘和取模； 4. 记录n/i，避免多次除法。

Code

 

  1 #include<iostream>  
  2 #include<string>  
  3 #include<algorithm>  
  4 #include<cstdio>  
  5 #include<cstring>  
  6 #include<cstdlib>  
  7 #include<cmath>
  8 using namespace std; 
  9 typedef long long s64;
 10  
 11 const int ONE = 1e6+5;
 12 const int MOD = 1e9+7;
 13 const int PHI = 1e9+6;
 14  
 15 int T;
 16 int n,m;
 17 int prime[ONE],p_num,miu[ONE];
 18 int F[ONE];
 19 bool isp[ONE];
 20 s64 Ans;
 21  
 22 int get() 
 23 {
 24         int res=1,Q=1;  char c;
 25         while( (c=getchar())<48 || c>57)
 26         if(c=='-')Q=-1;
 27         if(Q) res=c-48; 
 28         while((c=getchar())>=48 && c<=57) 
 29         res=res*10+c-48; 
 30         return res*Q; 
 31 }
 32  
 33 int Quickpow(int a,int b)
 34 {
 35         int res = 1;
 36         while(b)
 37         {
 38             if(b&1) res = (s64)res*a%MOD;
 39             a = (s64)a*a%MOD;
 40             b>>=1;
 41         }
 42         return res;
 43 }
 44  
 45 void Deal_first(int MaxN)
 46 {
 47         F[1]=1;
 48         F[0]=0; for(int i=2; i<=MaxN; i++) F[i] = ((s64)F[i-1]+F[i-2]) % MOD;
 49         F[0]=1; for(int i=2; i<=MaxN; i++) F[i] = (s64)F[i]*F[i-1] % MOD;
 50          
 51         miu[1] = 1;
 52         for(int i=2; i<=MaxN; i++)
 53         {
 54             if(!isp[i])
 55                 prime[++p_num] = i, miu[i] = -1;
 56             for(int j=1; j<=p_num, i*prime[j]<=MaxN; j++)
 57             {
 58                 isp[i * prime[j]] = 1;
 59                 if(i % prime[j] == 0)
 60                 {
 61                     miu[i * prime[j]] = 0;
 62                     break;
 63                 }
 64                 miu[i * prime[j]] = -miu[i];
 65             }
 66             miu[i] += miu[i-1];
 67         }
 68 }
 69  
 70 int f(int n,int m)
 71 {
 72         if(n > m) swap(n,m);
 73         s64 Ans = 0;
 74         for(int i=1,j=0; i<=n; i=j+1)
 75         {
 76             j = min(n/(n/i), m/(m/i));
 77             Ans += (s64)(n/i) * (m/i)%PHI * ((s64)(miu[j] - miu[i-1] + PHI)%PHI) % PHI;
 78             Ans %= PHI;
 79         }
 80         return Ans;
 81 }
 82  
 83 void Solve()
 84 {
 85         n=get();    m=get();
 86         if(n > m) swap(n,m);
 87         Ans = 1;
 88         for(int i=1,j=0; i<=n; i=j+1)
 89         {
 90             j = min(n/(n/i), m/(m/i));
 91             Ans = Ans * Quickpow( (s64)F[j] * Quickpow(F[i-1],MOD-2) % MOD , f(n/i,m/i) % PHI) % MOD;
 92         }
 93         printf("%lld\n",Ans);
 94 }
 95  
 96 int main()
 97 {
 98         Deal_first(ONE-1);
 99         T = get();
100         while(T--)
101             Solve();
102         return 0;
103 }
非卡常版

 

  1 #include<iostream>  
  2 #include<string>  
  3 #include<algorithm>  
  4 #include<cstdio>  
  5 #include<cstring>  
  6 #include<cstdlib>  
  7 #include<cmath>
  8 using namespace std; 
  9 typedef long long s64;
 10  
 11 const int ONE = 1e6+2;
 12 const int MOD = 1e9+7;
 13 const int PHI = 1e9+6;
 14  
 15 int T;
 16 int n,m;
 17 int prime[ONE],p_num,miu[ONE];
 18 int Niyu[ONE];
 19 int F[ONE];
 20 bool isp[ONE];
 21 int Ans;
 22  
 23 inline int get() 
 24 {
 25         int res=1,Q=1;  char c;
 26         while( (c=getchar())<48 || c>57)
 27         if(c=='-')Q=-1;
 28         if(Q) res=c-48; 
 29         while((c=getchar())>=48 && c<=57) 
 30         res=res*10+c-48; 
 31         return res*Q; 
 32 }
 33 
 34 inline int modmul(const int &a, const int &b,const int &M)
 35 {
 36         int ret;
 37         __asm__ __volatile__("\tmull %%ebx\n\tdivl %%ecx\n" : "=d"(ret) : "a"(a), "b"(b), "c"(M));
 38         return ret;
 39 }
 40 
 41 inline int Quickpow(int a,int b)
 42 {
 43         int res = 1;
 44         while(b)
 45         {
 46             if(b&1) res = modmul(res,a,MOD);
 47             a = modmul(a,a,MOD);
 48             b>>=1;
 49         }
 50         return res;
 51 }
 52  
 53 inline void Deal_first(int MaxN)
 54 {
 55         F[0]=0; F[1]=1;
 56         int val=1;
 57         for(int i=2; i<=MaxN; i++)
 58         {
 59             F[i] = F[i-1]+F[i-2];
 60             if(F[i] >= MOD) F[i] -= MOD;
 61             val = modmul(val,F[i],MOD);
 62         }
 63         Niyu[MaxN] = Quickpow(val, MOD-2);
 64         for(int i=MaxN-1;i>=1;i--) Niyu[i] = modmul(Niyu[i+1],F[i+1],MOD);
 65         Niyu[0] = Niyu[1]; F[0]=1;
 66         for(int i=1; i<=MaxN; i++) F[i] = modmul(F[i],F[i-1],MOD);
 67          
 68         miu[1] = 1;
 69         for(int i=2; i<=MaxN; i++)
 70         {
 71             if(!isp[i])
 72                 prime[++p_num] = i, miu[i] = -1;
 73             for(int j=1; j<=p_num, i*prime[j]<=MaxN; j++)
 74             {
 75                 isp[i * prime[j]] = 1;
 76                 if(i % prime[j] == 0)
 77                 {
 78                     miu[i * prime[j]] = 0;
 79                     break;
 80                 }
 81                 miu[i * prime[j]] = -miu[i];
 82             }
 83             miu[i] += miu[i-1];
 84         }
 85 }
 86  
 87 inline int f(int n,int m)
 88 {
 89         if(n > m) swap(n,m);
 90         int Ans = 0;
 91         for(int i=1,j=0; i<=n; i=j+1)
 92         {
 93             int x=n/i, y=m/i;
 94             j = min(n/x, m/y);
 95             Ans = ((s64)Ans + modmul(modmul(x,y,PHI) , ((s64)miu[j] - miu[i-1] + PHI), PHI) )%PHI;
 96         }
 97         return Ans;
 98 }
 99  
100 inline void Solve()
101 {
102         n=get();    m=get();
103         if(n > m) swap(n,m);
104         Ans = 1;
105         for(int i=1,j=0; i<=n; i=j+1)
106         {
107             int x=n/i, y=m/i;
108             j = min(n/x, m/y);
109             Ans = (s64)modmul(Ans , Quickpow( modmul(F[j],Niyu[i-1],MOD) , f(x,y)), MOD);
110         }
111         printf("%d\n",Ans);
112 }
113  
114 int main()
115 {
116         Deal_first(ONE-1);
117         T = get();
118         while(T--)
119             Solve();
120         return 0;
121 }
卡常版

 

posted @ 2017-04-11 21:46  BearChild  阅读(635)  评论(4)  编辑  收藏  举报