卢卡斯定理

$C_n^m\ mod\ p=C_{n/p}^{m/p}\times C_{n\%p}^{m\%p}\ mod\ p$

可以看作是将 $n$ 和 $m$ 用 $p$ 进制分解，对于每一位分别求组合数再乘起来。

扩展卢卡斯定理

由于求组合数时需要用到除法，如果 $p$ 不为质数可能不存在逆元，这时卢卡斯定理不成立。

扩展卢卡斯定理将 $p$ 分解成 $p=p_1^{q_1}\times p_2^{q_2}\times \cdots \times p_m^{q_m}$ 的形式，计算 $C_n^m$ 对每个 $p_i^{q_i}$ 分别取模的值再用中国剩余定理合并。

对于模 $p_i^{q_i}$ 时，将阶乘中的 $p_i$ 的倍数提出来用快速幂单独计算，一个阶乘提出 $p_i$ 的倍数再除掉 $p_i$ 后仍然是一个阶乘，于是可用三个循环得到 $p_i$ 的指数。

for(ll i = n;i > 0;i /= p)c += i / p;

for(ll i = m;i > 0;i /= p)c -= i / p;

for(ll i = n - m;i > 0;i /= p)c -= i / p;

去掉所有 $p$ 的倍数后与模数互质，可以求逆元，在求逆元时由于每 $p_i^{q_i}$ 循环一次，所以可以只计算出 $\begin{align}(\Pi_{i=2}^{p_i^{q_i}}i)^{\lfloor\frac n p\rfloor}\times\Pi_{i=2}^{n\%p_i^{q_i}}i\ mod\ p_i^{q_i}\end{align}$ ，由于在前一步中提出来了一个 $p_i$ 的指数，而这一步计算了 $p_i$ 的倍数， $p_i$ 的倍数除去指数剩下的还是一个阶乘，范围为 $1$ ~ $\lfloor\frac n p\rfloor$ ，可以递归计算。

ll tot = 0,divi[100000],times[100000];

ll sumw = 0;

void dev()

{

    ll k = MOD;

    for(ll i = 2;i <= k;++i)

    {

        if(k % i == 0)

        {

            divi[++tot] = i;

            times[tot] = 1;

            while(k % i == 0)

            {

                times[tot] *= i;

                k /= i;

            }

        }

    }

    return;

}

ll power(ll a,ll b,ll mod)//a ^ b % mod

{

    ll res = 1;

    while(b > 0)

    {

        if(b & 1)res = res * a % mod;

        a = a * a % mod;

        b = b >> 1;

    }

    return res;

}

ll fac(ll n,ll p,ll pr)//n! % pr

{

    if(n == 0)return 1;

    ll res = 1;

    for(ll i = 2;i <= pr;++i)

    {

        if(i % p != 0)

        {

            res = res * i % pr;

        }

    }

    res = power(res,n / pr,pr);

    ll r = n % pr;

    for(ll i = 2;i <= r;++i)

    {

        if(i % p != 0)

        {

            res = res * i % pr;

        }

    }

    res = res * fac(n / p,p,pr) % pr;

    return res;

}

ll exgcd(ll a,ll b,ll &x,ll &y)//ax + by = gcd(a,b)

{

    if(b == 0)

    {

        x = 1;y = 0;

        return a;

    }

    ll t = exgcd(b,a % b,x,y);

    ll k = x;

    x = y;

    y = k - (a / b) * y;

    return t;

}

ll inv(ll a,ll n)//a * a^{-1} = 1 (mod n)

{

    ll x,y;

    ll t = exgcd(a,n,x,y);

    x = (x % n + n) % n;

    return x;

}

ll lucas(ll n,ll m,ll p,ll pr)

{

    if(n < m)return 0;

    ll c = 0;

    for(ll i = n;i > 0;i /= p)c += i / p;

    for(ll i = m;i > 0;i /= p)c -= i / p;

    for(ll i = n - m;i > 0;i /= p)c -= i / p;

    ll res = fac(n,p,pr) * inv(fac(m,p,pr),pr) % pr * inv(fac(n - m,p,pr),pr) % pr * power(p,c,pr) % pr;

    return res;

}

ll CRT(ll n,ll m)

{

    ll ans = 0;

    for(int i = 1;i <= tot;++i)

    {

        ans = (ans + lucas(n,m,divi[i],times[i]) * inv(MOD / times[i],times[i]) % MOD * (MOD / times[i]) % MOD) % MOD;

    }

    return ans;

}