「SDOI2017」数字表格 - 莫比乌斯反演

题目大意

    Doris刚刚学习了fibonacci数列。用$f[i]$表示数列的第$i$项，那么
      $f[0]=0$
      $f[1]=1$
      $f[n]=f[n-1]+f[n-2],n\ge2$
    Doris用老师的超级计算机生成了一个$n\times m$的表格，第$i$行第$j$列的格子中的数是$f[\gcd(i,j)]$，其中$\gcd(i,j)$表示$i,j$的最大公约数。Doris的表格中共有$n\times m$个数，她想知道这些数的乘积是多少。答案对$10^9+7$取模。

---

题目分析

题目要求的是：

$$\begin{aligned} Ans &=\prod_{i=1}^n\prod_{j=1}^mf[\gcd(i,j)] \\ &=\prod_kf[k]^{\sum_{i=1}^n\sum_{j=1}^m [\gcd(i,j)==k]} \\ &=\prod_kf[k]^{\sum_{i=1}^{\lfloor\frac nk\rfloor}\sum_{j=1}^{\lfloor\frac mk\rfloor} [\gcd(i,j)==1]} \\ &=\prod_kf[k]^{\sum_{i=1}^{\lfloor\frac nk\rfloor}\sum_{j=1}^{\lfloor\frac mk\rfloor}\sum_{d\mid i,d\mid j}\mu(d)} \\ &=\prod_kf[k]^{\sum_d^{\lfloor\frac nk\rfloor}\mu(d)\lfloor\frac n{kd}\rfloor\lfloor\frac m{kd}\rfloor} \\ &=\prod_T^{\min(n,m)}\prod_{k\mid T}f(k)^{\mu(\frac Tk)\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor}\\ &=\prod_T^{\min(n,m)}(\prod_{k\mid T}f(k)^{\mu(\frac Tk)})^{\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor} \\ &=\prod_T^{\min(n,m)}g(T)^{\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor} \\ \end{aligned}$$

因为$g(T)=\prod_{k\mid T}f(k)^{\mu(\frac Tk)}$与$n,m$无关，所以可以考虑对其进行预处理。

筛出莫比乌斯函数后按照约数累加$f(k)$即可。
剩下的可以分块跳跃，达到$O(\sqrt n)$的复杂度。

虽然一个数的约数个数最多$\sqrt n$个，但是对于$1\rightarrow n$的所有约数个数求和，平均是$O(\log n)$的。
因此，总时间复杂度为$O(n\log n+T\sqrt n)$。

---

代码

#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<climits>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;

typedef long long LL;

inline const int Get_Int() {
	int num=0,bj=1;
	char x=getchar();
	while(x<'0'||x>'9') {
		if(x=='-')bj=-1;
		x=getchar();
	}
	while(x>='0'&&x<='9') {
		num=num*10+x-'0';
		x=getchar();
	}
	return num*bj;
}

const int maxn=1000005;
const LL mod=1e9+7;

int vst[maxn],Prime[maxn],cnt=0; //prime
int Mobius[maxn]; //mobius
int t;
LL n,m,f[maxn],invf[maxn],g[maxn],sum[maxn],invsum[maxn];

void Mobius_Table(int n) {
	Mobius[1]=1;
	for(int i=2; i<=n; i++) {
		if(!vst[i]) {
			Prime[++cnt]=i;
			Mobius[i]=-1;
		}
		for(int j=1; j<=cnt&&i*Prime[j]<=n; j++) {
			vst[i*Prime[j]]=1;
			if(i%Prime[j]==0) {
				Mobius[i*Prime[j]]=0;
				break;
			}
			Mobius[i*Prime[j]]=-Mobius[i];
		}
	}
}

LL Quick_Pow(LL a,LL b) {
	LL sum=1;
	while(b) {
		if(b&1)sum=sum*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return sum;
}

void initialize(int n) {
	Mobius_Table(n);
	f[1]=invf[1]=g[1]=1;
	for(int i=2; i<=n; i++) {
		g[i]=1;
		f[i]=(f[i-1]+f[i-2])%mod;
		invf[i]=Quick_Pow(f[i],mod-2);
	}
	for(int i=1; i<=n; i++)
		for(int j=i; j<=n; j+=i) {
			if(Mobius[j/i]==1)g[j]=g[j]*f[i]%mod;
			if(Mobius[j/i]==-1)g[j]=g[j]*invf[i]%mod;
		}
	sum[0]=invsum[0]=1;
	for(int i=1; i<=n; i++)sum[i]=sum[i-1]*g[i]%mod,invsum[i]=Quick_Pow(sum[i],mod-2);
}

int main() {
	t=Get_Int();
	initialize(1000000);
	while(t--) {
		n=Get_Int();
		m=Get_Int();
		int next;
		LL ans=1;
		for(int i=1; i<=min(n,m); i=next+1) {
			next=min(n/(n/i),m/(m/i));
			ans=ans*Quick_Pow(sum[next]*invsum[i-1]%mod,(n/i)*(m/i)%(mod-1))%mod;
		}
		printf("%lld\n",ans);
	}
	return 0;
}