组合数学 - 求组合数(递推 + 定义 + Lucas定理 + 高精度)

组合数学 - 求组合数(递推 + 定义 + Lucas定理)

1、递推法 - $O(n^2)$

$给定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

---

分析：

$a,b<=2000，可以用递推直接求出来，利用公式：$

$$C_n^m=C_{n-1}^m+C_{n-1}^{m-1}$$

代码：

#include<iostream>
#include<cstdio>

using namespace std;

const int N = 2010 , mod = 1e9+7;

int n,C[N][N];

void cal()
{
    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;
}

int main()
{
    cal();
    cin>>n;
    
    while(n--)
    {
        int a,b;
        scanf("%d%d",&a,&b);
        printf("%d\n",C[a][b]);
    }
    
    return 0;
}

---

2、定义法 - $O(nlogn)$

$给定n组询问，每组询问给定两个整数a，b，请你输出C_a^{b}mod(10^9+7)的值。$

输入格式
第一行包含整数n。

接下来n行，每行包含一组a和b。

输出格式
共n行，每行输出一个询问的解。

数据范围
1≤n≤10000,
1≤b≤a≤10⁵

输入样例：
3
3 1
5 3
2 2
输出样例：
3
10
1

---

分析：

$a,b<=100000，直接利用定义：$
$$C_n^m=\frac{n!}{m!(n-m)!}$$

$先预处理1$~$100000的所有阶乘的结果，直接带入定义式。$

$除数取模问题，由于10^9+7是质数，可以使用快速幂+费马小定理求逆元。$

$$\begin{gathered} 共预处理两个数组 \\ ①、fact[i]=i!=fact[i-1]\times i\%mod， \\ ②、infact[i]=fact[i{]}^{-1}=infact[i-1]\times i^{-1}\%mod=infact[i-1]\times i^{mod-2}\%mod。 \end{gathered}$$

代码：

#include<iostream>
#include<cstdio>

#define ll long long

using namespace std;

const int N = 100010, mod = 1e9+7;

ll fact[N],infact[N];
int n;

ll quick_pow(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}

int main()
{
    fact[0]=infact[0]=1;
    for(int i=1;i<N;i++)
    {
        fact[i]=fact[i-1]*i%mod;
        infact[i]=infact[i-1]*quick_pow(i,mod-2)%mod;
    }

    cin>>n;
    while(n--)
    {
        int a,b;
        scanf("%d%d",&a,&b);
        printf("%d\n",fact[a]*infact[b]%mod*infact[a-b]%mod);
    }
    
    return 0;
}

---

3、Lucas定理 - $O(lo{g}_{p}N\cdot (p+lo{g}_{2}p))$

$给定n组询问，每组询问给定三个整数a,b,p，其中p是质数，请你输出C_a^{b}modp的值。$

输入格式
第一行包含整数n。

接下来n行，每行包含一组a,b,p。

输出格式
共n行，每行输出一个询问的解。

数据范围
1≤n≤20,
1≤b≤a≤10¹⁸,
1≤p≤10⁵,

输入样例：
3
5 3 7
3 1 5
6 4 13
输出样例：
3
3
2

---

分析：

$a,b<=10^{18}，p<=10^5，采用卢卡斯定理：$

$$\begin{gathered} 设a=a_k{p}^k+a_{k-1}{p}^{k-1}+...+a_0{p}^0，b=b_k{p}^k+b_{k-1}{p}^{k-1}+...+b_0{p}^0, \\ 则有： \end{gathered}$$

$$C_a^b\equiv C_{a_k}^{b_k}\cdot C_{a_{k-1}}^{b_{k-1}}\cdot ...\cdot C_{a_0}^{b_0}\cdot (modp)$$

$$\begin{gathered} 证明： \\ 由p为素数，可知对1<=j<=p-1，都有C_p^j=\frac{p}{j}{C}_{p-1}^{j-1}\equiv 0(modp)， \end{gathered}$$

$于是，$
$$(1+x{)}^p=1+C_p^{1}x+...+C_p^{p-1}{x}^{p-1}+x^p\equiv 1+x^p(modp).$$

$利用上述结果，可知$
$$(1+x{)}^p=(1+x{)}^{a_0}((1+x{)}^p{)}^{a_1}...((1+x{)}^{p^k}{)}^{a_k}\equiv (1+x{)}^{a_0}(1+x^p{)}^{a_1}...(1+x^{p^k}{)}^{a_k}(modp).$$

$对比两边x^b项的系数，可得$
$$C_a^b\equiv C_{a_k}^{b_k}\cdot C_{a_{k-1}}^{b_{k-1}}\cdot ...\cdot C_{a_0}^{b_0}\cdot (modp)$$

具体实现可以通过递归：

$由于a=a_k{p}^k+a_{k-1}{p}^{k-1}+...+a_0{p}^0，b=b_k{p}^k+b_{k-1}{p}^{k-1}+...+b_0{p}^0,$

$只有最后项的系数a_0和b_0不含p的正整数次幂，故C_a^b对p取模后只剩下C_{a_0}^{b_0}，此时我们将a和b都除去一个p，得到：$

$a=a_k{p}^k+a_{k-1}{p}^{k-1}+...+a_1{p}^0，b=b_k{p}^k+b_{k-1}{p}^{k-1}+...+b_1{p}^0,$

$此时C_a^b对p取模能够得到C_{a_1}^{b_1}，这样递归k+1次就能够求得C_{a_k}^{b_k}\cdot C_{a_{k-1}}^{b_{k-1}}\cdot ...\cdot C_{a_0}^{b_0}\cdot (modp).$

$关系式：C_a^b\equiv C_{a\%p}^{b\%p}\cdot C_{a/p}^{b/p}(modp)$

代码：

#include<iostream>

#define ll long long 

using namespace std;

const int N=100010;

int quick_pow(ll a,ll b,int p)
{
    int ans=1;
    while(b)
    {
        if(b&1) ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}

ll C(int a,int b,int p)
{
    if(b>a) return 0;
    
    ll up=1, down=1;
    for(int i=1,j=a;i<=b;i++,j--)
    {
        up = up*j%p;
        down = down*i%p;
    }
    return up*quick_pow(down,p-2,p)%p;
}

int lucas(ll a,ll b,int p)
{
    if(a<p&&b<p) return C(a,b,p);
    return C(a%p,b%p,p)*lucas(a/p,b/p,p)%p;
}

int main()
{
    int n;
    cin>>n;
    
    while(n--)
    {
        ll a,b;
        int p;
        cin>>a>>b>>p;
        cout<<lucas(a,b,p)<<endl;
    }
    
    return 0;
}

---

4、定义 + 高精度 - $O(n)$

$输入a,b，求C_a^b的值。$

注意结果可能很大，需要使用高精度计算。

输入格式
共一行，包含两个整数a和b。

输出格式
共一行，输出$C_a^b$的值。

数据范围
1≤b≤a≤5000

输入样例：
5 3
输出样例：
10

---

分析：

$根据定义：$
$$C_n^m=\frac{n!}{m!(n-m)!}$$

$$\begin{gathered} 若直接计算，需要同时实现高精度乘法、除法，效率较低，我们可以将C_n^m先分解质因数， \\ 即C_n^m=P_1^{a_1}{P}_2^{a_2}...P_k^{a_k}，这样就可以转化为乘法取模，且仅需用高精度乘法即可实现。 \end{gathered}$$

$$\begin{gathered} 对阶乘a!分解质因数p_i： \\ 阶乘a!中p_i的幂次su{m}_i=\lfloor \frac{a}{p_i}\rfloor +\lfloor \frac{a}{p_i^2}\rfloor +\lfloor \frac{a}{p_i^3}\rfloor +... \end{gathered}$$

$则组合数C_n^m=\frac{n!}{m!(n-m)!}的质因数的幂为分子的质因数的幂减去分母的对应的质因数的幂。$

说明：

$a!表示1$~$$\begin{gathered} a中所有数的乘积， \\ 区间[1,a]中，包含p_i^1的数共有\lfloor \frac{a}{p_i}\rfloor ， \\ 同时这当中可能有些数包含的是p_i^2，这些数对幂次的贡献是2，但是我们只算了1次，故我们再加上\lfloor \frac{a}{p_i^2}\rfloor ， \\ 以此类推再加上\lfloor \frac{a}{p_i^3}\rfloor ,\lfloor \frac{a}{p_i^4}\rfloor ... \end{gathered}$$

具体落实：

$①、筛出5000以内的所有质数。$

$②、预处理计算a和b的每个质因子的次数。$

$③、高精度计算C_n^m=P_1^{a_1}{P}_2^{a_2}...P_k^{a_k}。$

代码：

#include<iostream>
#include<vector>

using namespace std;

const int N=5010;

int prime[N],cnt;
int sum[N];
bool st[N];

void get_prime(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!st[i]) prime[cnt++]=i;
        for(int j=0;prime[j]<=n/i;j++)
        {
            st[i*prime[j]]=true;
            if(i%prime[j]==0) break;
        }
    }
}

int cal(int x,int p)
{
    int res=0;
    while(x)
    {
        res+=x/p;
        x/=p;
    }
    return res;
}

vector<int> mul(vector<int> a,int b)
{
    vector<int> c;
    int t=0;
    for(int i=0;i<a.size();i++)
    {
        t+=a[i]*b;
        c.push_back(t%10);
        t/=10;
    }
    while(t)
    {
        c.push_back(t%10);
        t/=10;
    }
    return c;
}

int main()
{
    int a,b;
    cin>>a>>b;
    
    get_prime(a);
    
    for(int i=0;i<cnt;i++)
    {
        int p=prime[i];
        sum[i]=cal(a,p)-cal(a-b,p)-cal(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,prime[i]);
            
    for(int i=res.size()-1;i>=0;i--) cout<<res[i];
    puts("");
    
    return 0;
}