Lucas 定理 / CRT 中国剩余定理
Lucas 定理
背结论最简单。
这是一个很天才的结论，证明也是极具人类智慧。
我们知道 \(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 中国剩余定理
求解同余方程组：
所有 \(m\) 两两互素。
老祖宗的智慧。
在 \(\bmod \prod_{i=1}^nm_i\) 意义下，只有唯一解。
我们可以从一个任意解得出唯一解。所以我们要构造一组特殊解。
构造一组 \(y\)：
我们想构造一个解：\(x=\sum\limits_{i=1}^n y_i\)
对于第一个数 \(y_1\)，从它的角度来看，它满足：
令 \(M=\prod_{i=1}^nm_i\),\(M_i=\dfrac{M}{m_i}\) 所以得到
考虑到 \(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)\)。
于是我们得到了通项公式：
时间复杂度 \(O(n\log m)\)。
Upd:
为什么要对 \(M\) 取模？不会不合法吗？
考虑若有两个解 \(x_1,x_2\)，满足
则
所以 \(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 关系不大。
求解同余方程组：
不保证所有 \(m\) 两两互素。
考虑合并方程组。
我们看这两个方程组：
相当于求解 \(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\)，求：
Part 1.
对 \(p\) 进行质因数分解，得到 \(p_1^{k_1}p_2^{k_2}\cdots p_n^{k_n}\)
化为同余方程组，求：
最后使用 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\)
所以
前面一部分容易发现可以递归解决，后面一部分肯定不能硬算。
我们发现后面一部分是有循环节的。\(p^k\) 以后就开始重复了。
这样我们只需算到 \(p^k\) 就行了。
递归解决，时间复杂度 \(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;
}
