Luogu4720 【模板】扩展卢卡斯

Luogu4720 【模板】扩展卢卡斯

设\(p=p_1^{a_1}p_2^{a_2}\cdots\)。

求出每一个\({n \choose m} \mod p_i^{a_i}\)的值，进行\(CRT\)合并。

求解\(\frac{n!}{m!(n-m)!} \mod p^k\)。

然而\(\frac{n!}{m!(n-m)!}\)不一定与\(p^k\)互质，无法直接求逆元，因此我们可以尝试去转化，使它们互质。

\[\large \frac{\frac{n!}{p^x}}{\frac{m!}{p^y} \frac{(n-m)!}{p^z}} p^{x-y-z} \mod p^k \]

\[n!\\ =p^{\lfloor \frac{n}{p} \rfloor} (\lfloor \frac{n}{p} \rfloor)! (\prod_{i=1,i \equiv k(k \ne 0) \bmod p}^{p^k-1} i)^{\lfloor \frac{n}{p^k} \rfloor } \prod_{i=\lfloor \frac{n}{p^k} \rfloor p^k + 1,i \equiv k(k \ne 0)}^n i \]

做完哩！
\((\prod_{i=1,i \equiv k(k \ne 0) \bmod p}^{p^k-1} i)^{\lfloor \frac{n}{p^k} \rfloor } \prod_{i=\lfloor \frac{n}{p^k} \rfloor p^k + 1,i \equiv k(k \ne 0)}^n i\)可以\(O(p^k)\)求解，但是我们需要除去\((\lfloor \frac{n}{p} \rfloor)!\)中的\(p\)因子，递归求解。

根据上面的式子，当然也可以轻松统计\(p\)因子个数。

\(Code:\)

#include<iostream>
#include<cstdio>
#include<algorithm>
#define ll long long
#define ll long long
using namespace std;
ll n,m,zs;
int cnt,p,x,y,a[15],b[15],c[15],d[15];
ll fcnt=0;
void exgcd(int &x,int &y,int a,int b)
{
    if (!b)
    {
        x=1,y=0;
        return;
    }
    exgcd(x,y,b,a%b);
    ll z=x;
    x=y;
    y=z-a/b*y;
}
ll ksm(ll x,ll y,ll p)
{
    ll ans=1;
    while (y)
    {
        if (y & 1)
            ans=1LL*ans*x%p;
        x=1LL*x*x%p;
        y >>=1;
    }
    return ans;
}
int cutp(ll n,int p,int pk)
{
    int ans=1;
    if (n<p)
    {
        for (int i=1;i<=n;++i)
            ans=1LL*ans*i%pk;
        return ans;
    }
    fcnt+=n/p;
    ans=1LL*ans*ksm(zs,n/pk,pk)%pk;
    for (ll i=n/pk*pk+1;i<=n;++i)
        if (i % p)
            ans=1LL*ans*(i%pk)%pk;
    return 1LL*ans*cutp(n/p,p,pk)%pk;
}
int solve(ll p,ll k,ll pk)
{
    int x,y,z;
    ll xx,yy,zz;
    zs=1;
    for (int i=1;i<pk;++i)
        if (i % p)
            zs=1LL*zs*i%pk;
    fcnt=0,x=cutp(n,p,pk),xx=fcnt;
    fcnt=0,y=cutp(m,p,pk),yy=fcnt;
    fcnt=0,z=cutp(n-m,p,pk),zz=fcnt;
    ll mi=xx-yy-zz;
    if (mi>=k)
        return 0;
    int e,f,e1,e2;
    exgcd(e,f,y,pk);
    e1=(e%pk+pk)%pk;
    exgcd(e,f,z,pk);
    e2=(e%pk+pk)%pk;
    return 1LL*x*e1%pk*e2%pk*ksm(p,mi,pk)%pk;
}
int main()
{
    scanf("%lld%lld%d",&n,&m,&p);
    int zp=p;
    for (int i=2;i*i<=zp;++i)
        if (zp % i == 0)
        {
            a[++cnt]=i,c[cnt]=1;
            while (zp % i == 0)
                zp/=i,++b[cnt],c[cnt]*=i;
        }
    if (zp!=1)
        a[++cnt]=zp,b[cnt]=1,c[cnt]=zp;
    for (int i=1;i<=cnt;++i)
        d[i]=solve(a[i],b[i],c[i]);
    int ans=0;
    for (int i=1;i<=cnt;++i)
    {
        int T=p/c[i];
        exgcd(x,y,T,c[i]);
        x=(x%p+p)%p;
        ans=(ans+1LL*x*T%p*d[i]%p)%p;
    }
    ans=(ans%p+p)%p;
    printf("%d\n",ans);
    return 0;
}