題目描述

  \(n\) 件禮物，送給 \(m\) 個人，其中送給第 \(i\) 個人禮物數量為 \(w_i\)。計算出送禮物的方案數（兩個方案被認為是不同的，當且僅當存在某個人在這兩種方案中收到的禮物不同），答案對 \(p\) 取模。

  資料範圍：\(1\leq n\leq 10^9,1\leq m\leq 5\)，\(p\) 不一定是質數。

exLucas定理

  \(\text{Lucas}\) 定理中對於模數 \(p\) 要求必須為質數，對於 \(p\) 不是質數的情況，需要用的 \(\text{exLucas}\) 定理。

原問題

  首先對於 \(p\) 進行質因數分解：\(p=p_1^{k_1}p_2^{k_2}\cdots p_{n}^{k_n}\)

，如果可以求出每個 \(\text{C}_{n}^{m}\equiv a_i\pmod {p_i^{k_i}}\)，對於同餘方程組：

\[\begin{cases}\text{C}_{n}^{m}\equiv a_1\pmod {p_1^{k_1}}\\\text{C}_{n}^{m}\equiv a_2\pmod {p_2^{k_2}}\\\cdots\\\text{C}_{n}^{m}\equiv a_n\pmod {p_n^{k_n}}\end{cases} \]

  由於 \(p_i^{k_i}\) 也不一定是質數，接下來介紹如何計算 \(\text{C}_{n}^{m}\mod p^t\)

。

組合數模質數冪

  首先由求組合數的公式 \(\text{C}_{n}^{m}=\frac{n!}{m!(n-m)!}\)，如果可以分別計算出 \(n!,m!,(n-m)!\) 在模 \(p^t\) 意義下的值，就可以得到答案。由於 \(m!\) 和 \((n-m)!\) 可能包含質因子 \(p\)，所以不能直接求它們對於 \(p^t\) 的逆元。此時我們可以將 \(n!,m!,(n-m)!\) 中的質因子 \(p\) 全部提出來，最後乘回去即可。

  即變為下式：

\[\text{C}_{n}^{m}=\frac{\frac{n!}{p^{k_1}}}{\frac{m!}{p^{k_2}}\times\frac{(n-m)!}{p^{k_3}}}\times p^{k_1-k_2-k_3} \]

  \(\frac{m!}{p^{k_2}},\frac{(n-m)!}{p^{k_3}}\) 和 \(p^t\) 是互質的，可以直接求逆元。

階乘除去質因子後模質數冪

  要計算 \(\frac{n!}{p^k}\mod p^t\)，先考慮如何計算 \(n!\mod p^t\)。

  當 \(p=3,t=2,n=19\) 時，有：

\[\begin{aligned}n!=&1·2·3\cdots19\\=&(1·2·4·5·7·8·19·11·13·14·16·17·19)·(3·6·8·12·15·18)\\=&(1·2·4·5·7·8·19·11·13·14·16·17·19)·3^6·(1·2·3·4·5·6)\end{aligned} \]

  後面的部分在模意義下相當於 \((\frac{n}{p})!\)，可以遞迴進行計算；\(p^{\lfloor\frac{n}{p}\rfloor}\) 部分直接快速冪。

  前面一部分是以 \(p^t\) 為週期的，也就是 \((1\times2\times4\times5\times7\times8)\equiv(10\times 11\times13\times14\times16\times17)\pmod {3^2}\)，每一個迴圈節長度為 $p^t $，所以只需要計算最後不滿足一個週期的數是哪些就可以了（這個例子中只需要計算 \(19\)）。顯然，不滿足一個週期的數的個數為 \(n-\lfloor\frac{n}{p^t}\rfloor\times p^t\)，不超過 \(p^t\) 個。

long long fac(long long n,long long p,long long pt)//處理階乘,n! mod p^t,pt=p^t
{
	if(n==0)
		return 1;
	long long ans=1;
	for(int i=1;i<pt;i++)//計算出迴圈節然後快速冪 
		if(i%p!=0)
			ans=ans*i%pt;
	ans=quick_pow(ans,n/pt,pt);
	for(int i=1;i<=n%pt;i++)//不滿足週期的數 
		if(i%p!=0)
			ans=ans*i%pt;
	return ans*fac(n/p,p,pt)%pt;//遞迴 
}

組合數模質數冪

\[\text{C}_{n}^{m}=\frac{\frac{n!}{p^{k_1}}}{\frac{m!}{p^{k_2}}\times\frac{(n-m)!}{p^{k_3}}}\times p^{k_1-k_2-k_3} \]

long long C(long long n,long long m,long long p,long long pt)//C(n,m) mod p^t 
{
    long long sum=0;
    for(long long i=n;i;i=i/p)//+k1
        sum=sum+i/p;
    for(long long i=m;i;i=i/p)//-k2
        sum=sum-i/p;
    for(long long i=n-m;i;i=i/p)//-k3
        sum=sum-i/p;
    return fac(n,p,pt)*quick_pow(p,sum,pt)%pt*inv(fac(m,p,pt),pt)%pt*inv(fac(n-m,p,pt),pt)%pt;
}

分析

  答案為：

\[\dbinom{n}{sum}·\dbinom{sum}{w_1}·\dbinom{sum-w_1}{w_2}\cdots\dbinom{sum-w_1-w_2\cdots-w_{m-1}}{w_m} \]

  用 \(\text{exLucas}\) 定理處理即可。

程式碼

#include<bits/stdc++.h>
using namespace std;
void exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    long long temp=x;x=y;y=temp-(a/b)*y;
}
long long inv(long long a,long long b)//ax=1(mod b),即ax+by=1,求x
{
    long long x,y;
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}
long long quick_pow(long long a,long long b,long long mod)
{
    long long ans=1;
    while(b)
    {
        if(b&1)
            ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans%mod;
}
long long fac(long long n,long long p,long long pt)//處理階乘,n! mod p^t,pt=p^t
{
    if(n==0)
        return 1;
    long long ans=1;
    for(int i=1;i<pt;i++)//計算出迴圈節以及快速冪
        if(i%p!=0)
            ans=ans*i%pt;
    ans=quick_pow(ans,n/pt,pt);
    for(int i=1;i<=n%pt;i++)//不滿足週期的數
        if(i%p!=0)
            ans=ans*i%pt;
    return ans*fac(n/p,p,pt)%pt;//遞迴
}
long long C(long long n,long long m,long long p,long long pt)//C(n,m) mod p^t
{
    long long sum=0;
    for(long long i=n;i;i=i/p)
        sum=sum+i/p;
    for(long long i=m;i;i=i/p)
        sum=sum-i/p;
    for(long long i=n-m;i;i=i/p)
        sum=sum-i/p;
    return fac(n,p,pt)*quick_pow(p,sum,pt)%pt*inv(fac(m,p,pt),pt)%pt*inv(fac(n-m,p,pt),pt)%pt;
}
long long a[1000010],b[1000010];
int cnt;
long long CRT()
{
    long long ans=0,lcm=1,x,y;
    for(int i=1;i<=cnt;i++)
        lcm=lcm*b[i];
    for(int i=1;i<=cnt;i++)
        ans=(ans+a[i]*(lcm/b[i])%lcm*inv(lcm/b[i],b[i])%lcm)%lcm;
    return (ans+lcm)%lcm;
}
long long exLucas(long long n,long long m,long long p)
{
    cnt=0;
    long long temp=p;
    for(long long i=2;i*i<=p;i++)//列舉p的質因子
    {
        long long k=1;
        while(temp%i==0)
            k=k*i,temp=temp/i;//p^k
        a[++cnt]=C(n,m,i,k);
        b[cnt]=k;
    }
    if(temp>1)
    {
        a[++cnt]=C(n,m,temp,temp);
        b[cnt]=temp;
    }
    return CRT();
}
long long w[10];
int main()
{
    long long n,m,p;
    cin>>p;
    cin>>n>>m;
    long long sum=0;
    for(int i=1;i<=m;i++)
    {
        scanf("%lld",&w[i]);
        sum=sum+w[i];
    }
    if(n<sum)
    {
        puts("Impossible");
        return 0;
    }
    long long ans=exLucas(n,sum,p);
    for(int i=1;i<=m;i++)
    {
        ans=ans*exLucas(sum,w[i],p)%p;
        sum=sum-w[i];
    }
    cout<<ans<<endl;
    return 0;
}