HDU 4767 Bell (贝尔数中国剩余定理构造矩阵) ★ ★

最新推荐文章于 2020-04-02 23:44:09 发布

Ezereal

最新推荐文章于 2020-04-02 23:44:09 发布

阅读量908

点赞数

分类专栏：数学杂七杂八—矩阵快速幂

数学同时被 2 个专栏收录

37 篇文章

订阅专栏

杂七杂八—矩阵快速幂

10 篇文章

订阅专栏

斯特灵数[编辑]

维基百科，自由的百科全书

在组合数学，Stirling数可指两类数，都是由18世纪数学家James Stirling提出的。

第一类[编辑]

s(4,2)=11

第一类Stirling数是有正负的，其绝对值是 $n$ 个元素的项目分作 $k$ 个环排列的方法数目。常用的表示方法有 $s(n,k) , \left[\begin{matrix} n \\ k \end{matrix}\right]$ 。

换个较生活化的说法，就是有 $n$ 个人分成 $k$ 组，每组内再按特定顺序围圈的分组方法的数目。例如 $s(4,2)$ ：

{A,B},{C,D}
{A,C},{B,D}
{A,D},{B,C}
{A},{B,C,D}
{A},{B,D,C}
{B},{A,C,D}
{B},{A,D,C}
{C},{A,B,D}
{C},{A,D,B}
{D},{A,B,C}
{D},{A,C,B}

这可以用有向图来表示。

给定 $s(n,0)=0,s(1,1)=1$ ，有递归关系 $s(n+1,k)=s(n,k-1) + n \; s(n,k)$

递推关系的说明：考虑第n+1个物品，n+1可以单独构成一个非空循环排列，这样前n种物品构成k-1个非空循环排列，方法数为s(n,k-1)；也可以前n种物品构成k个非空循环排列，而第n+1个物品插入第i个物品的左边，这有n*s(n,k)种方法。

$| s(n,1) | =(n-1)!$
$s(n,k) = (-1)^{n+k} | s(n,k) |$
$s(n,n-1) = - C(n,2)$
$s(n,2) = (-1)^n (n-1)!\; H_{n-1}$
$s(n,3) = \frac{1}{2} (-1)^{n-1} (n-1)! [ (H_{n-1})^2 - H_{n-1}^{(2)} ]$

$H_n^{(m)}$ 是调和数的推广。

$s(n,k)$ 是递降阶乘多项式的系数：

$x^{\underline{n}}= x(x-1)(x-2)\ldots(x-n+1) = \sum_{k=1}^n s(n,k)x^k$

第二类[编辑]

第二类Stirling数是 $n$ 个元素的集定义k个等价类的方法数目。常用的表示方法有 $S(n,k) , S_n^{(k)} , \left\{\begin{matrix} n \\ k \end{matrix}\right\}$ 。

换个较生活化的说法，就是有 $n$ 个人分成 $k$ 组的分组方法的数目。例如有甲、乙、丙、丁四人，若所有人分成1组，只有所有人在同一组这个方法，因此 $S(4,1)=1$ ；若所有人分成4组，只可以人人独立一组，因此 $S(4,4)=1$ ；若分成2组，可以是甲乙一组、丙丁一组，或甲丙一组、乙丁一组，或甲丁一组、乙丙一组，或其中三人同一组另一人独立一组，即是：

{A,B},{C,D}
{A,C},{B,D}
{A,D},{B,C}
{A},{B,C,D}
{B},{A,C,D}
{C},{A,B,D}
{D},{A,B,C}

因此 $S(4,2)=7$ 。

给定 $S(n,n)=S(n,1)=1$ ，有递归关系 $S(n,k) = S(n-1,k-1) + k S(n-1,k)$
递推关系的说明：考虑第n个物品，n可以单独构成一个非空集合，此时前n-1个物品构成k-1个非空的不可辨别的集合，方法数为S(n-1,k-1)；也可以前n-1种物品构成k个非空的不可辨别的集合，第n个物品放入任意一个中，这样有k*S(n-1,k)种方法。
$S(n,n-1)=C(n,2)=n(n-1)/2$
$S(n,2)=2^{n-1} - 1$
$S(n,k) =\frac{1}{k!}\sum_{j=1}^{k}(-1)^{k-j} C(k,j) j^n$
$B_n=\sum_{k=1}^n S(n,k)$

$C(k,j)$ 是二项式系数，B_n是贝尔数。

两者关系[编辑]

$\sum_{n=0}^{\max\{j,k \}} S(n,j) s(k,n) = \sum_{n=0}^{\max\{j,k \}} s(n,j) S(k,n) = \delta_{jk}$

$\delta_{jk}$ 是克罗内克尔δ。

3个分类：

HDU 4767 Bell (中国剩余定理)

分类： ACM_数论and数学类 ACM_HDU 2013-09-28 20:48 84人阅读评论(0) 收藏举报

题意：

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