Yet Another Number Sequence——[矩阵快速幂]

Description

  Everyone knows what the Fibonacci sequence is. This sequence can be defined by the recurrence relation:

F1 = 1, F2 = 2, Fi = Fi - 1 + Fi - 2 (i > 2).

We'll define a new number sequence Ai(k) by the formula:

Ai(k) = Fi × ik (i ≥ 1).

In this problem, your task is to calculate the following sum: A1(k) + A2(k) + ... + An(k). The answer can be very large, so print it modulo 1000000007(109 + 7).

Input

  The first line contains two space-separated integers nk(1 ≤ n ≤ 1017; 1 ≤ k ≤ 40).

Output

  Print a single integer — the sum of the first n elements of the sequence Ai(k) modulo 1000000007(109 + 7).

Sample Input

Input
1 1
Output
1
Input
4 1
Output
34
Input
5 2
Output
316
Input
7 4
Output
73825

解题思路:
  这是一道矩阵快速幂求多项式的应用,思路很清晰:找到各个量之间的递推关系,用矩阵求幂转移状态。问题的关键在于推导A(n,k),sumA(n),
F(n)之间的关系。过程如下:
  1.令 U(n+1,k)=F(n+1)*(n+1)^k;
     V(n+1,k)=F(n)*(n+1)^k;
    Sum(n)=∑[i=1~n]A(i,k);
  2.利用二项式定理将(1+i)^n展开;
  则 U(n+1,k)=∑[i=0~k]C(i,k)*F(n)*(n)^i+∑[i=0~k]C(i,k)*F(n-1)*(n)^i

            =∑[i=0~k]C(i,k)*U(n,i)+∑[i=0~k]C(i,k)*V(n,i);
  同理 V(n+1,k)=
∑[i=0~k]C(i,k)*U(n,i);
  Sum(n)=Sum(n-1)+U(n,k);
  
  3.状态转移用矩阵向量表示:
sum(n-1) sum(n)
U(n,k)   U(n+1,k)

...

 ...
U(n,0)     =>U(n+1,0)
V(n,k)     V(n+1,k)

...

 ...
V(n,0) V(n+1,0)
 
 
 






代码如下:
 1 #include <iostream>
 2 #include <cstdio>
 3 #include <cstring>
 4 #include <time.h>
 5 #include <assert.h>
 6 #define time_ printf("time : %f\n",double(clock())/CLOCKS_PER_SEC)
 7 using namespace std;
 8 #define maxk 40
 9 #define MOD 1000000007
10 #define mod_(x) (x%MOD)
11 
12 typedef long long LL;
13 LL n;
14 int k;
15 
16 struct Matrix{
17     int row,col;
18     int m[2*maxk+5][2*maxk+5];
19     Matrix(int r=83,int c=83){
20         row=r;
21         col=c;
22         memset(m, 0, sizeof m);
23     }
24     void operator=(const Matrix& m2){
25         memcpy(m, m2.m, sizeof m);
26     }
27 };
28 Matrix operator*(const Matrix& a,const Matrix& b){
29     Matrix c(a.row,b.col);
30     for(int i=0;i<c.row;i++)
31         for(int j=0;j<c.col;j++){
32             for(int p=0;p<a.col;p++){
33                 c.m[i][j]=(c.m[i][j]+LL(a.m[i][p])*b.m[p][j])%MOD;
34                 //assert(c.m[i][j]>=0);
35             }
36         }
37     return c;
38 }
39 Matrix C;
40 void set_C(){
41     for(int i=k;i>=0;i--)
42         for(int j=k;j>=i;j--){
43             if(j==k||j==i)
44                 C.m[i][j]=1;
45             else
46                 C.m[i][j]=(C.m[i+1][j]+C.m[i+1][j+1])%MOD;
47             }
48 }
49 void cpy_C(Matrix& a,int pi,int pj){
50     int len=k+1;
51     for(int i=0;i<len;i++)
52         for(int j=0;j<len;j++)
53             a.m[pi+i][pj+j]=C.m[i][j];
54 }
55 void set_M(Matrix& a){
56     a.m[0][0]=1,a.m[0][1]=1;
57     cpy_C(a, 1, 1);
58     cpy_C(a, k+2, 1);
59     cpy_C(a, 1, k+2);
60 }
61 Matrix pow_mod(Matrix& a,LL n){
62     if(n==1) return a;
63     Matrix temp=pow_mod(a, n/2);
64     temp=temp*temp;
65     if(n%2){
66         temp=temp*a;
67     }
68     return temp;
69 }
70 int pow_2(int n){
71     if(n==0) return 1;
72     int temp=pow_2(n/2)%MOD;
73     temp=int((LL(temp)*temp)%MOD);
74     if(n%2)
75         temp=(temp*2)%MOD;
76     return temp;
77 }
78 int main(int argc, const char * argv[]) {
79     while(scanf("%lld%d",&n,&k)==2){
80         if(n==1){
81             printf("%d\n",1);
82             continue;
83         }
84         set_C();
85         Matrix I(2*k+3,2*k+3);
86         set_M(I);
87         Matrix A(2*k+3,1);
88         A.m[0][0]=1;
89         for(int i=1;i<k+2;i++)
90             A.m[i][0]=pow_2(k+1-i+1);
91         for(int i=k+2;i<2*k+3;i++)
92             A.m[i][0]=pow_2(2*k+2-i);
93         Matrix temp=pow_mod(I, n-1);
94         A=temp*A;
95         printf("%d\n",A.m[0][0]);
96         //time_;
97     }
98     return 0;
99 }

 







转载于:https://www.cnblogs.com/Kiraa/p/5317175.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值