HDU 3306 矩阵快速幂

题解: S(n) = S(n - 1) + A(n)^2 = S(n-1) + x^2*A(n-1)^2 + 2*x*y*A(n-1)A(n-2) +y^2*A(n-2)^2构造矩阵

code : 

#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;

const int p = 10007;

int n, x, y;

struct Matrix{
    int mat[4][4];

    void clr(){memset(mat, 0, sizeof mat);}
    void E(){
        clr();
        for(int i = 0; i < 4; ++i) mat[i][i] = 1;
    }
    Matrix operator * (const Matrix &rhs) const{
        Matrix res;
        res.clr();
        for(int i = 0; i < 4; ++i)
            for(int j = 0; j < 4; ++j)
                for(int k = 0; k < 4; ++k)
                       res.mat[i][j] = (res.mat[i][j] + mat[i][k] * rhs.mat[k][j]) % p;
        return res;
    }
    Matrix operator ^ (int b){
        Matrix res, tmp;
        res.E();
        memcpy(tmp.mat, mat, sizeof mat);
        while(b){
            if(b & 1) res = res * tmp;
            tmp = tmp * tmp;
            b >>= 1;
        }
        return res;
    }
};
int main(){
    //freopen("in.txt", "r", stdin);
    while(cin >> n >> x >> y){
        x %= p;
        y %= p;

        Matrix base;
        base.clr();

        base.mat[0][0] = 1; 
        /**row 1*/

        base.mat[1][0] = base.mat[1][1] = x * x % p;
        base.mat[1][2] = 1;
        base.mat[1][3] = x;
        /**row 2*/

        base.mat[2][0] = base.mat[2][1] = y * y % p;
        /**row 3*/

        base.mat[3][0] = base.mat[3][1] = 2 * x * y % p;
        base.mat[3][3] = y;
        /**row 4*/

        Matrix data;
        data.clr();
        data.mat[0][0] = 2;
        data.mat[0][1] = data.mat[0][2] = data.mat[0][3] = 1;

        Matrix res = data * (base ^ (n - 1));

        cout << res.mat[0][0] << endl;
    }
    return 0;
}