ZS and The Birthday Paradox
题目链接：http://codeforces.com/contest/711/problem/E
数学题（Legendre's formula）
这题是以生日悖论（如果有23个或23个以上的人，那么至少有两个人的生日相同的概率要大于50%）为背景的数学题。公式本身不难推，主要是对大数的处理。
首先，我们需要找到分子和分母的公因数，这里用到了Legendre's formula（举个例子：10!的因数中质数3的个数为[10/3+10/(3^2)]）。因为若2^n-x与2^n有公因数，那么x与2^n有相同的公因数，所以求(2^n)(2^n-1)(2^n-2)*...*(2^n-k+1)与(2^n)^k的公因数，就转化为了求(k-1)!与(2^n)^k的公因数。
之后就是分别求分子和分母：
对于分母来说，只需要用快速幂（也可以用费马小定理）就可以很容易的求出，再乘上公因数的逆元即可；
而对于分子来说，稍微有点麻烦：
1.如果k>=M，即(2^n)(2^n-1)(2^n-2)*...*(2^n-k+1)中至少有连续的M个整数，
那么(2^n)(2^n-1)(2^n-2)*...*(2^n-k+1)一定为M的倍数，所以它被M求模后余0；
2.如果k<M，因为M很小，所以枚举一下，就可以求出分子，再乘上公因数的逆元即可。
总的时间复杂度为O（M+lgk+lgn）
（感谢游少半夜教我证明Orz）
证明：(A/B)modM=(A*(BmodM)')modM，其中B'为B在M下的逆元
令B=b1*b2*b3*...*bn，则((BmodM)')modM
=((b1*b2*b3*...*bn mod M)')modM
=(b1modM*b2modM*...*bn mod M）'modM
=(b1*b2*...*bn)modM*(b1modM*b2modM*...*bn mod M）'modM*(b1*b2*...*bn)'modM
=(b1modM*b2modM*...*bn mod M)*(b1modM*b2modM*...*bn mod M）'modM*(b1*b2*...*bn)'modM
=1*(b1*b2*...*bn)'modM=(b1*b2*...*bn)'modM
=b1*b1'modM*b2*b2'modM...bn*bn'mod M*(b1*b2*...*bn)'modM
=(b1*b2*...*bn)modM*(b1'modM*b2'modM...bn'mod M)*(b1*b2*...*bn)'modM
=b1'modM*b2'modM...bn'mod M
代码如下：
1 #include<cstdio> 2 #include<iostream> 3 #define M (long long)(1e6+3) 4 using namespace std; 5 typedef long long LL; 6 LL n,k,cnt,molecular,numerator,x,y; 7 LL exGCD(LL a,LL b){ 8 if(b==0){ 9 x=1,y=0; 10 return a; 11 } 12 LL r=exGCD(b,a%b); 13 LL temp=x; 14 x=y; 15 y=temp-(a/b)*y; 16 return r; 17 } 18 LL mod(LL a,LL b){ 19 LL base=a,temp=1; 20 while(b){ 21 if(b&1)temp=(temp*base)%M; 22 base=(base*base)%M; 23 b>>=1; 24 } 25 return temp; 26 } 27 int main(void){ 28 cin>>n>>k; 29 if(n<=63&&k>(((LL)1)<<n)){//总人数大于总天数 30 cout<<"1"<<" "<<"1"<<endl; 31 return 0; 32 } 33 LL temp=2; 34 while(k-1>=temp){//根据Legendre's formula,求出(k-1)!的因数中质数2的个数 35 cnt+=((k-1)/temp); 36 temp<<=1; 37 } 38 if(cnt){//若有公因数，则求出公因数的逆元 39 LL gcd=mod(2,cnt); 40 exGCD(gcd,M); 41 x=(x+M)%M; 42 }else x=1;//若没有公因数，则令x=1，来消除对后面计算的影响 43 temp=mod(2,n);//计算2^n 44 numerator=(mod(temp,k-1)*x)%M;//计算(2^n)^(k-1) 45 if(k>=M)molecular=0;//若分子出现连续的M个整数，则分子一定为M的倍数 46 else{ 47 molecular=1; 48 for(LL i=1;i<k;++i){ 49 molecular=((temp-i+M)%M*molecular)%M;//计算分子 50 } 51 molecular=(molecular*x)%M;//除以公因数 52 } 53 molecular=(numerator-molecular+M)%M; 54 cout<<molecular<<" "<<numerator<<endl; 55 }