组合数

最新推荐文章于 2021-03-10 02:29:56 发布

CrystalDavid

最新推荐文章于 2021-03-10 02:29:56 发布

阅读量283

点赞数

分类专栏：数论

本文链接：https://blog.csdn.net/CrystalDavid/article/details/77170721

版权

数论专栏收录该内容

5 篇文章 0 订阅

订阅专栏

杨辉三角，复杂度O( $n^2$ )，可以用来预处理除2000以内的组合数

const int maxn = 2000+10;
const long long mod = 1e9+7;
long long comb[maxn][maxn];
void init()
{
    for(int i=0;i<maxn;i++){
        comb[i][0] = comb[i][i] = 1;
        for(int j=1;j<i;j++){
            comb[i][j] = comb[i-1][j-1] + comb[i-1][j];
            comb[i][j] %= mod;
        }
    }
}

2.逆元求解，复杂度O(n)，可以处理1e5以内的数据

const ll maxn = 1e5+10;
const ll mod = 1e9+7;
ll F[maxn],Finv[maxn],inv[maxn];//F是阶乘，Finv是阶乘的逆
void init()
{
    inv[1] = 1;
    for(int i=2;i<maxn;i++)
        inv[i] = (mod - mod/i)*inv[mod%i]%mod;
    F[0] = Finv[0] = 1;
    for(int i=1;i<maxn;i++){
        F[i] = F[i-1]*i%mod;
        Finv[i] = Finv[i-1]*inv[i]%mod;
    }
}
ll comb(ll n,ll m)
{
    if(m<0 || m>n)
        return 0;
    return F[n]*Finv[n-m]%mod*Finv[m]%mod;
}

3.卢卡斯定理，用于求解1< n ,m <1e18，1< p <1e5且p为质数的情况
C(n, m) % p = C(n / p, m / p) * C(n%p, m%p) % p

ll Lucas(ll n,ll m,ll p)
{
    if(m==0)
        return 1;
    return Lucas(n/p,m/p,p)*comb(n%p,m%p)%p;
}

Lucas定理在不同的题目中有不同的运用方法，唯本质如一。如果题目中的模数p不变，就可以预处理1~p的所有组合数。如果每组样例会提供模数，那可以只预处理阶乘，再用快速幂取模求逆元。

例题：http://acm.hdu.edu.cn/showproblem.php?pid=3037

题意大略是在n棵树上储存不超过m棵豆子的种类数。
过程就不说，结果是C(n+m,m),直接套Lucas的板子

#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
const ll maxn = 1e5+10;
int F[maxn];

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

void init(ll p)
{
    F[0] =1;
    for(ll i=1;i<=p;i++) //注意，中间结果可能爆int
        F[i] = F[i-1]*i%p;
}

ll inv(ll a,ll p)
{
    return pow_mod(a,p-2,p);
}

ll comb(ll n,ll m,ll p)
{
    if(m<0||m>n)
        return 0;
    return F[n]*inv(F[n-m],p)%p*inv(F[m],p)%p;
}

ll Lucas(ll n,ll m,ll p)
{
    if(m==0)
        return 1;
    return Lucas(n/p,m/p,p)*comb(n%p,m%p,p)%p;
}

int main()
{
    int T;
    ll n,m,p;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%I64d %I64d %I64d",&n,&m,&p);
        init(p);
        printf("%I64d\n",Lucas(n+m,m,p));
    }
}

4.特殊情况：1< m,n< 1e6，1< p <1e5，且p有可能是合数。
根据惟一分解定理暴力分解，再使用快速幂取模。
由于C(x,y)=x!/(y!* (x-y)!)，这里我们可以将x!分解素因子，并保存记录下来，同样的方法记录后面两个，由于x!必然能够整除(y!*(x-y)!)，所以后面两个数有的因子，x!比然有，只需要将他们的因子的指数相加减，就能得到最后结果的素因子分解的情况，然后最后使用快速幂取模，就能得到最后的结果。

注意：如何进行素因子分解？
首先要打表将所有的素因子求出来，这里有是将n！进行素因子分解，假设想要求出其中有多少个5，这里是有技巧的。
假设n=200，那么因子5的个数=200/5+40/5+8/5=49，怎么得到的呢？200中5的倍数有40个，这40个数中其中是25的倍数的有8个，所以还能分解出8个5，这8个数中还有一个是125的倍数，还能分解出一个5，就这样一直循环下去，就能求出指数的值。

//贴上大佬的模板
#include <iostream>
#include <string.h>
#include <stdio.h>
using namespace std;
typedef long long LL;
const int N = 200005;

bool prime[N];
int p[N];
int cnt;
void isprime()
{
    cnt = 0;
    memset(prime,true,sizeof(prime));
    for(int i=2; i<N; i++)
        if(prime[i])
        {
            p[cnt++] = i;
            for(int j=i+i; j<N; j+=i)
                prime[j] = false;
        }
}
LL quick_mod(LL a,LL b,LL m)
{
    LL ans = 1;
    a %= m;
    while(b)
    {
        if(b & 1)
        {
            ans = ans * a % m;
            b--;
        }
        b >>= 1;
        a = a * a % m;
    }
    return ans;
}

LL Work(LL n,LL p)
{
    LL ans = 0;
    while(n)
    {
        ans += n / p;
        n /= p;
    }
    return ans;
}

LL Solve(LL n,LL m,LL P)
{
    LL ans = 1;
    for(int i=0; i<cnt && p[i]<=n; i++)
    {
        LL x = Work(n, p[i]);
        LL y = Work(n - m, p[i]);
        LL z = Work(m, p[i]);
        x -= (y + z);
        ans *= quick_mod(p[i],x,P);
        ans %= P;
    }
    return ans;
}

int main()
{
    int T;
    isprime();
    cin>>T;
    while(T--)
    {
        LL n,m,P;
        cin>>n>>m>>P;
        n += m - 2;
        m--;
        cout<<Solve(n,m,P)<<endl;
    }
    return 0;
}