数论三连2·LUCAS和EXLUCAS

LUCAS：C(n,m)=C(n%P,m%P)*C(n/P,m/P)；

m比较小或P比较小，而且P是质数的时候可以用。

EXLUCAS：和LUCAS其实关系不大啊。。我一开始被名字误导了，分解质因数以后在非常傻地LUCAS求组合数，

然后理所当然地WA了。。正确做法是分解质因数以后用神奇的方法求出阶乘，再用CRT求解，解就是要求的组合数了。

神奇的求阶乘方法：设分解出的第i项是pi^ki,n=k*(pi^ki)+t，则

n!=(n/pi)! ->去掉了所有含pi的项，剩下的项和pi一定是互质的，可以求逆元，而那些含pi的项除以pi^1以后又形成了这个阶乘，可以递归实现。

*((pi^ki)!)^k->显然，在模（pi^ki）域下，k*(pi^ki)+t=q*(pi^ki+t）=t。于是对于q=0->k-1，形成了循环。

*t!->这是多出来的q=k的那些产生的影响。

(*pw(pi^(n/pi)))->这一项在预处理阶乘时不需要计算在内，而要分开统计，在算组合数时加减pi的幂次即可统计出pi的影响。

于是就成功避开了没有逆元的情况。

bzoj2142是一道板子。然而我还是调了很多天。。然后最后的代码还无比丑陋。。果然是菜的啊orz。

#include<bits/stdc++.h>
using namespace std;
#define rep(x,y,z) for(int x=y; x<=z; x++)
#define downrep(x,y,z) for(int x=y; x>=z; x--)
#define ms(x,y,z) memset(x,y,sizeof(z))
#define LL long long
#define repedge(x,y) for(int x=hed[y]; ~x; x=edge[x].nex)
inline int read(){
	int x=0; int w=0; char ch=0;
	while (ch<'0' || ch>'9') w|=ch=='-',ch=getchar();
	while (ch>='0' && ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	return w? -x:x;
}
const int PR=65;
const int N=100005;
int tot;
LL w[15],a[3][PR][N],pr[PR],P;
struct NODE{ LL a,b; }o[PR];
struct node{ LL res1,res2; };
LL pw(LL a,LL b,LL c){
	LL res=1; a%=c; while (b){ if (b&1) res=(res*a)%c; a=(a*a)%c; b>>=1; }
	return (res%c+c)%c;
}
void exgcd(LL a,LL b,LL c,LL &x,LL &y){
	if (!a){ x=0; y=c/b; return; }
	if (!b){ y=0; x=c/a; return; }
	exgcd(b,a%b,c,y,x);
	y-=(a/b)*x;
}
LL gcd(LL a,LL b){ return b? gcd(b,a%b):a; }
LL lcm(LL a,LL b){ return (a/gcd(a,b))*b; }
LL inv(LL u,LL c){
	LL x,y; u%=c; exgcd(u,c,1,x,y); x=(x%c+c)%c; return x;
}
LL excrt(NODE *o,int n){
	//rep(i,1,n) cout<<o[i].a<<' '<<o[i].b<<endl;
	LL m=o[1].b; LL x=((o[1].a%m)+m)%m;
	rep(i,2,n){
		//x'=x+k*m=a mod b; m*k=(a-x) mod b;
		LL newm=lcm(m,o[i].b); LL a=m; LL b=o[i].b; LL c=((o[i].a-x)%o[i].b+o[i].b)%o[i].b;  
		LL d=gcd(a,b); if (c%d) return -1; a/=d; b/=d; c/=d; LL X,Y; exgcd(a,b,1,X,Y); 
		X=((X%o[i].b)*c)%o[i].b; x=((x+X*m)%newm+newm)%newm; m=newm;
	}   
	return (x%m+m)%m;
}
void pre(){
	LL d=floor(sqrt(P)); tot=0; LL x=P;
	rep(i,2,d)if (!(x%i)){
		LL s=1; while (!(x%i)){ s*=i; x/=i;	}
		pr[++tot]=i; o[tot].b=s;
	}
	if (x>1){ pr[++tot]=x; o[tot].b=x; }
	rep(i,1,tot){ a[1][i][0]=1; rep(j,1,o[i].b){ a[1][i][j]=a[1][i][j-1]; 
	if (j%pr[i]) a[1][i][j]=(a[1][i][j]*j)%o[i].b; } }
	rep(i,1,tot){ a[2][i][0]=1; rep(j,1,o[i].b) a[2][i][j]=inv(a[1][i][j],o[i].b); }
} 
node jc(LL n,int t,int kd){
	node res=(node){1,0}; if (n<=1) return res;
	node x=jc(n/pr[t],t,kd); res.res2=x.res2+(n/pr[t]);
    res.res1=(pw(a[kd][t][o[t].b],n/o[t].b,o[t].b)*a[kd][t][n%o[t].b])%o[t].b;
    res.res1=(res.res1*x.res1)%o[t].b; return res;
}
LL C(LL n,LL m){
	rep(i,1,tot){
		node t1=jc(n,i,1); node t2=jc(m,i,2); node t3=jc(n-m,i,2);
		o[i].a=(((t1.res1*t2.res1)%o[i].b*t3.res1)%o[i].b*pw(pr[i],t1.res2-t2.res2-t3.res2,o[i].b))%o[i].b;
		o[i].a=(o[i].a%o[i].b+o[i].b)%o[i].b;
	}   
	return excrt(o,tot);
}
int main(){
	P=read(); LL n=read(); int m=read(); pre(); 
	rep(i,1,m) w[i]=read(); LL ans=1;
	rep(i,1,m){ if (n<w[i]){ puts("Impossible"); return 0; }
	ans=(ans*C(n,w[i]))%P; n-=w[i]; } ans=(ans%P+P)%P;
	printf("%lld\n",ans);
	return 0; 
}