今天下午学习了二项式反演，做了一道错排的题，开始了苦逼的经历。
显然答案是C(︀n,k)︀*H(n − k).
其中H(i)为长度为i的错排序列
然后经过课件上一番二项式反演的推导
我就写了个扩展卢卡斯然后交上去了。
一直t啊.....
我算了算复杂度差不多是O(T*P*log^3P)
后来剪了剪枝，应该低了点。
还是t啊.....
我搜了搜题解发现没有我这么写的。
看了一下错排是有规律的，果然还是打表大法吼啊。
发个正解
#include<bits/stdc++.h> using namespace std; typedef long long ll; ll qmod(ll n,ll m,ll p) { ll ans=1; while(m) { if(m&1)ans=ans*n%p; n=n*n%p;m>>=1; } return ans; } void exgcd(ll a,ll b,ll &x,ll &y) { if(!b){ x=1;y=0;return; } exgcd(b,a%b,y,x);y-=a/b*x; } ll inv(ll n,ll p) { if(!n)return 0; ll a=n,b=p,x=0,y=0; exgcd(a,b,x,y); x=(x%b+b)%b; if(!x)x+=b; return x; } ll mul(ll n,ll pi,ll pk) { if(!n)return 1ll; ll ans=1; if(n/pk) { for(ll i=2;i<=pk;++i) if(i%pi)ans=ans*i%pk; ans=qmod(ans,n/pk,pk)%pk; } for(ll i=2;i<=n%pk;++i) if(i%pi)ans=ans*i%pk; return ans*mul(n/pi,pi,pk)%pk; } int C(ll n,ll m,ll p,ll pi,ll pk) { if(m>n)return 0; ll a=mul(n,pi,pk),b=mul(m,pi,pk),c=mul(n-m,pi,pk); ll k=0,ans; for(ll i=n;i;i/=pi)k+=i/pi; for(ll i=m;i;i/=pi)k-=i/pi; for(ll i=n-m;i;i/=pi)k-=i/pi; ans=a*inv(b,pk)%pk*inv(c,pk)%pk*qmod(pi,k,pk)%pk; return ans*(p/pk)%p*inv(p/pk,pk)%p;//CRT } bool v[100005]; int pp[100005],cnt; void pri() { for(int i=2;i<=100000;++i) { if(!v[i]) { pp[++cnt]=i; } for(int j=1;j<=cnt&&i*pp[j]<=100000;++j) { v[i*pp[j]]=1; if(i%pp[j]==0)break; } } } ll calc(ll n,ll m,ll p) { ll ans=0; for(ll x=p,i=1;i<=cnt&&x;++i) { if(x==1)break; if(x%pp[i]==0) { ll num=1; while(x%pp[i]==0)x/=pp[i],num*=pp[i]; ans=(ans+C(n,m,p,pp[i],num))%p; } } return ans; } ll F(ll x,ll p) { ll ans=0; if(x==0)return 1; x=x%(2*p); if(x==0)x=2*p; for(int i=2;i<=x;++i) ans=(ans*i+(i%2==0?1:-1))%p; return (ans+p)%p; } int main() { ll n,m,p,t,ans=0; scanf("%I64d",&t);pri();v[1]=1; for(int ii=1;ii<=t;++ii) { scanf("%I64d%I64d%I64d",&n,&m,&p); ans=calc(n,m,p)%p; ans=ans*F(n-m,p)%p; printf("Case %d: %I64d ",ii,ans); } return 0; }
再补个我的辣鸡程序，路过的dalao帮忙看看也中啊、
#include<bits/stdc++.h> using namespace std; typedef long long ll; ll qmod(ll n,ll m,ll p) { ll ans=1; while(m) { if(m&1)ans=ans*n%p; n=n*n%p;m>>=1; } return ans; } void exgcd(ll a,ll b,ll &x,ll &y) { if(!b){ x=1;y=0;return; } exgcd(b,a%b,y,x);y-=a/b*x; } ll inv(ll n,ll p) { if(!n)return 0; ll a=n,b=p,x=0,y=0; exgcd(a,b,x,y); x=(x%b+b)%b; if(!x)x+=b; return x; } ll mul(ll n,ll pi,ll pk) { if(!n)return 1ll; ll ans=1; if(n/pk) { for(ll i=2;i<=pk;++i) if(i%pi)ans=ans*i%pk; ans=qmod(ans,n/pk,pk)%pk; } for(ll i=2;i<=n%pk;++i) if(i%pi)ans=ans*i%pk; return ans*mul(n/pi,pi,pk)%pk; } int C(ll n,ll m,ll p,ll pi,ll pk) { if(m>n)return 0; ll a=mul(n,pi,pk),b=mul(m,pi,pk),c=mul(n-m,pi,pk); ll k=0,ans; for(ll i=n;i;i/=pi)k+=i/pi; for(ll i=m;i;i/=pi)k-=i/pi; for(ll i=n-m;i;i/=pi)k-=i/pi; ans=a*inv(b,pk)%pk*inv(c,pk)%pk*qmod(pi,k,pk)%pk; return ans*(p/pk)%p*inv(p/pk,pk)%p;//CRT } bool v[100005]; int pp[100005],cnt; void pri() { for(int i=2;i<=100000;++i) { if(!v[i]) { pp[++cnt]=i; } for(int j=1;j<=cnt&&i*pp[j]<=100000;++j) { v[i*pp[j]]=1; if(i%pp[j]==0)break; } } } ll calc(ll n,ll m,ll p) { ll ans=0; for(ll x=p,i=1;i<=cnt&&x;++i) { if(x==1)break; if(x%pp[i]==0) { ll num=1; while(x%pp[i]==0)x/=pp[i],num*=pp[i]; ans=(ans+C(n,m,p,pp[i],num))%p; } } return ans; } int main() { ll n,m,p,t,ans=0; scanf("%I64d",&t);pri();v[1]=1; for(int ii=1;ii<=t;++ii) { scanf("%I64d%I64d%I64d",&n,&m,&p); ans=calc(n,m,p)%p;n-=m;ll pre=0;ll num=1,pos=max(0ll,n-p); if(!ans) { printf("Case %d: %I64d ",ii,pre*ans%p);continue; } for(ll k=n;k>=pos;--k) { if(n!=k)num=num*(n-k)%p; if(k&1ll)pre=(pre-calc(n,k,p)%p*num%p+p)%p; else pre=(pre+calc(n,k,p)%p*num%p)%p; if(!num)break; } printf("Case %d: %I64d ",ii,pre*ans%p); } return 0; }
好吧，蒟蒻苦逼的一下午。
同时纪念衡水人民224起义。