码迷,mamicode.com
首页 > 其他好文 > 详细

组合数

时间:2018-12-09 16:41:15      阅读:363      评论:0      收藏:0      [点我收藏+]

标签:fine   clu   include   exgcd   质因数   ||   pre   char   col   

我们要求的呢,就是这个东西\[C_n^m\% p\]

这里P不一定是质数

1.将P质因数分解

\[\prod\limits_{\rm{i}}^r {P_{}^{{c_i}}} \]

对下面这个进行中国剩余定理合并,便是答案

\[\begin{array}{l}
\begin{array}{*{20}{l}}
{{x_1} \equiv C_n^m(\,\bmod \,p_1^{{c_1}})}\\
{{x_2} \equiv C_n^m(\,\bmod \,p_2^{{c_2}})}\\
{{x_3} \equiv C_n^m(\,\bmod \,p_3^{{c_3}})}
\end{array}\\
\cdot \cdot \cdot \\
{x_r} \equiv C_n^m(\,\bmod \,p_r^{{c_r}})
\end{array}\]

2.求$C_n^m\% {p^k}$

这个该怎么求呢,很明显这个求不了逆元,但是因为分母的两个数和模数不互质,那么只要把所有的模数因子全部干掉就好了

也就是这样\[C_n^m\% {p^k} = \frac{{\frac{{n!}}{{{p^{{a_1}}}}}}}{{\frac{{m!}}{{{p^{{a_2}}}}}*\frac{{(n - m)!}}{{{p^{{a_3}}}}}}}*{p^{{a_1} - {a_2} - {a_3}}}\% {p^k}\]

但是这里的a1a2a3是不知道的,没事我们看下一步就知道了

3.求$n!\% {p^k}$

我们来举个栗子

\[\begin{array}{l}
22!\% {3^2} = (1*2*3*4*5*6*7*8*9*10*11*12*13*14*14*16*17*18*19*20*21*22)\% {3^2}\\
= (3*6*9*12*15*18*21)*(1*2*4*5*7*8)*(10*11*13*14*16*17)*(19*20*22)\% {3^2}\\
= {3^7}*7!*{(1*2*4*5*7*8)^2}\% {3^2}
\end{array}\]

我们推广一下

\[n!\% {p^k} = {p^{n/p}}*(n/p)!*(\prod\limits_1^{ <  = {p^k}} {nu{m_i}} )\% {p^k}\]

看到$(n/p)!$,这个很明显可以递归来解决(注意边界是n==1||2时候返回值是1)

看到那个${p^{n/p}}$这不就是第二步里没解决的a值吗,对,就是n/p;

但是n/p很明显是递归的,这个应该怎么解决,写个循环每次都除就好了

就是这样

 1 lt fac(const lt n,const lt pi,const lt pk)//n! % pi^ki pk=pi^ki
 2 {
 3     if(n==1||n==0)return 1;
 4     lt ans=1;
 5     for(rt i=2;i<pk;i++)
 6         if(i%pi)
 7             ans=ans*i%pk;
 8     ans=ksm(ans,n/pk,pk);        
 9     for(rt i=2;i<=n%pk;i++)
10         if(i%pi)ans=ans*i%pk;
11     return ksc(ans,fac(n/pi,pi,pk),pk);
12 } 

那么第二步的a解决之后,第二步也很好办了

1 lt C(const lt n,const lt m,const lt pi,const lt pk)////C(n,m)%pk  pk=pi^ki
2 {
3     lt up=fac(n,pi,pk),dl=fac(m,pi,pk),dr=fac(n-m,pi,pk),kk=0;
4     for(rt i=n;i;i/=pi)kk+=i/pi;
5     for(rt i=m;i;i/=pi)kk-=i/pi;
6     for(rt i=n-m;i;i/=pi)kk-=i/pi;
7     return up*inv(dl,pk)%pk*inv(dr,pk)%pk*ksm(pi,kk,pk)%pk;
8 }

然后我进行了一些基本函数丧心病狂的封装23333

 1 namespace basic
 2 {
 3     inline lt read()
 4     {
 5         rt x = 0; char zf = 1; char ch;
 6         while (ch != - && !isdigit(ch)) ch = getchar();
 7         if (ch == -) zf = -1, ch = getchar();
 8         while (isdigit(ch)) x = x * 10 + ch - 0, ch = getchar(); return x * zf;
 9     }
10     void write(lt x)
11     {
12         if (x<0)putchar(-),x=-x;
13         if (x>9) write(x/10);
14         putchar(x%10+0);
15     }
16     inline void writeln(const lt x){write(x);putchar(\n);}
17     lt gcd(lt a,lt b){return b?gcd(b,a%b):a;}
18     lt ksc(lt a,lt b,lt mod)
19     {
20         lt fina=0,kk=1;
21         if(a<0)a=-a,kk=-kk;if(b<0)b=-b,kk=-kk;
22         while(b)
23         {
24             if(b%2)fina=(fina+a)%mod;
25             a=(a+a)%mod,b>>=1;
26         }
27         return fina%mod*kk;
28     }
29     lt ksm(lt a,lt k,lt mod)
30     {
31         if(!k)return 1;
32         lt fina=1;
33         while(k)
34         {
35             if(k%2)fina=ksc(fina,a,mod);
36             a=a*a%mod,k>>=1;
37         }
38         return fina%mod;
39     }
40     void ex_gcd(lt a,lt b,lt &x,lt &y)
41     {
42         if(!b)x=1,y=0;
43         else
44         {
45             ex_gcd(b,a%b,x,y);lt tt=x;x=y,y=tt-a/b*x;
46         }
47     }
48     lt exgcd(lt a,lt b,lt c)//ax=c(mod b)
49     {
50         lt gc=gcd(a,b),x=0,y=0;;
51         if(c%gc)return -1;
52         a/=gc,b/=gc,c/=gc;
53         ex_gcd(a,b,x,y);
54         return (x*c%b+b)%b;
55     }
56     inline lt inv(lt a,lt p)
57     {
58         if(!a)return 0;
59         return exgcd(a,p,1);
60     }//ax=1(mod p)
61 }

然后lucas部分

 1 namespace lucas
 2 {
 3     using namespace basic;
 4     lt fac(const lt n,const lt pi,const lt pk)//n! % pi^ki pk=pi^ki
 5     {
 6         if(n==1||n==0)return 1;
 7         lt ans=1;
 8         for(rt i=2;i<pk;i++)
 9             if(i%pi)
10                 ans=ans*i%pk;
11         ans=ksm(ans,n/pk,pk);        
12         for(rt i=2;i<=n%pk;i++)
13             if(i%pi)ans=ans*i%pk;
14         return ksc(ans,fac(n/pi,pi,pk),pk);
15     } 
16     lt C(const lt n,const lt m,const lt pi,const lt pk)////C(n,m)%pk  pk=pi^ki
17     {
18         lt up=fac(n,pi,pk),dl=fac(m,pi,pk),dr=fac(n-m,pi,pk),kk=0;
19         for(rt i=n;i;i/=pi)kk+=i/pi;
20         for(rt i=m;i;i/=pi)kk-=i/pi;
21         for(rt i=n-m;i;i/=pi)kk-=i/pi;
22         return up*inv(dl,pk)%pk*inv(dr,pk)%pk*ksm(pi,kk,pk)%pk;
23     }
24     void exlucas()
25     {
26         n=read(),m=read(),P=read();    
27         lt ans=0,p=P;
28         for(rt i=2;i<=p;i++)
29             if(p%i==0)
30             {
31                 lt pi=i,pk=1;
32                 while(!(p%i))p/=i,pk*=i;
33                 ans=(ans+C(n,m,pi,pk)*(P/pk)%P*inv(P/pk,pk)%P)%P;
34             }
35         writeln(ans%P);
36     }
37 }

完整代码

  1 #include<cstdio>
  2 #include<algorithm>
  3 #include<cstring>
  4 #include<cctype>
  5 typedef long long lt;
  6 #define rt register lt
  7 using namespace std;
  8 lt n,m,p,P;
  9 namespace basic
 10 {
 11     inline lt read()
 12     {
 13         rt x = 0; char zf = 1; char ch;
 14         while (ch != - && !isdigit(ch)) ch = getchar();
 15         if (ch == -) zf = -1, ch = getchar();
 16         while (isdigit(ch)) x = x * 10 + ch - 0, ch = getchar(); return x * zf;
 17     }
 18     void write(lt x)
 19     {
 20         if (x<0)putchar(-),x=-x;
 21         if (x>9) write(x/10);
 22         putchar(x%10+0);
 23     }
 24     inline void writeln(const lt x){write(x);putchar(\n);}
 25     lt gcd(lt a,lt b){return b?gcd(b,a%b):a;}
 26     lt ksc(lt a,lt b,lt mod)
 27     {
 28         lt fina=0,kk=1;
 29         if(a<0)a=-a,kk=-kk;if(b<0)b=-b,kk=-kk;
 30         while(b)
 31         {
 32             if(b%2)fina=(fina+a)%mod;
 33             a=(a+a)%mod,b>>=1;
 34         }
 35         return fina%mod*kk;
 36     }
 37     lt ksm(lt a,lt k,lt mod)
 38     {
 39         if(!k)return 1;
 40         lt fina=1;
 41         while(k)
 42         {
 43             if(k%2)fina=ksc(fina,a,mod);
 44             a=a*a%mod,k>>=1;
 45         }
 46         return fina%mod;
 47     }
 48     void ex_gcd(lt a,lt b,lt &x,lt &y)
 49     {
 50         if(!b)x=1,y=0;
 51         else
 52         {
 53             ex_gcd(b,a%b,x,y);lt tt=x;x=y,y=tt-a/b*x;
 54         }
 55     }
 56     lt exgcd(lt a,lt b,lt c)//ax=c(mod b)
 57     {
 58         lt gc=gcd(a,b),x=0,y=0;;
 59         if(c%gc)return -1;
 60         a/=gc,b/=gc,c/=gc;
 61         ex_gcd(a,b,x,y);
 62         return (x*c%b+b)%b;
 63     }
 64     inline lt inv(lt a,lt p)
 65     {
 66         if(!a)return 0;
 67         return exgcd(a,p,1);
 68     }//ax=1(mod p)
 69 }
 70 namespace lucas
 71 {
 72     using namespace basic;
 73     lt fac(const lt n,const lt pi,const lt pk)//n! % pi^ki pk=pi^ki
 74     {
 75         if(n==1||n==0)return 1;
 76         lt ans=1;
 77         for(rt i=2;i<pk;i++)
 78             if(i%pi)
 79                 ans=ans*i%pk;
 80         ans=ksm(ans,n/pk,pk);        
 81         for(rt i=2;i<=n%pk;i++)
 82             if(i%pi)ans=ans*i%pk;
 83         return ksc(ans,fac(n/pi,pi,pk),pk);
 84     } 
 85     lt C(const lt n,const lt m,const lt pi,const lt pk)////C(n,m)%pk  pk=pi^ki
 86     {
 87         lt up=fac(n,pi,pk),dl=fac(m,pi,pk),dr=fac(n-m,pi,pk),kk=0;
 88         for(rt i=n;i;i/=pi)kk+=i/pi;
 89         for(rt i=m;i;i/=pi)kk-=i/pi;
 90         for(rt i=n-m;i;i/=pi)kk-=i/pi;
 91         return up*inv(dl,pk)%pk*inv(dr,pk)%pk*ksm(pi,kk,pk)%pk;
 92     }
 93     void exlucas()
 94     {
 95         n=read(),m=read(),P=read();    
 96         lt ans=0,p=P;
 97         for(rt i=2;i<=p;i++)
 98             if(p%i==0)
 99             {
100                 lt pi=i,pk=1;
101                 while(!(p%i))p/=i,pk*=i;
102                 ans=(ans+C(n,m,pi,pk)*(P/pk)%P*inv(P/pk,pk)%P)%P;
103             }
104         writeln(ans%P);
105     }
106 }
107 int main(){lucas::exlucas();return 0;}

 

组合数

标签:fine   clu   include   exgcd   质因数   ||   pre   char   col   

原文地址:https://www.cnblogs.com/Qin-Wei-Kai/p/10090972.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!