Lucas 定理 / CRT 中国剩余定理

Lucas 定理

背结论最简单。

$$C_n^m\equiv C_{n\bmod p}^{m\bmod p}\times C_{\lfloor\frac{n}{p}\rfloor} ^{\lfloor\frac{m}{p} \rfloor}$$

这是一个很天才的结论，证明也是极具人类智慧。

我们知道 $C_n^m$ 实际上是 $(1+x)^n$ 中 $x^m$ 的系数。

以下所有运算在 $\bmod p$ 情况下计算。
我们知道，在 $(1+x)^p$ 中，我们可以拆成 $\sum\limits_{i=0}^pC_p^ix^i$，其中，$C_p^i\equiv [i=0\or i=p]$。不难证明。

因此 $(1+x)^p=(1+x^p)$。

得到这个结论，我们就去往降幂的方向去思考。令 $n=k_1p+b_1\;\;m=k_2p+b_2$ 。

因此 $(1+x)^n\equiv (1+x)^{k_1p+b_1}\equiv (1+x^p)^{k_1}\times (1+x)^{b_1}$

我们要得到 $x^m$ 这个项，观察式子，$k_2p$ 这个系数只有在 $(1+x^p)^{k_1}$ 才有，$b_2$ 也只有后面的一项才有

根据二项式定理，得到 $x^m=x^{k_2p+b_2}$，系数为 $C_{k_1}^{k_2}\times C_{b_1}^{b_2}$，证毕。

时间复杂度 $O(p+\log_p n)$。

点击查看代码

#include<bits/stdc++.h>
#define ll long long
#define N 100005
using namespace std;
int inv[N];
int fac[N],ifac[N];
int T;
int n,m,p;
void init()
{
	ifac[0]=fac[0]=1;inv[1]=1;
	const int mod=p;
	for(int i=2;i<mod;i++) inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod; 
	for(int i=1;i<mod;i++) fac[i]=1ll*fac[i-1]*i%mod,ifac[i]=1ll*ifac[i-1]*inv[i]%mod;
}
int C(int n,int m)
{
	if(m>n) return 0;
	return 1ll*fac[n]*ifac[m]%p*ifac[n-m]%p;
}
int Lucas(int n,int m){return m?1ll*Lucas(n/p,m/p)*C(n%p,m%p)%p:1;}
int main()
{
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d%d%d",&n,&m,&p);
		init();
		printf("%d\n",Lucas(n+m,m));
	}
	return 0;
}

CRT 中国剩余定理

求解同余方程组：

$$\begin{cases}x≡a_1\;\;(mod\;\;m_1)\\x≡a_2\;\;(mod\;\;m_2)\\ \cdots \cdots\\x≡a_n\;\;(mod\;\;m_n)\\\end{cases}$$

所有 $m$ 两两互素。

老祖宗的智慧。

在 $\bmod \prod_{i=1}^nm_i$ 意义下，只有唯一解。

我们可以从一个任意解得出唯一解。所以我们要构造一组特殊解。

构造一组 $y$：

$$\begin{cases}y_1≡a_1\;\;(mod\;\;m_1)\\y_2≡a_2\;\;(mod\;\;m_2)\\ \cdots \cdots\\y_n≡a_n\;\;(mod\;\;m_n)\\\end{cases}$$

我们想构造一个解：$x=\sum\limits_{i=1}^n y_i$

对于第一个数 $y_1$，从它的角度来看，它满足：

$$\begin{cases}y_1≡a_1\;\;(\bmod m_1)\\\sum\limits_{i=1}^n y_i\equiv a_1(\bmod m_1) \end{cases}$$

令 $M=\prod_{i=1}^nm_i$,$M_i=\dfrac{M}{m_i}$ 所以得到

$$\begin{cases}y_1≡a_1 \;\;(\bmod m_1)\\y_1\equiv 0\;\;(\bmod M_1) \end{cases}$$

考虑到 $M_1$ 很大，令 $y_1=kM_1$，得到 $kM_1\equiv a_1\;\;(\bmod m_1)$。求出一个 $k\equiv a_1M_1^{-1}\;(\bmod m_1)$，回带，得到 $y_1=a_1M_1M_1^{-1}(\bmod m_1)$。

于是我们得到了通项公式：

$$x\equiv\sum\limits_{i=1}^n a_iM_iM_i^{-1}\;\;\;(\bmod M)$$

时间复杂度 $O(n\log m)$。

Upd:
为什么要对 $M$ 取模？不会不合法吗？
考虑若有两个解 $x_1,x_2$，满足

$$\begin{cases}x≡a_1\;\;(mod\;\;m_1)\\x≡a_2\;\;(mod\;\;m_2)\\ \cdots \cdots\\x≡a_n\;\;(mod\;\;m_n)\\\end{cases}$$

则

$$\begin{cases}x_1-x_2\equiv 0\;\;(\bmod m_1)\\x_1-x_2\equiv 0\;\;(\bmod m_2)\\ \cdots\cdots\\ x_1-x_2\equiv 0\;\;(\bmod m_n)\\ \end{cases}$$

所以 $M\mid x_1-x_2$。即 $x_1\equiv x_2\;\;(\bmod M)$。 证毕。

点击查看代码

#include<bits/stdc++.h>
#define ll long long
#define N 15
using namespace std;
int n;
ll a[N],b[N];
ll sum=1,res;
void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0){x=1;y=0;return;}
	exgcd(b,a%b,x,y);
	ll tmp=x;
	x=y;
	y=tmp-a/b*y;
}
inline ll inv(ll x,ll mod)
{
	ll a,b;
	exgcd(x,mod,a,b);
	return (a%mod+mod)%mod;
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
		scanf("%lld%lld",&a[i],&b[i]),sum*=a[i];
	for(int i=1;i<=n;i++)
		res+=__int128(b[i])*inv(sum/a[i],a[i])%sum*(sum/a[i])%sum,res%=sum;
	printf("%lld",res);
	return 0;
}

ExCRT

其实和 CRT 关系不大。

求解同余方程组：

$$\begin{cases}x≡a_1\;\;(mod\;\;m_1)\\x≡a_2\;\;(mod\;\;m_2)\\ \cdots \cdots\\x≡a_n\;\;(mod\;\;m_n)\\\end{cases}$$

不保证所有 $m$ 两两互素。

考虑合并方程组。
我们看这两个方程组：

$$\begin{cases}x≡a_1\;\;(mod\;\;m_1)\\x≡a_2\;\;(mod\;\;m_2)\end{cases}$$

相当于求解 $x=k_1m_1+a_1=k_2m_2+a_2$

移项，得到 $k_1m_1-k_2m_2=a_2-a_1$。
我们可以使用 Exgcd 得到一组 $(k_1,k_2)$ 或者判断无解。得到解 $x'$。

由上面的证明，我们知道，在 $\bmod lcm(m_1,m_2)$ 下有唯一解 $x$。这说明，对于以后所有的方程组，我们只需满足 $x\equiv x'(\bmod lcm(m_1,m_2))$ 即可，这样就完成了方程组合并。合并到最后就是答案。

注意，同余中最好不要出现负数，否则会出现奇奇怪怪的错误，特别是在 Exgcd 中，所以要对解 $x'$ 进行修正。

事件复杂度同样是 $O(n\log m)$。

点击查看代码

#include<bits/stdc++.h>
#define ll long long
#define N 100005
#define LL __int128
using namespace std;
int n;
ll a[N],b[N];
LL x,y,d;
void exgcd(LL a,LL b,LL &x,LL &y)
{
	if(b==0){x=1,y=0,d=a;return;}
	exgcd(b,a%b,x,y);
	ll tmp=x;
	x=y;
	y=tmp-a/b*y;
}
void EXCRT()
{
	LL A=a[1],m=b[1];
	for(int i=2;i<=n;i++)
	{
		LL k1,k2;
		exgcd(m,b[i],k1,k2);
		if(LL(a[i]-A)%d) return;
		k1*=LL(a[i]-A)/d;
		x=k1*m+A;
		m=m/d*b[i];
		A=(x%m+m)%m;
	}
	ll res=A;
	printf("%lld\n",res);
}
int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
		scanf("%lld%lld",&b[i],&a[i]);
	EXCRT();
	return 0;
}

ExLucas

要说 ExCRT 还可以拆拆式子回到 CRT，那么 ExLucas 和 Lucas 定理没有半毛钱关系。你只需要会 Exgcd 和 CRT 即可。
给出任意整数 $A,B,p$，求：

$$C_A^B\bmod p$$

Part 1.

对 $p$ 进行质因数分解，得到 $p_1^{k_1}p_2^{k_2}\cdots p_n^{k_n}$

化为同余方程组，求：

$$\begin{cases}x≡C_A^B\;\;(\bmod\;\;p_1^{k_1})\\x≡C_A^B\;\;(\bmod\;\;p_2^{k_2})\\\cdots \cdots\\\end{cases}$$

最后使用 CRT 合并。

Part 2.

如何求解 $x≡C_A^B\;\;(\bmod\;\;p_1^{k_1})$ ？
显然的，$C_n^m\equiv \dfrac{n!}{m!(n-m)!} \bmod p^k$

问题来了，$n!$ 这个东西有 $p$ 这样一类的数不存在逆元。

我们自然而然就是想把所有的 $p$ 提出来。
即 $\dfrac{\dfrac{n!}{p^x}}{\dfrac{m!}{p^y}\dfrac{(n-m)!}{p^z}}\times p^{x-y-z} \bmod p^k$

Part 3.1

如何求解 $n!\bmod p^k$ ？
我们定义函数 $f(n)=n!\bmod p^k$
把所有带 $p$ 的系数提出来，即 $p$,$2p$,$3p\cdots$
所以

$$n!\;=\prod_{i=1}^{\lfloor\frac{n}{p}\rfloor}i\times \prod_{i=1,i\not\equiv 0}^ni$$

前面一部分容易发现可以递归解决，后面一部分肯定不能硬算。

我们发现后面一部分是有循环节的。$p^k$ 以后就开始重复了。
这样我们只需算到 $p^k$ 就行了。

$$f(n)=n!\;=f(\lfloor\frac{n}{p}\rfloor)\times (\prod_{i=1,i\not\equiv 0}^{p^k}i)^{\lfloor\frac{n}{p^k}\rfloor}\times \prod_{i=\lfloor\frac{n}{p^k}\rfloor\times p^k}^ni$$

递归解决，时间复杂度 $O(p\log_p n)$。
边界为 $f(0)=1$。

Part 3.2

对于 $n$，如何求解 $p^x$？
其实就是把 $x$ 转为 $p$ 进制后每一位乘上对应指数。
最广泛流传的就是转成递归式，也就是 $g(n)=g(\lfloor\dfrac{n}{p}\rfloor)+\lfloor\dfrac{n}{p}\rfloor$

Part 4.

这还用说，把上面的合并了，分数部分 exgcd 求逆元解，右边使用快速幂，常理来说你也可以特判返回 $0$，最后使用 CRT 合并即可。

conclude

我认为 exLucas 完全不能算是一个算法，更多的算一道数论练习题，给我们提供拆模数用 CRT 合并的思想，这个思想会非常有用。

#include<bits/stdc++.h>
#define ll long long
#define N 1000005
using namespace std;
ll n,m,P;
int L;
ll p[30],pk[30];
ll res[30]; 
void Case1()//分解质因数 
{
	for(int i=2;i<=sqrt(P);i++)
	{
		if(P%i) continue;
		p[++L]=i,pk[L]=1;
		while(!(P%i)) P/=i,pk[L]*=i;
	}
	if(P!=1) p[++L]=P,pk[L]=P;
}
ll Qpow(ll a,ll x,ll p)
{
	ll res=1;
	while(x)
	{
		if(x&1) res=res*a%p;
		a=a*a%p;
		x>>=1;
	}
	return res;
}
void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0){x=1;y=0;return;}
	exgcd(b,a%b,x,y);
	ll tmp=x;
	x=y;
	y=tmp-a/b*y;
}
ll inv(ll a,ll p)
{
	ll x,y;
	exgcd(a,p,x,y);
	return (x%p+p)%p;
}
ll g(ll x,ll p){return x?g(x/p,p)+x/p:0;}
ll tmp;
ll f(ll x,ll p,ll pk)
{
	if(x==0) return 1;
	ll cur=1;
	for(int i=1;i<=x%pk;i++)
		if(i%p) cur=cur*i%pk;
	return cur*Qpow(tmp,x/pk,pk)%pk*f(x/p,p,pk)%pk;
}
void calc(ll &res,ll p,ll pk)
{
	tmp=1;
	for(int i=1;i<=pk;i++)
		if(i%p) tmp=tmp*i%pk;
	res=f(n,p,pk)*inv(f(m,p,pk)*f(n-m,p,pk)%pk,pk)%pk*Qpow(p,g(n,p)-g(m,p)-g(n-m,p),pk)%pk;
	return;
}
void Case2()// 对每个求答案 
{
	for(int i=1;i<=L;i++)
		calc(res[i],p[i],pk[i]);
}
ll ans,M;
void Case3()// CRT merges.
{
	for(int i=1;i<=L;i++)
		ans=(ans+res[i]*M/pk[i]%M*inv(M/pk[i],pk[i])%M)%M;
}
int main()
{
	scanf("%lld%lld%lld",&n,&m,&P);M=P;
	Case1();
	Case2();
	Case3();
	printf("%lld",ans);
	return 0;
}