HDU 2294 Pendant（DP+矩阵）

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=2294

题意：n种珠子，串成一个长度不超过L的链。要求每种珠子至少使用一个。求有多少种方案？

思路：设f[i][j]表示长度为i的链子使用了j种珠子，则f[i][j]=f[i-1][j]*j+f[i-1][j-1]*(n-j+1)。直接DP的复杂度太大。根据该转移方程定义输入矩阵f[1]为1*(n+1)的矩阵,该矩阵只有第二个元素为n其余为0表示f[1][1]=n，转移矩阵a为a[i-1][i]=n-i+1,a[i][i]=i。则f[i]=f[i-1]*a即f[i]=f[1]*a^(i-1)。则答案为f[1]+f[2]+……+f[L]=f[1]*(a^0+……+a^(L-1))的(1,n)号元素和n的乘积。

int n;

class Matrix
{
public:
    int a[N][N];

    void init(int x)
    {
        int i;
        if(!x) clr(a,0);
        else
        {
            clr(a,0);
            FOR0(i,N) a[i][i]=1;
        }
    }


    Matrix operator+(Matrix b)
    {
        int i,j;
        Matrix c;
        FOR0(i,n) FOR0(j,n) c.a[i][j]=((i64)a[i][j]+b.a[i][j])%mod;
        return c;
    }

    Matrix operator+(int x)
    {
        int i;
        Matrix c=*this;
        FOR0(i,n) (c.a[i][i]+=x)%=mod;
        return c;
    }


    Matrix operator*(Matrix b)
    {
        int i,j,k,t;
        Matrix p;
        p.init(0);
        FOR0(i,n) FOR0(j,n) FOR0(k,n)
        {
            p.a[i][j]=(p.a[i][j]+(i64)a[i][k]*b.a[k][j]%mod)%mod;
        }
        return p;
    }

    Matrix power(int t)
    {
        Matrix ans,p=*this;
        ans.init(1);
        while(t)
        {
            if(t&1) ans=ans*p;
            p=p*p;
            t>>=1;
        }
        return ans;
    }
};

Matrix a;

int C;
int L,m;

struct node
{
    Matrix a,b;
};

stack<node> st;

Matrix cal(Matrix a,int n)
{
    Matrix b;
    b.init(1);
    if(n==0) return b;
    if(n==1) return a+b;
    if(n==2) return a+a*a+b;
    if(n==3) return a+a*a+a*a*a+b;
    while(!st.empty()) st.pop();
    node p;
    while(n)
    {
        if(n%2==0)
        {
            p.a=a.power(n);
            p.b=a.power(n/2)+1;
            n=n/2-1;
        }
        else
        {
            p.a.init(0);
            p.b=a.power(n/2+1)+1;
            n/=2;
        }
        st.push(p);
    }
    while(!st.empty())
    {
        p=st.top();
        st.pop();
        b=p.a+p.b*b;
    }
    return b;
}

int main()
{
    RD(C);
    while(C--)
    {
        RD(L,m);
        if(L<m)
        {
            puts("0");
            continue;
        }
        n=m+1;
        a.init(0);
        int i;
        FOR1(i,m)
        {
            a.a[i-1][i]=m-i+1;
            a.a[i][i]=i;
        }
        a=cal(a,L-1);
        printf("%I64d\n",(i64)a.a[1][m]*m%mod);
    }
    return 0;
}