我们都知道组合数的公式是:
看起来直接写程序也是挺简单的。
但是我们能想出两个问题,溢出怎么办? 速度慢怎么办?
第一种情形
如果我直接告诉你,我要询问k次,每次给我求C(n,m), n,m都小于2000,我们怎么算呢?
对于每次询问,我们都取计算阶乘,这确实可以,但是太慢。
我们采用预处理的办法,讲预处理之前,我们要了解一个公式。
这个公式也很好理解,是一个动态规划的思想,对于第n个物品,我们选或者不选分别对应于式子右边的第一项和第二项,选第n个物品,那就从剩下的n-1个物品里面选m-1个;不选第n个物品,那就从剩下n-1个物品里面选m个。
由这个式子我们得出了动态规划的方法。
c[i][j] = c[i-1][j-1] + c[i-1][j]
例题为:
题目描述:
给定n组询问,每组询问给定两个整数a,b,请你输出C(a,b) mod (10^9+7)的值。
输入格式
第一行包含整数n。接下来n行,每行包含一组a和b。
输出格式
共n行,每行输出一个询问的解。
数据范围
1≤n≤10000,1≤b≤a≤2000
输入样例:
3
3 1
5 3
2 2
输出样例:
3
10
1
代码如下:
#include <iostream>
using namespace std;
const int maxn = 2005,mod = 1e9 + 7;
int c[maxn][maxn];
int main(){
int n,a,b;
scanf("%d",&n);
for(int i = 0;i <= 2000;i++){
for(int j = 0;j <= i;j++){
if(!j) c[i][j] = 1;
else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
}
}
while(n--){
scanf("%d%d",&a,&b);
printf("%d\n",c[a][b]);
}
return 0;
}
第二种情况
数据范围
1≤n≤10000,
1≤b≤a≤10^5
mod = 1e9 + 7
我们发现a和b的范围扩大到了1e5,所以不能使用动态规划的方法,时间复杂度与空间复杂度均超过10^9,然而发现题目中要求模为一个素数,但是除法是没有取模的(只有加,减,乘有类似,那么我们必须将除法改为乘法逆元的形式。由于模的是一个质数,所以我们可以用费马小定理去求这个逆元。
在这种情况下,我们保存阶乘的结果,和逆元的结果。
阶乘和逆元的更新公式为
int fact[N],infact[N] ;
fact[0] = infact[0] = 1;
for(int i = 1 ; i < N ; i++){
fact[i] = (ll)fact[i-1] * i % mod ; //预处理前5000的阶乘和逆元
infact[i] = (ll)infact[i-1] * qmi(i,mod - 2,mod) % mod ;
}
为什么逆元可以这样算呢?
依照上述定义,infact[i] = qmi(i的阶乘,m - 2,m) % m;
但是i的阶乘=(i-1)! * i;
所以infact[i]=((i - 1)! * i)^(m - 2) % m;把式子展开,则等价于:
[( ((i - 1)!) ^ (m - 2) % m ) * (i ^ (m - 2) ) % m] % m;
而依照infact数组的定义:infact[i]表示i!的逆元
可以知道( ((i - 1)!) ^ (m - 2) % m ) = infact[i - 1];
所以有infact[i] = (LL)infact[i - 1] * qmi(i,mod - 2,mod) % mod;
其中(LL)是为了防止爆int,最后在%m是因为 (a * b) % m = (a % m * b % m ) % m;
题目:
题目描述:
给定n组询问,每组询问给定两个整数a,b,请你输出C(a,b) mod (10^9+7)的值。
输入格式
第一行包含整数n。接下来n行,每行包含一组a和b。
输出格式
共n行,每行输出一个询问的解。
数据范围
1≤n≤100000,1≤b≤a≤10^5
输入样例:
3
3 1
5 3
2 2
输出样例:
3
10
1
代码:
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
const int maxn = 100005,mod = 1e9 + 7;
int fact[maxn],infact[maxn];
int qmi(int a,int k,int p){
int res = 1;
while(k){
if(k & 1) res = (ll) res * a % p;
k >>= 1;
a = (ll)a * a % p;
}
return res;
}
int main(){
fact[0] = infact[0] = 1;
for(int i = 1;i <= 100000;i++){
fact[i] = (ll)fact[i - 1] * i % mod;
infact[i] = (ll)infact[i - 1] * qmi(i,mod - 2,mod) % mod;
}
int n;
scanf("%d",&n);
while(n--){
int a,b;
scanf("%d%d",&a,&b);
printf("%d\n",(ll)fact[a] * infact[b] % mod * infact[a - b] % mod);
}
return 0;
}
第三种情况
数据范围
1≤n≤20,
1≤b≤a≤10^18,
1≤p≤10^5,且p均为素数
我们可以看到这里a,b的范围太大了,完全无法按照我们之前预处理的思路来。
这就意味着我们必须直接算。
解决办法是,Lucas定理:
用代码表示是:
ll lucas(ll a,ll b){ // lucas定理递归处理
if(a < p && b < p) return C(a,b) ;
return lucas(a/p,b/p) * C(a%p,b%p) % p ;
}
这些数论的东西我就不证明了,一般也没人看,会用就行。
题目:
题目描述:
给定n组询问,每组询问给定三个整数a,b,p,其中p是质数,请你输出C(a,b) mod p的值。
输入格式
第一行包含整数n。接下来n行,每行包含一组a,b,p。
输出格式
共n行,每行输出一个询问的解。
数据范围
1≤n≤20,1≤b≤a≤10^18,1≤p≤10^5,
输入样例:
3
5 3 7
3 1 5
6 4 13
输出样例:
3
3
2
代码:
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
int qmi(int a,int k,int p){
int res = 1;
while(k){
if(k & 1) res = (ll)res * a % p;
a = (ll)a * a % p;
k >>= 1;
}
return res;
}
int C(int a,int b,int p){
if(b > a) return 0;
int res = 1;
for(int i = 1,j = a;i <= b;i++,j--){
res = (ll)res * j % p;
res = (ll)res * qmi(i,p - 2,p) % p;
}
return res;
}
int lucas(ll a,ll b,int p){
if(a < p && b < p) return C(a,b,p);
return (ll)C(a % p,b % p,p) * lucas(a / p,b / p,p) % p;
}
int main(){
int n,p;
ll a,b;
cin>>n;
while(n--){
cin>>a>>b>>p;
cout<<lucas(a,b,p)<<endl;
}
return 0;
}
特别地,如果p不是质数,那么我们就用第四种办法来算。
第四种情况
数据范围
1≤b≤a≤5000 (无取模,需要高精度)
注意这里没有取模,没有取模就要考虑是不是高精度的问题。
但是直接算显然不是特别好的办法, 因为你要不只要实现高精度乘高精度,还要实现高精度除高精度。
我们知道任何数都能分解成质因数的乘积。
特别地,如果这个数是n!,我们还有另外一种巧妙办法。
我们用n依次除以p的k次幂就能得到质因数p出现了多少次
代码如下:
int get(int n, int p)
{
int res = 0;
while (n)
{
res += n / p;
n /= p;
}
return res;
}
我们怎么得到小于n的质数呢?用线性筛
void init(int n){ //线性筛筛出素数
for(int i = 2; i <= n ; i++){
if(!st[i]) primes[cnt++] = i ;
for(int j = 0 ; primes[j] * i <= n ; j++){
st[primes[j] * i] = 1 ;
if(i % primes[j] == 0) break ;
}
}
}
接下来就是高精度乘低精度了
vector<int> mul(vector<int> &A, int b) // C = A * b, A >= 0, b >= 0
{
vector<int> C;
int t = 0;
for (int i = 0; i < A.size() || t; i ++ )
{
if (i < A.size()) t += A[i] * b;
C.push_back(t % 10);
t /= 10;
}
while (C.size() > 1 && C.back() == 0) C.pop_back();
return C;
}
题目:
输入a , b 求c(a,b)需要用高精度。
输入格式:
共一行,包含两个整数a 和b 。
输出格式:
共一行,输出c(a,b)的值。
数据范围:
1 ≤ b ≤ a ≤ 5000
代码:
#include <iostream>
#include <cstring>
#include <algorithm>
#include <vector>
using namespace std;
const int N = 5010;
int primes[N], cnt, sum[N];
bool st[N];
void get_primes(int n) // 线性筛质数
{
for (int i = 2; i <= n; i ++ )
{
if (!st[i]) primes[cnt ++ ] = i;
for (int j = 0; primes[j] <= n / i; j ++ )
{
st[primes[j] * i] = true;
if (i % primes[j] == 0) break;
}
}
}
int get(int n, int p)
{
int res = 0;
while (n)
{
res += n / p;
n /= p;
}
return res;
}
vector<int> mul(vector<int> &A, int b) // C = A * b, A >= 0, b >= 0
{
vector<int> C;
int t = 0;
for (int i = 0; i < A.size() || t; i ++ )
{
if (i < A.size()) t += A[i] * b;
C.push_back(t % 10);
t /= 10;
}
while (C.size() > 1 && C.back() == 0) C.pop_back();
return C;
}
int main()
{
int a, b;
cin >> a >> b;
get_primes(a);
for (int i = 0; i < cnt; i ++ )
{
int p = primes[i];
sum[i] = get(a, p) - get(b, p) - get(a - b, p);
}
vector<int> res;
res.push_back(1);
for (int i = 0; i < cnt; i ++ )
for (int j = 0; j < sum[i]; j ++ )
res = mul(res, primes[i]);
for (int i = res.size() - 1; i >= 0; i--) printf("%d",res[i]);
puts("");
return 0;
}