Lucas定理用于求C(n,m)%p,(p为质数),而拓展Lucas定理中,p不一定为质数。
Lucas定理有两种形式,在线算法和离线打表,在线算法适用于p经常改变 或 查询次数较少 的情况,其他情况下用离线打表。
首先是一个在线算法。
1 ll pow(ll a, ll b, ll m) 2 { 3 ll ans = 1; 4 a %= m; 5 while(b) 6 { 7 if(b & 1)ans = (ans % m) * (a % m) % m; 8 b /= 2; 9 a = (a % m) * (a % m) % m; 10 } 11 ans %= m; 12 return ans; 13 } 14 ll inv(ll x, ll p)//x关于p的逆元,p为素数 15 { 16 return pow(x, p - 2, p); 17 } 18 ll C(ll n, ll m, ll p)//组合数C(n, m) % p 19 { 20 if(m > n)return 0; 21 ll up = 1, down = 1;//分子分母; 22 for(int i = n - m + 1; i <= n; i++)up = up * i % p; 23 for(int i = 1; i <= m; i++)down = down * i % p; 24 return up * inv(down, p) % p; 25 } 26 ll Lucas(ll n, ll m, ll p) 27 { 28 if(m == 0)return 1; 29 return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p; 30 }
然后放一个离线打表的。
1 const int maxn = 1e5 + 10; 2 ll fac[maxn];//阶乘打表 3 void init(ll p)//此处的p应该小于1e5,这样Lucas定理才适用 4 { 5 fac[0] = 1; 6 for(int i = 1; i <= p; i++) 7 fac[i] = fac[i - 1] * i % p; 8 } 9 ll pow(ll a, ll b, ll m) 10 { 11 ll ans = 1; 12 a %= m; 13 while(b) 14 { 15 if(b & 1)ans = (ans % m) * (a % m) % m; 16 b /= 2; 17 a = (a % m) * (a % m) % m; 18 } 19 ans %= m; 20 return ans; 21 } 22 ll inv(ll x, ll p)//x关于p的逆元,p为素数 23 { 24 return pow(x, p - 2, p); 25 } 26 ll C(ll n, ll m, ll p)//组合数C(n, m) % p 27 { 28 if(m > n)return 0; 29 return fac[n] * inv(fac[m] * fac[n - m], p) % p; 30 } 31 ll Lucas(ll n, ll m, ll p) 32 { 33 if(m == 0)return 1; 34 return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p; 35 }
但是除了这个我还有一个板子,感觉上应该更快,也没发现什么错的,也放上来吧。
1 #include <bits/stdc++.h> 2 #include<stdint.h> 3 using namespace std; 4 #define int long long 5 #define scan(n) scanf("%lld", &(n)) 6 #define scann(n, m) scanf("%lld%lld", &(n), &(m)) 7 #define scannn(a, b, c) scanf("%lld%lld%lld", &(a), &(b), &(c)) 8 #define prin(n) printf("%lld", (n)) 9 #define pb push_back 10 #define mp make_pair 11 #define ms(a) memset(a, 0, sizeof(a)) 12 #define fo(i, a, b) for (int i = (a); i <= (b); i++) 13 #define ro(i, a, b) for (int i = (a); i >= (b); i--) 14 #define dbg(args...) do {cout << #args << " : "<< args << endl;} 15 const int inf = 0x3f3f3f3f; 16 const int maxn = 1e5+100; 17 const int mod=1e9+7; 18 int f[maxn],inv[maxn]; 19 int qpow(int a,int b){ 20 int res=1; 21 while(b){ 22 if(b&1)res=res*a%mod; 23 a=a*a%mod;; 24 b>>=1; 25 } 26 return res; 27 } 28 void init(){ 29 f[0]=1; 30 for(int i=1;i<=maxn-1;i++)f[i]=f[i-1]*i%mod; 31 inv[maxn-1]=qpow(f[maxn-1],mod-2); 32 for(int i=maxn-2;i>=0;i--)inv[i]=inv[i+1]*(i+1)%mod; 33 } 34 int C(int n,int m){ 35 return f[n]*inv[m]%mod*inv[n-m]%mod; 36 } 37 int32_t main() { 38 int f,w,h; 39 cin>>f>>w>>h; 40 int fm=0,fz=0; 41 init(); 42 for(int i=1;i<=w;i++){ 43 fm+=(C(w-1,i-1)*C(f+1,i))%mod; 44 fm%=mod; 45 } 46 for(int i=1;i<=w/(h+1);i++){ 47 fz=fz+C(w-i*h-1,i-1)*C(f+1,i)%mod; 48 fz%=mod; 49 } 50 int ans=fz*qpow(fm,mod-2)%mod; 51 cout<<ans<<endl; 52 return 0; 53 }
拓展Lucas定理,只有在线算法。
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 const int maxn = 1e6 + 10; 5 const int mod = 1e9 + 7; 6 ll pow(ll a, ll b, ll m) 7 { 8 ll ans = 1; 9 a %= m; 10 while(b) 11 { 12 if(b & 1)ans = (ans % m) * (a % m) % m; 13 b /= 2; 14 a = (a % m) * (a % m) % m; 15 } 16 ans %= m; 17 return ans; 18 } 19 ll extgcd(ll a, ll b, ll& x, ll& y) 20 //求解ax+by=gcd(a, b) 21 //返回值为gcd(a, b) 22 { 23 ll d = a; 24 if(b) 25 { 26 d = extgcd(b, a % b, y, x); 27 y -= (a / b) * x; 28 } 29 else x = 1, y = 0; 30 return d; 31 } 32 ll mod_inverse(ll a, ll m) 33 //求解a关于模上m的逆元 34 //返回-1表示逆元不存在 35 { 36 ll x, y; 37 ll d = extgcd(a, m, x, y); 38 return d == 1 ? (m + x % m) % m : -1; 39 } 40 41 ll Mul(ll n, ll pi, ll pk)//计算n! mod pk的部分值 pk为pi的ki次方 42 //算出的答案不包括pi的幂的那一部分 43 { 44 if(!n)return 1; 45 ll ans = 1; 46 if(n / pk) 47 { 48 for(ll i = 2; i <= pk; i++) //求出循环节乘积 49 if(i % pi)ans = ans * i % pk; 50 ans = pow(ans, n / pk, pk); //循环节次数为n / pk 51 } 52 for(ll i = 2; i <= n % pk; i++) 53 if(i % pi)ans = ans * i % pk; 54 return ans * Mul(n / pi, pi, pk) % pk;//递归求解 55 } 56 57 ll C(ll n, ll m, ll p, ll pi, ll pk)//计算组合数C(n, m) mod pk的值 pk为pi的ki次方 58 { 59 if(m > n)return 0; 60 ll a = Mul(n, pi, pk), b = Mul(m, pi, pk), c = Mul(n - m, pi, pk); 61 ll k = 0, ans;//k为pi的幂值 62 for(ll i = n; i; i /= pi)k += i / pi; 63 for(ll i = m; i; i /= pi)k -= i / pi; 64 for(ll i = n - m; i; i /= pi)k -= i / pi; 65 ans = a * mod_inverse(b, pk) % pk * mod_inverse(c, pk) % pk * pow(pi, k, pk) % pk;//ans就是n! mod pk的值 66 ans = ans * (p / pk) % p * mod_inverse(p / pk, pk) % p;//此时用剩余定理合并解 67 return ans; 68 } 69 70 ll Lucas(ll n, ll m, ll p) 71 { 72 ll x = p; 73 ll ans = 0; 74 for(ll i = 2; i <= p; i++) 75 { 76 if(x % i == 0) 77 { 78 ll pk = 1; 79 while(x % i == 0)pk *= i, x /= i; 80 ans = (ans + C(n, m, p, i, pk)) % p; 81 } 82 } 83 return ans; 84 } 85 86 int main() 87 { 88 ll n, m, p; 89 while(cin >> n >> m >> p) 90 { 91 cout<<Lucas(n, m, p)<<endl; 92 } 93 return 0; 94 }