[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
2 3
4 5
6 7
Sample Output
1
6
960
6
960
$$\prod_{d=1}^n\prod_{i=1}^n\prod_{j=1}^mif(gcd(i,j)==d)f[gcd(i,j)]$$
直接把$f[d]$提出来
$$\prod_{d=1}^{n}f[d]^{\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]}$$
上面那个东西用莫比乌斯反演+数论分块可以$O(\sqrt n)$求
外面套的这一层也可以数论分块求
于是，我们就得到了一个$O(n)$的做法
但是显然还不够
把上面那坨东西拎出来看
$$\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]$$
太熟悉了
$$\sum_{i=1}^{n/d}\mu(i)[\frac{n}{id}][\frac{m}{id}]$$
还是老套路，
令$T=id$
直接把$T$在整个式子里面提出来
$$\prod_{T=1}^{n}\prod_{d|T}f[d]^{[n/T][m/T]\mu(T/d)}$$
有一些一样的东西
$$\prod_{T=1}^{n}(\prod_{d|T}f[d]^{\mu(T/d)})^{[n/T][m/T]}$$
然后怎么办。。。。
很明显，已经可以对$[n/T][m/T]$分块了
里面的东西不好处理，那么就暴力预处理F(T)数组
枚举i,j(i*j<=N)那么复杂度：
$\frac{n}{1}+\frac{n}{2}+.....\frac{n}{10^6}$
复杂度O(nlogn+Q√n)
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #include<cmath> 6 using namespace std; 7 typedef long long lol; 8 int Q[1001][2]; 9 const int NN=1000000; 10 int N=0; 11 int Mod=1e9+7; 12 int f[NN+5],F[NN+5],inv[NN+5],Inv[NN+5],ans; 13 int mu[NN+5],tot,prime[NN],n,m; 14 bool vis[NN+5]; 15 int qpow(int x,int y) 16 { 17 int res=1; 18 while (y) 19 { 20 if (y&1) res=1ll*res*x%Mod; 21 x=1ll*x*x%Mod; 22 y>>=1; 23 } 24 return res; 25 } 26 int pow(int x1,int x2,int y) 27 { 28 if (y<0) 29 { 30 return x2; 31 } 32 if (y==0) return 1; 33 if (y>0) 34 { 35 return x1; 36 } 37 } 38 void pre() 39 {int i,j; 40 mu[1]=1; 41 for (i=2;i<=N;i++) 42 { 43 if (vis[i]==0) 44 { 45 tot++; 46 prime[tot]=i; 47 mu[i]=-1; 48 } 49 for (j=1;j<=tot;j++) 50 { 51 if (1ll*i*prime[j]>N) break; 52 vis[i*prime[j]]=1; 53 if (i%prime[j]==0) 54 { 55 mu[i*prime[j]]=0; 56 break; 57 } 58 else mu[i*prime[j]]=-mu[i]; 59 } 60 } 61 for (i=1;i<=N;i++) 62 F[i]=1; 63 for (i=1;i<=N;i++) 64 { 65 for (j=1;j<=N;j++) 66 { 67 if (1ll*i*j>N) break; 68 F[i*j]=(1ll*F[i*j]*pow(f[i],inv[i],mu[j]))%Mod; 69 } 70 } 71 } 72 int main() 73 {int T; 74 int i,j; 75 cin>>T; 76 for (i=1;i<=T;i++) 77 { 78 scanf("%d%d",&Q[i][0],&Q[i][1]); 79 N=max(N,max(Q[i][0],Q[i][1])); 80 } 81 f[1]=1;f[0]=0; 82 for (i=2;i<=N;i++) 83 { 84 f[i]=f[i-1]+f[i-2]; 85 if (f[i]>Mod) f[i]-=Mod; 86 } 87 for (i=1;i<=N;i++) 88 inv[i]=qpow(f[i],Mod-2); 89 pre(); 90 F[0]=1; 91 for (i=1;i<=N;i++) 92 F[i]=(1ll*F[i]*F[i-1])%Mod; 93 for (j=1;j<=T;j++) 94 { 95 n=Q[j][0];m=Q[j][1]; 96 if (n>m) swap(n,m); 97 int pos=0; 98 ans=1; 99 for (i=1;i<=n;i=pos+1) 100 { 101 if (n/(n/i)<m/(m/i)) pos=n/(n/i); 102 else pos=m/(m/i); 103 ans=(1ll*ans*qpow((1ll*F[pos]*qpow(F[i-1],Mod-2))%Mod,1ll*(n/i)*(m/i)%(Mod-1)))%Mod; 104 } 105 printf("%d\n",ans); 106 } 107 }
