【例题】【矩阵乘法】NKOJ 3124 珍珠吊坠

NKOJ 3124 珍珠吊坠
时间限制 : 10000 MS 空间限制 : 65536 KB
问题描述
有k种珍珠，每种珍珠有N颗，问能组合出多少种含有k种珍珠且长度<=N的吊坠?(注意：吊坠为链状，非环状)

输入格式
一行，两个整数N和K

输出格式
一行，一个整数，表示所求方案总数,结果可能很大，mod 1234567891再输出

样例输入
样例输入1：
2 1
样例输入2：
3 2

样例输出
样例输出1：
2
样例输出2：
8

提示
1 ≤ N ≤ 1,000,000,000
1 ≤ K ≤ 30
样例2说明：
设两种珍珠编号为1和2
1 2
2 1
1 1 2
1 2 1
2 1 1
2 2 1
2 1 2
1 2 2

来源 HDU2294

思路：
设f[i][j]表示前i颗珠子共含j种的方案数。
可得到递推式f[i][j]=f[i-1][j-i]*(k-(j-1))+f[i-1][j]*j
1<=i<=n,1<=j<=k
而所求结果ans[i]=f[1][k]+f[2][k]+…+f[i][k]=ans[i-1]+f[i][k];
再考虑用矩阵乘法加速递推：
考虑构造F、A使得F*A=F’且F’中包含ans[i]；
要使F’中包含ans[i]，则F中需包含f[i][k]、ans[i-1]；
要包含f[i][k]，根据递推式，需包含f[i-1][k]、f[i-1][k-1]
…..
所以构造的F应为f[i][1],f[i][2],…,f[i][k],ans[i-1];
再考虑构造A：
由：
f[i+1][1]=k=f[i][1];
f[i+1][2]=f[i][1]*(k-1)+f[i][2]*2;
….
f[i+1][k]=f[i][k-1]*1+f[i][k]*k;
ans[i]=ans[i-1]+f[i][k];
可得A：

代码：

#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
#define ll long long
const int p=1234567891;
const int need=31;
typedef int int_[need][need];

int_ a,res,c;
int n,k,m;

void matrix_multi(int_ a,int_ b)
{
    memset(c,0,sizeof(c));
    for(int i=1;i<=m;i++)
     for(int j=1;j<=m;j++)
      for(int t=1;t<=m;t++)
       c[i][j]=((ll)a[i][t]*b[t][j]%p+c[i][j])%p;
    memcpy(a,c,sizeof(c));
}

void matrix_power(int b)
{
    for(int i=1;i<=m;i++) res[i][i]=1;
    while(b)
    {
        if(b&1) matrix_multi(res,a);
        b>>=1;
        matrix_multi(a,a);
    }
}

int main()
{
    scanf("%d%d",&n,&k);
    for(int i=1;i<=k;i++) 
    {
        a[i][i]=i;
        a[k-i][k-i+1]=i;
    }
    a[k][k+1]=1,a[k+1][k+1]=1;
    m=k+1;
    matrix_power(n); //ans[]下标从0开始，所以乘n次 
    printf("%d",((ll)k*res[1][m])%p);
}