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

 

正解:莫比乌斯函数。

水水的一道题,不过卡常数。。推导一波吧。。

$Ans=\prod_{i=1}^{n}\prod_{j=1}^{m}f(\gcd(i,j))$

$Ans=\prod_{d=1}^{min(n,m)}f(d)^{\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)==d]}$

直接跳过吧。。因为中间的都是老套路了。。

$Ans=\prod_{d=1}^{min(n,m)}f(d)^{\sum_{p=1}^{min(\left \lfloor \frac{n}{d} \right \rfloor,\left \lfloor \frac{m}{d} \right \rfloor)} \mu(p)\left \lfloor \frac{n}{dp} \right \rfloor\left \lfloor \frac{m}{dp} \right \rfloor}$

然后好像没办法往下化简了,不过这个式子用数论分块就能过了。。因为复杂度不是满的,好像是$O(Tn^{\frac{3}{4}})$吧。。

反正这就能$AC$了。。

 

 1 //It is made by wfj_2048~
 2 #include <algorithm>
 3 #include <iostream>
 4 #include <complex>
 5 #include <cstring>
 6 #include <cstdlib>
 7 #include <cstdio>
 8 #include <vector>
 9 #include <cmath>
10 #include <queue>
11 #include <stack>
12 #include <map>
13 #include <set>
14 #define rhl (1000000007)
15 #define inf (1<<30)
16 #define N (1000010)
17 #define il inline
18 #define RG register
19 #define ll long long
20 #define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
21 
22 using namespace std;
23 
24 int vis[N],inv[N],mu[N],prime[N],n,m,cnt,pos1,pos2;
25 ll f[N],ans;
26 
27 il int gi(){
28     RG int x=0,q=1; RG char ch=getchar();
29     while ((ch<'0' || ch>'9') && ch!='-') ch=getchar();
30     if (ch=='-') q=-1,ch=getchar();
31     while (ch>='0' && ch<='9') x=x*10+ch-48,ch=getchar();
32     return q*x;
33 }
34 
35 il ll qpow(RG ll a,RG ll b){
36     RG ll ans=1;
37     while (b){
38     if (b&1) ans=ans*a%rhl;
39     a=a*a%rhl,b>>=1;
40     }
41     return ans;
42 }
43 
44 il void pre(){
45     f[1]=vis[1]=mu[1]=1;
46     for (RG int i=2;i<N;++i){
47     if (!vis[i]) prime[++cnt]=i,mu[i]=-1;
48     for (RG int j=1,k;j<=cnt;++j){
49         k=i*prime[j]; if (k>=N) break; vis[k]=1;
50         if (i%prime[j]) mu[k]=-mu[i]; else break;
51     }
52     f[i]=f[i-1]+f[i-2]; if (f[i]>=rhl) f[i]-=rhl;
53     }
54     inv[0]=inv[1]=1,f[0]=1;
55     for (RG int i=2;i<N;++i)
56     mu[i]+=mu[i-1],f[i]*=f[i-1],f[i]%=rhl,inv[i]=qpow(f[i],rhl-2);
57     return;
58 }
59 
60 il void work(){
61     n=gi(),m=gi(),ans=1; if (n>m) swap(n,m);
62     for (RG int i=1;i<=n;i=pos1+1){
63     pos1=min(n/(n/i),m/(m/i)); RG ll res=0;
64     for (RG int j=1;j<=n/i;j=pos2+1){
65         pos2=min(n/i/(n/i/j),m/i/(m/i/j));
66         res+=(ll)(mu[pos2]-mu[j-1])*(n/i/j)*(m/i/j);
67     }
68     ans*=qpow(f[pos1]*(ll)inv[i-1]%rhl,res),ans%=rhl;
69     }
70     printf("%lld\n",ans); return;
71 }
72 
73 int main(){
74     File("product");
75     pre(); RG int T=gi();
76     while (T--) work();
77     return 0;
78 }

 

posted @ 2017-04-12 21:31  wfj_2048  阅读(175)  评论(0编辑  收藏