【spoj SEQN】【hdu 3439】Sequence
题意：
给出n、m、k 求C(n,k)*H(n-k)%m的值 H(n-k)为错排公式
题解：
先算H(n-k)
计算H(n)有个通式：
H(n)=(-1)^n+((-1)^(n-1))n+((-1)^(n-2))n(n-1)+...+n(n-1)(n-2)...3
证明详见维基百科：
因为我们是要算H(n-k)mod m的值 显然他的前不超过m项是>0的 而其它都为0 枚举求解即可
_________________________________________________________________________________
再说C(n,k)%m
我们知道世界上有个叫 Lucas 定理的东西：
C(n, m) % p = C(n / p, m / p) * C(n % p, m % p) % p
但是它要求p是质数 所以它和本题的正解没什么关系- -
_________________________________________________________________________________
由于除的数可能不与m互质 所以不能用乘法逆元
不妨将m分解 m=p1^a1*p2^a2*...*pn^an
分别计算C(n,k)%(pi^ai) 的值 最后用中国剩余定理计算答案
于是问题转换为求C(n,k)%(p^a)的值
C(n,k)=n!/((n-k)!*k!)
为了让乘法逆元能能使用 只要将n!、(n-k)!、k! 所有p的倍数 提出来单独计算即可
如果n比较小 枚举计算不为p的倍数的数的积即可 但是这里n<=10^9 我们就需要用一些比较高端的方法
设A(n)=(1*2*..*(p-1)*(p+1)*...*(2p-1)*(2p+1)*...*n)
n!=A(n)*(p*2p*...*(n/p)p)
=A(n)*(n/p)!*p^(n/p)
=A(n)*A(n/p)*(n/p/p)!*p^(n/p/p)*p^(n/p)
=...
就是对n!不断的分解为A(n)*(n/p)!*p^(n/p) 再递归处理(n/p)! 直到n/p为0
当然在计算的同时 要累加/减 p的指数sum
这样就能求出C(n,k)不包含p的乘积了 最后再乘p^sum 即为C(n,k)%(p^a)
_________________________________________________________________________________
问题结束了吗- - 显然木有 对于A(n)还需要用O(n)的时间求出 TLE！
设r=p^a
这里A(n)的值是要模r的 我们可以这样转换
A(n)=1*2*..*(p-1)*(p+1)*...*(2p-1)*(2p+1)*...*n
=(1*2*...*(r-1))*((r+1)*(r+2)*...*(2r-1))*(...*n)
=(1*2*...*(r-1))*(1*2*...*(r-1))*...*(1*2*...*n%r) (mod r)
我们就可以预处理save[i]存下1*2*...*i （p的倍数不乘，i<=r-1）
A(n)=save[r-1]^(n/r)*save[n%r] 快速幂就可用log(n)的时间求出A
加上预处理O(r) r<=m<=10^5
最后乘m的质因数的个数blabla- -具体多少我也不知道了 反正能过
_________________________________________________________________________________
代码：
1 #include <cstdio> 2 typedef long long ll; 3 const ll M=100001; 4 ll n,k,m,t,a[M],r[M],ans,p[M],pri[M],bo[M],zs[M],save[M],tot; 5 void makepri(){ 6 for (ll i=2;i<M;i++){ 7 if (!bo[i]) pri[++pri[0]]=i; 8 for (ll j=1;j<=pri[0] && i*pri[j]<M;j++){ 9 bo[i*pri[j]]=1; 10 if (!(i%pri[j])) break; 11 } 12 } 13 } 14 ll mi(ll a,ll b,ll mo){ 15 ll res=1%mo; 16 for (;b;b>>=1){ 17 if (b&1) res=res*a%mo; 18 a=a*a%mo; 19 } 20 return res; 21 } 22 ll extgcd(ll &x,ll &y,ll a,ll b){ 23 if (!b){ 24 x=1,y=0; 25 return a; 26 }else{ 27 ll res=extgcd(x,y,b,a%b); 28 ll t=x; 29 x=y; 30 y=t-a/b*y; 31 return res; 32 } 33 } 34 ll getny(ll a,ll b){ 35 ll x,y,gc=extgcd(x,y,a,b),mo=b/gc; 36 x=(x%mo+mo)%mo; 37 return x; 38 } 39 40 void makep(){ 41 p[0]=r[0]=0; 42 for (ll i=1,x=m;i<=pri[0] && x>1;i++) 43 if (!(x%pri[i])){ 44 p[++p[0]]=pri[i]; 45 r[++r[0]]=1; 46 while (!(x%pri[i])) x/=pri[i],r[r[0]]*=pri[i]; 47 } 48 } 49 void makesave(ll t){ 50 save[0]=1; 51 for (ll i=1,x=1;i<=r[t];i++){ 52 if (i%p[t]) x=x*i%r[t]; 53 save[i]=x; 54 } 55 } 56 ll jc(ll n,ll t,ll &z,ll bo){ 57 if (!n) return 1; 58 z+=bo*(n/p[t]); 59 return mi(save[r[t]],n/r[t],r[t])*save[n%r[t]]%r[t]*jc(n/p[t],t,z,bo)%r[t]; 60 } 61 ll getc(){ 62 makep(); 63 for (ll i=1;i<=p[0];i++){ 64 zs[i]=0; 65 makesave(i); 66 a[i]=jc(n,i,zs[i],1)*getny(jc(n-k,i,zs[i],-1)*jc(k,i,zs[i],-1)%r[i],r[i])%r[i]; 67 a[i]=a[i]*mi(p[i],zs[i],r[i])%r[i]; 68 } 69 ll res=0; 70 for (ll i=1;i<=p[0];i++) 71 res=(res+m/r[i]*getny(m/r[i],r[i])%m*a[i]%m)%m; 72 return res; 73 } 74 75 ll geth(ll n){ 76 if (n==1) return 0; 77 ll res=(n&1) ? -1 : 1; 78 for (ll i=n,x=n*res*-1%m;i>=3 && x;x=x*(--i)*-1%m) 79 res=((res+x)%m+m)%m; 80 return res; 81 } 82 83 int main(){ 84 scanf("%I64d",&t); 85 makepri(); 86 for (ll i=1;i<=t;i++){ 87 tot=i; 88 scanf("%I64d%I64d%I64d",&n,&k,&m); 89 ans=(getc()*geth(n-k)%m+m)%m; 90 printf("Case %I64d: %I64d\n",i,ans); 91 } 92 }
注：这个代码能过hdu的数据但是不能过spoj的... 可能有点小bug- -？（我说spoj → →）
