嵌套递推——矩阵快速幂

题目描述：

令f(n)=2*f(n-1)+3*f(n-2)+n,f(1)=1,f(2)=2
令g(n)=g(n-1)+f(n)+n*n,g(1)=2
告诉你n，输出g(n)的结果，结果对1e9+7取模

输入描述:

多组输入，每行一个整数n（1<=n<=1e9)，如果输入为0，停止程序。

输出描述:

在一行中输出对应g(n)的值，结果对1e9+7取模。

示例1

输入

输出

说明

多组输入，输入为0时，终止程序

备注:

项数极大，朴素算法无法在规定时间内得出结果


解题思路：
根据公式得出递推公式（含有6项），通过前一项和后一项的关系得出关系矩阵。然后求解。

 1 //嵌套公式 
 2 #include<bits/stdc++.h>
 3 const int MOD = 1e9+7; 
 4 #define mod(x) ((x)%MOD)
 5 const int maxn = 20;
 6 typedef long long ll;
 7 using namespace std;
 8 struct martix
 9 {
10     ll a[maxn][maxn];
11     friend martix operator * (martix a, martix b)
12     {
13         martix res;
14         for(int i = 1; i <= 6;i++)
15         {
16             for(int j = 1;j <= 6;j++)
17             {
18                 ll ans = 0;
19                 for(int k = 1;k <= 6;k++)
20                     ans += mod((ll)a.a[i][k] * b.a[k][j]);
21                 res.a[i][j] = mod(ans);
22             }
23         }
24         return res;
25     }
26     friend martix operator ^ (martix a ,ll n)
27     {
28         martix unit;
29         memset(unit.a,0,sizeof(unit.a));
30         for(int i = 1;i <= 6;i++) unit.a[i][i] = 1;
31         while(n)
32         {
33             if(n&1) unit = unit * a;
34             a = a * a;
35             n >>= 1;
36         }
37         return unit;
38     }
39 }; 
40 int main()
41 {
42     martix res,ans;
43     
44     ll n,m;
45     while(cin>>n)
46     {
47         if(n == 0) return 0;
48         else if(n == 1) cout<<2<<endl;
49         else if(n == 2) cout<<8<<endl;
50         else 
51         {
52             memset(res.a,0,sizeof(res.a));
53             memset(ans.a,0,sizeof(ans.a));
54             
55             res.a[1][1] = 8;res.a[1][2] = 2;res.a[1][3] = 1;res.a[1][4] = 4;res.a[1][5] = 2;res.a[1][6] = 1;
56             
57             ans.a[1][1] = 1;ans.a[1][2] = 0;ans.a[1][3] = 0;ans.a[1][4] = 0;ans.a[1][5] = 0;ans.a[1][6] = 0;
58             ans.a[2][1] = 2;ans.a[2][2] = 2;ans.a[2][3] = 1;ans.a[2][4] = 0;ans.a[2][5] = 0;ans.a[2][6] = 0;
59             ans.a[3][1] = 3;ans.a[3][2] = 3;ans.a[3][3] = 0;ans.a[3][4] = 0;ans.a[3][5] = 0;ans.a[3][6] = 0;
60             ans.a[4][1] = 1;ans.a[4][2] = 0;ans.a[4][3] = 0;ans.a[4][4] = 1;ans.a[4][5] = 0;ans.a[4][6] = 0;
61             ans.a[5][1] = 3;ans.a[5][2] = 1;ans.a[5][3] = 0;ans.a[5][4] = 2;ans.a[5][5] = 1;ans.a[5][6] = 0;
62             ans.a[6][1] = 2;ans.a[6][2] = 1;ans.a[6][3] = 0;ans.a[6][4] = 1;ans.a[6][5] = 1;ans.a[6][6] = 1;
63             ans = ans ^ (n-2);
64             ans = res * ans;
65             cout<<mod(ans.a[1][1])<<endl; 
66         }
67     }
68     return 0;
69 }