HDU 3306 Another kind of Fibonacci（矩阵快速幂）

Another kind of Fibonacci

As we all known , the Fibonacci series : F(0) = 1, F(1) = 1, F(N) = F(N - 1) + F(N - 2) (N >= 2).Now we define another kind of Fibonacci : A(0) = 1 , A(1) = 1 , A(N) = X * A(N - 1) + Y * A(N - 2) (N >= 2).And we want to Calculate S(N) , S(N) = A(0) 2 +A(1) 2+……+A(n) 2.

Input
There are several test cases.
Each test case will contain three integers , N, X , Y .
N : 2<= N <= 2 31 – 1
X : 2<= X <= 2 31– 1
Y : 2<= Y <= 2 31 – 1

Output
For each test case , output the answer of S(n).If the answer is too big , divide it by 10007 and give me the reminder.

Sample Input

2 1 1
3 2 3

Sample Output

6
196

题意
给出$A_0=1,A_1=1,A_n=X{A}_{n-1}+Y{A}_{n-2}$,求$S_n,(S_n=A_0^2+A_1^2+\cdots +A_n^2)$

题解
矩阵快速幂，唯一一点比较麻烦的就是这个递推公式了。
$$\left(\begin{array}{c}{A}_n\\ A_{n-1}\\ A_{n-2}\end{array}\right)=\left(\begin{array}{ccc}X& Y& 0\\ 1& 0& 0\\ 0& 1& 0\end{array}\right)\left(\begin{array}{c}{A}_{n-1}\\ A_{n-2}\\ A_{n-3}\end{array}\right)$$
$$\left(\begin{array}{c}{S}_n\\ A_n^2\\ A_{n-1}^2\\ A_n{A}_{n-1}\end{array}\right)=\left(\begin{array}{cccc}1& X^2& Y^2& 2XY\\ 0& X^2& Y^2& 2XY\\ 0& 1& 0& 0\\ 0& X& 0& Y\end{array}\right)\left(\begin{array}{c}{S}_{n-1}\\ A_{n-1}^2\\ A_{n-2}^2\\ A_{n-1}{A}_{n-2}\end{array}\right)=T^{n-1}\left(\begin{array}{c}{S}_1\\ A_1^2\\ A_0^2\\ A_1{A}_0\end{array}\right)$$
这是两个矩阵递推式，第一个没啥大用。。。主要是第二个。$T$是转换矩阵。其中$A_n,S_n,A_n^2,A_{n-1}^2$显而易见，而$A_n{A}_{n-1}=(X{A}_{n-1}+Y{A}_{n-2})A_{n-1}=X{A}_{n-1}^2+Y{A}_{n-1}{A}_{n-2}$

实在不记得自己写这题的时候脑子哪里抽了，死活找不到bug，看了半天实在受不了看了博客，发现博客写$S_1=1$我瞬间怀疑了人生，我没觉得自己的公式推错了，而且按着题中的要求，$S_1=2$才对啊，但是一直WA，怀疑人生一下午，问了实验室三个人，最后又把博客的$S_1=1$的公式又改回$S_1=2$又交一次，突然就AC了？？？我自己瞬间爆炸

代码

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <queue>
#include <cmath>
#include <string>
#include <cstring>
#include <map>
#include <set>
#include <cmath>
#include <vector>

//#include <bits/stdc++.h>
using namespace std;
#define me(x,y) memset(x,y,sizeof x)
#define MIN(x,y) x < y ? x : y
#define MAX(x,y) x > y ? x : y
typedef long long ll;
const int maxn = 20;
const int INF = 0x3f3f3f3f;
const int MOD = 10007;
const int eps = 1e-11;
const double PI = acos(-1.0);

struct  Matrix
{
    ll m[maxn][maxn];
};

Matrix Mul(Matrix A,Matrix B,int n){
    Matrix tmp;
    for(int i = 1; i <= n; ++i){
        for(int j = 1; j <= n; ++j){
            tmp.m[i][j] = 0;
        }
    }
    for(int i = 1; i <= n; ++i){
        for(int j = 1; j <= n; ++j){
            for(int k = 1; k <= n; ++k){
                tmp.m[i][j] += A.m[i][k]*B.m[k][j];
                tmp.m[i][j] %= MOD;
            }
        }
    }
    return tmp;
}
Matrix Qpow(Matrix A,int N,int n){ //求A的N次幂
    Matrix ans;
    for(int i = 1; i <= n; ++i){
        for(int j = 1; j <= n; ++j){
            if(i == j) ans.m[i][j] = 1;
            else ans.m[i][j] = 0;
        }
    }
        
    while(N){
        if(N&1) ans = Mul(ans,A,n);
        A = Mul(A,A,n);
        N >>=1;
    }
    return ans;
}

int main(){
    ll n,x,y;
    while(cin>>n>>x>>y){
        Matrix b;
        b.m[1][1] = 1,b.m[1][2] = x*x%MOD,b.m[1][3] = y*y%MOD,b.m[1][4] = 2*x*y%MOD;
        b.m[2][1] = 0,b.m[2][2] = x*x%MOD,b.m[2][3] = y*y%MOD,b.m[2][4] = 2*x*y%MOD;
        b.m[3][1] = 0,b.m[3][2] = 1,b.m[3][3] = 0,b.m[3][4] = 0;
        b.m[4][1] = 0,b.m[4][2] = x%MOD,b.m[4][3] = 0,b.m[4][4] = y%MOD;
        Matrix c;
        c = Qpow(b,n-1,4);
        ll ans = (c.m[1][1]*2+c.m[1][2]+c.m[1][3]+c.m[1][4])%MOD;
        cout<<ans<<endl;
    }
    
    return 0;
}

/*

*/