『Lucas定理以及拓展Lucas』
Lucas定理
在『組合數學基礎』中,我們已經提出了\(Lucas\)定理,并給出了\(Lucas\)定理的證明,本文僅將簡單回顧,并給出代碼。
\(Lucas\)定理:當\(p\)為質數時,\(C_n^m\equiv C_{n\ mod\ p}^{m\ mod\ p}*C_{n/p}^{m/p}(mod\ p)\)。
在計算模域組合數時,如果模數較小,那么就可以嘗試使用\(Lucas\)定理來遞歸求解,其時間復雜度為\(O(plog_p\min(n,m))\)。
\(Code:\)
inline long long Lucas(long long u,long long d) {if (u<=Mod&&d<=Mod)return C(u,d) % Mod;else return Lucas(u/Mod,d/Mod) * C(u%Mod,d%Mod) % Mod; }本文給出了一種遞歸寫法,由于取模操作的存在,在組合數計算函數\(C\)內也需要進行一些特判,避免一些無意義式子的計算,影響答案。
\(Code:\)
inline long long C(long long u,long long d) {if ( u>d || u<0 || d<0 )return 0LL;//... }關于組合數的計算,我們也有多種方法,常用的幾種如下:
\(1.\) 當\(n\),\(m\)都比較小的時候,可以根據組合數的階乘計算公式來預處理\(n\),\(m\)范圍以內的階乘以及階乘的逆元,然后根據詢問\(O(1)\)地回答組合數。
\(2.\) 當\(n\)的范圍較大,\(m\)的范圍較小時,由于不適合預處理,所以我們可以直接利用有關\(m\)組合數計算式\(C_n^m=\frac{n*(n-1)*...*(n-m+1)}{m*(m-1)*...*1}\)來模擬計算,同理,分母利用逆元操作即可,時間復雜度\(O(m)\)。
Exlucas算法
與之對應的,\(Exlucas\)算法是一種用來解決與\(Lucas\)定理形式很像但更具有一般性的問題的算法。為什么不叫\(Exlucas\)定理而叫\(Exlucas\)算法呢,是因為這個算法和\(Lucas\)定理沒什么關系,只是一個數論算法罷了。
對于求解模域組合數\(C_n^m\%p\),當\(p\)不一定是質數時,可以使用\(Exlucas\)算法求解。
對于這樣一個問題,我們考慮對模數進行分解,設\(p=\prod_{i=1}^kp_i^{a_i}\),答案\(x=C_n^m\%p\),則可以得到\[x\equiv C_n^ m(mod\ p_i^{a_i})\]
這是一個線性同余方程組,可以使用中國剩余定理求解。
那么現在我們的問題就轉換成了求解\(C_n^m\%p_i^{a_i}\)。將組合數寫為階乘的計算式,即\(C_n^m=\frac{n!}{m!(n-m)!}\)。
因為\(\frac{n!}{m!(n-m)!}\)與\(p_i^{a_i}\)不可避免的會存在公約數,所以不能直接使用\(Exeuclid\)等算法來求解逆元,那么我們就要考慮把\(\frac{n!}{m!(n-m)!}\)中有關\(p_i\)的項都提取掉,利用\(\frac{n!}{m!(n-m)!}\)是整數這一特點將分子該項的指數直接減掉分母該項的指數,再求解其余部分以及逆元,最后就能避免不互質的情況從而求解答案。
先考慮一個簡單的問題,如何計算一個形如\(n!\)的數中含有多少個因子\(p\),則和階乘分解一題類似,若設\(f(n)\)表示\(n!\)中含有因子\(p\)的個數,則有:
\[f(n)=\begin{cases}0\ (n<p)\\\lfloor \frac{n}{p} \rfloor +f(\lfloor \frac{n}{p} \rfloor)\ (n\geq p)\end{cases}\]
簡單地利用\(for\)循環求解即可。
求得階乘中因子\(p_i\)的數量后,我們的問題就變成為如何求解\(n!\)中除去所有\(p_i\)后取模\(p_i^{a_i}\)的值。對于每一項,這又需要我們分兩種情況考慮:
\(1.\) 該項是除去\(p_i\)后得到的,我們將所有這樣的項放在一起發現,發現除去\(p_i\)后,剩下的系數乘積也是一個階乘,使用遞歸求解。
\(2.\) 該項是階乘中原本的一項,除去情況\(1.\),我們發現階乘中若干個連續的項在模\(p_i^{a_i}\)意義下恰好構成了循環節,其長度不超過\(p_i^{a_i}\),先暴力計算一個循環節,然后快速冪即可。對于不包括在完整循環節中的項,由于其長度小于循環節,可以暴力計算。
然后就能將求得的答案通過計算逆元求解,最后,乘上提取出的\(p_i\)剩余的若干次方即可。
考慮一個簡單的例子就能加深理解:
求解\(n!\ mod\ p_i^{p_i}\),此時,\(n=19\),\(p_i=3\),\(a_i=2\)。
\[ n!=1*2*3*...*19\\=(1*2*4*5*7*8*10*11*13*14*16*17*19)*(3^6*6!) \\=(1*2*4*5*7*8)^2*19*(3^6*6!) \]
可以證明,其時間復雜度與\(O(plog_2p)\)同級。
\(Code:\)
#include <bits/stdc++.h> using namespace std; const int SIZE=300; long long A,B,Mod,m[SIZE],r[SIZE]; inline void input(void) {scanf("%lld%lld%lld",&A,&B,&Mod);//計算C(A,B)%Mod } inline long long power(long long a,long long b,long long p) {long long res = 1;while (b){if (1&b)res = res * a % p;b >>= 1;a = a * a % p;}return res; } inline long long Exeuclid(long long a,long long &x,long long b,long long &y,long long c) {if (!b){x=c/a,y=0;return a;}else{long long p = Exeuclid(b,x,a%b,y,c);long long x_ = x , y_ = y;x = y_ , y = x_ - a/b * y_;return p;} } inline long long inv(long long a,long long p) {long long x_,y_;Exeuclid(a,x_,p,y_,1);return (x_+p)%p; } inline long long calc(long long x,long long val,long long p) {//計算x!中除去val后模p意義下的值if (!x)return 1;long long res = 1 , last = x%p;for (long long i=1;i<=p;i++)if (i%val)res = res * i % p;res = power(res,x/p,p);for (long long i=1;i<=last;i++)if (i%val)res = res * i % p;return res * calc(x/val,val,p) % p; } inline long long C(long long d,long long u,long long val,long long Pow) {long long mulup=calc(d,val,Pow) , k=0;long long muldown1=calc(u,val,Pow) , muldown2=calc(d-u,val,Pow);for (long long i=d;i;i/=val)k+=i/val;for (long long i=u;i;i/=val)k-=i/val;for (long long i=d-u;i;i/=val)k-=i/val;//計算剩余的val的指數return mulup * inv(muldown1,Pow) % Pow * inv(muldown2,Pow) % Pow * power(val,k,Pow) % Pow; } inline long long CRT(int cnt) {long long m_ = Mod , M[SIZE] = {} , t[SIZE] = {} , res = 0;for (int i=1;i<=cnt;i++)M[i] = m_ / m[i];for (int i=1;i<=cnt;i++){long long y;Exeuclid(M[i],t[i],m[i],y,1);res = (res + r[i] % m_ * M[i] % m_ * t[i] % m_) % m_;}return ( res % m_ + m_ ) % m_; } inline long long Exlucas(void) {int temp = Mod , cnt = 0 ;for (long long i=2;i*i<=Mod;i++){if (temp%i==0){m[++cnt] = 1;while (temp%i==0)temp /= i , m[cnt] *= i;r[cnt] = C(A,B,i,m[cnt]);}}if (temp>1) m[++cnt] = temp , r[cnt] = C(A,B,temp,m[cnt]);return CRT(cnt); } int main(void) {input();printf("%lld\n",Exlucas());return 0; }轉載于:https://www.cnblogs.com/Parsnip/p/10764776.html
總結
以上是生活随笔為你收集整理的『Lucas定理以及拓展Lucas』的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 关于页面无法实现高度100%的原因及实现
- 下一篇: 32.C#--方法中使用out参数做登录