HDU 3944

最新推荐文章于 2019-07-09 23:45:45 发布

RaceBug2010

最新推荐文章于 2019-07-09 23:45:45 发布

阅读量1k

点赞数

文章标签： c 算法

本文链接：https://blog.csdn.net/RaceBug2010/article/details/6754751

版权

原来剑哥讲组合数取模的时候说过Lucas定理，但是不是很明白。总之C(n,k)%p = C(n/p,k/p)*C(n%p,k%p)%p。这里要注意如果k>n, C(n,k)=0。C(n,k)%p可以通过求逆元的方法计算。不过C(n,k)%p = n!/(k!*(n-k)!)%p = n!*((k!*(n-k)!)^(p-2))%p是更好的方法，编码简单一点，时间效率高一点。

还要注意一点，这题比较恶心，每次都初始化阶乘数组肯定会超时。考虑到10000以内只有1229个素数，所以可以先打表。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <memory.h>
#include <string>
#include <vector>
#include <map>
#include <set>
#include <queue>
#include <algorithm>
#include <iostream>
#include <sstream>

using namespace std;

// Lucas定理求解C(n,k)%p, p是素数
// C(n,k)%p = C(n%p,k%p)*C(n/p,k/p)%p
// C(n,k)%p = n!/(k!*(n-k)!)%p = n!*((k!*(n-k)!)^(p-2))%p
// pow_mod(a,b,p)求a^b%p的值, 快速幂取模算法
// com_mod(n,k,p)求C(n,k)%p的值
// lucas(n,k,p)使用lucas定理求C(n,k)%p的值

int pow_mod(int a, int b, int p) {
    int r = 1, t = a;
    while (b) {
        if (b & 1) r = r * t % p;
        t = t * t % p;
        b >>= 1;
    }
    return r;
}

int com_mod(int n, int k, int p, int *fac) {
    int a, b;
    if (k > n) return 0;
    a = fac[n];
    b = fac[n-k] * fac[k] % p;
    return a * pow_mod(b, p-2, p) % p;
}

int lucas(int n, int k, int p, int *fac) {
    int r = 1, a, b;
    while (n && k && r) {
        a = n % p;
        b = k % p;
        n /= p;
        k /= p;
        
        r = r * com_mod(a, b, p, fac) % p;
    }
    return r;
}

int pri[1250], cnt;
int fac[1250][10000];
bool not_prime[10005];

void init() {
    int i, j, p;
    cnt = 0;
    for (i = 2; i < 10000; i++) {
        if (not_prime[i]) continue;
        pri[cnt++] = i;
        for (j = i + i; j < 10000; j += i)
            not_prime[j] = 1;
    }
    for (i = 0; i < cnt; i++) {
        fac[i][0] = 1;
        p = pri[i];
        for (j = 1; j < p; j++)
            fac[i][j] = j * fac[i][j-1] % p;
    }
}

int *find(int p) {
    int l = 0, r = cnt - 1;
    while (l <= r) {
        int m = (l + r) / 2;
        if (pri[m] == p) return fac[m];
        else if (pri[m] < p) l = m + 1;
        else r = m - 1;
    }
    return fac[0];
}

int main() {
    int n, k, p, t = 1;
    init();
    while (scanf("%d %d %d", &n, &k, &p) != EOF) {
        if (k > n - k) k = n - k;
        printf("Case #%d: ", t++);
        printf("%d\n", (lucas(n+1, k, p, find(p)) + n - k) % p);
    }
    return 0;
}