HDU 4767 Bell (中国剩余定理)

题意：

Bell(n)是基数为n的集合的划分方法数。例如n = 3

{1, 2, 3}

{1}  {2, 3}

{1, 2} {3}

{1, 3} {2}

{1} {2} {3}

所以Bell(3) = 5

给你一个n，求Bell(n) % 95041567的值

解题思路：

首先知道Bell(n) = S(n, 1) + S(n, 2) + ... S(n, n)，S为第二类斯特灵数

然后google到结论。。 ：若p为任意质数，则有Bell(n+p) = Bell(n) + Bell(n+1) (mod  p)。

然后可以发现95041567分解质因子后为 31 * 37 * 41 * 43 * 47...都为质数。

利用矩阵二分幂可以很容易得到 Bell(n) % p (p = 31, 37, 41, 43, 47)的各个值，现在问题就变成了这种模型

X = a[i] (mod m[i])  (m[i]为质数)，求X mod (m[1]*m[2]...m[k])的值。也就是中国剩余定理的应用了。具体见代码

#include <stdio.h>

int Md[] = {31, 37, 41, 43, 47};
int s[6][55][55], w[6][55];

// 第二类斯特灵数的预处理
void init() {
    for(int i = 0;i < 5;i ++)
        s[i][0][0] = 1;
    for(int i = 1;i <= 50; i++) {
        for(int j = 0;j < 5;j ++)
            s[j][i][1] = 1;
        for(int j = 1;j <= i; j++) {
            for(int k = 0;k < 5; k++)
                s[k][i][j] = (s[k][i-1][j-1] + j*s[k][i-1][j])%Md[k];
        }
        for(int j = 0;j < 5; j++) {
            w[j][i] = 0;
            for(int k = 1;k <= i; k++)
               w[j][i] = (w[j][i] + s[j][i][k])%Md[j];
        }
    }
}

// 矩阵二分幂
int pow_mod(int id, int n, int mod) {
    int a[55], aa[55], q[55][55], qq[55][55];
    int sz = Md[id] + 1;
    for(int i = 1;i <= sz+1; i++)
        for(int j = 1;j <= sz+1; j++)
            q[i][j] = 0;
    for(int i = 1;i < sz; i++)
        q[i+1][i] = 1;
    q[2][sz] = q[3][sz] = 1;
    for(int i = 1;i <= sz; i++)
        a[i] = w[id][i];
    if(n <= Md[id]) return a[n];
    n --;
    while(n) {
        if(n & 1) {
            for(int i = 1;i <= sz; i++) {
                aa[i] = 0;
                for(int j = 1;j <= sz; j++) {
                    aa[i] += a[j]*q[j][i];
                }
            }
            for(int i = 1;i <= sz; i++)
                a[i] = aa[i]%mod;
        }
        for(int i = 1;i <= sz; i++) {
            for(int j = 1;j <= sz; j++) {
                qq[i][j] = 0;
                for(int k = 1;k <= sz; k++) {
                    qq[i][j] += q[i][k]*q[k][j];
                }
                qq[i][j] %= mod;
            }
        }
        for(int i = 1;i <= sz; i++)
            for(int j = 1;j <= sz; j++)
            q[i][j] = qq[i][j];
        n /= 2;
    }
    return a[1];
}
// 扩展欧几里得
int exgcd(int a, int b, int &x, int &y) {
    if(!b) {
        x = 1; y = 0;
        return a;
    }
    int ret = exgcd(b, a%b, y, x);
    y -= a/b*x;
    return ret;
}

// 中国剩余定理
int china(int n, int a[], int m[]) {
    int M = 1;
    for(int i = 0;i < n; i++) M *= m[i];
    int ret = 0;
    for(int i = 0;i < n; i++) {
        int w = M/m[i], x, y;
        int d = exgcd(w, m[i], x, y);
        ret = (ret + x*w*a[i])%M;
    }
    return (ret + M)%M;
}

int x[55];
int main() {
    init();
    int t, n;
    scanf("%d", &t);
    while(t--) {
        scanf("%d", &n);
        for(int i = 0;i < 5;i ++) {
            // 求Bell(n) % 各个质数的值
            x[i] = pow_mod(i, n, Md[i]);
        }
        int ans = china(5, x, Md);
        printf("%d\n", ans);
    }
}