数字表格
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
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 #include2 #include 3 #include 4 #include 5 #include 6 #include 7 #include 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 #include2 #include 3 #include 4 #include 5 #include 6 #include 7 #include 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 }