【BZOJ 4816】 4816: [Sdoi2017]數字表格 （莫比烏斯）

4816: [Sdoi2017]數字表格

Time Limit: 50 Sec  Memory Limit: 128 MB
Submit: 666  Solved: 312

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

HINT

Source

鳴謝infinityedge上傳

 

 

【分析】

　　額。。。加法變乘法有時候仍是忍不住求和。。

　　推式子。。

　　$$Ans=\Pi \Pi f(gcd(i,j))$$

　　$$=\Pi f(d)^{\sum 1 [gcd(i,j)==d]}$$

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

　　$$=\Pi f(d)^{\sum \mu(d')*\lfloor\dfrac{n}{d*d'}\rfloor \lfloor\dfrac{m}{d*d'}\rfloor}$$

　　【啊好辛苦

　　設D=d*d'

　　$$Ans=\Pi_{D}\Pi_{d|D} f(d)^{\lfloor\dfrac{n}{D}\rfloor \lfloor\dfrac{m}{D}\rfloor*\mu(\dfrac{D}{d})}$$

　　$$=\Pi \Pi (f(d)^{\mu(D/d)})^{\lfloor\dfrac{n}{D}\rfloor \lfloor\dfrac{m}{D}\rfloor}$$

　　設$g(D)=\Pi_{d|D}f(d)^{\mu(D/d)}$

　　這個$O(nlogn)$預處理，後面的D在詢問時$\sqrt n $分塊便可。

　　由於有快速冪，因此是$O(nlogn)+T*\sqrt n*logn$

 

 1 #include<cstdio>
 2 #include<cstdlib>
 3 #include<cstring>
 4 #include<iostream>
 5 #include<algorithm>
 6 using namespace std;
 7 #define Maxn 1000010
 8 #define Mod 1000000007
 9 
10 int g[Maxn],pri[Maxn],pl,mu[Maxn],f[Maxn];
11 int inv[Maxn],sm[Maxn];
12 bool vis[Maxn];
13 
14 int qpow(int x,int b)
15 {
16     int ans=1;
17     while(b)
18     {
19         if(b&1) ans=1LL*ans*x%Mod;
20         x=1LL*x*x%Mod;
21         b>>=1;
22     }
23     return ans;
24 }
25 
26 void init()
27 {
28     f[0]=0;f[1]=1;for(int i=2;i<=Maxn-10;i++) f[i]=(f[i-1]+f[i-2])%Mod;
29     for(int i=1;i<=Maxn-10;i++) inv[i]=qpow(f[i],Mod-2);
30     memset(vis,0,sizeof(vis));
31     mu[1]=1;
32     for(int i=2;i<=Maxn-10;i++)
33     {
34         if(!vis[i]) pri[++pl]=i,mu[i]=-1;
35         for(int j=1;j<=pl;j++)
36         {
37             if(i*pri[j]>Maxn-10) break;
38             vis[i*pri[j]]=1;
39             if(i%pri[j]==0) {mu[i*pri[j]]=0;break;}
40             mu[i*pri[j]]=-mu[i];
41         }
42     }
43     for(int i=0;i<=Maxn-10;i++) g[i]=1;
44     sm[0]=1;
45     int i;
46     for(i=1;i<=Maxn-10;i++)
47     {
48         for(int j=i;j<=Maxn-10;j+=i) if(mu[j/i]!=0)
49         {
50             if(mu[j/i]==-1) g[j]=1LL*g[j]*inv[i]%Mod;
51             else g[j]=1LL*g[j]*f[i]%Mod;
52         }
53         sm[i]=1LL*sm[i-1]*g[i]%Mod;
54     }
55 }
56 
57 
58 int main()
59 {
60     init();
61     int T;
62     scanf("%d",&T);
63     while(T--)
64     {
65         int n,m;
66         scanf("%d%d",&n,&m);
67         if(n>m) swap(n,m);
68         int ans=1;
69         for(int i=1;i<=n;)
70         {
71             int x=n/i,y=m/i,r1=n/x,r2=m/y;
72             r1=min(r1,r2);
73             ans=1LL*ans*qpow(1LL*sm[r1]*qpow(sm[i-1],Mod-2)%Mod,1LL*x*y%(Mod-1))%Mod;
74             i=r1+1;
75         }
76         printf("%d\n",ans);
77     }
78     return 0;
79 }

View Code

2017-04-26 20:21:03