组合狮子很好推，推出来是

$$(C_n^{sum}C_{sum}^{a_1}C_{sum-a_1}^{a_2}...C_{sum-a_1-a_2-..-a_{m-1}}^{a_m}) \mod k $$

现在只需要求组合数了，但是 $k$ 不是质数，而且 $n$ 巨大，暴力算阶乘会炸。

这时候我们再观察数据范围，发现 $k$ 能分解为一些小于 $10^5$ 的质因子。

那只要能求 $C_n^m \mod p^c$ (其中 $p$ 是质数) 就行了。

由于 $C_n^m = \frac{n!}{m!(n-m)!}$，其中分母部分可能含有因子 $p$。

此时 $\gcd(m!(n-m)!,p^c) \neq 1$，分母不存在逆元。

想想怎么把分母上的 $p$ 因子去掉，显然只要约分就行了。

我们设 $f_n$ 为 $n!$ 去掉所有 $p$ 因子的结果，$g_n$ 为 $n!$ 含有 $p$ 因子的个数。

则

$$C_n^m = \frac{n!}{m!(n-m)!} = \frac{f_n}{f_m f_{n-m}}p^{g_n-g_m-g_{n-m}} $$

这时候 $f_m$ 与 $p^c$ 互质，存在逆元，可以直接算。

现在我们只需要考虑怎么快速算 $f_n$ 了。

考虑把 $n!$ 展开，把所有形如 $kp$ 的拿到一边。

$$n! = (1 \times 2 \times 3 \times ... \times n) =(1 \times 2 \times ... \times (p-1) \times (p+1) \times ....)(p \times 2p \times 3p \times ...) $$$$= p^{\lfloor \frac{n}{p} \rfloor}\lfloor \frac{n}{p} \rfloor !(1 \times 2 \times 3 \times ... \times (p-1) \times (p+1) \times ...) $$

只考虑最右边那一坨取模 $p^c$ 的结果，则右边有许多个形如

$$((1+kp) \times (2+kp) \times ... \times (p-1+kp))$$

的东西的乘积。

发现当 $kp = p^c$ 时，这个狮子 $\mod p^c$ 的结果等于 $kp = 0$ 的结果，也就是这一坨狮子有循环节。

令这个循环节为 $s$，后面没循环完的一坨为 $t$，不然我真的写不下去了。

则:

$$n! = p^{\lfloor \frac{n}{p} \rfloor}\lfloor \frac{n}{p} \rfloor ! s^{\lfloor \frac{n}{p^c} \rfloor} t $$$$f_n = f_{\lfloor \frac{n}{p} \rfloor} s^{\lfloor \frac{n}{p^c} \rfloor} t $$

这个是个递归狮子，最多递归 $\log_p n$ 次，边界是 $f_0 = 1$。

然后 $g_n$ 也很好求，$g_n = g_{\lfloor \frac{n}{p} \rfloor} + {\lfloor \frac{n}{p} \rfloor}$。

最后 $CRT$ 合并一下就行了。

发现我们刚才发明了 $exLucas$ 算法.....

复杂度大概是基于模数的一只老哥，具体的自己想吧。

然后就理论可做了，实践见代码吧，我好累啊。

CODE

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define fi first
#define se second
int T,n,m;
int qpow(int a,int b,int mod){
	int r=1;
	while(b){
		if(b&1) r=r*a%mod;
		a=a*a%mod,b>>=1;
	}
	return r;
}
int calc_f(int x,int p,int pk){
	if(x==0) return 1;
	int res=1;
	for(int i=1;i<pk;i++){
		if(i%p) res=res*i%pk;
	}
	res=qpow(res,x/pk,pk);
	for(int i=1;i<=x%pk;i++){
		if(i%p) res=res*i%pk;
	}
	return calc_f(x/p,p,pk)*res%pk;
}
int calc_g(int x,int p){
	if(x<p) return 0;
	return calc_g(x/p,p)+x/p;
}
void exgcd(int a,int b,int &x,int &y){
	if(b==0) x=1,y=0;
	else{
		exgcd(b,a%b,x,y);
		int z=x;
		x=y,y=z-a/b*y;
	}
}
int getinv(int t,int p){
	int x,y;
	exgcd(t,p,x,y);
	return (x%p+p)%p;
}
int calc(int a,int b,int p,int pk){
	return calc_f(a,p,pk)*getinv(calc_f(b,p,pk),pk)%pk*getinv(calc_f(a-b,p,pk),pk)%pk*qpow(p,calc_g(a,p)-calc_g(b,p)-calc_g(a-b,p),pk)%pk;
}
int mod[64],res[64],Mod[64],tot;
int CRT(){
	int M=1;
	for(int i=1;i<=tot;i++){
		M*=mod[i];
	}
	int ans=0;
	for(int i=1;i<=tot;i++){
		Mod[i]=M/mod[i];
		ans=(ans+res[i]*Mod[i]%M*getinv(Mod[i],mod[i])%M)%M;
	}
	return ans;
}
int exlucas(int n,int m,int p){
	tot=0;
	for(int i=2;i<=101000;i++){
		if(p%i==0){
			int pk=1;
			while(p%i==0){
				p/=i,pk*=i;
			}
			mod[++tot]=pk,res[tot]=calc(n,m,i,pk);
		}
	}
	return CRT();
}
int a[7];
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0),cout.tie(0);
	int p;
	cin>>p>>n>>m;
	int sum=0;
	for(int i=1;i<=m;i++){
		cin>>a[i];
		sum+=a[i];
	}
	if(sum>n) cout<<"Impossible\n";
	else{
		int ans=exlucas(n,sum,p);
		for(int i=1;i<=m;i++){
			ans=ans*exlucas(sum,a[i],p)%p;
			sum-=a[i];
		}
		cout<<ans<<'\n';
	}
	return 0;
}

洛谷P2183 [国家集训队] 礼物