组合数,卢卡斯定理

有$q$个询问，每个询问两个整数 $n$, $m$ (1$<=m<=n<=2000$) 问 $C^{m}_{n} mod (10^{9} + 7)$

若不选则就是从剩余 n-1 选 m 个若选就是从剩下n-1中选m-1个由加法原理得。 $C^{m}_{n} = C^{m}_{n-1} + C^{m-1}_{n-1}$

通过打表不难发现就是杨辉三角数

时间复杂度$O(n)$

1
2
3
4
5
6
7
8


const int mod = 1e9 + 7;
const int N = 2010;
void get_c(){
for(int i = 0;i<N;++i)
    for(int j = 0;j<=i;++j)
        if(j == 0) c[i][j] = 1;//如果不选 也是一种方案数
        else c[i][j] =(c[i-1][j] + c[i-1][j-1]) % mod;
}

当数量级来到 10^5 , 上述递推法会TLE

考虑用高中组合公式直接计算

$C^{m}_{n} = \frac{n!}{(n-m)!m!}$

lucas

若给定$m,n,p$的值求 $C^{m}_{n} mod p 的值$ 其中 $1<=m<=n<=10^{8} , 1<=p<=10^{5}$ 保证$p$为质数。

可以发现由于给定的数值非常大,之前的方法会失效因为在求逆元的时候有可能 $m$ 为 $p$的倍数而不是互质的情况.

这时候就要用到 lucas定理

$C^{m}_{n} \equiv C^{\frac{m}{p}}_{\frac{n}{p}}C^{m\ mod p}_{n\ mod\ p}(mod \ p)$

$n\mod p , m \mod p$一定是小于 $p$的数可以直接求解

$C^{\frac{m}{p}}_{\frac{n}{p}}$可以继续用Lucas递归求解，而 m % p，n % p 一定是小于 $p$的数，既互质

当$m = 0$ 时, 返回1 .不选的方案数为1

证明：引理1： $C^{x}_p \equiv 0 (mod\ p)$ ($0<x<p$) 因$C^{x}_{p}=\frac{p!}{x!(p-x)!} = \frac{p(p-1)!}{x(x-1)!(p-x)!} = \frac{p}{x}C^{x-1}_{p-1}$

因为gcd(m,p) = gcd(n,p) = 1;

∴ $C^{x}_{p} \equiv p*x^{-1}C^{x-1}_{p-1} \equiv 0 (mod\ p)$ 其中 $x^{-1}为x的逆元$

引理2

$(1+x)^{p} \equiv 1+x^{p}\ (mod\ p)$

显然，由二项式定理展开$(1+x)^{p} = {\textstyle \sum_{p}^{i=0}} C^{i}_{p}X^{i} $ 又由引理1 显然成立

令$n = ap+b , m = cp+d$

$(1+x)^{n} \equiv {\textstyle \sum_{i=0}^{n}} C^{i}_{p}X^{i} (mod\ p)$

$$ (1+x)^{n} \equiv (1+x)^{ap+b} \quad\quad① \ \equiv ((1+x)^{p})^{a}\times (1+x)^{b}\ $$

通过引理2 展开 $$ \equiv (1+x^{p})^{a}\times (1+x)^{b}\ $$

$$ \equiv {\textstyle \sum_{i=0}^{a}} C^{i}_{a}X^{ip} {\textstyle \sum_{j=0}^{b}} C^{j}_{b}X^{j} (mod\ p)\quad \quad ②\ $$

观察一式带入系数$m$ 则有$x^{m}$的系数为 $C^{m}_{n}$

观察二式带入系数$m$ ,又∵$m = cp+d$ 则存在 $x^{cp}\times x^{d}$ ，其系数为$C^{c}_{a}C^{d}_{b}$，由于是同类项，显然存在 $C^{m}_{n} \equiv C^{c}_{a}C^{d}_{b}$

又∵$n = ap+b , m = cp+d$ , 显然$a = \lfloor \frac{n}{p}\rfloor $ , $c = \lfloor \frac{m}{p}\rfloor $

$b = n \ mod\ p , d = m \ mod \ p$

带入系数得证

 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


const int N = 1e5+10;
const int P = 1e9+7;
typedef long long ll;
ll f[N];ll g[N];
ll qp(ll a,ll b){// ll, long long 
	ll ans = 1;
	while(b){
		if(b&1) ans = ans*a % P;
		a = a * a % P;
		b>>=1;
	}
	return ans;
}
//O(nlogn)
void gen(){
	f[0] = g[0] = 1; // 0 ! = 1;
	for(int i = 1;i<N;++i){
		f[i] = f[i-1]*i % P;
		g[i] = g[i-1]*qp(i,P-2)% P;
	}
}

int main() {
	//fast IO
 	gen();
    int m , n ; cin>>n>>m;
    cout<< f[n]%P*g[n-m] % P *g[m] % P <<'\n';
    return 0;
}

高精度+筛

给定两个整数 $n$ ,$m$ ($1\le m\le n \le 10000$) 求$C^{m}_{n}$

组合数增长是非常快 $C^{50}_{1000}$， 30位数 $C^{5000}_{1000}$，300位数$C^{5000}_{10000}$ 3009位数

所以这个范围只能配合高精度

把 1~n质数筛选出
枚举所有质数

a. 计算$C^{n}_{m} = \frac{n!}{(n-m)!m!}中的质数个数$

b. 高精度计算$C[i] \times p$

由勒让德定理得

$n!$中 $p$的个数为 $s = \frac{n}{p}+\frac{n}{p^{2}}+\cdots$

计算完每个单元则整体的质数个数满足？

显然由于上下会约掉一些数，所以质数也会相应的减少，且n>=m>=1

用分母的个数减去分子的个数即可

 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


//不开long long 见祖宗 ,..血的教训!
const int N = 1e5+10;
ll f[N];
ll g[N];
ll qp(ll a,ll b,int p){
	ll ans = 1; 
	while(b){
		if(b&1) ans = ans*a%p;
		a = a*a%p;
		b>>=1;
	}
	return ans;
}
void init(int p){
	f[0] = g[0] = 1;
	for(int i =1;i<N;++i){
		f[i] = f[i-1] * i%p;
		g[i] = g[i-1] * qp(i,p-2,p) %p ;
	}
}
ll get_c(int n, int m, int p){
	return (f[n]%p * g[m]%p * g[n-m]%p)%p;
}
ll lucas(ll n, ll m,int p){
	if(m == 0) return 1;
	return lucas(n / p , m / p, p) * get_c(n % p , m % p , p) % p;
}
int lucas2(LL a, LL b, int p){
    if (a < p && b < p) return get_c(a, b, p);
    return (LL)get_c(a % p, b % p, p) * lucas(a / p, b / p, p) % p;
}

void solve(){
    int _ ;
    cin>>_;
    while(_--){
        ll n , m , p;
        cin>>n>>m>>p;
        init(p);
        cout<<lucas(n,m,p)<<'\n';
    }
}

待完成

扩展Lucas