三月校赛1006 wuli通通和Fibonacci （a[n]=f[n]*(n^m)的前k项和）

题目链接：
校赛 1006 wuli通通和Fibonacci
题意：
$计算a[n]=f[n]*(n^m)的前k项和mod 23333。f[n]是Fibonacci数列，f[1]=f[2]=1,k<=1e9,m<=40$
分析：
$这样明显的复杂递推公式题而且要求mod23333一定和矩阵快速幂有关，所以问题就在构造矩阵上了。 再读题目发现m<=40，m比较小，所以构造的矩阵可能和m有关。而m在递推式中又是和n挂在一起的,那就要联想到如何 构造矩阵来由(n-1)^m得到n^m，由组合数的知识我们知道$

$n^m=((n-1)+1)^m=C[m][m]*(n-1)^m+C[m][m-1]*(n-1)^(m-1)+C[m][m-2]*(n-1)^(m-2)+...+C[m][1]*(n-1)^(1)+C[m][0]*(n-1)^(0);其中C[i][j]是组合数$

那么就可以构造出如下矩阵了：

| C[m][0]    0         0        0    ...    0   |   | (n-1)^(0)  |   | n^0 |
| C[m][0]  C[m][1]     0        0    ...    0   |   | (n-1)^(1)  |   | n^1 |
| C[m][0]  C[m][1]  C[m][2]     0    ...    0   | * | (n-1)^(2)  | = | n^2 |
| C[m][0]  C[m][1]  C[m][2]  C[m][3] ...    0   |   | (n-1)^(3)  |   | n^3 |
|    .   .   .      .   .  .    .    .    . .   |   | .  .   .  .|   | . . |
| C[m][0]  C[m][1]  C[m][2]  C[m][3] .. C[m][m] |   | (n-1)^(m)  |   | n^m |

再想办法将Fibonacci数列添加进去，一开始我是想在构造一个右乘矩阵，如下：

| C[m][0]     0          0         0    ...       0   |    | f[n-1](n-1)^(0)  f[n-1](n-2)^(0) |
| C[m][0]   C[m][1]      0         0    ...       0   |    | f[n-1](n-1)^(1)  f[n-1](n-2)^(1) |
| C[m][0]   C[m][1]   C[m][2]      0    ...       0   |  * | f[n-1](n-1)^(2)  f[n-1](n-2)^(2) |
| C[m][0]   C[m][1]   C[m][2]   C[m][3] ...       0   |    | f[n-1](n-1)^(3)  f[n-1](n-2)^(3) |
|    .   .    .    .   .   .   .    .    .    .   .   |    | .  .   .  .       .    .   .  .  |
| C[m][0]   C[m][1]   C[m][2]   C[m][3] ...   C[m][m] |    | f[n-1](n-1)^(m)  f[n-1](n-2)^(m) |
|         |  
      |  1  1   |  
   *  |         | 
      |         |  = 
      |  1   0  |   
      |         |    
| f[n](n-1)^(0)  f[n-1](n-1)^(0) |
| f[n](n-1)^(1)  f[n-1](n-1)^(1) |
| f[n](n-1)^(2)  f[n-1](n-1)^(2) |
| f[n](n-1)^(3)  f[n-1](n-1)^(3) |
| .  .  .  .  .  .  .  .  .  . . |
| f[n](n-1)^(m)  f[n-1](n-1)^(m) |

但是这样就发现有一个问题那就是我无法再通过添加一维的方法将a[n]的和这一项添加进中间矩阵。
后来又想了到了一个方法，可以将上面的中间(m+1)2阶矩阵拉伸成(2(m+1))*阶矩阵，也就是这样:

| f[n-1](n-1)^(0) |
| f[n-1](n-2)^(0) |
| f[n-1](n-1)^(1) |
| f[n-1](n-2)^(1) |
| f[n-1](n-1)^(2) |
| f[n-1](n-2)^(2) |
| f[n-1](n-1)^(3) | 
| f[n-1](n-2)^(3) |
|  .  .   .  .    |  
|  .    .   .  .  |
| f[n-1](n-1)^(m) |
| f[n-1](n-2)^(m) |

然后再将左面的矩阵也跟着变一下：

| C[m][0]   C[m][0]      0         0    ...       0       0    |
| C[m][0]      0         0         0    ...       0       0    |
| C[m][0]   C[m][0]   C[m][1]   C[m][1] ...       0       0    |
| C[m][0]   C[m][1]      0         0    ...       0       0    |
|    .   .    .    .   .   .   .    .    .    .   .       0    |   
| C[m][0]   C[m][0]   C[m][1]   C[m][1] ...   C[m][m]  C[m][m] |
| C[m][0]      0      C[m][1]      0    ...   C[m][m]     0    |

最后再在最后一维各添加一行（注意维数都变了）：

|   S[n-2]  | 
和
| C[m][0]   C[m][0]   C[m][1]   C[m][1] ...   C[m][m]  C[m][m]  1 | 

就行了。【好累。。。o(╯□╰)o】

#include <cstdio>
#include <cstring>
using namespace std;

const long long mod=23333;
const int maxn=45;
int T,m,k;
long long c[maxn][maxn];

struct Matrix{
    int row,col;
    long long data[maxn*2][maxn*2];
};

inline void calc()
{
    for(int i=0;i<=42;i++){
        for(int j=0;j<=i;j++){
            if(j==0) c[i][j]=1;
            else c[i][j]=c[i][j-1]*(i-j+1)/j;
        }
    }
}

inline Matrix mul(Matrix a,Matrix b)
{
    Matrix ans;
    ans.row=a.row,ans.col=b.col;
    memset(ans.data,0,sizeof(ans.data));
    for(int i=1;i<=ans.row;i++){
        for(int j=1;j<=ans.col;j++){
            for(int k=1;k<=a.col;k++){
                ans.data[i][j]+=a.data[i][k]*b.data[k][j]%mod;
                ans.data[i][j]%=mod;
            }
        }
    }
    return ans;
}

inline Matrix quick_power(Matrix a,int time)
{
    Matrix ans,tmp=a;
    ans.row=ans.col=a.row;
    memset(ans.data,0,sizeof(ans.data));
    for(int i=1;i<=ans.row;i++) ans.data[i][i]=1;
    while(time){
        if(time&1) ans=mul(ans,tmp);
        tmp=mul(tmp,tmp);
        time>>=1;
    }
    return ans;
}

int main()
{
    //freopen("1006in.txt","r",stdin);
    //freopen("1006out.txt","w",stdout);
    calc();
    scanf("%d",&T);
    while(T--){
        scanf("%d%d",&k,&m);
        if(k==1){
            printf("1\n");
            continue;
        }
        Matrix ans,tmp;
        ans.row=ans.col=tmp.row=2*m+3,tmp.col=1;
        memset(ans.data,0,sizeof(ans.data));
        for(int i=1;i<=m+1;i++){
            for(int j=1;j<=i;j++){
                ans.data[2*i-1][2*j-1]=ans.data[2*i-1][2*j]=ans.data[2*i][2*j-1]=c[i-1][j-1];
                ans.data[2*i][2*j]=0;
            }
        }
        for(int i=1;i<=m+1;i++) ans.data[2*m+3][2*i-1]=ans.data[2*m+3][2*i]=c[m][i-1];
        ans.data[2*m+3][2*m+3]=1;
        for(int i=1;i<=m*2+2;i+=2) {
            tmp.data[i][1]=1;
            tmp.data[i+1][1]=0;
        }
        tmp.data[2*m+3][1]=1;
        ans=quick_power(ans,k-1);
        tmp=mul(ans,tmp);
        printf("%I64d\n",tmp.data[2*m+3][1]%mod);
    }
    return 0;
}