<更新提示>
<第一次更新>
<正文>
Lucas定理
在『组合数学基础』中，我们已经提出了(Lucas)定理，并给出了(Lucas)定理的证明，本文仅将简单回顾，并给出代码。
(Lucas)定理：当(p)为质数时，(C_n^mequiv C_{n mod p}^{m mod p}*C_{n/p}^{m/p}(mod p))。
在计算模域组合数时，如果模数较小，那么就可以尝试使用(Lucas)定理来递归求解，其时间复杂度为(O(plog_pmin(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)，则可以得到$$xequiv 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)的个数，则有：
简单地利用(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)。
可以证明，其时间复杂度与(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
",Exlucas());
return 0;
}
<后记>