题目链接：C. Yet Another Number Sequence
题目大意：令 (F_1=1,F_2=2,F_n=F_{n-1}+F_{n-2}(nge 3),A_i(k)=F_i imes i^k) ，求 (sum_{i=1}^{n} A_i(k)) 。
(nleq 10^{17},kleq 40)
Bonus: (kleq 200)
题解：我们令 (f(n,k)=sum_{i=1}^{n} F_i imes i^k) ，接下来划一波式子：
为了方便起见，令 (F_0=1,F_{-1}=0)
[f(n,j)=sum_{i=1}^{n} (F_{i-1}+F_{i-2}) imes i^k
\
=sum_{i=1}^{n} F_{i-1} imes i^k +sum_{i=1}^{n} F_{i-2} imes i^k
\
=sum_{i=1}^{n} F_{i-1}sum_{j=0}^{k} inom{k}{j} (i-1)^j+sum_{i=1}^{n} F_{i-2} sum_{j=0}^{k} inom{k}{j} (i-2)^j 2^{k-j}
\
=sum_{j=0}^{k} inom{k}{j}sum_{i=0}^{n-1} F_{i} imes i^j + 1 +sum_{j=0}^{k}inom{k}{j} 2^{k-j}sum_{i=1}^{n-2}F_i imes i^j +2^k
\
=sum_{j=0}^{k} inom{k}{j} f(n-1,j)+sum_{j=0}^{k} inom{k}{j}2^{k-j}f(n-2,j) +2^k+1]
划到这一步我们就可以直接矩阵快速幂了，矩阵大小为 ((2k+3) imes (2k+3)) ，时间复杂度 (O(k^3log n)) 。
接下来我们考虑 Bonus 中的数据范围怎么做。
我们令 (G) 为斐波那契数列的转移矩阵， (f(n,k)=sum_{i=1}^{n} G^{i} imes i^k)。
那么我们考虑从 (f(n,k)) 拓展到 (f(2n,k)) ：
[f(2n,k)=sum_{i=0}^{2n}G^{i} imes i^k
\
=sum_{i=0}^{n}G^i imes i^k + G^nsum_{i=0}^{n}G^i imes (i+n)^k
\
=f(n,k)+G^nsum_{i=0}^{n}G^isum_{j=0}^{k}inom{k}{j}i^jn^{k-j}
\
=f(n,k)+G^nsum_{j=0}^{k}inom{k}{j}n^{k-j}sum_{i=0}^{n}G^ii^j
\
=f(n,k)+G^nsum_{j=0}^{k}inom{k}{j}n^{k-j}f(n,j)]
直接递归下去处理即可，时间复杂度 (O(k^2log n + 8klog n)) 。
接下来是代码：
(O(k^3log n))
#include <cstdio>
typedef long long ll;
const int Maxk=40;
const int Mod=1000000007;
int C[Maxk+5][Maxk+5];
int pow_2[Maxk+5];
ll n;
int k;
int len;
void init(){
pow_2[0]=1;
C[0][0]=1;
for(int i=1;i<=k;i++){
C[i][0]=C[i][i]=1;
for(int j=1;j<i;j++){
C[i][j]=(C[i-1][j]+C[i-1][j-1])%Mod;
}
}
for(int i=1;i<=k+1;i++){
pow_2[i]=(pow_2[i-1]<<1)%Mod;
}
}
struct Matrix{
int a[2*Maxk+5][2*Maxk+5];
void init(){
for(int i=0;i<len;i++){
a[i][i]=1;
}
}
friend Matrix operator *(Matrix a,Matrix b){
Matrix ans;
for(int i=0;i<len;i++){
for(int j=0;j<len;j++){
ans.a[i][j]=0;
for(int k=0;k<len;k++){
ans.a[i][j]=(ans.a[i][j]+1ll*a.a[i][k]*b.a[k][j])%Mod;
}
}
}
return ans;
}
}trans,ans;
Matrix quick_power(Matrix a,ll b){
Matrix ans;
ans.init();
while(b){
if(b&1){
ans=ans*a;
}
a=a*a;
b>>=1;
}
return ans;
}
int main(){
scanf("%lld%d",&n,&k);
init();
len=(k<<1)+3;
ans.a[0][0]=1;
trans.a[0][0]=1;
for(int i=0;i<=k;i++){
ans.a[0][i+1]=1;
}
for(int i=k+2;i<len;i++){
trans.a[i-(k+1)][i]=1;
}
for(int i=0;i<=k;i++){
for(int j=0;j<=i;j++){
trans.a[j+1][i+1]=C[i][j];
}
for(int j=0;j<=i;j++){
trans.a[j+(k+2)][i+1]=1ll*C[i][j]*pow_2[i-j]%Mod;
}
trans.a[0][i+1]=(pow_2[i]+1)%Mod;
}
if(n==1){
printf("%d
",ans.a[0][k+1]);
return 0;
}
ans=ans*quick_power(trans,n-1);
printf("%d
",ans.a[0][k+1]);
return 0;
}
(O(k^2log n + 8klog n))
#include <cstdio>
typedef long long ll;
const int Maxk=40;
const int Mod=1000000007;
int C[Maxk+5][Maxk+5];
ll n;
int k;
void init(){
C[0][0]=1;
for(int i=1;i<=Maxk;i++){
C[i][0]=C[i][i]=1;
for(int j=1;j<i;j++){
C[i][j]=(C[i-1][j]+C[i-1][j-1])%Mod;
}
}
}
struct Matrix{
int a[2][2];
void init(){
a[0][0]=a[1][1]=1;
a[0][1]=a[1][0]=0;
}
void clear(){
a[0][0]=a[0][1]=a[1][0]=a[1][1]=0;
}
friend Matrix operator *(Matrix a,Matrix b){
Matrix ans;
for(int i=0;i<2;i++){
for(int j=0;j<2;j++){
ans.a[i][j]=0;
for(int k=0;k<2;k++){
ans.a[i][j]=(ans.a[i][j]+1ll*a.a[i][k]*b.a[k][j])%Mod;
}
}
}
return ans;
}
friend Matrix operator +(Matrix a,Matrix b){
Matrix ans;
for(int i=0;i<2;i++){
for(int j=0;j<2;j++){
ans.a[i][j]=(a.a[i][j]+b.a[i][j])%Mod;
}
}
return ans;
}
friend Matrix operator *(Matrix a,int b){
Matrix ans;
for(int i=0;i<2;i++){
for(int j=0;j<2;j++){
ans.a[i][j]=1ll*a.a[i][j]*b%Mod;
}
}
return ans;
}
}fib,f[Maxk+5],t[Maxk+5],tmp;
int pow_m[Maxk+5];
void solve(ll n){
if(n==1){
for(int i=0;i<=k;i++){
f[i]=fib;
}
tmp=fib;
return;
}
ll m=(n>>1);
solve(m);
for(int i=0;i<=k;i++){
t[i]=f[i];
f[i].clear();
}
pow_m[0]=1;
for(int i=1,tmp_m=m%Mod;i<=k;i++){
pow_m[i]=1ll*pow_m[i-1]*tmp_m%Mod;
}
for(int i=0;i<=k;i++){
for(int j=0;j<=i;j++){
f[i]=f[i]+t[j]*(1ll*C[i][j]*pow_m[i-j]%Mod);
}
f[i]=f[i]*tmp;
}
for(int i=0;i<=k;i++){
f[i]=f[i]+t[i];
}
tmp=tmp*tmp;
if(n&1){
tmp=tmp*fib;
pow_m[0]=1;
for(int i=1,tmp_m=n%Mod;i<=k;i++){
pow_m[i]=1ll*pow_m[i-1]*tmp_m%Mod;
}
for(int i=0;i<=k;i++){
f[i]=f[i]+tmp*pow_m[i];
}
}
}
int main(){
fib.a[0][0]=fib.a[0][1]=fib.a[1][0]=1;
fib.a[1][1]=0;
scanf("%lld%d",&n,&k);
init();
solve(n);
printf("%d
",f[k].a[0][0]);
return 0;
}