Lucas定理

在『组合数学基础』中，我们已经提出了\(Lucas\)定理，并给出了\(Lucas\)定理的证明，本文仅将简单回顾，并给出代码。

\(Lucas\)定理：当\(p\)为质数时，\(C_n^m\equiv C_{n\ mod\ p}^{m\ mod\ p}*C_{n/p}^{m/p}(mod\ p)\)。

在计算模域组合数时，如果模数较小，那么就可以尝试使用\(Lucas\)定理来递归求解，其时间复杂度为\(O(plog_p\min(n,m))\)。

\(Code:\)

inline long long Lucas(long long u,long long d)
{if (u<=Mod&&d<=Mod)return C(u,d) % Mod;else return Lucas(u/Mod,d/Mod) * C(u%Mod,d%Mod) % Mod;
}

本文给出了一种递归写法，由于取模操作的存在，在组合数计算函数\(C\)内也需要进行一些特判，避免一些无意义式子的计算，影响答案。

\(Code:\)

inline long long C(long long u,long long d)
{if ( u>d || u<0 || d<0 )return 0LL;//...
}

关于组合数的计算，我们也有多种方法，常用的几种如下：

\(1.\) 当\(n\)，\(m\)都比较小的时候，可以根据组合数的阶乘计算公式来预处理\(n\)，\(m\)范围以内的阶乘以及阶乘的逆元，然后根据询问\(O(1)\)地回答组合数。

\(2.\) 当\(n\)的范围较大，\(m\)的范围较小时，由于不适合预处理，所以我们可以直接利用有关\(m\)组合数计算式\(C_n^m=\frac{n*(n-1)*...*(n-m+1)}{m*(m-1)*...*1}\)来模拟计算，同理，分母利用逆元操作即可，时间复杂度\(O(m)\)。

Exlucas算法

与之对应的，\(Exlucas\)算法是一种用来解决与\(Lucas\)定理形式很像但更具有一般性的问题的算法。为什么不叫\(Exlucas\)定理而叫\(Exlucas\)算法呢，是因为这个算法和\(Lucas\)定理没什么关系，只是一个数论算法罢了。

对于求解模域组合数\(C_n^m\%p\)，当\(p\)不一定是质数时，可以使用\(Exlucas\)算法求解。

对于这样一个问题，我们考虑对模数进行分解，设\(p=\prod_{i=1}^kp_i^{a_i}\)，答案\(x=C_n^m\%p\)，则可以得到\[x\equiv C_n^ m(mod\ p_i^{a_i})\]
这是一个线性同余方程组，可以使用中国剩余定理求解。

那么现在我们的问题就转换成了求解\(C_n^m\%p_i^{a_i}\)。将组合数写为阶乘的计算式，即\(C_n^m=\frac{n!}{m!(n-m)!}\)。

因为\(\frac{n!}{m!(n-m)!}\)与\(p_i^{a_i}\)不可避免的会存在公约数，所以不能直接使用\(Exeuclid\)等算法来求解逆元，那么我们就要考虑把\(\frac{n!}{m!(n-m)!}\)中有关\(p_i\)的项都提取掉，利用\(\frac{n!}{m!(n-m)!}\)是整数这一特点将分子该项的指数直接减掉分母该项的指数，再求解其余部分以及逆元，最后就能避免不互质的情况从而求解答案。

先考虑一个简单的问题，如何计算一个形如\(n!\)的数中含有多少个因子\(p\)，则和阶乘分解一题类似，若设\(f(n)\)表示\(n!\)中含有因子\(p\)的个数，则有：
\[f(n)=\begin{cases}0\ (n<p)\\\lfloor \frac{n}{p} \rfloor +f(\lfloor \frac{n}{p} \rfloor)\ (n\geq p)\end{cases}\]

简单地利用\(for\)循环求解即可。

求得阶乘中因子\(p_i\)的数量后，我们的问题就变成为如何求解\(n!\)中除去所有\(p_i\)后取模\(p_i^{a_i}\)的值。对于每一项，这又需要我们分两种情况考虑：

\(1.\) 该项是除去\(p_i\)后得到的，我们将所有这样的项放在一起发现，发现除去\(p_i\)后，剩下的系数乘积也是一个阶乘，使用递归求解。
\(2.\) 该项是阶乘中原本的一项，除去情况\(1.\)，我们发现阶乘中若干个连续的项在模\(p_i^{a_i}\)意义下恰好构成了循环节，其长度不超过\(p_i^{a_i}\)，先暴力计算一个循环节，然后快速幂即可。对于不包括在完整循环节中的项，由于其长度小于循环节，可以暴力计算。

然后就能将求得的答案通过计算逆元求解，最后，乘上提取出的\(p_i\)剩余的若干次方即可。

考虑一个简单的例子就能加深理解：

求解\(n!\ mod\ p_i^{p_i}\)，此时，\(n=19\)，\(p_i=3\)，\(a_i=2\)。
\[ n!=1*2*3*...*19\\=(1*2*4*5*7*8*10*11*13*14*16*17*19)*(3^6*6!) \\=(1*2*4*5*7*8)^2*19*(3^6*6!) \]

可以证明，其时间复杂度与\(O(plog_2p)\)同级。

\(Code:\)

#include <bits/stdc++.h>
using namespace std;
const int SIZE=300;
long long A,B,Mod,m[SIZE],r[SIZE];
inline void input(void)
{scanf("%lld%lld%lld",&A,&B,&Mod);//计算C(A,B)%Mod
}
inline long long power(long long a,long long b,long long p)
{long long res = 1;while (b){if (1&b)res = res * a % p;b >>= 1;a = a * a % p;}return res;
}
inline long long Exeuclid(long long a,long long &x,long long b,long long &y,long long c)
{if (!b){x=c/a,y=0;return a;}else{long long p = Exeuclid(b,x,a%b,y,c);long long x_ = x , y_ = y;x = y_ , y = x_ - a/b * y_;return p;}
}
inline long long inv(long long a,long long p)
{long long x_,y_;Exeuclid(a,x_,p,y_,1);return (x_+p)%p;
}
inline long long calc(long long x,long long val,long long p)
{//计算x!中除去val后模p意义下的值if (!x)return 1;long long res = 1 , last = x%p;for (long long i=1;i<=p;i++)if (i%val)res = res * i % p;res = power(res,x/p,p);for (long long i=1;i<=last;i++)if (i%val)res = res * i % p;return res * calc(x/val,val,p) % p;
}
inline long long C(long long d,long long u,long long val,long long Pow)
{long long mulup=calc(d,val,Pow) , k=0;long long muldown1=calc(u,val,Pow) , muldown2=calc(d-u,val,Pow);for (long long i=d;i;i/=val)k+=i/val;for (long long i=u;i;i/=val)k-=i/val;for (long long i=d-u;i;i/=val)k-=i/val;//计算剩余的val的指数return mulup * inv(muldown1,Pow) % Pow * inv(muldown2,Pow) % Pow * power(val,k,Pow) % Pow;
}
inline long long CRT(int cnt)
{long long m_ = Mod , M[SIZE] = {} , t[SIZE] = {} , res = 0;for (int i=1;i<=cnt;i++)M[i] = m_ / m[i];for (int i=1;i<=cnt;i++){long long y;Exeuclid(M[i],t[i],m[i],y,1);res = (res + r[i] % m_ * M[i] % m_ * t[i] % m_) % m_;}return ( res % m_ + m_ ) % m_;
}
inline long long Exlucas(void)
{int temp = Mod , cnt = 0 ;for (long long i=2;i*i<=Mod;i++){if (temp%i==0){m[++cnt] = 1;while (temp%i==0)temp /= i , m[cnt] *= i;r[cnt] = C(A,B,i,m[cnt]);}}if (temp>1) m[++cnt] = temp , r[cnt] = C(A,B,temp,m[cnt]);return CRT(cnt);
}
int main(void)
{input();printf("%lld\n",Exlucas());return 0;
}