注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

一个蒟蒻的代码回收站

最后一次省选求rp

 
 
 

日志

 
 

数论专题之二:组合数取模  

2014-11-04 16:58:00|  分类: 算法 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

本文主要口胡组合数的取模(请要来D我的神犇绕道)

---------基本转自AC的空间

我们都知道


表示从n个里面选m个的方案数

当n,m比较大时,这东西可能会很大(比如C(10^9,5*10^8)=?),当然正常情况下我们都是会取模的,不然输出就爆炸了(“占了我整整十个硬盘~”)

但是取过模之后这个问题也是非常的蛋疼= =


一:n,m<=1000(我会暴力!)

大约就是n^2把杨辉三角求出来,然后就没有了

(组合数和杨辉三角的关系)


二:n,m<=10^6,p为质数(我会逆元!)

由于p为质数,所以有逆元。

逆元的话可以用扩展欧几里得算,也可以用费马小定理得到a*a^(p-2)=1(mod p),然后上快速幂

直接按公式推就好了

复杂度O(n)

三:n,m<=10^10,p为质数,且<=10^5(我会lucas定理!)

根据lucas定理,

把m,n转换成p进制,其中第i位为mi,ni

数论专题之二:组合数取模 - mxh1999 - 一个蒟蒻的代码回收站

然后转二

如果mi>ni则为0

四:n,m<=10^10,p=p1*p2*……*pk(p1,p2,……,pk为不同的质数)(我会中国剩余定理!)

把p质因数分解,得到p1,p2,p3,……,pk

然后令ansi=C(n,m)%pi(ansi怎么算?转三)

然后我们令mi=p/pi,g_mi=mi的逆元

则ans=∑ansi*mi*g_mi(mod p)

五:n,m<=10^10,p=p1^c1*p2^c2*……*pk^ck(pk^ck<=10^5)(我不会数论!)

lucas只对c=1的情况适用= =

孙子定理只有分解出来的东西互质时才可以= =

所以我们只能搞到C(n,m)%(pi^ci),然后呢?!

滚去百度

发现了扩展lucas定理= =

我们去算一下n!(mod p^c)

比如n=19,n!=1*2*3*4*5*6*7*8*9*……*19=(1*2*4*5*7*8)*(10*11*13*14*16*17)*3^6*(1*2*3*4*5*6)

我们发现这东西是按p^c循环的

然后暴力算出循环节,快速幂,再用中国剩余定理合并。

(好像口胡的太口胡了?不管了QAQ)

上代码(模拟赛的= =,求C(n,m+n)%p)

#include<cstdio>
#include<cmath>
#include<iostream>
using namespace std;

#define maxn 50
typedef long long ll;
int prime[maxn],c[maxn],cnt,sum[maxn];
int n,m,p,mo;
ll ans;
int quick_power(int a,int b,int mo)
{
if (b==0) return 1;
if (b==1) return a%mo;
ll ans=quick_power(a,b/2,mo);
ans*=ans;
ans%=mo;
if (b&1) ans*=a;
ans%=mo;
return ans;
}
int exGcd(int a,int b,int &x,int &y)
{
if (b==0){x=1;y=0;return a;}
int r=exGcd(b,a%b,x,y);
int t=x;x=y;y=t-a/b*y;
return r;
}
void prime_divide(int mo)
{
for (int i=2;i<=sqrt(mo);i++)
if (mo%i==0)
{
cnt++;prime[cnt]=i;sum[cnt]=1;c[cnt]=0;
while (mo%i==0) {sum[cnt]*=i;c[cnt]++;mo/=i;}
}
if (mo>1) {cnt++;prime[cnt]=mo;c[cnt]=1;sum[cnt]=mo;}
}
inline int cheng(int a,int b,int mo){return (int)(1LL*a*b%mo);}
int calc_fact(int q,int pos,int &num)
{
int tmp=q/sum[pos];
ll v=1;
for (int j=1;j<=sum[pos];j++)
if (j%prime[pos]) v=cheng(v,j,sum[pos]);
v=quick_power(v,tmp,sum[pos]);
for (int j=1;j<=q-tmp*sum[pos];j++)
if (j%prime[pos]) v=cheng(v,j,sum[pos]);
if (q/prime[pos]) v=cheng(v,calc_fact(q/prime[pos],pos,num),sum[pos]);
num+=q/prime[pos];
return v;
}
int calc(int x,int y,int pos)
{
int num1=0,num2=0,num3=0;
ll ans=calc_fact(y,pos,num1);
int tmp=calc_fact(x,pos,num2);
int tmp2=calc_fact(y-x,pos,num3);
int num=num1-num2-num3;
int p,q;
exGcd(tmp,sum[pos],p,q);
while (p<0) p+=sum[pos];
ans=cheng(ans,p,sum[pos]);
exGcd(tmp2,sum[pos],p,q);
while (p<0) p+=sum[pos];
ans=cheng(ans,p,sum[pos]);
ans=cheng(ans,quick_power(prime[pos],num,sum[pos]),sum[pos]);
return ans;
}
int main()
{
freopen("reward.in","r",stdin);
freopen("reward.out","w",stdout);
scanf("%d%d%d",&n,&m,&mo);
m+=n;
prime_divide(mo);
ans=0;
for (int i=1;i<=cnt;i++)
{
int tmp=calc(n,m,i);
int mi=mo/sum[i];
int tmp2=mi%sum[i];
int x,y;
exGcd(tmp2,sum[i],x,y);
while (x<0) x+=sum[i];
ans+=(ll)tmp*mi*x;
ans%=mo;
}
printf("%lld\n",ans);
fclose(stdin);
fclose(stdout);
return 0;
}

六,n,m,p<=10^10

我不会!!!!

求大神教导!

如果这个做出了……应该组合数取模就可以终结了吧

  评论这张
 
阅读(71)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018