BZOJ 4816 SDOI2017 数字表格 莫比乌斯反演
题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=4816
题意概述：
多组询问，给出一个n*m的表格，表格中的格子form[i][j]=f[gcd(i,j)]，f[i]表示fibonacci数列的第i项。
求表格中所有数字的乘积对 10^9+7 取模的结果。、
T<=1000,1<=n,m<=10^6.
分析：
一眼就可以认出来这是莫比乌斯反演。关键在于式子的推倒。
说一下推式子的要点，主要是能够灵活地利用反演公式，把式子踢来踢去！
啊这个博客实在是不好打公式......（至少我还没有研究出来），具体过程请看ATP的博客！(https://blog.csdn.net/fromatp/article/details/69944236)
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<cstdlib> 5 #include<algorithm> 6 #include<cmath> 7 #include<queue> 8 #include<set> 9 #include<map> 10 #include<vector> 11 #include<cctype> 12 using namespace std; 13 const int mo=1000000007; 14 const int maxp=80005; 15 const int maxn=1000005; 16 typedef long long LL; 17 18 int T,N; 19 int f[maxn],pri[maxp],mu[maxn],pcnt,g[maxn],_g[maxn]; 20 bool ntp[maxn]; 21 struct data{ int n,m; }Q[1005]; 22 23 void exgcd(LL a,LL b,LL &d,LL &x,LL &y){ 24 if(!b) d=a,x=1,y=0; 25 else exgcd(b,a%b,d,y,x),y-=a/b*x; 26 } 27 int inv(int a){ 28 LL d,x,y; 29 exgcd(a,mo,d,x,y); 30 return (x%mo+mo)%mo; 31 } 32 int qkpow(int x,LL y) 33 { 34 int re=1,tt=x; 35 for(int i=0;(1ll<<i)<=y;i++){ 36 if((1ll<<i)&y) re=1ll*re*tt%mo; 37 tt=1ll*tt*tt%mo; 38 } 39 return re; 40 } 41 void ready() 42 { 43 ntp[0]=ntp[1]=1,mu[1]=1; 44 for(int i=2;i<=N;i++){ 45 if(!ntp[i]) pri[++pcnt]=i,mu[i]=-1; 46 for(int j=1;j<=pcnt&&1ll*i*pri[j]<=N;j++){ 47 ntp[i*pri[j]]=1; 48 if(i%pri[j]==0){ mu[i*pri[j]]=0; break; } 49 mu[i*pri[j]]=-mu[i]; 50 } 51 } 52 f[0]=0,f[1]=1,g[0]=_g[0]=g[1]=_g[1]=1; 53 for(int i=2;i<=N;i++){ 54 f[i]=(f[i-1]+f[i-2])%mo; 55 g[i]=_g[i]=1; 56 } 57 for(int i=1;i<=N;i++) 58 for(int j=1;1ll*i*j<=N;j++){ 59 if(!mu[j]) continue; 60 if(mu[j]==1) g[i*j]=1ll*g[i*j]*f[i]%mo; 61 else _g[i*j]=1ll*_g[i*j]*f[i]%mo; 62 } 63 for(int i=1;i<=N;i++) 64 g[i]=1ll*g[i]*inv(_g[i])%mo; 65 for(int i=2;i<=N;i++) 66 g[i]=1ll*g[i]*g[i-1]%mo; 67 } 68 int solve(int n,int m) 69 { 70 int re=1,last=0; 71 for(int t=1;t<=n;t=last+1){ 72 last=min(n/(n/t),m/(m/t)); 73 re=1ll*re*qkpow(1ll*g[last]*inv(g[t-1])%mo,1ll*(n/t)*(m/t))%mo; 74 } 75 return re; 76 } 77 int main() 78 { 79 scanf("%d",&T); 80 for(int i=1;i<=T;i++){ 81 scanf("%d%d",&Q[i].n,&Q[i].m); 82 if(Q[i].n>Q[i].m) swap(Q[i].n,Q[i].m); 83 N=max(N,Q[i].n); 84 } 85 ready(); 86 for(int i=1;i<=T;i++) 87 printf("%d\n",solve(Q[i].n,Q[i].m)); 88 return 0; 89 }
