hdu2294---Pendant (矩阵)

Problem Description
On Saint Valentine’s Day, Alex imagined to present a special pendant to his girl friend made by K kind of pearls. The pendant is actually a string of pearls, and its length is defined as the number of pearls in it. As is known to all, Alex is very rich, and he has N pearls of each kind. Pendant can be told apart according to permutation of its pearls. Now he wants to know how many kind of pendant can he made, with length between 1 and N. Of course, to show his wealth, every kind of pendant must be made of K pearls.
Output the answer taken modulo 1234567891.

Input
The input consists of multiple test cases. The first line contains an integer T indicating the number of test cases. Each case is on one line, consisting of two integers N and K, separated by one space.
Technical Specification

1 ≤ T ≤ 10
1 ≤ N ≤ 1,000,000,000
1 ≤ K ≤ 30

Output
Output the answer on one line for each test case.

Sample Input

2 2 1 3 2

Sample Output

2 8

Source
The 4th Baidu Cup final

Recommend
lcy | We have carefully selected several similar problems for you: 1588 3117 2971 2256 1757

dp[i][j]表示长度为i,有j种珍珠的吊坠的数目
dp[i][j] = (k - j + 1) * dp[i - 1][j - 1] + j * dp[i - 1][j];
n太大,考虑用矩阵

由于每一个状态只和长度比他小1的状态有关,可以降到一维
设矩阵fn 表示 长度为n的吊坠的数目(0种珍珠,1种珍珠, ….. k种珍珠)
fn:
dp[0] dp[1] ….. dp[k]
0 0 ….. 0
… .. … 0
0 0 …. 0
构造这样一个(k+1)*(k+1)的矩阵

fn = fn-1 Ａ　＝　f1 A^(n-1)

sum = f1 + f2 + … + fn = (E + A + … + A ^(n - 1)) * f1
根据递推关系可得
Ａ：
　０ k …… …. 0
　 0 1 k - 1 … 0
　 0 0 2 k - 2 0
　 …………………..
　 0 0 ….. …… k

显然f1就等于Ａ

则sum(n) = (A + A^2 + …. + A^n) = (E + A^(n / 2)) * (A + A ^ 2 + … + A ^ (n / 2)) + (n & 1) ? A^n : 0
斜体加粗部分就是sum(n / 2)
二分就行

/*************************************************************************
    > File Name: hdu2294.cpp
    > Author: ALex
    > Mail: zchao1995@gmail.com 
    > Created Time: 2015年03月11日 星期三 15时52分41秒
 ************************************************************************/

#include <map>
#include <set>
#include <queue>
#include <stack>
#include <vector>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

const double pi = acos(-1);
const int inf = 0x3f3f3f3f;
const double eps = 1e-15;
typedef long long LL;
typedef pair <int, int> PLL;

const int mod = 1234567891;

int k;

struct MARTIX
{
    LL mat[32][32];
    int n;

    void init(int _n)
    {
        n =_n;
        for (int i = 0; i < n; ++i)
        {
            for (int j = 0; j < n; ++j)
            {
                mat[i][j] = 0;
            }
        }
    }

    MARTIX operator + (const MARTIX &b) const
    {
        MARTIX res;
        res.n = b.n;
        for (int i = 0; i < b.n; ++i)
        {
            for (int j = 0; j < b.n; ++j)
            {
                res.mat[i][j] = mat[i][j] + b.mat[i][j];
                res.mat[i][j] %= mod;
            }
        }
        return res;
    }

    MARTIX operator * (const MARTIX &b)const
    {
        MARTIX res;
        res.n = b.n;
        for (int i = 0; i < b.n; ++i)
        {
            for (int j = 0; j < b.n; ++j)
            {
                res.mat[i][j] = 0;
                for (int l = 0; l < b.n; ++l)
                {
                    res.mat[i][j] += mat[i][l] * b.mat[l][j];
                    res.mat[i][j] %= mod;
                }
            }
        }
        return res;
    }

    MARTIX & operator = (const MARTIX & b)
    {
        this -> n = b.n;
        for (int i = 0; i < b.n; ++i)
        {
            for (int j = 0; j < b.n; ++j)
            {
                this -> mat[i][j] = b.mat[i][j];
            }
        }
        return *this;
    }
}A, E;

MARTIX fastpow(MARTIX ret, int n)
{
    MARTIX ans;
    ans.init(k + 1);
    for (int i = 0; i < ans.n; ++i)
    {
        for (int j = 0; j < ans.n; ++j)
        {
            ans.mat[i][j] = (i == j);
        }
    }
    while(n)
    {
        if (n & 1)
        {
            ans = ans * ret;
        }
        n >>= 1;
        ret = ret * ret;
    } 
    return ans;
}

void Debug(MARTIX A)
{
    for (int i = 0; i < A.n; ++i)
    {
        for (int j = 0; j < A.n; ++j)
        {
            printf("%lld ", A.mat[i][j]);
        }
        printf("\n");
    }
}

MARTIX binsearch(int cnt)
{
    if (cnt == 1)
    {
        return A;
    }
    MARTIX tmp = binsearch(cnt >> 1);
    MARTIX c = fastpow(A, cnt >> 1);
    c = c + E;
    tmp = c * tmp;
    if (cnt & 1)
    {
        MARTIX d = fastpow(A, cnt);
        tmp = tmp + d;
    }
    return tmp;
}

int main()
{
    int t;
    scanf("%d", &t);
    while (t--)
    {
        int n;
        scanf("%d%d", &n, &k);
        A.init(k + 1);
        E.init(k + 1);
        for (int i = 0; i < E.n; ++i)
        {
            for (int j = 0; j < E.n; ++j)
            {
                E.mat[i][j] = (i == j);
            }
        }
        A.mat[0][0] = 0;
        A.mat[0][1] = k;
        for (int i = 1; i < A.n; ++i)
        {
            A.mat[i][i] = i;
            if (i + 1 < A.n)
            {
                A.mat[i][i + 1] = k - i;
            }
        }
        MARTIX ans;
        ans.init(k + 1);
        ans = binsearch(n);
//      Debug(ans);
        printf("%lld\n", ans.mat[0][k]);
    }
    return 0;
}