【模板】卢卡斯与拓展卢卡斯

xiaoxiao2025-11-08  2

LUCAS

( n m ) ≡ ( ⌊ n p ⌋ ⌊ m p ⌋ ) ⋅ ( n % p m % p )   m o d    p \dbinom n m\equiv \dbinom{\lfloor\dfrac{n}{p}\rfloor}{\lfloor\dfrac{m}{p}\rfloor}·\dbinom{n\% p}{m\%p} \ \mod p (mn)(pmpn)(m%pn%p) modp


EXLUCAS

设模数为 P ( P = p 1 k 1 p 2 k 2 . . . p n k n ) , p i P(P=p_1^{k_1}p_2^{k_2}...p_n^{k^n}),p_i P(P=p1k1p2k2...pnkn),pi为互不相同的质数。 首先分别求出 x i = ( n m ) m o d    p i k i x_i=\dbinom n m\mod p_i^{k_i} xi=(mn)modpiki。 问题便转换成了求解: x ≡ x i m o d    p i k i x\equiv x_i\mod p_i^{k_i} xximodpiki。套用CRT解决。

那么问题在于如何快速求出所有的 ( n m ) m o d    p i k i \dbinom n m\mod p_i^{k_i} (mn)modpiki

k i = 1 k_i=1 ki=1时,套用LUCAS。

k i > 1 k_i>1 ki>1时, ( n m ) = n ! m ! ( n − m ) ! \dbinom n m=\dfrac{n!}{m!(n-m)!} (mn)=m!(nm)!n!无法对质数的次幂快速取模,考虑分别对 n ! , m ! , ( n − m ) ! n!,m!,(n-m)! n!,m!,(nm)!进行拆分取模,最终答案即为 n ! × i n v ( m ! , p i k i ) × i n v ( ( n − m ) ! , p i k i ) n!\times inv(m!,p_i^{k_i})\times inv((n-m)!,p_i^{k_i}) n!×inv(m!,piki)×inv((nm)!,piki)( i n v ( x , p ) inv(x,p) inv(x,p)表示 x x x在模 p p p意义下的逆元,exgcd求解)。

拆分取模方式如下:

假设当前 p i = 3 , k i = 2 , P = 3 2 p_i=3,k_i=2,P=3^2 pi=3,ki=2,P=32,求解 19 ! m o d    P 19!\mod P 19!modP

19 ! = ( 1 × 2 × 4 × 5 × 7 × . . . × 14 × 16 × 17 × 19 ) × ( 3 6 ) × ( 1 × 2 × 3 × 4 × 5 × 6 ) 19!=(1\times2\times4\times5\times7\times...\times14\times16\times17\times 19)\times(3^{6})\times(1\times2\times3\times4\times5\times6) 19!=(1×2×4×5×7×...×14×16×17×19)×(36)×(1×2×3×4×5×6)

对于第一部分 1 × 2 × 4 × 5 × 7 × 8 ≡ 10 × 11 × 13 × 14 × 16 × 17 m o d    P 1\times2\times 4\times 5\times 7\times 8 \equiv10\times11\times13\times14\times16\times17\mod P 1×2×4×5×7×810×11×13×14×16×17modP,等价于 ( ∏ i = 1 , i ∤ p i P − 1 i ) ⌊ n P ⌋ m o d    P (\prod\limits_{i=1,i\nmid p_i}^{P-1}i)^{\lfloor{\frac{n}{P}}\rfloor}\mod P (i=1,ipiP1i)PnmodP,剩余的数必然少于 P P P个,而将右边的 ⌊ n p i ⌋ ! \lfloor{\frac{n}{p_i}}\rfloor! pin!递归下去求解即可。


代码

#include<bits/stdc++.h> #define debug(x) cerr<<#x<<":"<<x<<endl using namespace std; typedef long long ll; const int N=5010; ll n,m,p,q[N],a[N]; int tot; namespace exlucas{ inline int fp(ll x,ll y,ll p) { ll re=1; for(;y;y>>=1,x=x*x%p) if(y&1) re=re*x%p; return re; } ll exgcd(ll n,ll m,ll &x,ll &y) { if(!m) {x=1;y=0;return n;} ll re=exgcd(m,n%m,y,x);y-=(n/m)*x; return re; } inline ll inv(ll n,ll p) { ll x,y; exgcd(n,p,x,y); return (x%p+p)%p; } inline ll C(ll n,ll m,ll p) { if(m>n) return 0; ll re=1,x,y,z,i; for(i=1;i<=m;++i){ x=n-i+1;y=(int)inv(i,p); re=re*x%p*y%p; } return re; } ll lucas(ll n,ll m,ll p) { if(m==0 || n==m) return 1; return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p; } ll cal(ll n,ll a,ll b,ll p) { if(!n) return 1; ll i,j,re=1; for(i=1;i<p;++i) if(i%a) re=re*i%p; re=fp(re,n/p,p); for(i=n/p*p+1;i<=n;++i) if(i%a) re=re*i%p; return re*cal(n/a,a,b,p)%p; } inline ll multilucas(ll n,ll m,ll a,ll b,ll p) { ll i,s=0; for(i=n;i;i/=a) s+=i/a; for(i=m;i;i/=a) s-=i/a; for(i=n-m;i;i/=a) s-=i/a; return fp(a,s,p)*cal(n,a,b,p)%p*inv(cal(m,a,b,p),p)%p*inv(cal(n-m,a,b,p),p)%p; } inline ll cal(ll n,ll m,ll p) { ll i,j,b,c,d,ans,mod,x,y; for(i=2;i*i<=p;++i) if(p%i==0){ q[++tot]=1; for(j=0;p%i==0;p/=i) q[tot]*=i,j++; if(j==1) a[tot]=lucas(n,m,i); else a[tot]=multilucas(n,m,i,j,q[tot]); } if(p>1){q[++tot]=p;a[tot]=lucas(n,m,p);} ans=a[1];mod=q[1]; for(i=2;i<=tot;++i){ b=exgcd(mod,q[i],x,y); c=a[i]-ans; if(c%b) return -1LL; d=q[i]/b;x=((c/b)%d*x%d+d)%d; ans=ans+mod*x; mod=mod/b*q[i]; } return ans; } } int main(){ scanf("%lld%lld%lld",&n,&m,&p); printf("%lld\n",exlucas::cal(n,m,p)); return 0; }
转载请注明原文地址: https://www.6miu.com/read-5039294.html

最新回复(0)