bzoj2142 禮物 擴充Lucas+CRT

http://www.elijahqi.win/archives/3181

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

考慮分析題目 他所要求的方案數其實是

Cw1n×Cw2n−w1×Cw3n−w1−w2 C n w 1 × C n − w 1 w 2 × C n − w 1 − w 2 w 3

是以把組合數打開 發現是很多階乘形式

n!w1!×w2!×...(n−w1−w2−...)! n ! w 1 ! × w 2 ! × . . . ( n − w 1 − w 2 − . . . ) !

那麼考慮模數任意的情況需要使用擴充Lucas

即首先将模數拆分成

P=pq11×pq22... P = p 1 q 1 × p 2 q 2 . . . 這樣的形式 針對每個質因數算一個模數 然後使用CRT将他們組合起來 因為模數互質是以使用普通CRT即可

普通CRT: (∑i=1totc[i]×Ppi×inv(Ppi,pi))%P ( ∑ i = 1 t o t c [ i ] × P p i × i n v ( P p i , p i ) ) % P

考慮每個組合數怎麼去算

把每個階乘和目前模數有關的項提出來 單獨處理 其他和他沒關的項分成三部分處理

首先就是 pi p i 的幂部分 就利用公式

ans=∑i=1,i=i∗pinni a n s = ∑ i = 1 , i = i ∗ p i n n i 大概思想就是比如三次的時候算一次的時候有他二次的時候有他 三次時候有他 分三次累計

然後可以發現剩餘部分變成 ⌊npi⌋! ⌊ n p i ⌋ ! 的結果 遞歸做即可

然後還有一部分 分為 pqii p i q i 部分和 %pqii % p i q i 模的部分暴力做即可

除的部分 算出等于 npqii n p i q i 這麼多次 然後對每個塊内暴力算 算完用快速幂倍增即可 然後這樣就算完了

最後要記得把前面幂次沒有被消去的部分快速幂求一下

具體例子和方法 參見擴充Lucas求法

https://blog.csdn.net/clove_unique/article/details/54571216

#include<cstdio>
#include<cctype>
#include<algorithm>
#define ll long long
using namespace std;
inline char gc(){
    static char now[<<],*S,*T;
    if (T==S){T=(S=now)+fread(now,,<<,stdin);if (T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=,f=;char ch=gc();
    while(!isdigit(ch)) {if (ch=='-') f=-;ch=gc();}
    while(isdigit(ch)) x=x*10+ch-'0',ch=gc();
    return x*f; 
}
const int N=;
int prime[N],mo[N],p,n,m,tot,w[],c[N];
inline void get_frac(int md){static int tp;tp=md;
    for (int i=;i*i<=tp;++i){
        if (md%i) continue;prime[++tot]=i;mo[tot]=;
        while(md%i==) mo[tot]*=i,md/=i; 
    }if (md>) prime[++tot]=md,mo[tot]=md;
}
inline int ksm(ll b,int t,int md){static ll tmp;
    for (tmp=;t;b=b*b%md,t>>=) if(t&) tmp=tmp*b%md;return tmp;
}
inline int jc(int n,int p,int md){
    if (!n) return ;ll tmp=;
    for (int i=;i<md;++i) if (i%p) tmp=tmp*i%md;
    tmp=ksm(tmp,n/md,md);int r=n%md;
    for (int i=;i<=r;++i) if (i%p) tmp=tmp*i%md;
    tmp=tmp*jc(n/p,p,md)%md;return tmp;
}
inline int exgcd(int a,int b,int &x,int &y){
    if(!b) {x=;y=;return a;}int g=exgcd(b,a%b,x,y);
    int t=x;x=y;y=t-a/b*y;return g;
}
inline int inv(int a,int b){
    static int t1,t2;exgcd(a,b,t1,t2);t1=(t1%b+b)%b;return t1;
}
inline int calc(int p,int pk){
    static ll t1,nm1,t2;
    t1=jc(n,p,pk),nm1=;t2=;
    for (ll i=p;i<=n;i*=p) nm1+=n/i;
    for (int i=;i<=m;++i) {
        t2=t2*jc(w[i],p,pk)%pk;
        for (ll j=p;j<=w[i];j*=p) nm1-=w[i]/j;
    }t2=inv(t2,pk);t1=t1*t2%pk*ksm(p,nm1,pk)%pk;return t1;
}
inline void inc(int &x,int v){x=x+v>=p?x+v-p:x+v;}
int main(){
    freopen("bzoj2142.in","r",stdin);
    p=read();n=read();m=read();ll sum=;
    for (int i=;i<=m;++i) sum+=(w[i]=read());
    if (sum>n) {puts("Impossible");return ;}
    get_frac(p);w[++m]=n-sum;int ans=;
    for (int i=;i<=tot;++i) c[i]=calc(prime[i],mo[i]);
    for (int i=;i<=tot;++i)
        inc(ans,(ll)c[i]*(p/mo[i])%p*inv(p/mo[i],mo[i])%p);
    printf("%d\n",ans);
    return ;
}