hdoj 4767 Bell 【矩阵快速幂 + CRT】

最新推荐文章于 2018-09-07 00:03:38 发布

笑着走完自己的路

最新推荐文章于 2018-09-07 00:03:38 发布

阅读量497

点赞数

分类专栏： gcd && lcm && exgcd && CRT 矩阵

本文链接：https://blog.csdn.net/chenzhenyu123456/article/details/49232183

版权

矩阵同时被 2 个专栏收录

29 篇文章 0 订阅

订阅专栏

gcd && lcm && exgcd && CRT

20 篇文章 1 订阅

订阅专栏

Bell

Time Limit: 6000/3000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 666 Accepted Submission(s): 286

Problem Description

What? MMM is learning Combinatorics!?
Looks like she is playing with the bell sequence now:
bell[n] = number of ways to partition of the set {1, 2, 3, ..., n}

e.g. n = 3:
{1, 2, 3}
{1} {2 3}
{1, 2} {3}
{1, 3} {2}
{1} {2} {3}
so, bell[3] = 5

MMM wonders how to compute the bell sequence efficiently.

Input

T -- number of cases
for each case:
　　n (1 <= n < 2^31)

Output

for each case:
　　bell[n] % 95041567

Sample Input

Sample Output

Bell数：讲解特详细

看会上面这个，就可以KO这道题了。

先找到95041567的小于100的质因子：

#include <cstdio>
int main()
{
    int n = 95041567;
    while(n != 1)
    {
        for(int i = 2; i <= 100; i++)//没有必要判断是否是质数 就可以找到。。。
        {
            if(n % i == 0)
            {
                printf("%d\n", i);
                n /= i;
                break;
            }
        }
    }
    return 0;
}

我构建的矩阵如下：

AC代码：

#include <cstdio>
#include <cstring>
#include <algorithm>
#define LL long long
#define debug printf("1\n");
using namespace std;
int m[6] = {0, 31, 37, 41, 43, 47};
int rm[6];
int a[6];//方程组x % m = a
int B[6][60][60];//保存对m[i]取余后的 Bell数
int F[60];
void getBell(int id, int MOD)//预处理贝尔数
{
    B[id][0][1] = 1;
    for(int i = 1; i <= 50; i++)
    {
        for(int j = 1; j <= i+1; j++)
        {
            if(j == 1)
                B[id][i][j] = B[id][i-1][i];
            else
                B[id][i][j] = (B[id][i-1][j-1] + B[id][i][j-1]) % MOD;
        }
    }
}
struct MATRIX
{
    struct Matrix{
        int a[60][60];
        int r, c;
    };
    Matrix ori, res;
    void init(int id, int p)//建立p*p矩阵
    {
        for(int i = 1; i <= p; i++)
            F[i] = B[id][i][1] % p;
        //debug;
        ori.r = res.r = ori.c = res.c = p;
        memset(ori.a, 0, sizeof(ori.a));
        memset(res.a, 0, sizeof(res.a));
        for(int i = 1; i <= p; i++)
        {
            res.a[i][i] = 1;
            if(i < p)
                ori.a[i][i] = ori.a[i+1][i] = 1;
            else
                ori.a[i][i] = ori.a[1][i] = ori.a[2][i] = 1;
        }
    }
    Matrix multi(Matrix x, Matrix y, int MOD)
    {
        Matrix z;
        memset(z.a, 0, sizeof(z.a));
        z.r = x.r; z.c = y.c;
        for(int i = 1; i <= x.r; i++)
        {
            for(int k = 1; k <= x.c; k++)
            {
                if(x.a[i][k] == 0) continue;
                for(int j = 1; j <= y.c; j++)
                    z.a[i][j] = (z.a[i][j] + (x.a[i][k] * y.a[k][j]) % MOD) % MOD;
            }
        }
        return z;
    }
    void solve(int n, int MOD)
    {
        while(n)
        {
            if(n & 1)
                res = multi(ori, res, MOD);
            ori = multi(ori, ori, MOD);
            n >>= 1;
        }
    }
};
MATRIX MA;
struct CRT
{
    int gcd(int a, int b){
        return b == 0 ? a : gcd(b, a%b);
    }
    void exgcd(int a, int b, int &d, int &x, int &y)
    {
        if(b == 0) {d = a, x = 1, y = 0;}
        else
        {
            exgcd(b, a%b, d, y, x);
            y -= x * (a / b);
        }
    }
    int solve(int l, int r, int *m, int *a)
    {
        int LCM = 1;
        for(int i = l; i <= r; i++)
            LCM = LCM / gcd(LCM, m[i]) * m[i];
        for(int i = l+1; i <= r; i++)
        {
            int A = m[l], B = m[i], d, x, y, c = a[i] - a[l];
            exgcd(A, B, d, x, y);
            if(c % d) return -1;
            int mod = m[i] / d;
            int k = ((x * c / d) % mod + mod) % mod;
            a[l] = m[l] * k + a[l];
            m[l] = m[l] / d * m[i];
        }
        if(a[l] == 0) return LCM;
        return a[l];
    }
};
CRT crt;
int main()
{
    for(int i = 1; i <= 5; i++)
        getBell(i, m[i]);
    int t;
    int n;
    scanf("%d", &t);
    while(t--)
    {
        scanf("%d", &n);
        if(n <= 50)
        {
            for(int i = 1; i <= 5; i++)
                a[i] = B[i][n][1];
            memcpy(rm, m, sizeof(m));
            printf("%d\n", crt.solve(1, 5, rm, a));
            continue;
        }
        for(int i = 1; i <= 5; i++)
        {
            MA.init(i, m[i]);//构建矩阵
            int temp = n / m[i];
            if(n % m[i] == 0)
                temp -= 1;
            MA.solve(temp, m[i]);
            int remain = n % m[i];//找到位置 即第多少列
            if(remain == 0)
                remain = m[i];
            a[i] = 0;
            for(int j = 1; j <= m[i]; j++)//直接求解当前列的值就可以了
                a[i] = (a[i] + (MA.res.a[j][remain] * F[j]) % m[i]) % m[i];
        }
        memcpy(rm, m, sizeof(m));//一开始忘了这个了。。。
        printf("%d\n", crt.solve(1, 5, rm, a));//求解
    }
    return 0;
}