HDU 6172-Array Challenge

Array Challenge
Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 153428/153428 K (Java/Others)
Total Submission(s): 97 Accepted Submission(s): 38

Problem Description
There’s an array that is generated by following rule.
h0=2,h1=3,h2=6,hn=4hn−1+17hn−2−12hn−3−16
And let us define two arrays bnandan as below.
bn=3hn+1hn+9hn+1hn−1+9h2n+27hnhn−1−18hn+1−126hn−81hn−1+192(n>0)
an=bn+4n
Now, you have to print ⌊√(an)⌋ , n>1.
Your answer could be very large so print the answer modular 1000000007.

Input
The first line of input contains T (1 <= T <= 1000) , the number of test cases.
Each test case contains one integer n (1 < n <= 1015) in one line.

Output
For each test case print ⌊√(an)⌋ modular 1000000007.

Sample Input

3
4
7
9

Sample Output

1255
324725
13185773

Source
2017 Multi-University Training Contest - Team 10

题意：h0=2,h1=3,h2=6,hn=4hn−1+17hn−2−12hn−3−16
bn=3hn+1hn+9hn+1hn−1+9h2n+27hnhn−1−18hn+1−126hn−81hn−1+192(n>0)
an=bn+4n
求⌊√(an)⌋
一波打表后可以发现答案和hn十分接近，由此猜想答案和hn 有较接近的递推公式，实际上ansn=4ansn−1+17ansn−2−12ansn−3，
事实上进一步地有ansn=4ansn−1−7ansn−2
n≤1015,用矩阵快速幂做，[0−417][ansnansn+1]=[ansn+1ansn+2]

• AC代码

#include <bits/stdc++.h>using namespace std;#define mem(s,v) memset(s,v,sizeof(s))typedef long long ll;const double PI = acos(-1);const ll MOD = 1000000007;const int inf = 0x3f3f3f3f;const ll INF = 0x3f3f3f3f3f3f3f3f;const int N = 2;//方阵维数struct Mat{    ll mat[N][N];}E;//初始化单位矩阵Mat init(){    Mat E;    for(int i = 0; i < N; i++){        for(int j = 0; j < N; j++){            if(i == j)            E.mat[i][i] = 1;            else            E.mat[i][j] = 0;        }    }    return E;}//重载乘法Mat operator *(Mat a,Mat b){    Mat c;    memset(c.mat,0,sizeof(Mat));    for(int i = 0; i < N; i++){        for(int j = 0; j < N; j++){            for(int k = 0; k < N; k++){                if(a.mat[i][k] && b.mat[k][j]){                    c.mat[i][j] += a.mat[i][k] * b.mat[k][j];                    c.mat[i][j] %= MOD;                }            }        }    }    return c;}//重载加法Mat operator +(Mat a,Mat b){    Mat c;    memset(c.mat,0,sizeof(Mat));    for(int i = 0; i < N; i++){        for(int j = 0; j < N; j++){            c.mat[i][j] = a.mat[i][j] + b.mat[i][j];        }    }    return c;}//矩阵快速幂Mat operator ^(Mat A,ll x){    Mat c;    c = init();    for(; x ; x >>= 1){        if(x&1){            c = c*A;        }        A = A*A;    }    return c;}ll get(ll x){    if(x == 2)return 31;    if(x == 3)return 197;    Mat EE;    EE.mat[0][0] = 0;    EE.mat[0][1] = 1;    EE.mat[1][0] = -4;    EE.mat[1][1] = 7;    ll beg[2] = {31,197};    EE = EE^(x - 3);    ll s1 = EE.mat[0][0]*beg[0] + EE.mat[0][1]*beg[1];    ll s2 = EE.mat[1][0]*beg[0] + EE.mat[1][1]*beg[1];    while(s2<0) s2+=MOD;    s2%=MOD;    return s2;}int main(){    int T;    scanf("%d",&T);    ll n;    while(T--){        scanf("%lld",&n);        printf("%lld\n",get(n));    }}