题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3306

题目大意：A(0) = 1 , A(1) = 1 , A(N) = X * A(N - 1) + Y * A(N - 2) (N >= 2)；给定三个值N，X，Y求S(N):S(N) = A(0)^2 +A(1)^2+……+A(n)^2。

 学了这几题，还是不太很懂，后来看题解，渐渐也是懂了一点。

题目的意思是求出A(0)^2 +A(1)^2+……+A(n)^2

考虑1*4 的矩阵【s[n-2],a[n-1]^2,a[n-2]^2,a[n-1]*a[n-2]】

我们需要找到一个4×4的矩阵A，使得它乘以A得到1×4的矩阵

【s[n-1],a[n]^2,a[n-1]^2,a[n]*a[n-1]】

即：【s[n-2],a[n-1]^2,a[n-2]^2,a[n-1]*a[n-2]】* A = 【s[n-1],a[n]^2,a[n-1]^2,a[n]*a[n-1]】

= 【s[n-2]+a[n-1]^2 , x^2 * a[n-1]^2 + y^2 * a[n-2]^2 + 2*x*y*a[n-1]*a[n-2] ,

a[n-1]^2 , x*a[n-1]^2 + y*a[n-2]a[n-1]】

可以构造矩阵A为：

1     0    0    0

1    x^2   1    x

0    y^2   0    0

0    2xy   0    y

故：【S[0],a[1]^2,a[0]^2,a[1]*a[0]】 * A^(n-1) = 【s[n-1],a[n]^2,a[n-1]^2,a[n]*a[n-1]】

所以：【S[0],a[1]^2,a[0]^2,a[1]*a[0]】 * A^(n) = 【s[n],a[n+1]^2,a[n]^2,a[n+1]*a[n]】

若A = (B * C ) 则AT = ( B * C )T = CT * BT

#include<stdio.h>
#include<string.h>
#define N 4
#define M 10007
struct Matrix
{
    __int64 a[N][N];
}res,tmp,ans,origin;
Matrix A={1,0,0,0,
          1,0,0,0,
          1,0,0,0,
          1,0,0,0};   //相当于那个1*4的矩阵的转置,即[s0,a1^2,a0^2,a1*a0]T
Matrix mul(Matrix x,Matrix y)
{
    int i,j,k;
    memset(tmp.a,0,sizeof(tmp.a));
    for(i=0;i<4;i++)
        for(j=0;j<4;j++)
            for(k=0;k<4;k++)
            {
                tmp.a[i][j]+=(x.a[i][k]*y.a[k][j])%M;
                tmp.a[i][j]%=M;
            }
    return tmp;
}
Matrix quickpow(int k)
{
    int i;
    memset(res.a,0,sizeof(res.a));
    for(i=0;i<4;i++)
        res.a[i][i]=1;
    while(k)
    {
        if(k&1)
            res=mul(res,origin);
        origin=mul(origin,origin);
        k>>=1;
    }
    return res;
}
int main()
{
    int n,x,y;
    while(scanf("%d%d%d",&n,&x,&y)!=EOF)
    {
        x%=M;
        y%=M;
        memset(origin.a,0,sizeof(origin.a));
        origin.a[0][0]=origin.a[0][1]=origin.a[2][1]=1;
        origin.a[1][1]=(x*x)%M;
        origin.a[1][2]=(y*y)%M;
        origin.a[1][3]=(2*x*y)%M;
        origin.a[3][1]=x;
        origin.a[3][3]=y;
        res=quickpow(n);  //求构造出的矩阵A^n
        ans=mul(res,A);   //A^n * [s0,a1^2,a0^2,a1*a0]T
        printf("%I64d\n",ans.a[0][0]%M);
    }
    return 0;
}