快速求排列组合C（m，n）%mod

写在前面：
1. 为防止产生n和m的歧义，本博文一律默认n >= m
2. 本博文默认mod = 10^6+3
3. 本博文假设读者已知排列组合公式

$C （ m ， n ） = n ！（ n - m ）！ * m ！$

4. 普通的小数据就不用多说了，直接用公式，当然别忘了取模

$C （ m ， n ） = C （ m - 1 ， n - 1 ） + C （ m ， n - 1 ）$

现在我们讨论当n可达10^9数量级大小时的算法。

步骤一：我们先把分子阶乘写成以下形式

$n ！ = X * m o d Y$

步骤二：对分母元素乘机求逆元。此时我们假设得到了以下方程式：

$n ！（ n - m ）！ * m ！ = A * m o d B G a m m a ( T ) * m o d D = A * C * m o d B m o d D$

其中Gamma（T）表示分母剩余数字的乘积，C为他的逆元。

步骤三：显然根据上式我们就可以得出结论了

如果B > D ，那么我们的最终答案为 C（m，n）% mod = 0

否则我们的答案为C（m，n） = （A * C） % mod

注意事项

求A的过程中，我们会发现最后的结果会变成：

$[1 * 2 * . . . * (m o d - 1) * (m o d + 1) * \dots] * m o d B * [1 * 2 * \dots * (n / m o d)]$

如果有多组数据，而mod的大小又不变，那么我们完全可以对k！%mod进行预处理

在对分母的每个元素求逆元时，我们可以由mod是素数直接用欧拉函数求出其逆元
1. 欧拉函数
  
  $p h i (m o d) = m o d - 1$
2. 求逆元，用快速幂，同时别忘了取模
  
  $i n v (x) = x p h i (m o d) - 1$

最后代码如下：

@Frosero#include 
     
      #include 
      
       using namespace std;const long long mod = 1000003;long long mul[1000100];long long pow_mod(long long a ,long long p){        //递归快速幂    if(p == 0) return 1;    long long ans = pow_mod(a,p/2);    ans = ans * ans % mod;    if(p % 2) ans = ans * a % mod;    return ans;}void ask(long long n,long long &x,long long y){     //将阶乘 n！ 拆分成 x * mod ^ y 的形式    y = n / mod;    x = (pow_mod(mul[mod - 1],y) * mul[n - y * mod] * mul[y]) % mod;}int main(){    mul[0] = mul[1] = 1;    for(long long i = 2 ;i < mod ;i++){             //预处理 0 至 mod-1 的阶乘取模值        mul[i] = mul[i - 1] * i % mod;    }    long long m,n,x,y;    long long A,B,C,D;    while(cin >> m >> n){       //对应上文中所讲的 即转换成 A * C * mod ^ B / mod ^ D 的形式        ask(n,A,B); ask(m,C,D); ask(n-m,x,y);        D += y;        C = pow_mod(C,mod - 2) * pow_mod(x,mod - 2) % mod;        if(B > D) cout<<0<

补充