Codeforces Gym 101138 G. LCM-er

Description

在 ([a,b]) 之间选择 (n) 个数 (可以重复) ,使这 (n) 个数的最小公倍数能被 (x) 整除,对 (10^9+7) 取膜.

(1leqslant n leqslant 100,1leqslant a,b,x leqslant 10^9) 

Sol

容斥原理+状态压缩.

这是容斥原理非常好的一道题啊QAQ.

转化成补集,有多少不能被 (x) 整除的序列.

首先质因数分解 (x) 因为最小的 (9) 个质数的积已经超过 (10^9) 所以 (x) 的质因子个数不会超过 (9)

那么 (x=p_1^{k_1}p_2^{k_2}p_3^{k_3}...p_m^{k_m})

另集合 (S={p_1^{k_1},p_2^{k_2},p_3^{k_3},...,p_m^{k_m}})

(f[S']) 表示 ([a,b]) 中不能被集合 (S') 的任意一个元素整除的个数.

这个可以通过枚举子集和容斥原理来实现.

那么答案还需要一下容斥就是 (ans=prod (-1)^{|S'|}egin{pmatrix}f[S']+n-1\ nend{pmatrix})

组合数可以最后乘上 (n) 的逆元.

枚举集合子集的复杂度是 (O(2^{|S|}))

枚举子集的子集复杂度是 (O(3^{|S|})) ,每个元素分为 (3) 种情况:不在子集中,在子集中不在子集的子集中,在子集的子集中.

所以总复杂度 (O(m3^m+n2^m)) .

Code

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<iostream>
using namespace std;

#define debug(a) cout<<#a<<"="<<a<<" "
typedef long long LL;
const int N = 1<<11;
const int M = 15;
const int W = 100005;
const LL p = 1000000007;

LL n,a,b,x,ans;
LL f[N];
LL pr[W],isp[W],s[M],c,cnt;
int pow2[M],lim;

void out(int S){ for(int i=0;i<c;i++) if(S & pow2[i]) putchar('1');else putchar('0');putchar('
'); }
LL gcd(LL a,LL b){ return a>b?(!b?a:gcd(a%b,b)):(!a?b:gcd(b%a,a)); }
LL Pow(LL a,LL b,LL res=1){ for(;b;b>>=1,a=a*a%p) if(b&1) res=res*a%p;return res; }
inline int count(int S,int res=0){ for(int i=0;i<c;i++) if(S & pow2[i]) res++;return res; }
void Pre(LL N){
	for(int i=2;i<=N;i++){
		if(!isp[i]) pr[++cnt]=i;
		for(int j=1;j<=cnt && (LL)pr[j]*i<=N;j++){
			isp[i*pr[j]]=1;
			if(i%pr[j]==0) break;
		}
	}
}
void work(LL x){
	for(int i=1;i<=cnt;i++) if(x%pr[i]==0){
		s[c++]=pr[i],x/=pr[i];
		while(x%pr[i]==0) s[c-1]*=pr[i],x/=pr[i];
	}if(x>1) s[c++]=x;
}
LL calc(int S){
	LL res=1;
	for(int i=0;i<c;i++) if(S & pow2[i]) res=res/gcd(res,s[i])*s[i];
	res=b/res-a/res;
	return res;
}
LL Sum(LL m){
	LL res=1;
	for(int i=1;i<=n;i++) res=res*(m-i+1)%p;
	return res;
}
int main(){
//	freopen("in.in","r",stdin);
	cin>>n>>a>>b>>x;a--;
	Pre(sqrt(x)+1);work(x);
	pow2[0]=1;for(int i=1;i<M;i++) pow2[i]=pow2[i-1]<<1;
	lim=1<<c;
	
//	for(int i=1;i<=cnt;i++) cout<<pr[i]<<" ";cout<<endl;
//	debug(c)<<endl;
//	for(int i=0;i<c;i++) cout<<s[i]<<" ";cout<<endl;
//	debug(gcd(2,3)),debug(gcd(2,6)),debug(gcd(3,6)),debug(gcd(6,6))<<endl;
	
	for(int S=0;S<lim;S++){
		f[S]=b-a;
		if(count(S)&1) f[S]=f[S]-calc(S);
		else if(count(S)) f[S]=f[S]+calc(S);
		for(int T=S&(S-1);T;T=(T-1)&S)
			if(count(T)&1) f[S]=f[S]-calc(T);
			else f[S]=f[S]+calc(T);
	}
//	cout<<lim<<endl;
//	for(int i=0;i<lim;i++) cout<<f[i]<<" ";cout<<endl;
	
	for(int S=0;S<lim;S++){
		if(count(S)&1) ans=(ans-Sum(f[S]+n-1)+p)%p;
		else ans=(ans+Sum(f[S]+n-1))%p;
	}
	LL tmp=1;
	for(int i=2;i<=n;i++) tmp=tmp*(LL)i%p;
	cout<<ans*Pow(tmp,p-2)%p<<endl;
	return 0;
}