gcdlcm[组合数学]

gcdlcm[组合数学]：

描述：
问有多少个 k 元组(a1,a2,a3,...,ak)(ai>=1)满足 gcd(a1,a2,...,ak)=D，
并且 lcm(a1,a2,a3,...,ak)=L。
输入：
第一行有三个整数 k，D,L 分别为题面所描述的。
输出：
输出答案 mod(1e9+7)
输入输出样例：
gcdlcm.in gcdlcm.out
2 1 6 4
数据范围：
对于 30%数据 2<=k<=100,D<=100,L<=100
对于 100%数据 2<=k<=1000000000 D<=1000000000 L<=1000000000

数学弱成渣了，毫无办法- -

---------------------------------------------------------------------------------------------------------------------------------------------------------

首先，D|L.  否则不合法，输出0.

然后，将D,L分解为D=p1^b1*p2^b2*p3^b3*```*pq^bq  ,L=p1^c1*p2^c2*p3^c3*````*pq^cq

然后选k个数出来，使得k个数分解出来的cnt[i][pj]满足：

cj>=(cnt[i][pj])>=bj        ,求方案数

考虑对于每个因子pi   我们要在[bi,ci]里面选k个，并且满足min(pij)=bi,max(pij)=ci.得到一个方案gi

然后把所有gi乘起来！

问题就转化为求gi了！

这里要用到数学的一个容斥：(数学太弱的我表示不会想到)

合法条件在pi中一定至少有一个选了bi,ci个.  不合法的就是 没有选bi个+没有选ci个-(没有选bi&&ci个)

我们考虑用总的减去不合法的

总的就是在[bi,ci]中选k个，为(ci-bi+1)^k

不选bi或者不选ci都是(ci-bi)^k

gi=(ci-bi+1)^k-2*(ci-bi)^k+(ci-bi-1)^k

这个用快速幂求一下就好了

复杂度其实不算很大，不过需要考虑分解出来的pq可能会非常大(如果>sqrt(n)就单独记录，因为数组开不下)

代码：

#include<cstdio>
#include<iostream>
#include<queue>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#define LL long long
using namespace std;
int k;
int D,L;
const LL M=1e9+7;
/*gcd(a1,a2,a3---an)=D=p1^b1 * p2^b2 * p3^b3
//lcm(a2,a2,a3---an)=L=p1^c1 * p2^c2 * p3^c3                  pk
//选法就是考虑每个因子pi 从[b1,c1]中选k个 得到方案数 gi   ans=π gi
                                                            //i=1 
对于每个gi,我们考虑从[b1,c1]中选k个，方案数 (c1-b1+1)^k
合法状态一定有至少一个b1,c1
那么不合法状态,  不选b1  ---(b1,c1]选k个，方案数  (c1-b1)^k
                 不选c1  ---[b1,c1)选k个，方案数  (c1-b1)^k
多减去了(b1,c1)的方案数  (c1-b1-1)^k
gi=(c1-b1+1)^k-2*(c1-b1)^k+(c1-b1-1)^k; 
复杂度 log max(D,L) * 
*/ 
int b[1000000],c[1000000];
int xx=-1;
int ml=0;
void fenjie(int x,int *a)//把x拆分到a数组 
{
	memset(a,0,sizeof(a));
	int t=(int)sqrt(1.0*x);
	for(int i=2;i<=t;i++)
	{
		if(x%i==0)
		{
			ml=max(ml,i);
			while(x%i==0)
			{
				x/=i;
				a[i]++;
			}
		}
	}
	if(x>1)
	{
		if(x>1000000)xx=x;
		else 
		{
			ml=max(ml,x);
			a[x]++;
		}
	}
	//
}
LL qk(LL x,LL y)//快速计算x^y mod M 
{
	LL res=1;
	while(y)
	{
		if(y%2==1)res=((res%M)*(x%M))%M;
		y/=2;
		x=((x%M)*(x%M))%M;
	}
	return res;
}
LL ans;
int main()
{
	freopen("gcdlcm.in","r",stdin);
	freopen("gcdlcm.out","w",stdout);
	scanf("%d%d%d",&k,&D,&L);
	if(L%D)printf("0\n");
	else
	{
		fenjie(D,b);
	    fenjie(L,c);
	    LL g=0;
	    ans=1;
	    for(int i=2;i<=ml;i++)
	 {
		if(!c[i])continue;
		LL d=(LL)c[i]-(LL)b[i];
		if(!d)g=1;
		else g=qk(d+1,k)%M-2*qk(d,k)%M+qk(d-1,k);
		g%=M;
		ans=((ans%M)*g)%M;
	 }
	    if(xx!=-1)
	    {
	    	if(D%xx!=0)
	    	{
	    		LL d=1;
	    	    g=qk(d+1,k)%M-2*qk(d,k)%M+qk(d-1,k);
	    	    ans=((ans%M)*(g%M))%M;
	    	}	    	
	    }
	    ans=(ans+M)%M;
	    printf("%I64d\n",ans);
	}
	
	return 0;
}