大组合数的模板

来源:互联网 发布:java web项目开发2017 编辑:程序博客网 时间:2024/05/18 04:12

//C(n,m) % MOD n, m<=10^6
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long LL;
const int MAXN = 2000001;
const int MAXP = 150000;
const int MOD = 364875103;
int n, m, num, prime[MAXP];
bool isprime[MAXN + 1];

int sieve(int n) {
    int p = 0;
    memset(isprime, 1, sizeof(isprime));
    for (int i = 2; i <= n; i++) {
        if (isprime[i])
            prime[p++] = i;
        for (int j = 0; j < p && i * prime[j] <= n; j++) {
            isprime[i * prime[j]] = 0;
            if (i % prime[j] == 0)
                break;
        }
    }
    isprime[0] = isprime[1] = 0;
    return p;
}

LL mul(LL a, LL b, LL c) { // a*b % c
    LL ret = 0, tmp = a % c;
    for (; b; b >>= 1) {
        if (b & 1)
            if((ret += tmp) >= c)
                ret -= c;
        if ((tmp <<= 1) >= c)
            tmp -= c;
    }
    return ret;
}

LL power(LL a, LL b, LL c) { // a^b % c
    LL res = 1;
    for (; b; b >>= 1) {
        if (b & 1) res = mul(res, a, c);
        a = mul(a, a, c);
    }
    return res;
}

int addfact(int n, int cnt[], int d) {
    int i;
    for (i = 0; i < num && prime[i] <= n; i++) {
        int x = n;
        while (x) {
            cnt[i] += d * (x / prime[i]);
            x /= prime[i];
        }
    }
    return i;
}
LL C(int n, int m) {
    int p, cnt[MAXP] = {0};
    p = addfact(n, cnt, 1);//分子n!
    addfact(m, cnt, -1);//分母m!
    addfact(n - m, cnt, -1);//分母(n-m)!
    LL ret = 1;
    for (int i = 0; i < p; i++)
        ret = mul(ret, power(prime[i], cnt[i], MOD), MOD);//mul(ret * power(prime[i],cnt[i],MOD)) % MOD;
    return ret % MOD;
}

int id[16];

int main() {
    //freopen("in.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
    num = sieve(MAXN);
    //cout << num << endl;
    /*while (scanf("%d%d", &n, &m) != EOF) {
        printf("%I64d\n", C(n, m));
    }*/
    int cas, n, m;
    scanf("%d", &cas);
    for (int T = 1; T <= cas; T++) {
        scanf("%d%d", &n, &m);
        int k;
        scanf("%d", &k);
        for (int i = 0; i < k; i++) {
            scanf("%d", &id[i]);
        }
        id[k] = m;
        int cnt;
        long long ans = 1;
        for (int i = 1; i <= k; i++) {
            cnt = id[i] - id[i - 1];
            ans = (ans * C(n + cnt - 1, cnt)) % MOD;
        }
        printf("Case #%d: %I64d\n", T, ans);
    }
}

分析:
首先观察全是12222……这种形态的for结构。
也就是如下结构:
for (int i = 0; i < n; i++) {
for (int j = i; j < n; j++) {
for (int k = j; k < n; k++) {
……
}
}
}
这样的是平时写程序经常用到的一种循环结构,这个简化的问题等价于n种不同颜色的球,可放回地取m个,问形成最后的结果集合不同的种类数。
为什么等价呢?因为由于是集合,我们可以用排序后的数列来表示一个集合,那么对于每个取球方案都可以对应到一个数列,如:
3种不同颜色的球,可放回地取3个,共有111,112,113,122,123,133,222,223,233,333这10种不同的方案,显然因为方案是有序的,所以他们跟简化题意是可以一一对应的,比如223就对应a[0]=2,a[1]=2,a[2]=3的一种合法方案。
那么现在问题就是求这个简化问题的公式了,有很多方法可以求这个公式,这里介绍一下我的做法:首先列n+m个空位,现在我们固定最右边的空位里放的是第n种颜色的球,现在我们只要在剩下的n+m-1个位置里放入n-1个球,然后从左到右标记为1~n-1种颜色,即能得到一种方案,这个方案里每种颜色的球取的个数为它左边的连续的空位数。比如我们还是取n=m=3为例,列6个空位:oooooo,最右边填上颜色3:ooooo3,现在在剩下的空位里放入颜色1,2的球,比如有方案:o12oo3,则表示1取了1个,2取了0个,3取了2个,也就是对应133这种方案。至于为什么最右边要填第n种颜色的球是因为我们要现在令的是每个球左边连续的空位数是它对应取的个数,所以最右边的球右边不能有空位,故位置固定。所以现在问题就十分简单了,公式为C(n+m-1,n-1)=C(n+m-1,m)。

接下来搞的是一般结构,由于一共有最多15个1-type的for-loop,所以只要每个以1type开头后面跟若干2type形成的组合结构单独处理,15个组合结构答案相乘即可,因为每个1-type重新开头后答案跟之前都是独立的。

最后的问题就是如何计算C(n+m-1,m)%MOD,其中0<=n+m-1<=1100000,0<=m<=100000,MOD是两个质数的乘积:97*3761599。
这里我的方法是处理任意MOD数的问题,由于上限是1000000,求出在1000000里面约有15w个质数,因为幂项相除等于指数相减,对于每个组合数我们算得它分解后的质数里面每一项的幂v,对于分子我们将幂+=v,对于组合数分母的两项,我们将幂-=v,所以一共最多15w项质数,最后就没有乘法了,可以二分快速幂取模得到答案。首先筛法得到质数后,每组Case计算量约为15w*15*log(n+m-1)。即可得到答案。
另一个更好的方法是用Lucas定理,由于364875103=97*3761599,于是我们可以分别对这两个质数进行组合数的计算 ,然后用中国剩余定理求出对于 MOD 的余数,Lucas的每次计算组合数的复杂度大概是 O(log n * log p),于是每 组Case计算量大概是15 * log n * log p。