CF - 392 C. Yet Another Number Sequence
这个题看了十几分钟直接看题解了，然后恍然大悟，发现纸笔难于描述于是乎用Tex把初始矩阵以及转移矩阵都敲了出来
(nle 1e17) 这个数量级求前缀和，发现递推关系之后矩阵快速幂是可以求出来的，所以就尝试把(A_i(k)) 的递推式求出来。
[A_{i-1}(k) = F_{i-1} * (i-1) ^ k\
A_{i-2}(k) = F_{i-2} * (i-2) ^ k
]
[egin{aligned}
A_i(k) =& F_i * i ^ k\
=&(F_{i-1} + F_{i-2}) * i ^ k\
=& F_{i-1} * [(i-1) + 1] ^ k + F_{i-2} * [(i-2) + 2] ^ k;\
=& sum_{j=0}^k{C_k^j *F_{i-1} * (i-1) ^ j} + sum_{j=0}^k{C_k^j * F_{i-2} * (i-2) ^ j * 2 ^ {k-j}}\
=& sum_{j=0}^{k}A_{i-1}(j)*C_k^j + sum_{j=0}^kA_{i-2}(j)*C_k^j*2^{k-j}
end{aligned}
]
到这里递推式就求出来了
[A_i(k) = sum_{j=0}^kA_{i-1}(j) *C_k^j + sum_{j=0}^kA_{i-2}(j)*C_k^j*2^{k-j}
]
由于最后求得是(sum_{i=1}^n A_i(k))
所以要把它放到矩阵中，然后矩阵中其他的元素也就理所当然的可以摆出来了
[egin{bmatrix}
sum_{j=1}^i A_j(k)&A_i(0)&A_i(1)&cdots&A_i(k)&A_{i-1}(0)&A_{i-1}(1)&cdots A_{i-1}(k)
end{bmatrix}
]
然后根据递推式以及原始矩阵设计转移矩阵
[egin{bmatrix}1&0&0&cdots&0&0&0&cdots&0\C_k^0&C_0^0&C_1^0&cdots&C_k^0&1&0&cdots&0\C_k^1&0&C_1^1&cdots&C_k^1&0&1&cdots&0\vdots&vdots&vdots&ddots&vdots&vdots&vdots&ddots&vdots\C_k^k&0&0&cdots&C_k^k&0&0&cdots&1\C_k^0*2^k&C_0^0*2^0&C_1^0*2^1&cdots&C_k^0*2^k&0&0&cdots&0\C_k^1*2^{k-1}&0&C_1^1*2^0&cdots&C_k^1*2^{k-1}&0&0&cdots&0\vdots&vdots&vdots&ddots&vdots&vdots&vdots&ddots&vdots\C_k^k*2^0&0&0&cdots&C_k^k*2^{0}&0&0&cdots&0end{bmatrix}
]
可能直接放出来不好理解，那么对这个矩阵划分一下几个区
左边一列是用来求和的，中间这两部分是用来计算(A_{i+1}^j (jin[0,k])) 的，右侧是用来转移(A_{i}^j(jin [0,k])) 的。
最后直接矩阵快速幂就好了，整个过程细节比较多。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 101;
const int mod = 1e9 + 7;
struct matrix{
ll mat[N][N];
int r,c;
matrix(){}
matrix(int rr,int cc){r = rr,c= cc;}
void clear(){
memset(mat,0,sizeof mat);
}
};
void MOD(ll &x){x=(x%mod+mod)%mod;}
matrix operator * (const matrix&a,const matrix&b){
matrix c(a.r,b.c);
c.clear();
for(int k=0;k<a.c;k++){
for(int i=0;i<c.r;i++){
for(int j=0;j<c.c;j++){
MOD(c.mat[i][j] += a.mat[i][k] * b.mat[k][j] % mod);
}
}
}
return c;
}
ll C[55][55],p[55];
void prework(){
p[0] = 1;
for(int i=1;i<=50;i++)p[i] = p[i-1] * 2 % mod;
for(int i=0;i<=50;i++)C[i][0] = 1;
for(int i=1;i<=50;i++){
for(int j=1;j<=i;j++)
C[i][j] = (C[i-1][j] + C[i-1][j-1]) % mod;
}
}
ll solve(ll n,ll k){
if(n == 1)return 1;
if(n == 2)return (p[k+1] + 1) % mod;
ll sum = 0;
int kk = k * 2 + 3;
matrix x(1,kk),y(kk,kk);
x.mat[0][0] = (p[k+1] + 1) % mod;
for(int i=0;i<=k;i++){
x.mat[0][i+1] = p[i+1];
x.mat[0][i+k+2] = 1;
}
y.mat[0][0] = 1;
for(int i=0;i<=k;i++){
y.mat[i+1][0] = C[k][i];
y.mat[i+k+2][0] = C[k][i] * p[k-i] % mod;
y.mat[i+1][i+k+2] = 1;
}
for(int i=0;i<=k;i++){
for(int j=0;j<=i;j++){
y.mat[j+1][i+1] = C[i][j];
y.mat[j+k+2][i+1] = C[i][j] * p[i-j] % mod;
}
}
n -= 2;
for(;n;n >>= 1){
if(n & 1)x = x * y;
y = y * y;
}
return x.mat[0][0];
}
int main(){
ll n,k;
prework();
scanf("%lld%lld",&n,&k);
printf("%lld
",solve(n,k));
return 0;
}