扩展Lucas定理

扩展 Lucas 定理 exLucas

其他题解讲的都很不清楚，而且很奇怪。。。

所以我来发一篇良心题解，求大家点赞~

问题

求

C_n^m \bmod{p}

其中p\leq 10^6

算法

设

p=p_1^{\alpha_1}p_2^{\alpha_2}...

求出

\left\{\begin{aligned}C_n^m \bmod\ {p_1^{\alpha_1}} \\C_n^m \bmod\ {p_2^{\alpha_2}}\\ ......\end{aligned} \right.

最后用中国剩余定理合并即可。

假设我们现在要求

C_n^m\bmod{P^k}([P\in {\texttt{prime}}]) \therefore\ =\frac {n!}{m!(n-m)!}\bmod{P^k}

我们无法求得m!,(n-m)!的逆元，理由很简单，不一定有解。

怎么办呢？发现a对p有逆元的充要条件为

\gcd(a,p)=1

所以我们换一个式子：

=\frac {\frac {n!}{P^{x}}}{\frac {m!}{P^{y}}\frac {(n-m)!}{P^{z}}}P^{x-y-z}\bmod{P^k}

其中

其他都是同理的。 所以如果我们对于一个$n$，可以求出$\frac {n!}{P^x}$，那不就完事了吗？（这样我们就可以求逆元了） 问题等价于求 $$\frac {n!}{P^x}\bmod{P^k}$$ 我们对$n!$进行变形： $$n!=1\cdot 2\cdot 3\cdot ...\cdot n$$ $$=(P\cdot 2P\cdot 3P...)(1\cdot 2\cdot ...)$$ 左边是$1\sim n$中所有$\leq n$的$P$的倍数。 右边反之。 因为$1\sim n$中有$\lfloor \frac nP\rfloor$个$P$的倍数， $$\therefore\ =P^{\lfloor \frac nP\rfloor}(1 \cdot 2 \cdot 3...)(1\cdot 2\cdot ...)$$ $$=P^{\lfloor \frac nP\rfloor}(\lfloor \frac nP\rfloor)!\prod_{i=1,i\not\equiv 0\pmod {P}}^ni$$ 显然后面这个$\bmod\ P$是有循环节的。 $$=P^{\lfloor \frac nP\rfloor}(\lfloor \frac nP\rfloor)!(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)$$ 其中 $$(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}$$ 是循环节$1\sim P$中所有无$P$因子数的乘积。 $$(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)$$ 是余项的乘积。 比如 $$22!\equiv(1\cdot 2\cdot 4\cdot 5\cdot 7\cdot 8)(10\cdot 11\cdot 13\cdot 14\cdot 16\cdot 17)(19\cdot 20\cdot 22)(3\cdot 6\cdot 9\cdot 12\cdot 15\cdot 18\cdot 21)\pmod {3^2}$$ $$\equiv(1\cdot 2\cdot 4\cdot 5\cdot 7\cdot 8)^2(19\cdot 20\cdot 22)3^7(1\cdot 2\cdot 3\cdot 4\cdot 5\cdot 6\cdot 7)\pmod {3^2}$$ $$\equiv3^77!(1\cdot 2\cdot 4\cdot 5\cdot 7\cdot 8)^2(19\cdot 20\cdot 22)\pmod {3^2}$$ 发现正好是一样的。$\lfloor \frac {22}{3^2}\rfloor=2

理解了吧...

=P^{\lfloor \frac nP\rfloor}(\lfloor \frac nP\rfloor)!(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)

发现前面这个P^{\lfloor \frac nP\rfloor}是要除掉的。

然而(\lfloor \frac nP\rfloor)!里可能还包含P。

所以，我们定义函数：

f(n)=\frac {n!}{P^x} \therefore\ f(n)=f(\lfloor \frac nP\rfloor)(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)

这样就好了。时间复杂度是O(\log_{P}n)。

边界f(0)=1

看看原式

=\frac {\frac {n!}{P^{x}}}{\frac {m!}{P^{y}}\frac {(n-m)!}{P^{z}}}P^{x-y-z}\bmod{P^k} =\frac {f(n)}{f(m)f(n-m)}P^{x-y-z}\bmod{P^k}

由于f(m)一定与P^k互质，所以我们可以直接用exgcd求逆元啦。

最后我们还有一个问题，P^{x-y-z}咋求?

其实很好求。

比如说，我们要求f(n)=\frac {n!}{P^x}中的x

设g(n)=x(上述的x)

再看看阶乘的式子

n!=P^{\lfloor \frac nP\rfloor}(\lfloor \frac nP\rfloor)!(\prod_{i=1,i\not\equiv 0\pmod {P}}^{P^k}i)^{\lfloor \frac n{P^k}\rfloor}(\prod_{i=P^k\lfloor \frac n{P^k}\rfloor,i\not\equiv 0\pmod {P}}^ni)

很显然我们要的就是P^{\lfloor \frac nP\rfloor}

而且(\lfloor \frac nP\rfloor)!可能还有P因子。

所以我们可以得到以下递推式

g(n)=\lfloor \frac nP\rfloor+g(\lfloor \frac nP\rfloor)

时间复杂度是O(\log_{P}n)。

边界g(n)=0(n<P)

所以答案就是

=\frac {\frac {n!}{P^{x}}}{\frac {m!}{P^{y}}\frac {(n-m)!}{P^{z}}}P^{g(n)-g(m)-g(n-m)}\bmod{P^k}

最后用中国剩余定理合并答案即可得到C_n^m。

代码:

#include<bits/stdc++.h>
#define ll long long
using namespace std;
#ifndef Fading
inline char gc(){
    static char now[1<<16],*S,*T;
    if (T==S){T=(S=now)+fread(now,1,1<<16,stdin);if (T==S) return EOF;}
    return *S++;
}
#endif
#ifdef Fading
#define gc getchar
#endif
void exgcd(ll a,ll b,ll &x,ll &y){
    if (!b) return (void)(x=1,y=0);
    exgcd(b,a%b,x,y);
    ll tmp=x;x=y;y=tmp-a/b*y;
}
ll gcd(ll a,ll b){
    if (b==0) return a;
    return gcd(b,a%b); 
}
inline ll INV(ll a,ll p){
    ll x,y;
    exgcd(a,p,x,y);
    return (x+p)%p;
}
inline ll lcm(ll a,ll b){
    return a/gcd(a,b)*b;
}
inline ll mabs(ll x){
    return (x>0?x:-x);
}
inline ll fast_mul(ll a,ll b,ll p){
    ll t=0;a%=p;b%=p;
    while (b){
        if (b&1LL) t=(t+a)%p;
        b>>=1LL;a=(a+a)%p;
    }
    return t;
}
inline ll fast_pow(ll a,ll b,ll p){
    ll t=1;a%=p;
    while (b){
        if (b&1LL) t=(t*a)%p;
        b>>=1LL;a=(a*a)%p;
    }
    return t;
}
inline ll read(){
    ll x=0,f=1;char ch=gc();
    while (!isdigit(ch)) {if (ch=='-') f=-1;ch=gc();}
    while (isdigit(ch)) x=x*10+ch-'0',ch=gc();
    return x*f;
}
inline ll F(ll n,ll P,ll PK){
    if (n==0) return 1;
    ll rou=1;//循环节
    ll rem=1;//余项 
    for (ll i=1;i<=PK;i++){
        if (i%P) rou=rou*i%PK;
    }
    rou=fast_pow(rou,n/PK,PK);
    for (ll i=PK*(n/PK);i<=n;i++){
        if (i%P) rem=rem*(i%PK)%PK;
    }
    return F(n/P,P,PK)*rou%PK*rem%PK;
}
inline ll G(ll n,ll P){
    if (n<P) return 0;
    return G(n/P,P)+(n/P);
}
inline ll C_PK(ll n,ll m,ll P,ll PK){
    ll fz=F(n,P,PK),fm1=INV(F(m,P,PK),PK),fm2=INV(F(n-m,P,PK),PK);
    ll mi=fast_pow(P,G(n,P)-G(m,P)-G(n-m,P),PK);
    return fz*fm1%PK*fm2%PK*mi%PK;
}
ll A[1001],B[1001];
//x=B(mod A)
inline ll exLucas(ll n,ll m,ll P){
    ll ljc=P,tot=0;
    for (ll tmp=2;tmp*tmp<=P;tmp++){
        if (!(ljc%tmp)){
            ll PK=1;
            while (!(ljc%tmp)){
                PK*=tmp;ljc/=tmp;
            }
            A[++tot]=PK;B[tot]=C_PK(n,m,tmp,PK);
        }
    }
    if (ljc!=1){
        A[++tot]=ljc;B[tot]=C_PK(n,m,ljc,ljc);
    }
    ll ans=0;
    for (ll i=1;i<=tot;i++){
        ll M=P/A[i],T=INV(M,A[i]);
        ans=(ans+B[i]*M%P*T%P)%P;
    }
    return ans;
}
signed main(){
    ll n=read(),m=read(),P=read();
    printf("%lld\n",exLucas(n,m,P));
    return 0;
}