hdu3306：Another kind of Fibonacci

A(0)=A(1)=1，A(i)=X*A(i-1)+Y*A(i-2)，求S(n)=A(0)^2+A(1)^2+A(2)^2+A(3)^2+……+A(n)^2。

这个矩阵有点毒。。

#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<math.h>
//#include<iostream>
using namespace std;

#define LL long long
LL n,x,y;
typedef LL mat[6][6];
mat a,f,ans;
const int mod=10007;
void copy(mat &a,mat b)
{
    for (int i=1;i<=4;i++)
        for (int j=1;j<=4;j++)
            a[i][j]=b[i][j];
}
void mul(mat a,mat b,mat &ans)
{
    mat t;memset(t,0,sizeof(t));
    for (int i=1;i<=4;i++)
        for (int j=1;j<=4;j++)
            for (int k=1;k<=4;k++)
                t[i][j]=(t[i][j]+a[i][k]*b[k][j]%mod)%mod;
    copy(ans,t);
}
void init(mat &a)
{
    memset(a,0,sizeof(a));
    for (int i=1;i<=4;i++) a[i][i]=1;
}
void pow(mat a,LL b,mat &ans)
{
    mat tmp,last;copy(tmp,a);init(last);
    while (b)
    {
        if (b&1) mul(last,tmp,last);
        mul(tmp,tmp,tmp);
        b>>=1;
    }
    copy(ans,last);
}
int main()
{
    while (~scanf("%lld%lld%lld",&n,&x,&y))
    {
        memset(a,0,sizeof(a));
        a[1][1]=a[1][2]=a[3][2]=1;
        a[2][2]=x*x%mod;
        a[2][3]=y*y%mod;
        a[2][4]=2*x*y%mod;
        a[4][2]=x;a[4][4]=y;
        f[1][1]=f[2][1]=f[3][1]=f[4][1]=1;
        pow(a,n,a);
        mul(a,f,ans);
        printf("%lld
",ans[1][1]);
    }
    return 0;
}