SDOI2017数字表格

莫比乌斯反演,适当结合2个式子,换元T=gd,调和级数

题意:求(prod^n_{i=1}prod^m_{j=1}fib[gcd(i,j)],fib[1]=fib[2]=1)

(T)组数组(n,min[1,1e6],Tin[1,1e3])

SOL:

易得(ans=prod^n_{g=1}fib(g)^{S(n/g,m/g)})

(S(n,m)=sum^n_{i=1}sum^m_{j=1}[gcd(i,j)]==1=sum^n_{d=1}mu(d)(n/d)(m/d))

默认(n<=m)

注意:

  1. 结合2个式子一起
  2. 新建变量代替乘积,设(T=gd)

(ans=prod^n_{T=1}(prod_{g|T}fib(g)^{mu(frac{T}{g})})^{(n/T)(m/T)})

加上预处理即可

时间复杂度(O(nlog_n+Tsqrt nlog_n))

#include<bits/stdc++.h>
using namespace std;
inline int read(){
	int x=0,f=1;char c=getchar();
	while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}
	while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
	return f==1?x:-x;
}
#define ll long long
const int N=1e6+4,mod=1e9+7;
inline int ksm(int x,ll r){
	int ret=1;
	for(int i=0;(1ll<<i)<=r;i++){
		if((r>>i)&1)ret=(ll)ret*x%mod;
		x=(ll)x*x%mod; 
	}
	return ret;
}
int n,m,mu[N],pri[N],b[N],fib[N],s[N],cnt,ans;
inline void pre(){
	mu[1]=1;
	for(int i=2;i<N;i++){
		if(!b[i]){pri[++cnt]=i;mu[i]=-1;}
		for(int j=1;j<=cnt&&i*pri[j]<N;j++){
			b[i*pri[j]]=1;
			if(i%pri[j])mu[i*pri[j]]=-mu[i];
			else{
				mu[i*pri[j]]=0;
				break;
			}
		}
	}
	for(int i=0;i<N;i++)s[i]=1;
	fib[1]=fib[2]=1;
	for(int i=3;i<N;i++)fib[i]=(fib[i-1]+fib[i-2])%mod;
	for(int i=1,inv;i<N;i++){
		inv=ksm(fib[i],mod-2);
		for(int j=i;j<N;j+=i)
			if(mu[j/i]==1)s[j]=(ll)s[j]*fib[i]%mod;
			else if(mu[j/i]==-1)s[j]=(ll)s[j]*inv%mod; 
	}
	for(int i=2;i<N;i++)s[i]=(ll)s[i]*s[i-1]%mod;
}
int main(){
	pre();
	int T=read();
	while(T--){
		n=read();m=read();
		if(n>m)n^=m^=n^=m;
		ans=1;
		for(int l=1,r;l<=n;l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans=(ll)ans*ksm((ll)s[r]*ksm(s[l-1],mod-2)%mod,(ll)(n/l)*(m/l))%mod;
		}
		cout<<ans<<"
";
	} 
	return (0-0);
}
原文地址:https://www.cnblogs.com/aurora2004/p/12505840.html