基本数学

中国剩余定理(CRT)

中国剩余定理是一个很玄学的东西，记不住老是忘，所以今天来写一下。

其实就是让我们求解一类方程组，形如

\[ \begin{cases} x \equiv a_1 \ (\bmod\ b_1)\\ x \equiv a_2 \ (\bmod\ b_2)\\ ...\\ x \equiv a_n \ (\bmod\ b_n)\\ \end{cases} \]

其中\(b_1,b_2,b_3,...,b_n\)互质，这就是普通\(CRT\)成立的先决条件。那么我们如何得到\(CRT\)的通解呢。

答案为$$ans=\sum\limits_{i=1}^{n}a_ic_i$$

那么\(c_i\)是什么呢，我们设\(Mod=\prod_{i=1}^{n}b_i\)，设\(m_i=\frac{Mod}{b_i}\),\(c_i\)就等于\(m_i\times m_i^{-1}\)(\(m_i^-1\)是在模\(b_i\)意义下的逆元),那么由于\(b_i\)间互质，所以\(m_i\)与\(b_i\)也是互质的，所以我们可以证明上面那种构造方法是正确的，证明如下：

证明：
对于每一组方程，我们有

\[ \sum_{j=1}^{n}a_jc_j \equiv\ (\bmod\ b_i) \]

由于除\(m_i\)以外的所有\(m_j\)都是可以整除\(b_i\)的，所以

\[ \begin{aligned} & \sum_{j=1}^{n}a_jc_j \equiv\ (\bmod\ b_i)\\ \Rightarrow & a_ic_i \equiv\ (\bmod\ b_i)\\ \Rightarrow & a_im_im_i^{-1} \equiv\ (\bmod\ b_i)\\ \Rightarrow & a_i \equiv\ (\bmod\ b_i)\\ \end{aligned} \]

所以我们证明完了，这就是一组\(CRT\)的通解，所以我们求解\(CRT\)的步骤如下

\(\bullet\ 1\),求出\(Mod=\prod_{i=1}^{n}b_i\)

\(\bullet\ 2\),求出\(m_i=\frac{Mod}{b_i}\)

\(\bullet\ 3\),用\(exgcd\)求出\(m_i\)在模\(b_i\)意义下的逆元\(m_i^{-1}\)

\(\bullet\ 4\),\(ans=\sum_{i=1}^{n}a_im_im_i^{-1}\),注意此处模\(Mod\),而不是模\(b_i\),如果模的是\(b_i\),那么\(m_im_i^{-1}\ \bmod \ b_i=1\),就没有意义了。

\(code\)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=20;
int n;
long long a[N],b[N];
long long exgcd(long long a,long long b,long long &x,long long &y){
    if(b==0){x=1,y=0;return a;}
    long long ret=exgcd(b,a%b,x,y);
    long long t=x;
    x=y;
    y=t-a/b*y;
    return ret;
}
long long inv(long long a,long long mod){
    long long x,y;
    exgcd(a,mod,x,y);
    x=(x%mod+mod)%mod;
    return x;
}
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++)scanf("%lld%lld",&a[i],&b[i]);
    long long Mod=1;
    for(int i=1;i<=n;i++)Mod*=a[i];
    long long ans=0;
    for(int i=1;i<=n;i++){
        long long m=Mod/a[i];
        long long ny=inv(m,a[i]);
        ans=(ans+b[i]*m%Mod*ny%Mod)%Mod;
    }
    printf("%lld\n",ans);
    return 0;
}

扩展中国剩余定理(exCRT)

普通的中国剩余定理只能解决模数之间互质的线性同余方程，但是当模数之间不互质了怎么办，那就要请出我们的\(exCRT\)了。(虽然说我感觉这个东西与\(CRT\)没有关系，反而与\(exgcd\)有很大的关联)

对于如下这个同余方程组，我们可以一个一个合并，我们举个例子

\[\begin{aligned} & x \equiv a_1 \ (\bmod\ b_1)\\ & x \equiv a_2 \ (\bmod\ b_2)\\ \Rightarrow & x=k_1\times b_1+a_1=k_2\times b_2+a_2\\ \Rightarrow & k_1\times b_1+a_1=k_2\times b_2+a_2\\ \Rightarrow & k_1\times b_1-k_2\times b_2=a_2-a_1\\ \Rightarrow & k_1\times b_1+k_2\times b_2=a_2-a_1\\ \end{aligned} \]

我们发现最后那个式子很像\(exgcd\)的形式，我们可以通过\(exgcd\)求出一组\(k_1,k_2\),当\(gcd(b_1,b_2)\nmid a_2-a_1\)时，显然无解，当\(gcd(b_1,b_2)\mid a_2-a_1\)时，我们求出\(k_1\)的最小正整数解，将\(ans\)更新为\(k_1\times b_1+a_1\),将模数更新为\(lcm(b_1,b_2)\)即可递推,即为下列转化

\[ \begin{cases} x \equiv a_1 \ (\bmod\ b_1)\\ x \equiv a_2 \ (\bmod\ b_2)\\ \end{cases} \Rightarrow x\equiv(k_1\times b_1+a_1)\ (\bmod\ lcm(b_1,b_2)) \]

\(code\)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int N=2e5+10;
int n;
long long a[N],b[N];
long long gcd(long long a,long long b){long long r=a%b;while(r){a=b,b=r,r=a%b;}return b;}
long long mul(long long a,long long b,long long p){long long ans=0;while(b){if(b&1)ans=(ans+a)%p;a=(a+a)%p;b>>=1;}return ans;}
long long lcm(long long a,long long b){return a/gcd(a,b)*b;}
long long exgcd(long long a,long long b,long long& x,long long& y){
    if(b==0){x=1,y=0;return a;}
    long long ret=exgcd(b,a%b,x,y);
    long long t=x;
    x=y;
    y=t-a/b*y;
    return ret;
}
long long ans,mod;
inline long long read(){
    long long x=0;
    int f=0;
    char ch;
    while(!isdigit(ch=getchar()))if(ch=='-')f=1;
    do{
        x=(x<<1)+(x<<3)+(ch^48);
    }while(isdigit(ch=getchar()));
    return f?-x:x;
}
int stk[200],top;
inline void write(long long x){
    if(x<0)putchar('-'),x=-x;
    while(x){
        stk[++top]=x%10;
        x/=10;
    }
    while(top)putchar(stk[top--]|48);
    putchar('\n');
}
int main(){
    scanf("%d",&n);
    for(int i=1;i<=n;i++)b[i]=read(),a[i]=read();
    ans=a[1],mod=b[1];
    for(int i=2;i<=n;i++){
        long long c=a[i]-ans;
        long long x=0,y=0;
        long long d=exgcd(mod,b[i],x,y);
        if(c%d!=0){
            printf("No solution!\n");
            return 0;
        }
        x=x*c/d;
        x=(x%(b[i]/d)+b[i]/d)%(b[i]/d);
        long long M=lcm(mod,b[i]);
        ans=(ans%M+mul(x,mod,M)+M)%M;
        mod=M;
    }
    write(ans%mod);
    return 0;
}

卢卡斯定理(lucas)

定理如下：

\[{n\choose m}\equiv {n/p\choose m/p}\times {n\bmod p\choose m\bmod p}(\bmod\ p),p\in prime \]

首先对于一个质数\(p\),有

\[{p\choose i}= \begin{cases} 0,(i>1\land i<p)\\ 1,(i=1\lor i=p)\\ \end{cases} (1) \]

证明：

我们考虑对于一个质数\(p\)，\({p\choose i}\)的分子上时刻都有质因子\(p\)，而分母中只有当\(i=1\)或\(i=p\)的时候才有质因子\(p\),所以得证。

那么我们设\(n=sp+q(q < p),m=tp+r(r < p)\)

我们考虑一个式子\((1+x)^n\)的系数，我们开始证明:

\[\begin{aligned} & (1+x)^n\\ \equiv & (1+x)^{sp+q}\\ \equiv & (1+x)^{sp}\times (1+x)^q\\ \equiv & ((1+x)^p)^s\times (1+x)^q\\ \end{aligned} (\bmod\ p) \]

由\((1)\)得，\((1+x)^p\equiv (1+x^p)(\bmod\ p)\)

那么

\[\begin{aligned} & (1+x)^n\\ \equiv & (1+x^p)^s\times (1+x)^q\\ \equiv & \sum\limits_{i=0}^{s}{s\choose i}x^{ip}\times \sum\limits_{j=0}^{q}{q\choose j}x^j \end{aligned} (\bmod\ p)\ (2) \]

我们发现直接求\((1+x)^n\)的\(x^{tp+r}\)项的系数为\({sp+q\choose tp+r}\)

而对于\((2)\)式来说，只能找到一个次幂为\(tp+r\)的组合，其系数为\({s\choose t}\times {q\choose r}\)

则有

\[{sp+q\choose tp+r}\equiv {s\choose t}\times {q\choose r} \]

代换得到

\[{n\choose m}\equiv {n/p\choose m/p}\times {n\bmod p\choose m\bmod p} \]

所以卢卡斯定理得证。

我们就可以解决模数较小但\(n,m\)较大的组合数问题了。

\(code\)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=1e5+10;
long long n,m,p,fac[N],ny[N],nyy[N];
int T;
void init(){
    fac[0]=fac[1]=nyy[0]=nyy[1]=ny[0]=ny[1]=1;
    for(int i=2;i<p;i++){
        fac[i]=1ll*fac[i-1]*i%p;
        ny[i]=1ll*(p-p/i)*ny[p%i]%p;
        nyy[i]=1ll*nyy[i-1]*ny[i]%p;
    }
}
long long C(long long n,long long m){
    if(m>n)return 0;
    return 1ll*fac[n]*nyy[m]%p*nyy[n-m]%p;
}
long long lucas(long long n,long long m){
    if(n<p)return C(n,m);
    return lucas(n/p,m/p)%p*C(n%p,m%p)%p;
}
void work(){
    scanf("%lld%lld%lld",&n,&m,&p);
    init();
    printf("%lld\n",lucas(n+m,m));
}
int main(){
    scanf("%d",&T);
    while(T--)work();
    return 0;
}

扩展卢卡斯定理(exlucas)

我们知道卢卡斯定理只能用来求解模数较小且为质数的组合数问题，那么有没有一种方法能用来求解模数较小但不一定是质数的组合数问题呢，扩展卢卡斯就应运而生了(虽说与 \(lucas\) 没什么关系)。

我们要求解 \({n\choose m}\bmod\;p\) 这类问题，我们将\(p\)用唯一分解定理分解，即 \(p=\prod_{i=1}^{k}p_i^{c_i}\)，那么我们的答案即为求解一个同余方程组，即求解

\[\begin{cases} {n\choose m}\equiv a_1 (\bmod\;p_1^{c_1})\\ {n\choose m}\equiv a_2 (\bmod\;p_2^{c_2})\\ ...\\ {n\choose m}\equiv a_k (\bmod\;p_k^{c_k}) \end{cases} \]

因为 \(p_1^{c_1},p_2^{c_2},...,p_k^{c_k}\) 两两之间是互质的，那么我们就能用 \(CRT\) 将其合并，所以我们就可以知道答案了。

那么现在的问题是对于每一个 \(i\),求 \({n\choose m}\bmod\;p_i^{c_i}\)。

我们发现： \({n\choose m}=\frac{n!}{m!(n-m)!}\)，其中 \(m!\) 或 \((n-m)!\) 在模 \(p_i^{c_i}\) 意义下可能没有逆元，那么我们转化一下式子：

\[ {n\choose m}=\frac{n!}{m!(n-m)!}=\frac{\frac{n!}{p_i^x}}{\frac{m!}{p_i^y}\times \frac{(n-m)!}{p_i^z}}\times p_i^{x-y-z} \]

那么其实分母就与 \(p_i^{c_i}\) 互质了，我们就能用 \(exgcd\) 求解其逆元了。

那么现在的问题变为了求解形如 \(\frac{n!}{p_i^k}\bmod p_i^{c_i}\)

我们举一个例子，对于 \(n=19,p_i=3,c_i=2\) 时的求解

\[\begin{aligned} & 19!=1\times 2\times 3\times ...\times 19\\ = & (1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19)\times 3^6\times (1\times 2\times 3\times 4\times 5\times 6) \end{aligned} \]

我们发现 \((1\times 2\times 4\times 5\times 7\times 8)\equiv (10\times 11\times 13\times 14\times 16\times 17)(\bmod\; p_i^{c_i})\)，那么我们就能推出式子

\[ \frac{n!}{p_i^k} \equiv (\prod\limits_{i=1,p_i\nmid i}^{p_i^{c_i}} i)^{\lfloor \frac{n}{(p_i^{c_i})} \rfloor}\times \prod\limits_{i=1,p_i\nmid i}^{n \bmod p_i^{c_i}}i (\bmod \;p_i^{c_i}) \]

那么我们就可以求出 \({n\choose m}\bmod p\) 的答案了。

\(code\)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
long long exgcd(long long a,long long b,long long& x,long long& y){
    if(b==0){x=1,y=0;return a;}
    long long ret=exgcd(b,a%b,x,y);
    long long t=x;
    x=y;
    y=t-a/b*y;
    return ret;
}
long long inv(long long a,long long p){
    long long x=0,y=0;
    exgcd(a,p,x,y);
    x=(x%p+p)%p;
    return x;
}
long long ksm(long long b,long long k,long long p){
    long long ans=1;
    while(k){
        if(k&1)ans=(ans*b)%p;
        b=b*b%p;
        k>>=1;
    }
    return ans;
}
long long calc(long long n,long long p_i,long long p_k){
    if(n==0)return 1ll;
    long long res=1;
    for(long long i=2;i<=p_k;i++){
        if(i%p_i)res=(res*i)%p_k;
    }
    res=ksm(res,n/p_k,p_k);
    for(long long i=2;i<=n%p_k;i++){
        if(i%p_i)res=(res*i)%p_k;
    }
    return res*calc(n/p_i,p_i,p_k)%p_k;
}
long long C(long long n,long long m,long long p_i,long long p_k){
    long long x=calc(n,p_i,p_k),y=calc(m,p_i,p_k),z=calc(n-m,p_i,p_k);
    long long cnt=0;
    for(long long i=n;i;)cnt+=i/=p_i;
    for(long long i=m;i;)cnt-=i/=p_i;
    for(long long i=n-m;i;)cnt-=i/=p_i;
    return x*inv(y,p_k)%p_k*inv(z,p_k)%p_k*ksm(p_i,cnt,p_k)%p_k;
}
long long crt(long long a,long long b,long long M){
    return a*inv(M/b,b)%M*(M/b)%M;
}
long long exlucas(long long n,long long m,long long p){
    long long s=sqrt(p),tmp=p;
    long long res=0;
    for(int i=2;i<=s;i++){
        if(tmp%i==0){
            long long p_k=1;
            while(tmp%i==0){
                p_k*=i;
                tmp/=i;
            }
            res=(res+crt(C(n,m,i,p_k),p_k,p)%p)%p;
        }
    }
    if(tmp>1ll)res=(res+crt(C(n,m,tmp,tmp),tmp,p)%p)%p;
    return res;
}
long long n,m,p;
int main(){
    scanf("%lld%lld%lld",&n,&m,&p);
    printf("%lld\n",exlucas(n,m,p));
    return 0;
}