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);    }}