CodeForces 392 C.Yet Another Number Sequence（矩阵快速幂）

Description

定义$F_{1}=1,F_{2}=2,F_{i}=F_{i-1}+F_{i-2}(i\ge 3)$

定义$A_{i}(k)=F_{i}\cdot i^k(i\ge 1)$

求$A_{1}(k)+A_{2}(k)+...+A_{n}(k)$，结果模$10^9+7$

Input

两个正整数$n$和$k$ $(1\leq n \leq 10^17,1\le k\le 40)$

Output

输出$A_{1}(k)+A_{2}(k)+...+A_{n}(k)$，结果模$10^9+7$

Sample Input

1 1

Sample Output

1

Solution

$A_{i}(k)=F_{i}\cdot i^k=(F_{i-1}+F_{i-2})\cdot i^k=F_{i-1}\cdot (i-1+1)^k+F_{i-2}\cdot (i-2+2)^k=\sum\limits_{j=0}^iC_{i}^{j}A_{i-1}(j)+\sum\limits_{j=0}^iC_{i}^{j}2^{i-j}A_{i-2}(j)$

用矩阵快速幂加速该转移即可，时间复杂度$O((2k+3)^3log_{2}n)$

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
#define mod 1000000007
#define maxn 88
struct Mat
{
    ll mat[maxn][maxn];//矩阵 
    int row,col;//矩阵行列数 
};
Mat mod_mul(Mat a,Mat b,ll p)//矩阵乘法 
{
    Mat ans;
    ans.row=a.row;
    ans.col=b.col;
    memset(ans.mat,0,sizeof(ans.mat));
    for(int i=0;i<ans.row;i++)      
        for(int k=0;k<a.col;k++)
            if(a.mat[i][k])
                for(int j=0;j<ans.col;j++)
                {
                    ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                    ans.mat[i][j]%=p;
                }
    return ans;
}
Mat mod_pow(Mat a,ll k,ll p)//矩阵快速幂 
{
    Mat ans;
    ans.row=a.row;
    ans.col=a.col;
    for(int i=0;i<a.row;i++)
        for(int j=0;j<a.col;j++)
            ans.mat[i][j]=(i==j);
    while(k)
    {
        if(k&1)ans=mod_mul(ans,a,p);
        a=mod_mul(a,a,p);
        k>>=1;
    }
    return ans;
}
int k;
ll n,C[maxn][maxn];
void init()
{
    memset(C,0,sizeof(C));
    C[0][0]=1;
    for(int i=1;i<=40;i++)
    {
        C[i][0]=C[i][i]=1;
        for(int j=1;j<i;j++)C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
    }
}
Mat A,B;
int main()
{
    init();
    while(~scanf("%I64d%d",&n,&k))
    {
        if(n==1)printf("1\n");
        else
        {
            A.row=A.col=2*k+3;
            memset(A.mat,0,sizeof(A.mat));
            A.mat[0][0]=A.mat[0][k+1]=1;
            for(int i=0;i<=k;i++)
                for(int j=0;j<=i;j++)
                {
                    A.mat[i+1][j+1]=C[i][j];
                    A.mat[i+1][j+k+2]=(1ll<<(i-j))%mod*C[i][j]%mod;
                }
            for(int i=k+2;i<=2*k+2;i++)A.mat[i][i-(k+1)]=1;
            B.row=2*k+3,B.col=1;
            B.mat[0][0]=1;
            for(int i=1;i<=k+1;i++)B.mat[i][0]=(1ll<<i)%mod;
            for(int i=k+2;i<=2*k+2;i++)B.mat[i][0]=1;
            A=mod_pow(A,n-2,mod);
            A=mod_mul(A,B,mod);
            printf("%I64d\n",(A.mat[0][0]+A.mat[k+1][0])%mod); 
        }
    }
    return 0;
}