$bzoj4816-SDOI2016$ 数字表格 莫比乌斯反演

  • 题面描述

    • \(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取模。

  • 输入格式

    • 有多组测试数据。

      第一个一个数T,表示数据组数。

      接下来T行,每行两个数n,m

      T<=1000,1<=n,m<=10^6

  • 输出格式

    • 输出T行,第i行的数是第i组数据的结果
  • 题解

    • 根据题意,我们能列出初始计算公式\(ans=\prod_{i=1}^n\prod_{j=1}^mf_{gcd(i,j)}\),这个式子在\(n,m\leq10^6​\)时很明显不可做。

    • \[\prod_{i=1}^n\prod_{j=1}^mf_{gcd(i,j)}\ (n<m)\\ \prod_{d=1}^nf_{d}^{\sum_{x=1}^{\frac{n}{d}}\sum_{y=1}^{\frac{m}{d}}[gcd(x,y)=1]}\\ \prod_{d=1}^nf_{d}^{\sum_{x=1}^{\frac{n}{d}}\sum_{y=1}^{\frac{m}{d}}\sum_{d'|gcd(x,y)}\mu(d')}\\ \prod_{d=1}^nf_{d}^{\sum_{d'=1}^{\frac{n}{d}}\lfloor\frac{n}{d*d'}\rfloor\lfloor\frac{m}{d*d'}\rfloor\mu(d')}\\ \prod_{p=1}^n\prod_{k|p}f_k^{\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor\mu(\frac{p}{k})}\\ \prod_{p=1}^n\prod_{k|p}(f_k^{\mu(\frac{p}{k})})^{\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor}\\ \prod_{p=1}^npre_p^{\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor} \]

    • 得到这个式子后,利用\(\lfloor\frac{n}{p}\rfloor\)\(\lfloor\frac{m}{p}\rfloor\)\(o(n)\)个取值,分块+前缀积处理即可,时间复杂度\(o(10^6log(mod)+Tnlog(mod))\),其中\(log(mod)\)为求逆元的快速幂复杂度。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
const int MAXN=1e6+6;
int T,n,m;
int pri[MAXN],mu[MAXN];
bool use[MAXN];
int fb[MAXN],f[MAXN][3];
int pre[MAXN];
int fm[MAXN],fr[MAXN];
int mod_pow(int x,int n){
	int ret=1;
	while (n){
		if (n&1) ret=(ll)ret*x%mod;
		x=(ll)x*x%mod;
		n>>=1;
	}
	return ret;
}
int main(){
	use[1]=0; mu[1]=1;
	for (int i=2;i<=1e6;i++){
		if (!use[i]) pri[++pri[0]]=i,mu[i]=-1;
		for (int j=1;j<=pri[0]&&i*pri[j]<=1e6;j++){
			use[i*pri[j]]=1;
			if (i%pri[j]!=0) mu[i*pri[j]]=-mu[i];
			else break;
		}
	}
	fb[1]=1; fb[2]=1;
	for (int i=3;i<=1e6;i++) fb[i]=(ll)(fb[i-1]+fb[i-2])%mod;
	for (int i=1;i<=1e6;i++){
		f[i][0]=mod_pow(fb[i],mod-2);
		f[i][1]=1; f[i][2]=fb[i];
	}
	for (int i=1;i<=1e6;i++) pre[i]=1;
	for (int i=1;i<=1e6;i++){
		for (int j=i;j<=1e6;j+=i){
			pre[j]=(ll)pre[j]*f[i][mu[j/i]+1]%mod;
		}
	}
//	for (int i=1;i<=10;i++) printf("%d\n",pre[i]);
//	cout<<pre[(int)1e6]<<endl;
	fm[0]=1;
	for (int i=1;i<=1e6;i++) fm[i]=(ll)fm[i-1]*pre[i]%mod;
	fr[(int)1e6]=mod_pow(fm[(int)1e6],mod-2);
//	cout<<fr[(int)1e6]<<endl;
	for (int i=(int)1e6-1;i>=0;i--) fr[i]=(ll)fr[i+1]*pre[i+1]%mod;
//	for (int i=1;i<=10;i++) printf("%d %d\n",fm[i],fr[i]);
	scanf("%d",&T);
	while (T--){
		scanf("%d%d",&n,&m);
		if (n>m) swap(n,m);
		int j=0,ans=1;
		for (int i=1;i<=n;i=j+1){
			j=min(n/(n/i),m/(m/i));
			ans=(ll)ans*mod_pow((ll)fr[i-1]*fm[j]%mod,(ll)(n/i)*(m/i)%(mod-1))%mod;
//			cout<<i<<" "<<j<<" "<<ans<<endl;
		}
		printf("%d\n",ans);
	}
	return 0;
}

posted @ 2019-05-02 23:19  paul120090105  阅读(64)  评论(0编辑  收藏