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

时间：2021-05-26 20:01:38

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 n, k(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 }

 







标签：

相关文章