各种友(e)善(xin)数论总集，从入门到绝望5---组合数取模？(扩展)卢卡斯套上直接摁在地上锤

目录

有些只是我个人看着像数论，其实不一定是数论

最近数论越来越短，只是单纯因为这个专题竟然只能选三个，更好的体验果然还是要上博客园吧。

参考资料

说是资料，就是抄袭吧不过还是会做一些取舍，把一些我当初也不是很懂的东西做一些注解，可以说是加强版吧。

普通卢卡斯推导：https://www.luogu.com.cn/blog/28007/lucas

扩展卢卡斯推导：

普通卢卡斯定理

例题

万能的洛谷

定义

没有定义的话这篇证明还是挺难看的。

话说只是为了更好的抄题解把。。。

一下定义仅适用于普通卢卡斯。

设\(C(n,m,p)\)为\(C_{n}^{m}\mod p\)，如果右边的定义你都不会的话，你学什么卢卡斯，先去把组合数学过一过吧，所以本文章几乎所有内容均在同余的意义下。

同时注意，\(C(n,m,p)(n>m\)和\(n<1)\)在这篇证明中都算是\(0\)，即无意义，默认定理中的组合数都是有效的，证明中可能出现的无意义组合数请在代码中特判为\(0\)。

做法

许多人也许会说，为什么不能爆搞逆元？

注意！有时候答案不一定是\(0\)，但是由于\(p<n,m\)，所以在计算\(C(n,m,p)\)中分母的\(m!\)的逆元就会无解。

但是其实有个很有意思的事情，就是许多人都会发现其实是可以把分子分母的\(p\)给全部提取出来的，但是最后都是直接去看题解了，这种其实也是一种做法，就是卢卡斯定理！！！

先说结论吧：

对于\(C(n,m,p)\)

设\(n\)可以拆分为：\(a_{i}p^{i}+a_{i-1}p^{i-1}+...+a_{0}\)
\(m\)也是，不过把\(a\)改成\(b\)。（\(b_{i}\)为\(0\)也要补上）

那么\(C(n,m,p)=C(a_{i},b_{i},p)C(a_{i-1},b_{i-1},p)...\)

但是其实我最好奇的是为什么只要有一个上面的大于下面的就全部都为0了，但又感觉证明没什么问题。。。

定理1：\(C_{p}^{i}=\frac{p}{i}C_{p-1}^{i-1}\)

\(C_{p}^{i}=\frac{p!}{i!(p-1)!}=\frac{(p-1)!}{(i-1)!(p-i)!}\frac{p}{i}=\frac{p}{i}C_{p-1}^{i-1}\)

定理2：\((1+x)^p≡1+x^{p}\mod p\)

这个在我未来要写的斐波那契循环节中也是特别重要的。。。

\((1+x)^p≡1+C(p,1,p)x+C(p,2,p)x^2+...+x^p\mod p\)

那么关键就是：\(C_p^i(0<i<p)\)到底能不能被\(p\)整除？

首先我们观察这个式子：\(\frac{p}{i}*C_{p-1}^{i-1}\)他又大又圆

我们可以发现因为\(p\)是质数，所以不会被除掉，所以就可以整除啦。

所以

\((1+x)^p≡1+x^p\mod p\)

得证。

定理三：证明开头说的那段话。简单暴力的偷懒

设：\(n=ap+b,m=ip+j(0≤b,j<p)\)

我们现在只需要证明\(C(n,m,p)≡C(a,i,p)*C(b,j,p)\)，然后递归一下就可以得出结论。

模仿抄袭题解证明\(C(n,m,p)≡C(a,i,p)*C(b,j,p)\)：

\((1+x)^{n}=(1+x)^{ap}(1+x)^{b}\)

展开\((1+x)^{ap}\)

\((1+x)^{ap}≡((1+x)^p)^a≡(1+x^p)^a\mod p\)

观察\((1+x)^n≡(1+x)^{ap}(1+x)^{b}≡(1+x^p)^a(1+x)^b\)

观察第\(x^m\)项的系数（其实\(x^m\)的系数在模意义下就是\(C(n,m,p)\)）：\(C(n,m,p)x^m≡C(a,i,p)x^{pi}*C(b,j,p)x^j\)

\(C(n,m,p)x^m≡C(a,i,p)*C(b,j,p)x^m\)

所以：\(C(n,m,p)≡C(n\%p,m\%p,p)*C(n/p,m/p,p)\)

得证。

那么其实代码里面也是拿这个去递归的。

这个证明我觉得特别优美，为什么，因为在拆分成\((1+x)^n≡(1+x)^{ap}(1+x)^{b}≡(1+x^p)^a(1+x)^b\)，\(x^m\)的系数表示是两个组合数相乘，没有加号的，也就是利用\(ip+j\)的表现形式。

当然，卢卡斯是\(O(log_{p}^{n}+p)\)级别的时间复杂度，还是很快的

#include<cstdio>
#include<cstring>
#define  N  110000
using  namespace  std;
typedef  long  long  LL;
LL  inv[N],jie[N]/**/;
LL  n,m,p;
inline  LL  C(LL  x,LL  y)//计算组合数
{
	if(x>y)return  0;
	return  (jie[y]*inv[jie[x]]*inv[jie[y-x]])%p;
}
inline  LL  Lucas(LL  x,LL  y)//卢卡斯递归过程
{
	if(!y)return  1;//y=0的别忘了特判
	return  (Lucas(x/p,y/p)*C(x%p,y%p))%p;
}
int  main()
{
	int  T;scanf("%d",&T);
	for(int  i=1;i<=T;i++)
	{
		scanf("%lld%lld%lld",&n,&m,&p);
		inv[0]=inv[1]=1;for(int  i=2;i<p;i++)inv[i]=(p-p/i)*inv[p%i]%p;
		jie[1]=1;for(int  i=2;i<=p;i++)jie[i]=jie[i-1]*i%p;//预处理逆元和阶乘
		printf("%lld\n",Lucas(m,n+m));
	}
	return  0;
}

扩展卢卡斯

例题

又是万能的洛谷

做法

这个做法其实也是很暴力的。

首先，\(p\)不是质数，看起来很不友善，那我们就把他们化成质数吗。

不过其实也讲的不是很严谨，应该说是质数次幂。（为了方便后面，其实也有因为两个相同的模数无法合并的问题。）

把\(p\)拆分成\(p_{1}^{a_{1}}p_{2}^{a_{2}}...p_{i}^{a_{i}}(p_{j}\)都是质数，\(1≤j≤i)\)

分别求出在模\(p_{j}^{a_{j}}\)意义下的值，然后中国剩余定理合并即可。

我们根据之前在卢卡斯之前提到过的。

在这里，我们假设我们现在模的是\(p^k\)，那么我们就把所有的\(p\)移出来，然后就可以求逆元了。

前面提到的方便就是我们只需要把所有的\(p\)移出来，然后剩下的数字绝对与\(p^k\)互质。

So，我们把式子化成了这个样子（注意，后面的除以\(p\)的几次方表示全部提出来）：\(C(n,m,p^k)≡\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}*p^{x-y-z}\mod p^k\)

对于为什么\(x-y-z\)这个玩意没可能小于\(0\)的话，ZFY奆佬告诉我们，前面那个分数的分子与\(p\)互质，后面又全是\(p\)，所以如果呈现的形式是\(\frac{1}{p}\)的话，就会变成分数，但是这玩意结果是整数，所以就是大于等于\(0\)的。(大雾

在这里插入图片描述

然后关键就是怎么计算\(\frac{n!}{p^x}\)以及快速返回\(x\)。（分母也是这样干，然后逆元。）

是\(p\)的倍数的个数其实是很容易计算的，就是\(\frac{n}{p}\)个。

为了快速计算非\(p\)倍数的乘法，这里我们进行分组：

对于\(1*2*3*...*n\)

把他们\(p^k\)一组一组的分。

然后对于长度为\(p^k\)的一组（也就是满的组），在\(\mod p^k\)的以一下，其实他们的值是一样的，都是\(1*...*(p^k-1)\)（解释一下，模运算里面有个很神奇的性质，\(a*b*c*...≡(a\mod p)*b*c*...\mod p\)，然后就可以得到这个性质了）。

剩下的那一组就只能单独计算了。

那么其实仔细想想就是：

\((\prod\limits_{i=1,(i\mod p)>0}^{p^k}i)^{\frac{n}{p^k}}*(\prod\limits_{i=\left \lfloor \frac{n}{p^k} \right \rfloor * p^k+1,(i\mod p)>0)}^{n}i)\)

但是其实由于提出了\(p\)，我们前面还跟了一些东西，就是：

\(p^{\left \lfloor \frac{n}{p} \right \rfloor}\left \lfloor \frac{n}{p} \right \rfloor!\)

那么其实我们前面的东西可以直接返回去，但是后面怎么又有阶乘？

那就递归求解吗，但是也要把前面的\(p\)次方传上去。

So，分析时间复杂度：

首先，对于\(p^k\)，中间的求阶乘和快速幂的部分，时间复杂度为：\(O(p^k*log_2\frac{n}{p^k}\log_{p}{n})\)。

由于后面的递归总和最多是等于第一次求阶乘的，所以最多带个\(2\)的常数。

对于\(p\)而言，时间复杂度为所有的\(p_{j}^{a_{j}}\)的单词时间复杂度相加，那么记为\(A\)。

我们采用扩展欧几里得求逆元，用CRT合并，CRT的时间复杂度为\(O(\log^2{p})\)。

所以总的时间复杂度为：\(O(A+\log^2{p})\)。其实可以接近于认为是\(O(p\log{p})\)（看某位大佬说的）。

这种东西的时间复杂度就图一乐，太玄学了。

#include<cstdio>
#include<cstring>
using  namespace  std;
typedef  long  long  LL;
inline  LL  ksm(LL  x,LL  k,LL  mod)
{
	LL  now=1;
	while(k)
	{
		if(k&1)now=(now*x)%mod;
		x=(x*x)%mod;
		k>>=1;
	}
	return  now;
}
inline  LL  jc(LL  l,LL  r,LL  x,LL  mod)//求阶乘
{
	LL  ans=1;
	for(LL  i=l;i<=r;i++)
	{
		if(i%x!=0)ans=ans*(i%mod)%mod;
	}
	return  ans;
}
inline  void  exgcd(LL  a,LL  b,LL  &x,LL  &y)//求逆元
{
	if(!a)return  (void)(x=0,y=1);
	exgcd(b%a,a,x,y);
	LL  tmp=x;x=y-(b/a)*x;y=tmp;
}
LL  gcd(LL  x,LL  y)//这个函数好像一直没有用到。
{
	if(!x)return  y;
	return  gcd(y%x,x);
}
inline  LL  inv(LL  x,LL  mod)//求逆元 
{
	LL  a,b;
	exgcd(x,mod,a,b);
	/*b=-b,因为mod其实是负数形式*/
	a=(a%mod+mod)%mod;
	return  a;
}//返回逆元
LL  calc(LL  n,LL  p,LL  mod,LL  &now)//表示计算(x!/P^now)%P^y
{
	if(!n)return  1;
	now+=n/p;//加上P的次幂 
	LL  ans=calc(n/p,p,mod,now);//先得到阶乘；
	return  (ans*ksm(jc(1,mod,p,mod),n/mod,mod)%mod)*jc((n/mod)*mod+1,n,p,mod)%mod;
}
inline  LL  C(LL  x,LL  y,LL  p,LL  mod)//C_x^y
{
	LL  n=0,m=0;LL  ans=calc(x,p,mod,m);n+=m;m=0;
	ans*=inv(calc(y,p,mod,m),mod);ans%=mod;n-=m;m=0;
	ans*=inv(calc(x-y,p,mod,m),mod);ans%=mod;n-=m;m=0;
	return  ans*ksm(p,n,mod)%mod;
}
inline  LL  fen(LL  x,LL  a1,LL  a2)//分解质因数
{
	LL  y=2;
	LL  nmo=1,ans=0;
	while(y<=x)
	{
		if(x%y==0)
		{
			LL  z=1;
			while(x%y==0)z*=y,x/=y;
            //z为模数，y为质数，也是ZFY奆佬教给我的一个分解质因数的方法
			LL  now=C(a1,a2,y,z);//组合数
			LL  zong=z*nmo;
			ans=(ans*z*inv(z,nmo)+now*nmo*inv(nmo,z))%zong;
			nmo=zong;
            //两个同余式合并，因为彼此互质，所以可以用比较简单的CRT。
		}
		y++;
	}
	return  ans;
}
int  main()
{
//	freopen("std.in","r",stdin);
//	freopen("vio.out","w",stdout);
	LL  n,m,p;scanf("%lld%lld%lld",&n,&m,&p);
	printf("%lld\n",fen(p,n,m));
	return  0;
}

一种新式的奇技淫巧

最近刚刚学习的一种神仙的奇技淫巧。

就是叫你计算\(C_{x}^y\mod p^a\)

\(p\)一般比较小，但是\(p^a\)可以在long long范围内，\(x,y\)也是。

这种奇技淫巧放在扩展卢卡斯定理中就可以解决质因数比较小的long long范围的质数了。

那么到底是个什么神仙东西呢？

例题

没有链接，就是叫你求\(C_{x}^{y}\mod 5^{23}\)

这个是某个恶心集训队互测的一种恶心的子问题。那道题DP很好，非要搞一个组合数取模

思路

这种神奇的思想依旧是按照算出除\(5\)以外的阶乘与\(5\)的个数。

但是特别神奇的是他在算阶乘的时候采用的是树套树套树套树拆括号。

没错就是小学的拆括号，但是这个拆括号极其有技术含量。

我们选择一个\(2*5\)（后面就知道为什么要乘\(2\)了）的倍数作为基数。

题解里使用的是\(10\)（其实是10的倍数也可以，但是取最小的一般常数是最小的），我们现在要计算\(n!\)，那么我们设\(k=\left \lfloor \frac{n}{10} \right \rfloor *10\)，\(k+1\)~\(n\)的我们暴力算，也就在\(10\)以内。但是关键是\(k!\)呢？

一种也同样明显隐藏的思路就是把\(k\)分成两份。

\(1-(k/2)\)和\((k/2+1)-k\)，然后类似分治的思想。

这也就解释了为什么基数要选择偶数了。

假设我们知道了\(1-(k/2)\)，但是怎么求出\((k/2+1)-k\)，我们列出括号（设\(o=\frac{k}{2}\)，容易知道：\(5|o\)（\(o\)被\(5\)整除））：\((1+o)*(2+o)*(3+o)*...*(o+o)\)，由于是\(5^{23}\)次幂，在选的时候最多只能选\(22\)个\(o\)

仔细观察，发现\(o!\)我们是知道的，但是\(o^i\)的系数我们是不知道的。

注意到我说了\(o^i\)的系数，明白人都发现那不是生成函数（自行学习吧）吗！！！

也就是说对于\(o^i\)的系数我们就需要知道在\(1-o\)中选\(o-i\)个不等的数字的乘积，而\(0≤i≤22\)。

所以我们只要接着构造出\(1-o\)的生成函数就行了，也就是\((1+x)*(2+x)*(3+x)*...*(o+x)\)，所以问题就成了通过\(1-o\)的生成函数的值求\((o+1)-k\)的生成函数的值，而\((o+1)-k\)的生成函数为\((o+1+x)*(o+2+x)*...*(k+x)\)。

等会，生成函数？难道是FFT，不不不，由于\(0≤i≤22\)，所以我们只要记录第\(0-22\)项，所以直接暴力\(n^2\)多项式乘法，不就可以了吗。

所以我们只需要对于\(x^i\)的系数乘以\(o^i\)加到\((o+1)-k\)的常数项中（即阶乘值），等会，那其他的\((o+1)-k\)其他的项你不理他了？还是要理会的，对于\(x^i\)的系数我们要转化到\((o+1)-k\)的\(x^j(i≥j)\)的系数，我们就需要在\(i\)个\(1\)中挑选\(i-j\)个变成\(o\)加入到\(x^j\)的系数中
，所以不仅要乘\(o^{i-j}\)，还要乘\(C_{i}^{i-j}\)。

关于为什么挑选了\(o\)就是非\(x\)项？

把式子展开不难发现，\(o\)其实是在括号内的常数项的，只不过把\(x\)的\(1\)系数变成了常数项的系数的一部分（这个过程可以说是把多个值融合成结果）。

所以就成了。

题解中可能为了减少特判，还直接把\(1-10000\)的生成函数的值全部求了！！！

当然\(5\)的个数就在递归中慢慢算了，注意，上面的全部过程都是把\(5\)的倍数拿了出来，所以请自行递归算\(\left \lfloor \frac{n}{5} \right \rfloor!\)。

代码在这里：

LL  mod=(LL)11920928955078125;
inline  LL  mul(LL  x,LL  y){return  (x*y-(LL)((long  double)x*y/mod+1e-10)*mod);}
inline  LL  pow(LL  x,LL  k)//求逆元专用 
{
	LL  ans=1;
	while(k)
	{
		if(k&1)ans=mul(x,ans);
		x=mul(x,x);k>>=1;
	}
	return  ans;
}
namespace  Big_num//大整数组合数 
{
	LL  pw[S]={1},C[S][S],K/*表示选出几个集合*/;
	struct  poly//没错，这个向量，他又来了
	{
		LL  a[S];
		poly(LL  x=0,LL  y=0){memset(a,0,sizeof(a));a[0]=x;a[1]=y;}
		void  init(LL  k)
		{
			static  LL  ret[S];memset(ret,0,sizeof(ret)); 
			for(int  i=1;i<S;i++)pw[i]=mul(pw[i-1],k);
			for(int  i=0;i<S;i++)
			{
				for(int  j=0;j<=i;j++)ret[j]=(ret[j]+mul(a[i],mul(pw[i-j],C[i][j/*从不选的里面挑出几个*/])))%mod;
			}
			memcpy(a,ret,sizeof(ret));
		}
		poly  operator*(poly  x)
		{
			poly  z;
			for(int  i=0;i<S;i++)
			{
				if(x.a[i])//打不打无所谓
				{
					for(int  k=i;k<S;k++)z.a[k]=(z.a[k]+mul(x.a[i],a[k-i]))%mod;//就是这里的推导
				}
			}
			return  z;
		}
	}P[10005];
	//-----------poly
	poly  facpoly(LL  n)//求阶乘 
	{
		if(n<=10000)return  P[n];
		LL  k=n/10*10;
		poly  t1=facpoly(k>>1),t2=t1;
		t2.init(k>>1);//用一个生成函数求另外一个生成函数的过程
		t1=t1*t2;
		for(LL  i=k+1;i<=n;i++)
		{
			if(i%5!=0)t1=t1*poly(i,1);
		}
		return  t1;
	}
	LLp  solve(LL  n)//递归
	{
		LLp  ret=make_pair(facpoly(n).a[0],n/5);
		if(n>=5)
		{
			LLp  tmp=solve(n/5);
			ret.first=mul(ret.first,tmp.first);
			ret.second+=tmp.second;
		}
		return  ret;
	}
	LL  Combk(LL  n)
	{
		if(n<K)return  0;
		LLp  f1=solve(n),f2=solve(K),f3=solve(n-K);
		f1.second-=f2.second+f3.second;
		return  mul(mul(f1.first,pow(mul(f2.first,f3.first),mod/5*4-1)),pow(5,f1.second));
	}
	void  Init()
	{
		C[0][0]=1;
		for(int  i=1;i<S;i++)
		{
			C[i][0]=1;
			for(int  j=1;j<=i;j++)
			{
				C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod;
			}
		}
		P[0]=poly(1,0);
		for(int  i=1;i<=10000;i++)
		{
			if(i%5!=0)P[i]=P[i-1]*poly(i,1);
			else  P[i]=P[i-1];
		}
	}
};

然后这个玩意再拿扩展中国剩余定理合并，就可以处理\(p\)在long long范围的情况了，不过最大的质因子要小一点。

但对于\(p^a\)而言（\(p\)是质数），且基数选择为\(2p\)，时间复杂度为：\(O((a^2+p)\log^2{n})\)，这个复杂度\(p\)大了受不了的（虽然\(p\)大了其中一个\(log\)几乎就是常数了）。

小结

好像没什么好说的。

posted @ 2020-01-11 14:16  敌敌畏58  阅读(280)  评论(0)    收藏  举报