题目描述

  \(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;
}