CRT、EX-CRT & Lucas、Ex-Lucas

中国剩余定理

中国剩余定理，Chinese Remainder Theorem，又称孙子定理，给出了一元线性同余方程组的有解判定条件，并用构造法给出了通解的具体形式。
现 在 有 方 程 组 ： ( S ) : { x ≡ a 1 ( m o d   m 1 ) x ≡ a 2 ( m o d   m 2 )      .      .      . x ≡ a n ( m o d   m n ) 中 国 剩 余 定 理 指 出 ： 若 ∀ i , j ， 1 ≤ i &lt; j ≤ n ， i , j ∈ Z ， 都 有 g c d ( m i , m j ) = 1 ， 则 对 任 意 整 数 ： a 1 , a 2 , . . . , a n ， 方 程 组 ( S ) 均 有 解 ， 且 通 解 可 以 用 如 下 方 法 构 造 。 设 M = ∏ i = 1 n m i ， 并 设 M i = M m i ， ∀ i ∈ { 1 , 2 , . . . , n } ， 设 t i = M i − 1 ( m o d   m i ) ， 则 方 程 通 解 为 ： x = k M + ∑ i = 1 n a i t i M i ， k ∈ Z 。 \begin{aligned} &amp;现在有方程组：\\ &amp;(S):\begin{cases} x\equiv a_1(mod\space m_1)\\ x\equiv a_2(mod\space m_2)\\ \space\space\space\space. \\ \space\space\space\space. \\ \space\space\space\space. \\ x\equiv a_n(mod\space m_n)\\ \end{cases}\\ &amp;中国剩余定理指出：若\forall i,j，1\le i&lt;j\le n，i,j\in Z，都有gcd(m_i,m_j)=1，\\ &amp;则对任意整数：a_1,a_2,...,a_n，方程组(S)均有解，且通解可以用如下方法构造。\\ &amp;设M=\prod_{i=1}^n m_i，并设M_i={M\over m_i}，\forall i\in\{1,2,...,n\}，\\ &amp;设t_i=M_i^{-1}(mod\space m_i)，则方程通解为：x=kM+\sum_{i=1}^n a_it_iM_i，k\in Z。\\ \end{aligned} ​现在有方程组：(S):⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧​x≡a1​(mod m1​)x≡a2​(mod m2​)    .    .    .x≡an​(mod mn​)​中国剩余定理指出：若∀i,j，1≤i<j≤n，i,j∈Z，都有gcd(mi​,mj​)=1，则对任意整数：a1​,a2​,...,an​，方程组(S)均有解，且通解可以用如下方法构造。设M=i=1∏n​mi​，并设Mi​=mi​M​，∀i∈{1,2,...,n}，设ti​=Mi−1​(mod mi​)，则方程通解为：x=kM+i=1∑n​ai​ti​Mi​，k∈Z。​
下面给出解的正确性证明：
∵ 对 ∀ i , j ， 1 ≤ i &lt; j ≤ n ， i , j ∈ Z ， g c d ( m i , m j ) = 1 ∴ g c d ( m i , M i ) = 1 说 明 ∃ t i ， 使 得 t i M i ≡ 1 ( m o d   m i ) 又 a i t i M i ≡ a i ⋅ 1 ≡ a i ( m o d   m i ) ， 且 对 ∀ i , j ∈ { 1 , 2 , . . . , n } ， i ≠ j ， a i t i M i ≡ 0 ( m o d   m j ) ∴ 对 ∀ i ∈ { 1 , 2 , . . , n } ， x 0 = ∑ j = 1 n a j t j M j = a i t i M i + ∑ j = 1 , j ≠ i n a j t j M j ≡ a i t i M i + ∑ j = 1 , j ≠ i n 0 ( m o d   m i ) 说 明 x 0 是 方 程 组 的 一 个 解 假 设 x 1 , x 2 都 是 方 程 组 的 解 ， 那 么 ∀ i ∈ { 1 , 2 , . . , n } ， x 1 − x 2 ≡ 0 ( m o d   m i ) 又 对 ∀ i , j ， 1 ≤ i &lt; j ≤ n ， i , j ∈ Z ， g c d ( m i , m j ) = 1 ， ∴ ( x 1 − x 2 ) % M = ( x 1 − x 2 ) % ( ∏ i = 1 n m i ) = 0 ∴ ∃ k ∈ Z ， x 1 − x 2 = k M ， ∴ 方 程 通 解 为 ： x = k M + ∑ i = 1 n a i t i M i ， k ∈ Z \begin{aligned} &amp;∵对\forall i,j，1\le i&lt;j\le n，i,j\in Z，gcd(m_i,m_j)=1∴gcd(m_i,M_i)=1 \\ &amp;说明\exist t_i，使得t_iM_i\equiv1(mod\space m_i)\\ &amp;又a_it_iM_i\equiv a_i\cdot1\equiv a_i(mod\space m_i)，且对\forall i,j\in\{1,2,...,n\}，i\ne j，a_it_iM_i\equiv0(mod\space m_j)\\ &amp;∴对\forall i\in\{1,2,..,n\}，x_0=\sum_{j=1}^n a_jt_jM_j=a_it_iM_i+\sum_{j=1,j\ne i}^n a_jt_jM_j\equiv a_it_iM_i+\sum_{j=1,j\ne i}^n0(mod\space m_i)\\ &amp;说明x_0是方程组的一个解\\ &amp;假设x_1,x_2都是方程组的解，那么\forall i\in \{1,2,..,n\}，x_1-x_2\equiv0(mod\space m_i)\\ &amp;又对\forall i,j，1\le i&lt;j\le n，i,j\in Z，gcd(m_i,m_j)=1，∴(x_1-x_2)\%M=(x_1-x_2)\%(\prod_{i=1}^n m_i)=0\\ &amp;∴\exist k\in Z，x_1-x_2=kM，∴方程通解为：x=kM+\sum_{i=1}^n a_it_iM_i，k\in Z \end{aligned} ​∵对∀i,j，1≤i<j≤n，i,j∈Z，gcd(mi​,mj​)=1∴gcd(mi​,Mi​)=1说明∃ti​，使得ti​Mi​≡1(mod mi​)又ai​ti​Mi​≡ai​⋅1≡ai​(mod mi​)，且对∀i,j∈{1,2,...,n}，i̸​=j，ai​ti​Mi​≡0(mod mj​)∴对∀i∈{1,2,..,n}，x0​=j=1∑n​aj​tj​Mj​=ai​ti​Mi​+j=1,j̸​=i∑n​aj​tj​Mj​≡ai​ti​Mi​+j=1,j̸​=i∑n​0(mod mi​)说明x0​是方程组的一个解假设x1​,x2​都是方程组的解，那么∀i∈{1,2,..,n}，x1​−x2​≡0(mod mi​)又对∀i,j，1≤i<j≤n，i,j∈Z，gcd(mi​,mj​)=1，∴(x1​−x2​)%M=(x1​−x2​)%(i=1∏n​mi​)=0∴∃k∈Z，x1​−x2​=kM，∴方程通解为：x=kM+i=1∑n​ai​ti​Mi​，k∈Z​

扩展中国剩余定理

在一般情况下，要求任两个数互质这个条件太苛刻了，CRT派不上用场，我们需要一个更具普遍性的结论，这就是EX-CRT。虽然是称为EX-CRT，但这个定理并没有直接用到CRT的结论。
$$\begin{array}{rl}& 我们单独考虑方程组的前两个式子组成的同余方程组S^{\prime }：\\ & \{\begin{array}{l}x\equiv a_1(mod{m}_1)\\ x\equiv a_2(mod{m}_2)\end{array}\\ & 我们设整数K_1,K_2满足x=K_1{m}_1+a_1,x=K_2{m}_2+a_2，那么有K_1{m}_1+a_1=K_2{m}_2+a_2，\\ & 即m_1{K}_1-m_2{K}_2=a_2-a_1，根据EXGCD，方程在gcd(m_1,m_2)|(a_2-a_1)时有解\\ & 把由EXGCD求得的特解记为k_1,k_2，设S’的一个解为x_0，那么有：\\ & \{\begin{array}{l}{x}_0=a_1+k_1\times \frac{a_2-a_1}{gcd(m_1,m_2)}\times m_1\\ x_0=a_2+k_2\times \frac{a_2-a_1}{gcd(m_1,m_2)}\times m_2\end{array}\\ & 这样我们就得到了S^{\prime }的一个解x_0，怎么求出它的通解呢？\end{array}$$

---

$$\begin{array}{rl}& 假设S^{\prime }存在两个解x_1,x_2(x_1<x_2)，由于x_1,x_2均满足S^{\prime }，则：\\ & \exists {\sigma }_1,{\sigma }_2\in Z，满足：\{\begin{array}{l}{x}_2=x_1+{\sigma }_1\times m_1\\ x_2=x_1+{\sigma }_2\times m_2\end{array}\\ & 显然{\sigma }_1\times m_1={\sigma }_2\times m_2，设d=gcd(m_1,m_2),n_1=\frac{m_1}{d},n_2=\frac{m_2}{d}\\ & 那么上式\iff {\sigma }_1\times n_{1}d={\sigma }_2\times n_{2}d\iff {\sigma }_1\times n_1={\sigma }_2\times n_2\therefore n_1|({\sigma }_2\times n_2)\\ & 又由定义可知，gcd(n_1,n_2)=1\therefore n_1|{\sigma }_2，也即(n_1{n}_2)|({\sigma }_2\times n_2)\to (n_1\times n_{2}d)|({\sigma }_2\times n_{2}d)\\ & \to (n_1\times n_{2}d)|({\sigma }_2\times m_2)\to (n_1\times n_{2}d)|(x_2-x_1)(由x_2=x_1+{\sigma }_2\times m_2可得)\\ & n_1\times n_{2}d=\frac{n_{1}d\times n_{2}d}{d}=\frac{m_1{m}_2}{gcd(m_1,m_2)}=lcm(m_1,m_2)\\ & 由此我们得到了S’任两解之间的关系：lcm(m_1,m_2)|(x_2-x_1)\end{array}$$

---

$$\begin{array}{rl}& 由解的关系有x_2=x_1+k\cdot lcm(m_1,m_2)，k\in Z，结合之前得到的特解x_0，可以得出S^{\prime }的通解X：\\ & X=x_0+k\cdot lcm(m_1,m_2)，k\in Z\\ & 上式又可以写成X\equiv x_0(modlcm(m_1,m_2))，这样，我们就把S^{\prime }转化为了一条同余式\\ & 对于n条同余式，只需逐条合并即可。\end{array}$$

代码如下：

typedef long long ll;
const int maxn = 111;

// m为模数组，a为余数数组，0~n-1
ll m[maxn], a[maxn];

ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1; y = 0;
        return a;
    }
    ll ans = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return ans;
}

ll excrt() {
    ll lcm = m[0], last_a = a[0];
    for(int i = 1; i < n; i++) {
        ll lcm_a = ((a[i] - last_a) % m[i] + m[i]) % m[i];
        ll k = lcm, x, y;
        ll gcd = exgcd(lcm, m[i], x, y);
        ll mod = m[i] / gcd;
        x = (x * lcm_a / gcd % mod + mod) % mod;
        lcm = lcm / gcd * m[i], last_a = (last_a + k * x) % lcm;
    }
    return (last_a % lcm + lcm) % lcm;
}

卢卡斯定理

卢卡斯定理是关于组合数和同余的定理，它表明当p为素数时：
$$\begin{gathered} C_n^m=\prod _{i=0}^k{C}_{n_i}^{m_i}(modp) \\ 其中，m=m_k{p}^k+m_{k-1}{p}^{k-1}+...+m_{1}p+m_0 \\ n=n_k{p}^k+n_{k-1}{p}^{k-1}+...+n_{1}p+n_0 \\ 即m_i,n_i为m,n的p进制展开中对应次数为i的系数 \end{gathered}$$
因为当m>n时，二项式系数为0，那么二项式系数即组合数能被p整除等价于在p进制下，存在某一位m的数值大于对应的n的数值。

基于母函数可以简单证明这个定理。
$$\begin{array}{rl}& \because 无论p是质数还是合数，对\forall i\in N,0<i<p，都有C_p^i=\frac{p!}{i!(p-i)!}\\ & \therefore 在母函数中，对于素数p，(1+x{)}^p\equiv 1+x^p(modp)\\ & 运用数学归纳法可以得到，对\forall i\in N，(1+x{)}^{p^i}\equiv 1+x^{p^i}(modp)\\ & 又对\forall t\in N以及素数p，把t表示为p进制数有：t=\sum _{i=0}^k{t}_i{p}^i\\ & 那么，把m,n都用p进制表示，有\\ & \sum _{m=0}^n\left(C_n^m\cdot x^m\right)\\ & =(1+x{)}^n=\prod _{i=0}^k{\left[(1+x{)}^{p^i}\right]}^{m_i}\\ & \equiv \prod _{i=0}^k{\left(1+x^{p^i}\right)}^{m_i}=\prod _{i=0}^k\left(\sum _{m_i=0}^{n_i}{C}_{n_i}^{m_i}\cdot x^{m_i{p}^i}\right)\\ & =\prod _{i=0}^k\left(\sum _{m_i=0}^{p-1}{C}_{n_i}^{m_i}\cdot x^{m_i{p}^i}\right)=\sum _{m=0}^n\left[\left(\prod _{i=0}^k{C}_{n_i}^{m_i}\right)\cdot x^m\right](modp)\end{array}$$
代码如下：

typedef long long ll;
const int mod = 1e9 + 7;
const int maxn = 1e5 + 100;

void init() {
	F[0] = 1;
	for(int i = 2; i < maxn; i++)
		F[i] = i * F[i - 1] % mod;
}

ll qpow(ll a, ll b) {
	ll ans = 1;
	while(b) {
		if(b & 1) ans = ans * a % mod;
		b >>= 1; a = a * a % mod;
	}
	return ans;
}

ll lucas(ll N, ll M) {
	ll ans = 1;
	while(N & M) {
		ll n = N % mod, m = M % mod;
		if(n < m) return 0;
		ans = ans * F[a] % mod * qpow(F[m] * F[n - m] % mod, mod - 2) % mod;
		N /= p; M /= p;
	}
	return ans;
}

扩展卢卡斯定理

卢卡斯定理同样不能处理模数不是素数的情况，这时便需要扩展卢卡斯定理。我们一步步分析如何求解模数不是素数的组合数问题。

代码和推导部分参考了这篇博客。
$$\begin{array}{rl}& 首先，我们要解决的问题是求C_n^m\%p，其中p不一定是素数。\\ & 对于非素数，我们首先会联想到质因分解后结合CRT解决问题。\\ & 假设分解得到t个质数，质数p_i对应的个数为k_i，对p质因分解有p=\prod _{i=1}^t{p}_i^{k_i}。\\ & 显然对\forall i,j，1\le i<j\le t，gcd(p_i^{k_i},p_j^{k_j})=1，假设对\forall i\in [1,t]，我们求出了a_i=C_n^m\%p_i^{k_i}，\\ & 那么我们可以得到同余方程组S：\\ & S：\{\begin{array}{l}{C}_n^m\equiv a_1(mod{p}_1^{k_1})\\ C_n^m\equiv a_2(mod{p}_2^{k_2})\\ .\\ .\\ .\\ C_n^m\equiv a_t(mod{p}_t^{k_t})\end{array}\\ & 这时我们便可以套用CRT解决问题，那么问题便转化为如何求解C_n^m\%p_i^{k_i}。\end{array}$$

---

$$\begin{array}{rl}& 现在我们要求的是C_n^m\%p^k，其中p是素数。\\ & 又C_n^m=\frac{n!}{m!(n-m)!}，显然需要求出m!和(n-m)!关于模p^k的逆元，\\ & 但考虑到这些项中可能包含p（含有p则不互质，逆元不存在），所以需要先提取p，\\ & 得到：C_n^m=\frac{\frac{n!}{p^{k_n}}}{\frac{m!}{p^{k_m}}\cdot \frac{(n-m)!}{p^{k_{n-m}}}}\times p^{k_n-k_m-k_{n-m}}，这里的阶乘是指提取之后的结果。\\ & 这时就可以计算\frac{m!}{p^{k_m}}和\frac{(n-m)!}{p^{k_{n-m}}}关于p^k的逆元了。\\ & 这里，为了形式的统一，同时提取了n!中的p。那么，问题又转化为了如何求n!\%p^k。\end{array}$$

---

$$\begin{array}{rl}& 目标：计算n!\%p^k，p为质数。上一步中提到，我们需要先提取p。\\ & 提取结果为：n!=p^{\lfloor \frac{n}{p}\rfloor }\times \lfloor \frac{n}{p}\rfloor !\times \prod _{i=1,i\%p̸=0}^{n}i。\\ & 第一部分很好理解，对于每一个p的倍数，都可以提取出一个p，一共有\lfloor \frac{n}{p}\rfloor 个；\\ & 第二部分为p的倍数被提取p之后余下的，是一个阶乘的形式。显然在n!中，\\ & 对于p的幂，p的个数不止1个，也就是说第二部分仍然需要提取，这一部分可以递归解决。\\ & 第三部分是n!剔除了p的倍数之后余下的。\\ & \because 对\forall i<p^k，t\in N，都有i\equiv (i+t\cdot p^k)(mod{p}^k)\\ & \therefore 对\forall t\in N，\prod _{i=1,i\%p̸=0}^{p^k}i\equiv \prod _{i=1,i\%p̸=0}^{p^k}(i+t\cdot p^k)(mod{p}^k)\\ & 这也就是说，第三部分其实是存在循环节的，\prod _{i=1,i\%p̸=0}^{p^k}i一共循环了\lfloor \frac{n}{p^k}\rfloor 次。\\ & 除去循环节的余项长度在p^k之内，直接累乘即可。\end{array}$$

完整代码如下：

typedef long long ll;
const int N = 1e6 + 100;

ll n, m, p;

ll qpow(ll a, ll b, ll mod) {
	ll ans = 1;
	while(b) {
		if(b & 1) ans = ans * a % mod;
		b >>= 1; a = a * a % mod;
	}
	return ans;
}

ll fac(ll n, ll p, ll pk) {
    if (!n) return 1;
    ll ans = 1;
    for (int i = 1; i < pk; i++)
        if (i % p) ans = ans * i % pk;
    ans = qpow(ans, n / pk, pk);
    int npk = n % pk;
    for (int i = 1; i <= npk; i++)
        if (i % p) ans = ans * i % pk;
    return ans * fac(n / p, p, pk) % pk;
}

ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) {
        x = 1; y = 0;
        return a;
    }
    ll ans = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return ans;
}

ll inv(ll a, ll p) {
    if(!a) return 0;
    int x, y;
    exgcd(a, p, x, y);
    return (x % p + p) % p;
}

ll C(ll n, ll m, ll p, ll pk) {
    if (n < m) return 0;
    ll fn = fac(n, p, pk),
       fm = fac(m, p, pk),
       fn_m = fac(n - m, p, pk),
       cnt = 0;
    for (ll i = n; i; i /= p)
        cnt += i / p;
    for (ll i = m; i; i /= p)
        cnt -= i / p;
    for (ll i = n - m; i; i /= p)
        cnt -= i / p;
    return fn * inv(fm * fn_m % pk, pk) % pk * qpow(p, cnt, pk) % pk;
}

ll a[N], mod[N]; // a[]是通过卢卡斯分解出来的组合数值，m[]是对应的模数
int cnt; // 质因数的种数

ll CRT() {
    ll M = 1, ans = 0;
    for (int i = 0; i < cnt; i++)
        M *= mod[i];
    for (int i = 0; i < cnt; i++)
        ans = (ans + a[i] * (M / mod[i]) % M * inv(M / mod[i], mod[i]) % M) % M;
    return ans;
}

ll exlucas(ll n, ll m, ll p) {
    ll sqrtp = sqrt(p + 0.5);
    for (int i = 2; p > 1 && i <= sqrtp; i++) {
        ll pk = 1;
        while (p % i == 0)
            p /= i, pk *= i;
        if (pk > 1)
            a[cnt] = C(n, m, i, pk), mod[cnt++] = pk;
    }
    if (p > 1)
        a[cnt] = C(n, m, p, p), mod[cnt++] = p;
    return CRT();
}