BZOJ 2142 礼物 （扩展Lucas）

2142: 礼物

Time Limit: 10 Sec Memory Limit: 259 MB
Description

一年一度的圣诞节快要来到了。每年的圣诞节小E都会收到许多礼物，当然他也会送出许多礼物。不同的人物在小E
心目中的重要性不同，在小E心中分量越重的人，收到的礼物会越多。小E从商店中购买了n件礼物，打算送给m个人
，其中送给第i个人礼物数量为wi。请你帮忙计算出送礼物的方案数（两个方案被认为是不同的，当且仅当存在某
个人在这两种方案中收到的礼物不同）。由于方案数可能会很大，你只需要输出模P后的结果。

Input

输入的第一行包含一个正整数P，表示模；
第二行包含两个整整数n和m，分别表示小E从商店购买的礼物数和接受礼物的人数；
以下m行每行仅包含一个正整数wi，表示小E要送给第i个人的礼物数量。

Output

若不存在可行方案，则输出“Impossible”，否则输出一个整数，表示模P后的方案数。

Sample Input
1
00
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。

思路：
公式很容易得到C(sum,w[1])∗C(sum−w[1],w[2])…
接下来就是求组合数的问题了，发现p并不保证是一个质数，所以就要用到扩展Lucas了。
摘自 < Frods >链接
扩展Lucas：
若p不是素数，我们将p分解质因数，将Cnm分别按照（1）中的方法求对p的质因数的模，然后用中国剩余定理合并。
例如：
当我们需要计算C(n,m)mod(p)，其中p=p^q1×p^q2×…×p^qk，我们可以求出：
C(n,m)≡ai(modp^qi)(1 < i < k)
然后对于方程组：
x≡ai(modp^qi)(1 < i < k)
我们可以求出满足条件的最小的x,记为x0。
那么我们有：
C(n,m)≡x0(modp)
但是，我们发现，p^qi并不是一个素数，它是某个素数的某次方。
下面我们介绍如何计算C(n,m)modp^t(t≥2,p为素数)。
计算C(n,m)modp^t：
我们知道，C(n,m)=n!/(m!(n−m)!)，若我们可以计算出n! mod p^t，我们就能计算出m! mod p^t以及(n−m)! mod p^t。我们不妨设x=n! mod p^t, y=m! mod p^t, z=(n−m)! mod p^t,那么答案就是x * reverse(y, p^t) * reverse(z, p^t) (reverse(a,b)表示计算a对b的乘法逆元)。那么下面问题就转化成如何计算n! mod p^t。例如p=3,t=2,n=19
我们发现
n!=1×2×3×4×5×6×7×8×…×19=(1×2×4×5×7×8×…×16×17×19)×(3×6×9×12×15×18)=(1×2×4×5×7×8×…×16×17×19)×36×(1×2×3×4×5×6)
后面的部分恰好是(n/p)!，于是递归即可。前半部分是以p^t为周期的(1×2×4×5×7×8)≡(10×11×13×14×16×17)(mod 9)。下面是孤立的19，可以知道孤立出来的长度不超过 pt,于是直接计算即可。对于最后剩下的36这些数我们只要计算出n!,m!,(n−m)!里含有多少个p(不妨设分别有x , y ,z)，那么x−y−z就是C(n,m)中p的个数，直接计算就行。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#define LL long long  
using namespace std;

const int N=7;

LL MOD, n, m, w[N];

inline LL read(){
    LL x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

inline LL power(LL a, LL b, LL mod) {
    LL ret = 1;
    while(b) {
        if(b & 1) ret = ret * a % mod;
        a = a * a % mod;
        b >>= 1;
    }
    return ret;
} 

inline LL exgcd(LL a, LL b, LL &x, LL &y) {
    if( !b ) {
        x = 1; y = 0;
        return a;
    }
    LL d = exgcd(b, a%b, y, x);
    y -= a / b * x;
    return d;
}

inline LL inv(LL a, LL b){  
    LL x, y;  
    exgcd(a, b, x, y);  
    return (x % b + b) % b;  
}  

LL mud(LL n, LL p, LL pr){//求n！在mod(p)意义下的ans 
    if(n == 0) return 1;
    LL res = 1;
    for(LL i=2; i<=pr; i++)//以p^t为周期的阶乘 
        if(i % p) res = res * i % pr;
    res = power(res, n/pr, pr);//有多少个周期 
    LL r = n % pr;//有多少个单的 
    for(LL i=2; i<=r; i++) 
        if(i % p) res = res * i % pr;
    return res * mud(n/p, p, pr) % pr;
}

LL C(LL n, LL m, LL p, LL pr){
    if(n < m) return 0;
    LL x = mud(n, p, pr), y = mud(m, p, pr), z = mud(n-m, p, pr);
    //x=n! mod p^t, y=m! mod p^t, z=(n-m)! mod p^t (pr == p^t)
    LL c = 0;
    for(int i=n; i; i/=p) c += i / p;
    for(int i=m; i; i/=p) c -= i / p;
    for(int i=n-m; i; i/=p) c -= i / p;
    //计算出n!,m!,(n-m)!里含有多少个p
    LL a = x * inv(y, pr) % pr * inv(z, pr) % pr * power(p, c, pr) % pr;
    return a * (MOD / pr) % MOD * inv(MOD / pr, pr) % MOD;//中国剩余定理 
}

LL Lucas(LL n, LL m){
    LL x = MOD, res = 0;
    for(LL i=2; i<=x; i++) 
        if(x % i == 0){
            LL pr = 1;
            while(x % i == 0) x /= i, pr *= i;//分解出一个质因数 
            res = (res + C(n, m, i, pr)) % MOD;
        }
    return res;
}

int main(){
    MOD=read(); n=read(); m=read();
    LL sum = 0;
    for(int i=1; i<=m; i++) w[i] = read(), sum += w[i];
    if(sum > n){ puts("Impossible"); return 0;}
    LL ans = 1;
    for(int i=1; i<=m; i++) 
        ans = ans * Lucas(n, w[i]) % MOD, n -= w[i];
    printf("%lld",ans);
    return 0;
}