HDU2294_DP+矩阵加速(实在妙题)

/*
*State: 2294    609MS    944K    4236 B    C++
*题目大意：
*        有k种珍珠，每种有n个，然后要求组合成长度为1~n的项链的总数。
*        （项链的长度为珍珠的个数），并要求项链中至少含有k种珍珠。
*解题思路：
*        复杂的组合题，至今不理解其递推式，但是知道dp的状态转移表达式为
*        f[i][j] = f[i - 1][j - 1] * (k - j + 1) + f[i - 1][j] * j;
*        f[i][j]下标的意思是有i颗珍珠时，含有j种能组合的个数。然后有了
*        这个状态转移，可以得出结果为：f[1][k] + f[2][k] + …… + f[n][k],
*        然后发现这个状态转移够囧，n<1 000 000 000,那么一般根据二维dp要
*        两维来实现，结果复杂度就变成了O(nk),超级计算机也难以接受。
*        考虑用别的方式来实现dp,猥琐地搜了解题报告后发现，原来可以用矩阵
*        加速，有一步跳跃比较高，就是一个矩阵里面把其中一维全部包括了。
*        F[i] = F[i - 1] * A;其中A为转移矩阵。然后F[i]这个矩阵要设定为：
*        f[i-1][0] f[i-1][1] f[i-1][2] …… f[i-1][k]
*        0         0         0            0
*       0         0         0            0
*        …………
*        0         0         0            0
*        这个矩阵有k+1维，为什么是k+1维，自己想。
*        转移矩阵A为
*        0 k 0   0 …… 0 0 0 
*        0 1 k-1 0 …… 0 0 0 
*        0 0 2   0 …… 0 0 0
*        ……
*        0 0 0   0 …… 0 0 k
*        为什么转移矩阵是如此的？把上面的F[i]跟转移矩阵相乘可发现跟状态转移
*        方程一致。
*        
*        有了矩阵的表达式，那么要求结果必须先求F[1] + F[2] + F[3] + …… + F[k]
*        转化为F[1] + F[1]*A + F[1]*A^2 + …… + F[1]*A^(k-1).提取F[1],之后就是
*        二分求sum(A^n)了。
*题目困惑：
*        题意看了好久，感觉题目表示得不好，表示一条项链必须含有k钟珍珠的时候。
*        to show his wealth, every kind of pendant must be made of K pearls.
*        好像是在说需要k个珍珠一样，不是应该加个 K kind of pearls么？
*        最后理清思路了，还是wa了3次，因为最后的模溢出了，还是自己走投无路给自己
*        乱出数据得出来的，这个在实在没办法的时候，不失为一种方法。
*/

#include <stdio.h>
#include <string.h>
#include <iostream>

#define MAX_DIMENSION 35 
using namespace std;

typedef __int64 MATRIX_TYPE;
typedef __int64 MAX_INT_TYPE; 
typedef MATRIX_TYPE Matrix[MAX_DIMENSION][MAX_DIMENSION];
int ndim=MAX_DIMENSION;
int mod=1234567891;//mod

void m_zero(Matrix  x)
{
    memset(x, 0, sizeof(MATRIX_TYPE)*MAX_DIMENSION*MAX_DIMENSION);
}

void m_one(Matrix  x)
{
    int i;
    m_zero(x);
    for(i=0;i<ndim;++i)x[i][i]=1;
}

void m_copy(Matrix  src,Matrix  dest)
{
    memcpy(dest,src, sizeof(MATRIX_TYPE)*MAX_DIMENSION*MAX_DIMENSION);
}

//z=x+y;
void m_add(Matrix  x,Matrix  y,Matrix  z)
{
    int i,j;
    for(i=0;i<ndim;i++)
        for(j=0;j<ndim;j++)
            if(mod<=1)z[i][j]=x[i][j]+y[i][j];
            else z[i][j]=(x[i][j]+(MAX_INT_TYPE)y[i][j])%mod;//module
}


//c=a*b
void m_multiple(Matrix  a,Matrix b,Matrix c)
{
    int i,j,k;
    MAX_INT_TYPE t;

    for(i=0;i<ndim;i++)
        for(j=0;j<ndim;j++)
        {
            t=0;
            if(mod<=1)
                for(k=0;k<ndim;k++) t+=a[i][k]*b[k][j];//module
            else
                for(k=0;k<ndim;k++){
                    t+=(a[i][k]*(MAX_INT_TYPE)b[k][j])%mod;
                    t%=mod;
                }//module
            c[i][j]=t;
        }
}

void m_pow_with_saved(Matrix  x_p[],unsigned int n, Matrix y)
{
    Matrix temp;
    m_one(y);
    for(int k=1;n;++k,n>>=1)
    {
        if ((n & 1) != 0)
        {
            m_multiple(x_p[k],y,temp);
            m_copy(temp,y);
        }
    }
}

void m_pow_sum1(Matrix  x_p[],unsigned int n, Matrix y)
{
    m_zero(y);
    if(n==0)return;
    if(n==1) m_copy(x_p[1],y);
    else
    {
        Matrix temp;
        //calculate x^1+...+^(n/2)
        m_pow_sum1(x_p,n>>1,temp);
        //add to y
        m_add(temp,y,y);
        //calculate x^(1+n/2)+...+x^n
        Matrix temp2;
        m_pow_with_saved(x_p,n>>1,temp2);
        Matrix temp3;
        m_multiple(temp,temp2,temp3);
        //add to y
        m_add(temp3,y,y);
        if(n&1)
        {
            //calculate x^(n-1)
            m_multiple(temp2,temp2,temp3);
            //calculate x^n
            m_multiple(temp3,x_p[1],temp2);
            //add x^n
            m_add(temp2,y,y);
        }
    }

}

void m_pow_sum(Matrix x, unsigned int n, Matrix y)
{
    //calculate x^1 x^2 x^4 ... x^logn
    unsigned int i=0,logn,n2=n;
    for(logn=0,n2=n;n2;logn++,n2 >>= 1);
    Matrix *pow_arry=new Matrix[logn==0?2:(logn+1)];
    m_one(pow_arry[0]);
    m_copy(x,pow_arry[1]);
    for(i=1;i<logn;i++)
    {
        m_multiple(pow_arry[i],pow_arry[i],pow_arry[i+1]);
    }

    m_pow_sum1(pow_arry,n,y);
    delete []pow_arry;
}

void view_mat(Matrix a, int n)
{
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
            cout << a[i][j] << " ";
        cout << endl;
    }
}

int main(void)
{
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif

    int cas;
    scanf("%d", &cas);
    while(cas--)
    {
        int n, k;
        scanf("%d %d", &n, &k);
        Matrix a;
        ndim = k + 1;
        m_zero(a);
        for(int i = 0; i <= k; i++)
        {
            if(i == 0)
                continue;
            a[i][i] = i;
            a[i - 1][i] = k - i + 1;
        }    
        Matrix asum, res, addsum, f1, ans;
        m_one(res);
        m_pow_sum(a, n - 1, asum);
        m_add(res, asum, addsum);
        m_zero(f1);
        f1[0][1] = k;
        m_multiple(f1, addsum, ans);
        printf("%d\n", ((ans[0][k] % mod) + (__int64)mod) % mod);
    }
    return 0;
}

这道题目把dp的弱性跟矩阵的优美强悍地给表现出来了。