[BZOJ4816][SDOI2017]数字表格(莫比乌斯反演)

4816: [Sdoi2017]数字表格
Time Limit: 50 Sec Memory Limit: 128 MB
Submit: 1259 Solved: 625
[Submit][Status][Discuss]
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上传

[Submit][Status][Discuss]

不知道为什么要用fibonacci，感觉既没有用到矩乘又没有用到$gcd(f[i],f[j])=f[gcd(i,j)]$的性质。

首先列出连乘式，可以发现很像莫比乌斯反演，先试着推一下式子，从每个数出现的次数入手。

$$Ans(n,m)=prod_{d}f(d)^{sum_{i=1}^{n}sum_{j=1}^m[(i,j)=d]}=prod_{d=1}^{min(n,m)}f(d)^{sum_{p=1}^{frac{min(n,m)}{d}}mu(d)lfloor frac{n}{pd} floor lfloorfrac{m}{pd} floor}$$
到这里，有一种想法（可以拿60分）：
设$$g(n,m,d)=sum_{i=1}^{min(n,m)}mu(d)lfloor frac{n}{i} floor lfloor frac{m}{i} floor$$
这样$$Ans(n,m)=prod_{d=1}^{lfloor frac{min(n,m)}{d} floor}f(d)^{g(lfloor frac{n}{d} floor lfloor frac{m}{d} floor,d)}$$
这个看上去用分块优化可以做到$$egin{aligned}O(Tint_1^nsqrt{frac{n}{x}}dx)& =O(Tsqrt{n}int_1^nsqrt{frac{1}{x}}dx)\ & =O(2Tsqrt{n}sqrt{n})\ & =O(Tn)end{aligned}$$

于是就有60分了。
那么我们继续化简刚才的式子：
$$egin{aligned}Ans(n,m)& =prod_{d=1}^{min(n,m)}f(d)^{sum_{d|T}^{min(n,m)}mu(frac{T}{d})lfloor frac{n}{T} floor lfloor frac{m}{T} floor}\ & =prod_{T=1}^{min(n,m)}prod_{d|T}f(d)^{mu(frac{T}{d})lfloor frac{n}{T} floorlfloorfrac{m}{T} floor}\ & =prod_{T=1}^{min(n,m)}g(T)^{lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor}end{aligned}$$ $g(T)$可以预处理出来，这样总复杂度就是$O(nlog n+Tsqrt{min(n,m)})$了。
注意$mu$在指数上时会有$-1$，这个要求逆元，如果先预处理所有$f[i]$的逆元的话会快一倍。
这道题总体并不难，主要就是将莫比乌斯反演中的一些套路从加法变成乘法了，实现的时候有几个细节注意一下就好。

（以后再也不手写LaTeX莫比乌斯反演的题解了）

 1 #include<cstdio>
 2 #include<algorithm>
 3 #define rg register int
 4 #define rep(i,l,r) for (rg i=l; i<=r; i++)
 5 typedef long long ll;
 6 using namespace std;
 7 
 8 const int N=1000100,mod=1000000007;
 9 bool b[N];
10 int n,m,T,ans,tot,p[N],f[N],g[N],G[N],G1[N],miu[N],F[N][3];
11 
12 int ksm(int a,int b){
13     int res;
14     for (res=1; b; a=1ll*a*a%mod,b>>=1)
15         if (b & 1) res=1ll*res*a%mod;
16     return res;
17 }
18 
19 void pre(){
20     miu[1]=1;
21     for (rg i=2; i<N; i++){
22         if (!b[i]) p[++tot]=i,miu[i]=-1;
23         for (rg j=1; j<=tot && p[j]*i<N; j++){
24             rg t=p[j]*i; b[t]=1;
25             if (i%p[j]) miu[t]=-miu[i]; else break;
26         }
27     }
28     for (rg i=1; i<N; i++) g[i]=1,F[i][0]=ksm(f[i],mod-2),F[i][1]=1,F[i][2]=f[i];
29     for (rg i=1; i<N; i++)
30         for (rg j=i; j<N; j+=i)
31             g[j]=1ll*g[j]*F[i][miu[j/i]+1]%mod;
32     G[0]=1; G[1]=g[1]; for (int i=2; i<N; i++) G[i]=1ll*G[i-1]*g[i]%mod;
33     G1[N-1]=ksm(G[N-1],mod-2); for (int i=N-2; ~i; i--) G1[i]=1ll*G1[i+1]*g[i+1]%mod;
34 }
35 
36 int main(){
37     freopen("product.in","r",stdin);
38     freopen("product.out","w",stdout);
39     f[1]=1; for (rg i=2; i<N; i++) f[i]=(f[i-1]+f[i-2])%mod;
40     pre();
41     for (scanf("%d",&T); T--; ){
42         scanf("%d%d",&n,&m); rg lst=0;
43         if (n>m) swap(n,m); ans=1;
44         for (rg i=1; i<=n; i=lst+1){
45             lst=min(n/(n/i),m/(m/i));
46             ans=1ll*ans*ksm(1ll*G[lst]*G1[i-1]%mod,1ll*(n/i)*(m/i)%(mod-1))%mod;
47         }
48         printf("%d
",ans);
49     }
50     return 0;
51 }