BZOJ 2142 礼物 组合数学 CRT 中国剩余定理

2142: 礼物

Time Limit: 10 Sec  Memory Limit: 259 MB
Submit: 1450  Solved: 593
[Submit][Status][Discuss]

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 2

Sample 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

Solution

根据题目大意可以列出式子C(n,sum)*C(sum, w[1])*C(sum-w[1], w[2])··········

如果p是质数的话，显然就是Lucas的裸题。

但现在p不是质数了，要怎么办呢？把p分解成一个个pi^ci，这样的话就保证了两两互质，CRT求解。

那怎么算C(x,y)%(pi^ci)呢？把C(x,y)转化成x!、y!和(x-y)!。

由于需要求逆元，而n!不一定与pi^ci互质，我们就需要把n!分解为k*(pi^x)，其中，k与pi互质。

对于当前的n!，我们发现可以把它整理为一些模(pi^ci)下的循环串乘上(pi^(n/pi))再乘上(n/pi)!

举个栗子，17!%(3^2) = (1*2*4*5*7*8*10*11*13*14*16*17)*(3^5)*(1*2*3*4*5)

很显然，栗子最后的那部分就是5!%3这个显然可以递归下去。

对于最前的那部分，我们可以发现，它是由一些循环节（模3^2下）组成的，且均不含有质因数3。

仔细分析，就是[1,3^2]中与质因数3互质的数。求出循环节之后就可以直接用快速幂完成。

对于中间的部分，我们收集起来，算组合数的时候统一再做一次快速幂。

本题主要的问题是处理好逆元的问题，把n!分解成k*(pi^x)，含pi的部分单独算，k与pi^ci互质，可以直接求解逆元。

Code

  1 #include <cstdio>
  2 #include <cstdlib>
  3 #include <cstring>
  4 #include <string>
  5 #include <algorithm>
  6 
  7 using namespace std;
  8 
  9 #define REP(i, a, b) for (int i = (a), i##_end_ = (b); i <= i##_end_; ++i)
 10 #define fi first
 11 #define se second
 12 #define mp make_pair
 13 typedef long long LL;
 14 typedef pair<LL, LL> pa;
 15 const int maxn = 1e5;
 16 LL p, a[105];
 17 LL prime[maxn+10], p_cnt;
 18 bool isNotPrime[maxn+10];
 19 LL p_num[105], p_mod[105], p_sum[105], cnt;
 20 
 21 void prepare()
 22 {
 23     REP(i, 2, maxn)
 24     {
 25         if (!isNotPrime[i]) prime[++p_cnt] = i;
 26         REP(j, 1, p_cnt)
 27         {
 28             if (i*prime[j] > maxn) break ;
 29             isNotPrime[i*prime[j]] = 1;
 30             if (i%prime[j] == 0) break ; 
 31         }
 32     }
 33     LL temp = p;
 34     REP(i, 1, p_cnt)
 35     {
 36         if (prime[i] > temp || temp == 1) break ;
 37         if (temp%prime[i] == 0)
 38         {
 39             p_num[++cnt] = prime[i], p_mod[cnt] = 1, p_sum[cnt] = 0;
 40             while (temp%prime[i] == 0)
 41                 p_mod[cnt] *= prime[i], p_sum[cnt] ++, temp /= prime[i];
 42         }
 43     }
 44 }
 45 
 46 LL power(LL x, LL y, LL MOD)
 47 {
 48     LL ret = 1;
 49     while (y > 0)
 50     {
 51         if (y&1) ret = (ret*x)%MOD;
 52         x = (x*x)%MOD;
 53         y >>= 1;
 54     }
 55     return ret;
 56 }
 57 
 58 pa calc(int k, LL num)
 59 {
 60     if (num == 0) return mp(0, 1);
 61     LL x = num/p_num[k], y = num/p_mod[k];
 62     LL ans = 1;
 63     if (y)
 64     {
 65         LL t_num = 1;
 66         REP(i, 1, p_mod[k]-1)
 67             if (i%p_num[k] != 0) t_num = (t_num*i)%p_mod[k];
 68         ans = power(t_num, y, p_mod[k]);
 69     }
 70     if (num%p_mod[k] != 0)
 71     {
 72         REP(i, 1, num%p_mod[k])
 73             if (i%p_num[k] != 0) ans = (ans*i)%p_mod[k];
 74     }
 75     pa temp = calc(k, x);
 76     return mp(temp.fi+x, temp.se*ans%p_mod[k]);
 77 }
 78 
 79 LL ex_gcd(LL aa, LL bb, LL &x, LL &y)
 80 {
 81     if (bb == 0)
 82     {
 83         x = 1, y = 0;
 84         return aa;
 85     }
 86     LL ret = ex_gcd(bb, aa%bb, x, y);
 87     LL temp = x;
 88     x = y, y = temp-(aa/bb)*y;
 89     return ret;
 90 }
 91 
 92 LL inv(LL k, LL MOD)
 93 {
 94     LL x, y;
 95     ex_gcd(k, MOD, x, y);
 96     x = (x%MOD+MOD)%MOD;
 97     return x;
 98 }
 99 
100 LL CRT()
101 {
102     LL ret = 0;
103     REP(i, 1, cnt)
104     {
105         LL Mi = p/p_mod[i];
106         ret = (ret+((Mi*inv(Mi, p_mod[i]))%p*a[i])%p)%p;
107     }
108     ret = (ret%p+p)%p;
109     return ret;
110 }
111 
112 LL work(LL x, LL y)
113 {
114     REP(i, 1, cnt)
115     {
116         pa t1 = calc(i, x), t2 = calc(i, y), t3 = calc(i, x-y);
117         a[i] = power(p_num[i], t1.fi-t2.fi-t3.fi, p_mod[i]);
118         a[i] = (((a[i]*t1.se)%p_mod[i]*inv(t2.se, p_mod[i]))%p_mod[i]*inv(t3.se, p_mod[i]))%p_mod[i];
119     }
120     return CRT();
121 }
122 
123 int main()
124 {
125 //    freopen("a.in", "r", stdin);
126 //    freopen("a.out", "w", stdout);
127     LL n, m, w[15], sum = 0;
128     scanf("%lld %lld %lld", &p, &n, &m);
129     prepare();
130     REP(i, 1, m) scanf("%lld", &w[i]), sum += w[i];
131     if (sum > n) { puts("Impossible"); return 0; }
132     LL ans = work(n, sum);
133     REP(i, 1, m)
134         ans = (ans*work(sum, w[i]))%p, sum -= w[i];//, printf("%lld\n", ans);
135     ans = (ans+p)%p;
136     printf("%lld\n", ans);
137     return 0;
138 }
View Code

 

　

posted @ 2017-03-28 16:10  Splay  阅读(367)  评论(0)  编辑  收藏  举报