[bzoj4816]数字表格莫比乌斯反演

author author     2022-10-04     684

关键词:

第一次见到这种形式,推了两步就被卡住了。没想到还能这么迁移QAQ,真强!

题目设 $f(i)$为第$i$项斐波那契数列

要求计算

$$ans=prod_{i=1}^n prod_{j=1}^m f(gcd(i,j))$$

枚举$d=gcd$得到

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

然后反演得到

$$ans=prod_{d=1}^n f(d)^{sum_{T=1}^{lfloor frac{n}{d} floor}mu(T) lfloor frac{n}{dT} floor lfloor frac{m}{dT} floor}$$

然后呢,令$p=dT$,得到

$$ans=left( prod_{d=1}^nf(d)^{mu(frac{p}{d})} ight) ^{sum_{p=1}^n  lfloor frac{n}{p} floor lfloor frac{m}{p} floor}$$

这样我们就发现我们可以对$p$进行分块去搞了

$$g(p)=prod_{d|p} f(d)^{mu(frac{p}{d})}$$

然后我们就可以通过预处理$g(i)$和它的前缀积,然后枚举$p$分块去搞,用$last$前缀积除以$i-1$的前缀积就可以得到这一段的乘积,然后将它们对 $lfloor frac{n}{i} floor lfloor frac{m}{i} floor$做快速幂即可

我们怎么预处理$g(i)$呢?

我们可以枚举$d$,然后统计它对它的倍数的贡献

观察式子$f(d)^{mu(frac{p}{d})}$,$mu$的取值只有三种。为$0$时不用管,为$1$时需要乘$f(i)$,为$-1$时需要乘$f(i)^{-1}$

关键问题就是预处理出$f(i)^{-1}$

我们线性推就好了。设$$s[x]=prod_{i=1}^x f(x)$$ $$rs[x]=s[x]^{-1}$$ $$ni[x]=f(x)^{-1}$$

于是$$rs[x-1]=rs[x]*f[x]$$ $$ni[x]=rs[x]*s[x-1]$$我们只需要求出$rs[n]$然后倒着推一遍就行了(撒花!)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#define N 1000100
#define pos(i,a,b) for(LL i=(a);i<=(b);i++)
#define pos2(i,a,b) for(LL i=(a);i>=(b);i--)
using namespace std;
#define LL long long
const int p=(LL)1e9+7;
int t;
LL n,m;
int notprime[N],prime[N],mu[N];
LL x,y;
LL s[N],rs[N],f[N],ni[N];
LL g[N];
void kuogcd(LL a,LL b,LL &x,LL &y){
	if(!b){
		x=1;y=0;return;
	}
	kuogcd(b,a%b,x,y);
	LL temp=x;x=y;y=temp-a/b*y;
}
LL qpow(LL a,LL k){
	LL res=1;
	while(k>0){
		if(k&1){
			res=(res*a)%p;
		}
		a=(a*a)%p;k>>=1;
	}
	return res;
}
void get_pre(){
	notprime[1]=mu[1]=1;
	pos(i,2,N-10){
		if(!notprime[i]){
			prime[++prime[0]]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=prime[0]&&prime[j]*i<=N-10;j++){
			notprime[i*prime[j]]=1;
			if(i%prime[j]==0){
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
	f[1]=1;
	pos(i,2,N-10) f[i]=(f[i-1]+f[i-2])%p;
	s[0]=1;
	pos(i,1,N-10){
		s[i]=s[i-1]*f[i]%p;
	}
	kuogcd(s[N-10],p,x,y);
	while(x<0) x+=p;
	rs[N-10]=x%p;
	pos2(i,N-10,1){
		rs[i-1]=rs[i]*f[i]%p;
		ni[i]=rs[i]*s[i-1]%p;
	}
	pos(i,1,N-10) g[i]=1;
	pos(i,1,N-10){
		pos(j,1,N-10){
			if(i*j>N-10) break;
			int pp=i*j;
			if(mu[j]==1) g[pp]=(g[pp]*f[i])%p;	
			else if(mu[j]==-1) g[pp]=(g[pp]*ni[i])%p;
		}
	}
	g[0]=1;
	pos(i,1,N-10) g[i]=(g[i]*g[i-1])%p;
}
LL get_ans(){
	LL res(1);
	LL last(0);
	for(LL i=1;i<=n;i=last+1){
		last=min(n/(n/i),m/(m/i));
		kuogcd(g[i-1],p,x,y);
		while(x<0) x+=p;
		res=res*qpow(g[last]*x%p,(n/i)*(m/i))%p;
	}
	return res%p;
}
int main(){
	get_pre();
	scanf("%d",&t);
	while(t--){
		scanf("%lld%lld",&n,&m);if(n>m) n^=m,m^=n,n^=m;
		printf("%lld
",get_ans());
	}
	return 0;
}

  

[bzoj4816]数字表格莫比乌斯反演

第一次见到这种形式,推了两步就被卡住了。没想到还能这么迁移QAQ,真强!题目设$f(i)$为第$i$项斐波那契数列要求计算$$ans=prod_{i=1}^nprod_{j=1}^mf(gcd(i,j))$$枚举$d=gcd$得到$$ans=prod_{d=1}^nf(d)^{sum_{i=1}^nsum_{j=1}^m[gcd(i,j)==d]}$$然后反演得... 查看详情

bzoj4816[sdoi2017]数字表格莫比乌斯反演

【BZOJ4816】[Sdoi2017]数字表格DescriptionDoris刚刚学习了fibonacci数列。用f[i]表示数列的第i项,那么f[0]=0f[1]=1f[n]=f[n-1]+f[n-2],n>=2Doris用老师的超级计算机生成了一个n×m的表格,第i行第j列的格子中的数是f[gcd(i,j)],其中gcd(i,j)表示i,... 查看详情

bzoj.4816.[sdoi2017]数字表格(莫比乌斯反演)(代码片段)

题目链接总感觉博客园的\(Markdown\)很。。\(gouzhi\),可以看这的。这个好像简单些啊,只要不犯sb错误\(Description\)  用\(f[i]\)表示\(Fibonacci\)数列的第\(i\)项,求\[\prod_i=1^n\prod_j=1^mf[\gcd(i,j)]\mod(10^9+7)\]\(Solution\)\[\beginalignedAns& 查看详情

bzoj4816:[sdoi2017]数字表格

二次联通门: BZOJ4816:[Sdoi2017]数字表格    /*BZOJ4816:[Sdoi2017]数字表格莫比乌斯反演妈呀,我的μ没有啦*/#include<cstdio>#include<iostream>#definergregister#defineMax1000050#definemo1000000007 查看详情

bzoj2154crash的数字表格(莫比乌斯反演)

【BZOJ2154】Crash的数字表格(莫比乌斯反演)题面BZOJ简化题意:给定(n,m)求[sum_{i=1}^nsum_{j=1}^mlcm(i,j)]题解以下的一切都默认(n<m)我们都知道(lcm(i,j)=frac{ij}{gcd(i,j)})所以所求化简[sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)}]看到(g 查看详情

bzoj2154:crash的数字表格[莫比乌斯反演]

2154:Crash的数字表格TimeLimit: 20Sec  MemoryLimit: 259MBSubmit: 2924  Solved: 1091[Submit][Status][Discuss]Description今天的数学课上,Crash小朋友学习了最小公倍数(LeastCommonMultiple)。对 查看详情

bzoj4816[sdoi2017]数字表格(代码片段)

题面https://www.lydsy.com/JudgeOnline/problem.php?id=4816题解显然是莫比乌斯反演首先得出然后发现我们要把d提出去这样就好做了跟SDOI2015的一道题类似因为$leftlfloorfracnpightfloorleftlfloorfracmpightfloor$只有$(sqrtn+sqrtm)$种取 查看详情

bzoj2154crash的数字表格莫比乌斯反演

...最小公倍数,他画了一张N*M的表格。每个格子里写了一个数字,其中第i行第j列的那个格子里写着数为LCM(i,j)。一个4*5的表格如下:1&nbs 查看详情

bzoj2154crash的数字表格(莫比乌斯反演)

题意:       思路:如上From http://blog.csdn.net/regina8023/article/details/44243911最后的F(x,y)的推法和求gcd(x,y)=1的(x,y)对数差不多,只不过在推导过程中把原来1的地方换成x*y。那么我们预处理出i^2*u[i]的前缀和... 查看详情

crash的数字表格(莫比乌斯反演)

Crash的数字表格Description今天的数学课上,Crash小朋友学习了最小公倍数(LeastCommonMultiple)。对于两个正整数a和b,LCM(a,b)表示能同时被a和b整除的最小正整数。例如,LCM(6,8)=24。回到家后,Crash还在想着课上学的东西,为了研究最小... 查看详情

bzoj2301莫比乌斯反演

求$(i,j)=k$的一系列模板题之一。但是这里i,j是有下界的,注意用容斥去掉重复组,其他都一样了。  /**@Date:2017-09-0919:21:18*@FileName:bzoj2301莫比乌斯反演多组范围内GCD=k.cpp*@Platform:Windows*@Author:Lweleth([email protected])*@Link:https: 查看详情

bzoj3309dzylovesmath(莫比乌斯反演)

【BZOJ3309】DZYLovesMath(莫比乌斯反演)题面求\[\sum_i=1^a\sum_j=1^bf(gcd(a,b))\]其中,\(f(x)\)表示\(x\)分解质因数之后,最高的幂次题解完全不会莫比乌斯反演了。先来推式子\[\sum_d=1^a\sum_i=1^a/d\sum_j=1^b/d[gcd(i,j)=1]f(d)\]\[\sum_d=1^a 查看详情

bzoj4816数字表格

...函数,话说上次考莫比鸟斯就是去年吧,好像题目名也叫数字表格,只不过多了一个前缀"Crash的"。慢慢推吧, 查看详情

bzoj2045双亲数莫比乌斯反演

【BZOJ2045】双亲数Description小D是一名数学爱好者,他对数字的着迷到了疯狂的程度。我们以d=gcd(a,b)表示a、b的最大公约数,小D执著的认为,这样亲密的关系足可以用双亲来描述,此时,我们称有序数对(a,b)为d的双亲数。与正常双... 查看详情

bzoj2440:[中山市选2011]完全平方数二分答案+容斥原理+莫比乌斯反演

...挺好想的。二分一个答案val,需要[1,val]之间存在的合法数字个数>=k即可。怎么判断呢?可以容斥,开始的时候有ans=val个,但是这里明显有些数字不符合。ans-=([1...val]中有多少个2^2倍 +[1...val] 查看详情

[bzoj]1101zap||莫比乌斯反演

原题求x<=a,y<=b,并且gcd(x,y)=d的x,y对数(sum^{n}_{x=1}sum^{m}_{y=1}[gcd(x,y)==d])(=sum^{lfloorfrac{n}{d} floor}_{x=1}sum^{lfloorfrac{m}{d} floor}_{y=1}[gcd(x,y)==1])根据莫比乌斯反演,可以把式子化为(=sum^ 查看详情

[bzoj2301]problemb莫比乌斯反演

纪念莫比乌斯反演首题!同时感谢Antileaf学长的耐心讲解!我们可以用容斥原理搞令f(d)为1<=x<=n,1<=y<=m且gcd(x,y)=d的数对(x,y)的个数则可设: 可得F(i)为1<=x<=n,1<=y<=m且i|gcd(x,y)的数对(x,y)的个数化简可得莫比乌斯... 查看详情

bzoj-2818莫比乌斯反演初步

要使用分块的技巧#include<iostream>#include<algorithm>#include<cstdio>#include<cstring>#include<cstdlib>#include<cmath>#include<string>#include<vector>#include&l 查看详情