2142: 礼物
Time Limit: 10 Sec Memory Limit: 259 MBSubmit: 2286 Solved: 1009[][][]Description
一年一度的圣诞节快要来到了。每年的圣诞节小E都会收到许多礼物，当然他也会送出许多礼物。不同的人物在小E心目中的重要性不同，在小E心中分量越重的人，收到的礼物会越多。小E从商店中购买了n件礼物，打算送给m个人，其中送给第i个人礼物数量为wi。请你帮忙计算出送礼物的方案数（两个方案被认为是不同的，当且仅当存在某个人在这两种方案中收到的礼物不同）。由于方案数可能会很大，你只需要输出模P后的结果。Input
输入的第一行包含一个正整数P，表示模；第二行包含两个整整数n和m，分别表示小E从商店购买的礼物数和接受礼物的人数；以下m行每行仅包含一个正整数wi，表示小E要送给第i个人的礼物数量。Output
若不存在可行方案，则输出“Impossible”，否则输出一个整数，表示模P后的方案数。
Sample Input
100 4 2 1 2Sample Output
12 【样例说明】 下面是对样例1的说明。 以“/”分割，“/”前后分别表示送给第一个人和第二个人的礼物编号。12种方案详情如下： 1/23 1/24 1/34 2/13 2/14 2/34 3/12 3/14 3/24 4/12 4/13 4/23 【数据规模和约定】 设P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct，pi为质数。 对于100%的数据，1≤n≤109，1≤m≤5，1≤pi^ci≤10^5。HINT
Source
[ ][ ][ ]
数论大集合。
答案显然$C_n^{n - w[1]}C_{n - w[1]}^{w[2]}C_{n - w[1] - w[2]}^{w[3]}......$，模数不互质成为难点。
扩展Lucas：将模数质因数分解，再用CRT合并。
问题就只剩下求$n!\%p_i^{k_i}$了。
首先由于不互质无法找到逆元，把$n!$中的$p_i$全部取出来：$ans=\lfloor\frac{n}{p}\rfloor+\lfloor\frac{n}{p^2} \rfloor+\lfloor\frac{n}{p^3}\rfloor+...$
然后考虑剩下部分怎么做，用$F(n,p_i,p_i^{k_i})$表示答案，$f(n,p_i,p_i^{k_i})$表示$\prod_{j=1,j\perp p_i}^{n}j(mod\ p_i^{k_i})$，则：
$F(n,p_i,p_i^{k_i})=\lfloor\frac{n}{p_i}\rfloor!\times p_i^{\lfloor\frac{n}{p_i}\rfloor}\times f(p_i^{k_i},p_i,p_i^{k_i})^{\lfloor\frac{n}{p_i^{k_i}}\rfloor}\times f(n\%p_i^{k_i})\%p_i^{k_i}$
递归处理即可，要用到exgcd。
复杂度好像是，$O(\sqrt{P}+m\log_{2}^{2}p_i^{k_i}\log{p_i})$不过肯定跑不满。
1 #include2 #include 3 #define rep(i,l,r) for (int i=l; i<=r; i++) 4 typedef long long ll; 5 using namespace std; 6 7 const int N=10010; 8 ll mod,n,m,w[10],ans,x,y,Mod[N],st[N],r[N],num; 9 10 ll ksm(ll a,ll b,ll p){11 ll res;12 for (res=1; b; a=a*a%p,b>>=1)13 if (b&1) res=res*a%p;14 return res;15 }16 17 ll exgcd(ll a,ll b,ll &x,ll &y){18 ll d=a;19 if (b) d=exgcd(b,a%b,y,x),y-=a/b*x; else x=1,y=0;20 return d;21 }22 23 ll inv(ll t,ll p){ ll x,y; exgcd(t,p,x,y); return (x+p)%p; }24 25 ll F(ll n,ll pi,ll pk){26 if (!n) return 1;27 ll ans=1;28 rep(i,2,pk) if (i%pi) ans=ans*i%pk;29 ans=ksm(ans,n/pk,pk);30 rep(i,2,n%pk) if (i%pi) ans=ans*i%pk;31 return ans*F(n/pi,pi,pk)%pk;32 }33 34 ll exlucas(ll n,ll m,ll pi,ll pk){35 if (m>n) return 0;36 ll a=F(n,pi,pk),b=F(m,pi,pk),c=F(n-m,pi,pk);37 ll k=0;38 for (ll i=n; i; i/=pi) k+=i/pi;39 for (ll i=m; i; i/=pi) k-=i/pi;40 for (ll i=n-m; i; i/=pi) k-=i/pi;41 return a*inv(b,pk)%pk*inv(c,pk)%pk*ksm(pi,k,pk)%pk;42 }43 44 ll CRT(ll n,ll r[],ll m[]){45 ll M=1,res=0,w;46 rep(i,1,n) M*=m[i];47 rep(i,1,n) w=M/m[i],res=(res+w*inv(w,m[i])*r[i])%M;48 return (res+M)%M;49 }50 51 ll par(ll n,ll m[],ll st[]){52 ll num=0;53 for (ll i=2; i*i<=n; i++) if (n%i==0){54 ll pk=1;55 while (n%i==0) pk*=i,n/=i;56 m[++num]=pk; st[num]=i;57 }58 if (n>1) m[++num]=n,st[num]=n;59 return num;60 }61 62 ll excomb(ll n,ll m){63 rep(i,1,num) r[i]=exlucas(n,m,st[i],Mod[i]);64 return CRT(num,r,Mod);65 }66 67 int main(){68 freopen("bzoj2142.in","r",stdin);69 freopen("bzoj2142.out","w",stdout);70 scanf("%lld\n",&mod); scanf("%lld%lld",&n,&m);71 ll sum=0; rep(i,1,m) scanf("%lld",&w[i]),sum+=w[i];72 if (n