bzoj 2142 : 礼物

求$C_{n}^{m}\%p$。

把p拆成$p1^{q1}*p2^{q2}...$最后用CRT合并。

把每个阶乘拆成$x*p^y$的形式,因为x与$p^q$互质,可以直接用Euler定理求逆元,y就直接减。

拆的时候把每个p的倍数提出一个p,变为$tmp*(p^x*(1*2*3*4...))$,tmp为与p互质的数,在$mod p^q$下有循环节,后半部分递归解决。

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<algorithm>
  4 #include<cstring>
  5 #define N 100005
  6 #define int long long
  7 #define pa pair<int,int>
  8 #define mp make_pair
  9 using namespace std;
 10 int xx,yy;
 11 bool flag;
 12 int pw(int x,int y,int p)
 13 {
 14     if(!x)return 0;
 15     x%=p;
 16     int lst=1;
 17     while(y)
 18     {
 19         if(y&1)lst=lst*x%p;
 20         y>>=1;
 21         x=x*x%p;
 22     }
 23     return lst;
 24 }
 25 void exgcd(int a,int b)
 26 {
 27     if(b==0)
 28     {
 29         xx=1;yy=0;
 30         return ;
 31     }
 32     exgcd(b,a%b);
 33     int tmp=xx;
 34     xx=yy;
 35     yy=tmp-(a/b)*yy;
 36     return ;
 37 }
 38 int ans[100],p,n,m,w[10],s[100],su[100],tot;
 39 int jie[15][N];
 40 pa fact(int k,int n)
 41 {
 42     if(!n)return mp(1,0);
 43     int x=n/s[k];int ass=1;
 44     int y=n/su[k];
 45     if(y)
 46     {
 47         for(int i=2;i<su[k];i++)
 48         {
 49             if(i%s[k]!=0)ass=ass*i%su[k];
 50         }
 51         ass=pw(ass,y,su[k]);
 52     }
 53     for(int i=y*su[k]+1;i<=n;i++)
 54     {
 55         if(i%s[k]!=0)ass=ass*i%su[k];
 56     }
 57     pa tmp=fact(k,x);
 58     return mp(ass*tmp.first%su[k],x+tmp.second);
 59 }
 60 int qiu(int x,int y)
 61 {
 62     exgcd(x,y);
 63     return (x%y+y)%y;
 64 }    
 65 int lus(int n,int m,int x)
 66 {
 67     if(m>n)return 0;
 68     pa a=fact(x,n),b=fact(x,m),c=fact(x,n-m);
 69     int num=a.second-b.second-c.second;
 70     int q=a.first*pw(b.first,su[x]/s[x]*(s[x]-1)-1,su[x])%su[x]*pw(c.first,su[x]/s[x]*(s[x]-1)-1,su[x])%su[x];
 71     return pw(s[x],num,su[x])*q%su[x];
 72 }
 73 void div(int x)
 74 {
 75     for(int i=2;i*i<=x;i++)
 76     {
 77         if(x%i==0)
 78         {
 79             tot++;su[tot]=1;s[tot]=i;
 80             while(x%i==0)su[tot]*=i,x/=i;
 81         }
 82     }
 83     if(x!=1)su[++tot]=x,s[tot]=x;
 84 }
 85 void yu(int x)
 86 {
 87     jie[x][0]=1;
 88     for(int i=1;i<s[x];i++)jie[x][i]=jie[x][i-1]*i%su[x];
 89 }
 90 signed main()
 91 {
 92     scanf("%lld%lld%lld",&p,&n,&m);
 93     for(int i=1;i<=m;i++)scanf("%lld",&w[i]),w[i]+=w[i-1];
 94     if(n<w[m]){puts("Impossible");return 0;}
 95     div(p);
 96     for(int i=1;i<=tot;i++)yu(i);
 97     for(int i=1;i<=tot;i++)
 98     {
 99         ans[i]=1;
100         for(int j=1;j<=m;j++)ans[i]=ans[i]*lus(n-w[j-1],w[j]-w[j-1],i)%su[i];
101     }
102     if(tot==1)
103     {
104         printf("%lld\n",ans[1]);
105         return 0;
106     }
107     int as=0;
108     for(int i=1;i<=tot;i++)
109     {
110         int tmp=p/su[i];
111         exgcd(tmp,su[i]);
112         xx=(xx%su[i]+su[i])%su[i];
113         xx=xx*tmp%p;
114         as+=xx*ans[i]%p;
115         as%=p;
116     }
117     printf("%lld\n",as);
118     return 0;
119 }

 

转载于:https://www.cnblogs.com/ezyzy/p/6590660.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值