扩展Lucas求组合数取模

发表于 | 更新于

做法来自jianglangcaijin

Problem

求 是 质 数

Solution

实在没看出扩展在哪。。和Lucas有任何相似之处吗。。。 看了很久才看懂啊。。。太弱了。 首先，如果可以分别求出来 ，用中国剩余定理就可以 质 因 数 个 数 地解出来了。 现在的问题是要求 是 质 数 。 而Cnm=n!m!(n−m)!C_n^m=\dfrac {n!}{m!(n-m)!} 而xx有关于yy的逆元，当且仅当gcd(x,y)==1\gcd(x,y)==1 所以要把除以分母变成乘它的逆元，必须先把分母中的含有的pp的幂先处理掉。 现在单独看n!n!的话

n!=∏p∤ii×p[pn][np]!n!=\prod_{p\nmid i}i\times p^{[\dfrac pn]}[\frac np]!

以n=10,p=3n=10,p=3为例 10!=1×2×..×9×1010!=1\times 2\times ..\times 9\times 10 把33的倍数提出来，就是 10!=1×2×4×5×7×8×10×33(1×2×3)10!=1\times 2\times 4\times 5\times7\times8\times10\times 3^3(1\times2\times3) 分别考虑这些项 [np]![\dfrac np]!递归处理 的值是循环的，每一段[i+1,i×pt][i+1,i\times p^t]区间的乘积 是相等的(不会证明啊。。求证明啊。。)，于是就可以循环部分用快速幂，超出循环的部分直接累乘。 而p[pn]p^{[\dfrac pn]}就是需要处理掉的部分，每一层递归的p[pn]p^{[\dfrac pn]}累乘起来就是在n!n!出现的所有pp了。 可以直接整体地把n!n!的质因子分解中pp的次数reprep求出来，把prepp^{rep}从阶乘中提出来，就可以消除其影响了。 于是就得到了上面问题的解决方法： 把n!,m!,(n−m)!n!,m!,(n-m)!中pp的次数repn,repm,repn−mrep_n,rep_m,rep_{n-m}分别求出来，而这个是比较好求的。幂的相除相当于指数相减，所以CnmC_n^m中包含的pp的次数为repn−repm−repn−mrep_n-rep_m-rep_{n-m}。消除了pp的影响，对剩下的两项递归求解即可。 呃虽然好求我也是看了一会儿才懂。。 

1
for(int i=n;i;i/=p) res+=i/p;

直观地感受一下，当i==ni==n时，i/pi/p求的是ii之内有多少个数字的是pp的倍数，i/=pi/=p相当于把ii以内的数集体/p/p，这样只含p1p^1的数字就被干掉了，剩下的数字继续进行这个过程。

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int MAXN=100;
typedef long long LL;
int P;
int Pow(int a,int v)
{
	int res=1;
	while(v)
	{
		if(v&1)
			res*=a;
		a*=a;
		v>>=1;
	}
	return res;
}
int Pow(int a,int v,int P)
{
	int res=1;
	while(v)
	{
		if(v&1)
			res=(LL)res*a%P;
		a=(LL)a*a%P;
		v>>=1;
	}
	return res;
}
int exgcd(int a,int b,LL &x,LL &y)
{
	if(!b)
	{
		x=1,y=0;
		return a;
	}
	else
	{
		int d=exgcd(b,a%b,x,y);
		LL nx=y,ny=x-(a/b)*y;
		x=nx,y=ny;
		return d;
	}
}
int inv(int a,int n)
{
	LL x,y;
	int d=exgcd(a,n,x,y);
	assert(d==1);
	return (x%n+n)%n;
}
namespace iLucas
{
	int P,dc;
	int dv[MAXN],dt[MAXN],vt[MAXN];
	//A[]
	//M[i]=M/m[i]
	//t[i]*M[i]=1 (mod m[i])
	int A[MAXN],M[MAXN],t[MAXN],m[MAXN],tM[MAXN];
	void Add(int np,int &x)
	{
		dv[++dc]=np,dt[dc]=0;
		while(x%np==0) dt[dc]++,x/=np;
		m[dc]=vt[dc]=Pow(dv[dc],dt[dc]);
	}
	void Init(int x)
	{
		P=x;
		for(int i=2;i*i<=x;i++) if(x%i==0) Add(i,x);
		if(x!=1) Add(x,x);
		int pM=1;//pi M<x
		for(int i=1;i<=dc;i++) pM*=m[i];
		for(int i=1;i<=dc;i++)
		{
			M[i]=pM/m[i];
			t[i]=inv(M[i],m[i]);
			tM[i]=(LL)M[i]*t[i]%P;
		}
	}
	int Calc(int n,int v,int t)
	{
		if(n==0) return 1;
		int vt=Pow(v,t),prod=1;
		for(int i=1;i<=vt;i++)
			if(i%v) prod=(LL)prod*i%vt;
		int rpt=n/vt,res=Pow(prod,rpt,vt);
		for(int i=rpt*vt+1;i<=n;i++)
			if(i%v) res=(LL)res*i%vt;
		return (LL)res*Calc(n/v,v,t)%vt;
	}
	int C(int n,int m,int v,int t)//C(n,m) mod v^t
	{
		int vt=Pow(v,t),vtot=0;
		for(int i=n;i;i/=v) vtot+=i/v;
		for(int i=m;i;i/=v) vtot-=i/v;
		for(int i=n-m;i;i/=v) vtot-=i/v;
		int prod=Pow(v,vtot,vt),
			cn=Calc(n,v,t),
			cm=Calc(m,v,t),
			cnmm=Calc(n-m,v,t);
		return (LL)prod*cn*inv(cm,vt)*inv(cnmm,vt)%vt;
	}
	int CRT()
	{
		int res=0;
		for(int i=1;i<=dc;i++) res=((LL)A[i]*tM[i]%P+res)%P;
		return res;
	}
	int C(int n,int m)
	{
		for(int i=1;i<=dc;i++) A[i]=C(n,m,dv[i],dt[i]);
		return CRT();
	}
}
int C(int n,int m)
{
	if(n<m) return 0;
	return iLucas::C(n,m);
}
bool ud[MAXN];
int a[MAXN],ac,n,m;
//前ac个数<=a[i] 所有的数>=1的方案数
int calc()
{
	int sum=m,cnt=0;
	for(int i=1;i<=ac;i++)
		if(ud[i]) cnt++,sum-=a[i];
	int res=C(sum-1,n-1);
	if(cnt&1)//奇数
		res=((P-res)%P+P)%P;
	return res;
}
int DFS(int pos)
{
	if(pos==ac+1)
		return calc();
	else
	{
		int res=0;
		ud[pos]=0;res=((LL)res+DFS(pos+1))%P;
		ud[pos]=1;res=((LL)res+DFS(pos+1))%P;//+=
		return res;
	}
}
void WORK()
{
	int n1,n2,t;
	get(n),get(n1),get(n2),get(m);
	for(int i=1;i<=n1;i++) get(a[i]);
	ac=n1;
	for(int i=1;i<=n2;i++) get(t),m-=t-1;
	printf("%d\n",DFS(1));
}
int main()
{
	int T;get(T),get(P);iLucas::Init(P);
	while(T--) WORK();
	return 0;
}

貌似跑得很慢。应该还有更优的做法。待进一步研究。 学习了一下别人的写法，学到了这些

小小的优化

  1. 对 一 个 合 适 的 值 的 值打表。
  2. 如果算出来的repn−repm−repn−mrep_n-rep_m-rep_{n-m}比当前需要算的tt（即模数ptp^t的指数）还大，直接返回0。

彻底改进

对于每个数质因子pip_i，预处理table[i][]table[i][]

table[i][j]=∏pi∤k,k<jktable[i][j]=\prod_{p_i\nmid k,k<j}k

上面的很多东西就O(1)O(1)查表即可。见Code。

Modified Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
//Code by Lucida
#include<bits/stdc++.h>
#define get(x) scanf("%d",&x)
//#define debug(x) std::cout<<#x<<'='<<x<<std::endl
template <class T> inline bool chkmx(T &a,const T &b){return a<b?a=b,1:0;}
template <class T> inline bool chkmn(T &a,const T &b){return a>b?a=b,1:0;}
const int MAXN=10;
typedef long long LL;
int P;
inline int Pow(int a,int v)
{
	int res=1;
	while(v)
	{
		if(v&1)
			res*=a;
		a*=a;
		v>>=1;
	}
	return res;
}
inline int Pow(int a,int v,int P)
{
	int res=1;
	while(v)
	{
		if(v&1)
			res=(LL)res*a%P;
		a=(LL)a*a%P;
		v>>=1;
	}
	return res;
}
int exgcd(int a,int b,LL &x,LL &y)
{
	if(!b)
	{
		x=1,y=0;
		return a;
	}
	else
	{
		int d=exgcd(b,a%b,x,y);
		LL nx=y,ny=x-(a/b)*y;
		x=nx,y=ny;
		return d;
	}
}
inline int inv(int a,int n)
{
	LL x,y;exgcd(a,n,x,y);
	return (x%n+n)%n;
}
namespace iLucas
{
	int P,dc;
	int dv[MAXN],dt[MAXN],dvt[MAXN];
	//A[]
	//M[i]=M/m[i]
	//t[i]*M[i]=1 (mod m[i])
	int A[MAXN],M[MAXN],t[MAXN],m[MAXN],tM[MAXN];
	const int FN=1e5,MAXF=FN+10;
	int table[MAXN][MAXF];
	inline void Add(int np,int &x)
	{
		++dc;dv[dc]=np;dt[dc]=0,dvt[dc]=1;
		while(x%np==0) dt[dc]++,x/=np,dvt[dc]*=np;
		int *ptr=table[dc];ptr[0]=1;
		for(int i=1;i<=dvt[dc];i++)
			if(i%np==0) ptr[i]=ptr[i-1];
			else ptr[i]=(LL)ptr[i-1]*i%dvt[dc];
		m[dc]=dvt[dc];M[dc]=P/m[dc];t[dc]=inv(M[dc],m[dc]);
		tM[dc]=(LL)t[dc]*M[dc]%P;
	}
	inline void Init(int x)
	{
		P=x;
		for(int i=2;i*i<=x;i++)
			if(x%i==0) Add(i,x);
		if(x!=1) Add(x,x);
	}
	int Calc(int n,int idx)
	{
		int *ptr=table[idx],v=dv[idx],vt=dvt[idx];
		if(n<v) return ptr[n];
		return (LL)Calc(n/v,idx)*Pow(ptr[vt],n/vt,vt)*ptr[n%vt]%vt;
	}
	int C(int n,int m,int idx)
	{
		int times,v,t,vt;
		v=dv[idx],t=dt[idx],vt=dvt[idx];times=0;
		for(int i=n;i;i/=v) times+=i/v;
		for(int i=m;i;i/=v) times-=i/v;
		for(int i=n-m;i;i/=v) times-=i/v;
		if(times>=t) return 0;
		else
		{
			int pt,fn,fm,fnmm;
			pt=Pow(v,times);
			fn=Calc(n,idx),fm=Calc(m,idx),fnmm=Calc(n-m,idx);
			return (LL)pt*fn*inv(fm,vt)*inv(fnmm,vt)%vt;
		}
	}
	inline int CRT()
	{
		int res=0;
		for(int i=1;i<=dc;i++) res=((LL)A[i]*tM[i]+res)%P;
		return res;
	}
	inline int C(int n,int m)
	{
		if(n<m) return 0;
		for(int i=1;i<=dc;i++) A[i]=C(n,m,i);
		return CRT();
	}
}
int a[MAXN],ac,n,m;
//前ac个数<=a[i] 所有的数>=1的方案数
void WORK()
{
	int n1,n2,t;
	get(n),get(n1),get(n2),get(m);
	for(int i=1;i<=n1;i++) get(a[i]);
	ac=n1;
	for(int i=1;i<=n2;i++) get(t),m-=t-1;
	int sum,cnt,ANS=0,res;
	for(int i=(1<<n1)-1;~i;i--)
	{
		sum=m,cnt=0;
		for(int j=1;j<=ac;j++)
			if(i&(1<<j>>1)) cnt++,sum-=a[j];
		res=iLucas::C(sum-1,n-1);
		if(cnt&1) ANS-=res,((ANS%=P)+=P)%=P;
		else ANS+=res,ANS%=P;
	}
	printf("%d\n",ANS);
}
int main()
{
	int T;get(T),get(P);iLucas::Init(P);
	while(T--) WORK();
	return 0;
}