P3704 [SDOI2017] 数字表格题解-526互联

一、题目描述：

　　用 $f_i$ 表示斐波那契数列的第 $i$ 项，那么有：

　　$ f_0=0,f_1=1;f_n=f_{n-1}+f_{n-2},n\ge2 $

　　现在有一个 $n$ 行 $m$ 列的数字表格，第 $i$ 行第 $j$ 列的数字是 $f_{\gcd(i,j)}$ 。

　　求这个表格所有数的乘积。共有 $T$ 组数据，答案对 $10^9+7$ 取模。

　　数据范围：$1\le T\le 10^3,1\le n,m\le 10^6$ 。

二、解题思路：

　　我迄今位置做过的最难的莫比乌斯反演题！

　　话不多说，直接开始推式子。

　　$$\prod_{i=1}^{n}\prod_{j=1}^mf(\gcd(i,j))=$$

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

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

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

　　$$令g(x,y)=\sum_{i=1}^x\sum_{j=1}^y[\gcd(i,j)=1]$$

　　$$则原式=\prod_{d=1}^nf(d)^{g(\lfloor\frac{n}{d} \rfloor,\lfloor\frac{m}{d} \rfloor)}$$

　　$$考虑g(x,y)=\sum_{i=1}^x\sum_{j=1}^y[\gcd(i,j)=1]=$$

　　$$\sum_{i=1}^x\sum_{j=1}^y\sum_{d\mid \gcd(x,y)}\mu(d)=\sum_{d=1}^x\mu(d)\times \lfloor\frac{x}{d} \rfloor\lfloor\frac{y}{d} \rfloor$$

　　$$求出斐波那契数列的乘积前缀和，再求出逆元即可$$

　　$$时间复杂度O(T\times (n+\sqrt n\times log_2^n))$$

　　$$时间复杂度显然过不去，考虑优化$$

　　$$将g(x,y)带入原式，得到原式=\prod_{d=1}^nf(d)^{\sum_{k=1}^{n/d}\mu(k)\times \lfloor\frac{n}{dk} \rfloor\lfloor\frac{m}{dk} \rfloor}$$

　　$$令T=dk,则原式=\prod_{T=1}^n \left(\prod_{d\mid T}f(d)^{\mu(\frac{T}{d})}\right)^{ \lfloor\frac{n}{T} \rfloor\lfloor\frac{m}{T} \rfloor}$$

　　$$令g(x)=\prod_{d\mid x}f(d)^{\mu(\frac{x}{d})},原式=\prod_{T=1}^n g(T)^{ \lfloor\frac{n}{T} \rfloor\lfloor\frac{m}{T} \rfloor}$$

　　$$考，原来里面的g(x)可可以直接暴力求出来，调和级数！$$

　　$$时间复杂度O((T\sqrt n+n)\times log_2^n)$$

三、完整代码：

 1 #include<iostream>
 2 #define N 1000010
 3 #define lim 1000000
 4 #define ll long long
 5 #define M 1000000007
 6 using namespace std;
 7 ll ans,f[N],g[N],inf[N],ing[N];
 8 int T,n,m,cnt,mu[N],p[N],vis[N];
 9 ll ksm(ll base,ll q)
10 {
11     ll res=1;
12     while(q){
13         if(q&1) res*=base,res%=M;
14         base*=base,base%=M,q>>=1;
15     }
16     return res;
17 }
18 void pre_work()
19 {
20     mu[1]=f[1]=g[0]=ing[0]=1;
21     for(int i=2;i<=lim;i++)
22         f[i]=(f[i-1]+f[i-2])%M;
23     for(int i=1;i<=lim;i++)
24         inf[i]=ksm(f[i],M-2),g[i]=1;
25     for(int i=2;i<=lim;i++)
26     {
27         if(!vis[i]) p[++cnt]=i,mu[i]=-1;
28         for(int j=1;j<=cnt&&i*p[j]<=lim;j++)
29         {
30             vis[i*p[j]]=1;
31             if(i%p[j]==0)break;
32             mu[i*p[j]]=-mu[i];
33         }
34     }
35     for(int i=1;i<=lim;i++)
36         for(ll j=i;j<=lim;j+=i)
37         {
38             if(mu[j/i]==1) g[j]=g[j]*f[i]%M;
39             if(mu[j/i]==-1) g[j]=g[j]*inf[i]%M;
40         }
41     for(int i=1;i<=lim;i++)
42         g[i]=g[i]*g[i-1]%M,ing[i]=ksm(g[i],M-2);
43 }
44 void solve()
45 {
46     cin>>n>>m; ans=1;
47     if(n>m) swap(n,m);
48     for(int l=1,r;l<=n;l=r+1)
49     {
50         r=min(n/(n/l),m/(m/l));ans%=M;
51         ans*=ksm(ing[l-1]*g[r]%M,1ll*(n/l)*(m/l)%(M-1));
52     }
53     cout<<ans%M<<'\n';
54 }
55 int main()
56 {
57     ios::sync_with_stdio(false);
58     cin.tie(0);cout.tie(0);
59     cin>>T;pre_work();
60     for(int i=1;i<=T;i++)
61         solve();
62     return 0;
63 }