本身當模數為質數時，可以直接使用盧卡斯定理：

\[\begin{pmatrix}a\\ b\end{pmatrix} \equiv \begin{pmatrix}\left\lfloor \frac{a}{p} \right\rfloor \\ \left\lfloor\frac{b}{p}\right\rfloor\end{pmatrix} \times \begin{pmatrix}a\mod p\\ b\mod p\end{pmatrix} \pmod p \]

然後今天做到一道模數非質數的題，需要擴盧。

[AH2017/HNOI2017]拋硬幣

套娃

學了下發現並不是很難。

首先可以將模數進行質因數拆分，然後列出同餘方程組，將右邊項求得便可以用 CRT/EXCRT 求解。

\[\begin{cases} \begin{pmatrix}a\\ b\end{pmatrix}\mod {p_1}^{x_1} =s_1 \\ \begin{pmatrix}a\\ b\end{pmatrix}\mod {p_2}^{x_2} =s_2 \\ \cdots \end{cases} \]

然後我們就需要去求每一個 \(s_i\)。

式子過幾天自己推一下，防止遺忘，今天就咕咕咕。

（今天因為調模板題調了一晚上，已經把式子刻在 DNA 裡了）

調了一個晚上，結果發現快速冪寫成了這個鬼東西：

if(x>=pk){sum=sl[pk];for(LL i=x/pk,s=sl[pk];i;i>>=1){if(i&1)sum=sum*s%pk;s=s*s%pk;}}
 

wssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssbwssb

#include <stdio.h>
#define LL long long
using namespace std;
#define sort random_shuffle
const int P=1e6+3;
const int D=8e4+3;
const int G=23;
int p,k,pk;
LL n,m;
inline LL rin()
{
    LL s=0;
    bool bj=false;
    char c=getchar();
    for(;(c>'9'||c<'0')&&c!='-';c=getchar());
    if(c=='-')bj=true,c=getchar();
    for(;c>='0'&&c<='9';c=getchar())s=(s<<1)+(s<<3)+(c^'0');
    if(bj)s=-s;
    return s;
}

struct gyq
{
    int x,k;
    int pk;
    LL s;
}a[G];
int nam;

bool vit[P];
int d[D];
int tail;
inline void prime_()
{
    int i,j;
    for(i=2;i<=p;i++)
    {
        if(!vit[i])d[++tail]=i;
        for(j=1;j<=tail;j++)
        {
            int now=i*d[j];
            if(now>p)break;
            vit[now]=true;
            if(i%d[j]==0)break;
        }
    }
    for(i=1;i<=tail;i++)
    {
        if(p%d[i]==0)
        {
            a[++nam].x=d[i];
            for(;p%d[i]==0;p/=d[i])a[nam].k++;
            if(p==1)break;
        }
    }
    return;
}

LL sl[P];
LL cuts=0;
inline LL f(LL x)
{
    if(x<p)return sl[x];
    cuts+=x/p;
    LL sum=1;
    if(x>=pk){sum=1;for(LL i=x/pk,s=sl[pk];i;i>>=1){if(i&1)sum=sum*s%pk;s=s*s%pk;}}
    for(LL i=(x/pk)*pk+1;i<=x;i++)if(i%p)sum=sum*(i%pk)%pk;
    sum=sum*f(x/p)%pk;
    return sum;
}
inline void Ex_Gcd(LL a,LL b,LL &x,LL &y){(!b)?(x=y=1):(Ex_Gcd(b,a%b,y,x),y-=x*(a/b));return;}
inline void Ex_Lucus()
{
    sl[0]=sl[1]=1;
    for(int i=1;i<=nam;i++)
    {
        p=a[i].x;k=a[i].k;pk=1;
        for(int j=k;j;j--)pk*=p;a[i].pk=pk;
        for(int j=2;j<=pk;j++){sl[j]=sl[j-1];if(j%p)sl[j]=sl[j]*j%pk;}
        LL cutt=0;
        LL sr_1,sr_2;
        cuts=0;a[i].s=f(n);cutt+=cuts;
        cuts=0;sr_1=f(m);cutt-=cuts;
        cuts=0;sr_2=f(n-m);cutt-=cuts;
        Ex_Gcd(sr_1,pk,sr_1,cuts);(sr_1=sr_1%pk+pk)%pk;
        Ex_Gcd(sr_2,pk,sr_2,cuts);(sr_2=sr_2%pk+pk)%pk;
        a[i].s=(a[i].s*sr_1%pk)*sr_2%pk;
        a[i].s=(a[i].s%pk+pk)%pk;
        if(cutt>=k)a[i].s=0;
        else for(LL s=p;cutt;cutt>>=1){if(cutt&1)a[i].s=a[i].s*s%pk;s=s*s%pk;}
    }
    return;
}
inline void Ex_Crt()
{
    LL ans=a[1].s;
    LL sum=a[1].pk;
    for(int i=2;i<=nam;i++)
    {
        Ex_Gcd(sum,a[i].pk,m,n);m=m*sum;
        sum*=a[i].pk;
        ans+=m*(a[i].s-ans)%sum;
        ans=(ans%sum+sum)%sum;
    }
    printf("%lld\n",ans);
}
int main()
{
    int i,j;
    n=rin();m=rin();p=rin();
    prime_();
    Ex_Lucus();
    Ex_Crt();
    return 0;
}