关键词:
第一次见到这种形式,推了两步就被卡住了。没想到还能这么迁移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 查看详情