2015 ICL, Finals, Div. 1 Ceizenpok’s formula && Gym - 100633J

上来先感谢大佬的文章：https://blog.csdn.net/clove_unique/article/details/54571216

题目链接：http://codeforces.com/gym/100633/problem/J

题意：

求

$$C_n^k\%p$$ 的值，其中 $1\le n \le 10^{18}$, $0\le k \le n$, $2 \le p \le 1000000$.

分析：

若p为质数，可以用lucas定理直接得，但题目中p不一定为质数，则需要用拓展lucas定理+中国剩余定理。

结论：

拓展lucas定理：

可求方程：

$$ans\equiv c(mod\, p_1^{k_1})$$

中国剩余定理（又名孙子定理）

可求

$$\begin{cases} ans\equiv c_1(mod \, m_1)\\ ans\equiv c_2(mod\, m_2)\\ ....\\ ans\equiv c_q(mod\, m_q) \end{cases}$$
的最小正整数解，解为 $$(\sum_{i=1}^q*c_i*\frac{M}{m_i}*inv(\frac{M}{m_i},m_i))\%M$$ 其中 $$M=\prod_{1}^{q} m_i$$

题解：

对于p，将它表示成

$$p=p_1^{k_1}*p_2^{k_2}...p_q^{k_q}$$ 的形式，对于每个 $p_i$，用lucas定理求解，再用中国剩余定理求和。

唯一的问题在于在求

$$\frac{n!}{(n-m)!m!}\,\%p_i^{k_i}$$ 的时候,由于分母与 $p_i^{k_i}$不一定互质，没法求逆元，所以需要先把分子分母中 $p_i$全部除去，最后再乘回来。
代码实现：

long long x=0,ans;
for(long long i=n;i;i/=pi) x+=i/pi;
for(long long i=k;i;i/=pi) x-=i/pi;
for(long long i=n-k;i;i/=pi) x-=i/pi;
ans=a*mod_inverser(b,pk)%pk*mod_inverser(c,pk)*pow_mod(pi,x,pk)%pk;

解释一下$\frac{i}{p_i}$的含义：
以求$n!$为例：假设$n=100,p_i=2$
$100!=100*99*98*...*1$
$=(99*97*95*...*1)*2^{50}*(1*2*3*...*50)$
$=(99*97*95*...*1)*2^{50}*(1*3*...*49)*2^{25}*(1*2*...*25)$
……
所以，一直除$p_i$可将$n!$表示为

$$p_i^{x}*a$$ a为其余的项的积。通过这个方法可求得 $\frac{n!}{(n-m)!m!}\,\%p_i^{k_i}$

代码：

#include <iostream>
#include <bits/stdc++.h>
using namespace std;

long long exgcd(long long a,long long b,long long &x,long long &y)
{
    long long d=a;
    if(b!=0)
    {
        d=exgcd(b,a%b,y,x);
        y-=(a/b)*x;
    }else{
        x=1;
        y=0;
    }
    return d;
}
long long pow_mod(long long a,long long b,long long p)
{
    long long res=1;
    while(b)
    {
        if(b&1) res=res*a%p;
        a=a*a%p;
        b>>=1;
    }
    return res;
}
long long mod_inverser(long long a,long long m)
{
    long long x,y;
    exgcd(a,m,x,y);
    return (m+x%m)%m;
}
long long cal(long long n,long long pi,long long pk)
{
    if(n==0) return 1;
    long long ans=1;
    if(n/pk)
    {
         for(long long i=2;i<=pk;i++)
         {
             if(i%pi)
                ans=ans*i%pk;
         }
         ans=pow_mod(ans,n/pk,pk);
    }
    for(long long i=2;i<=n%pk;i++)
    {
        if(i%pi)
            ans=ans*i%pk;
    }
    return ans*cal(n/pi,pi,pk)%pk;
}
long long C(long long n,long long k,long long p,long long pi,long long pk)
{
    if(k>n) return 0;
    long long a=cal(n,pi,pk),b=cal(k,pi,pk),c=cal(n-k,pi,pk);
    long long x=0,ans;
    for(long long i=n;i;i/=pi) x+=i/pi;
    for(long long i=k;i;i/=pi) x-=i/pi;
    for(long long i=n-k;i;i/=pi) x-=i/pi;
    ans=a*mod_inverser(b,pk)%pk*mod_inverser(c,pk)*pow_mod(pi,x,pk)%pk;
    return ans*(p/pk)%p*mod_inverser(p/pk,pk)%p;
}

int main()
{
    long long n,k,p;
    cin>>n>>k>>p;
    long long ans=0;
    for(long long x=p,i=2;i<=p;i++)
    {
        if(x%i==0)
        {
            long long pk=1;
            while(x%i==0)
            {
                pk*=i;
                x/=i;
            }
            ans=(ans+C(n,k,p,i,pk))%p;
        }
    }
    cout<<ans<<endl;
    return 0;
}