本篇题解用于作者本人对于矩阵乘法的印象加深，也欢迎大家的阅读。

题目大意

众所周知，斐波那契数列为 \(f(0)=1\) ， \(f(1)=1\) ，\(f(n)=f(n-1)+f(n-2)~(n>=2)\) 。定义另一种斐波那契数列： \(A(0)=1\) ， \(A(1)=1\) ， \(A(n)=x*A(n-1)+y*A(n-2)~(n>=2)\) 。我们要计算 \(S(n)\) ， \(S(n)=A(0)^2+A(1)^2+...+A(n)^2\) 。

题解

我们可以很轻易的发现这是一道矩阵乘法的题，因为他是求关于一个递推式的平方和，而本题的难点就在于如何构建出合适的加速矩阵。

本题求的是 \(S(n)\) ，所以我们可以从 \(S(n)\) 递推式入手。

\[S(n)=S(n-1)+A^2(n) \]

所以， \(A^2(n)\) 肯定要列入我们的矩阵中。我们再来看看 \(A^2(n)\) 的递推式：

\[\begin{array}{} A^2(n)=(x*A(n-1)+y*A(n-2))^2\\ \\ =x^2A^2(n-1)+2xyA(n-1)A(n-2)+y^2A^2(n-2)\\ \end{array} \]

所以， \(A^2(n)\) ， \(A(n)A(n-1)\) 和 \(A^2(n-1)\) 也是需要加入矩阵的。因此我们的状态矩阵就是：

\[\left[\begin{matrix} S(n)&A^2(n)&A^2(n-1)&A(n)A(n-1)\\ \end{matrix}\right] \]

其中每一个元素的递推式如下：

\[\begin{array}{} S(n)=S(n-1)+x^2A^2(n-1)+2xyA(n-1)A(n-2)+y^2A^2(n-2)\\ \\ A^2(n)=x^2A^2(n-1)+2xyA(n-1)A(n-2)+y^2A^2(n-2)\\ \\ A^2(n-1)=A^2(n-1)\\ \\ A(n)A(n-1)=xA^2(n-1)+yA(n-1)A(n-2) \end{array} \]

我们再根据状态矩阵以及状态矩阵元素的递推式，我们可以求出加速矩阵：

\[\left[\begin{matrix} 1&0&0&0\\ x^2&x^2&1&x\\ y^2&y^2&0&0\\ 2xy&2xy&0&y\\ \end{matrix}\right] \]

即：

\[\begin{array}{} \left[\begin{matrix} S(n-1)&A^2(n-1)&A^2(n-2)&A(n-1)A(n-2)\\ \end{matrix}\right] * \left[\begin{matrix} 1&0&0&0\\ x^2&x^2&1&x\\ y^2&y^2&0&0\\ 2xy&2xy&0&y\\ \end{matrix}\right] = \left[\begin{matrix} S(n)&A^2(n)&A^2(n-1)&A(n)A(n-1)\\ \end{matrix}\right]\\ \end{array} \]

最后我们再套一个矩阵快速幂的模板就可以了。

代码如下：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int MOD=10007;
struct Matrix
{
    int n,m;
    ll h[5][5];
    Matrix()
    {
        memset(h,0,sizeof(h));
    }
    void Re1(int a)
    {
        n=m=a;
        memset(h,0,sizeof(h));
        for(int i=1;i<=a;++i)
        h[i][i]=1;
    }
};
Matrix operator * (const Matrix a,const Matrix b)
{
    Matrix ans;
    ans.n=a.n;
    ans.m=b.m;
    for(int i=1;i<=a.n;++i)
    {
        for(int j=1;j<=b.m;++j)
        {
            for(int k=1;k<=a.m;++k)
            {
                ans.h[i][j]+=a.h[i][k]*b.h[k][j]%MOD;
                ans.h[i][j]%=MOD;
            }
        }
    }
    return ans;
}
Matrix operator ^ (const Matrix xx,const ll kk)
{
    Matrix ans,x=xx;
    ll k=kk;
    ans.Re1(4);
    while(k>0)
    {
        if(k&1)
        ans=ans*x;
        x=x*x;
        k>>=1;
    }
    return ans;
}
int n;
ll x,y;
Matrix tmp,stp,ans;
int main()
{
    tmp.n=1;
    tmp.m=4;
    tmp.h[1][1]=2;
    tmp.h[1][2]=1;
    tmp.h[1][3]=1;
    tmp.h[1][4]=1;
    while(cin>>n>>x>>y)
    {
        stp.n=stp.m=4;
        stp.h[1][1]=stp.h[2][3]=1;
        stp.h[2][1]=stp.h[2][2]=x*x%MOD;
        stp.h[3][1]=stp.h[3][2]=y*y%MOD;
        stp.h[4][1]=stp.h[4][2]=2*x*y%MOD;
        stp.h[2][4]=x;
        stp.h[4][4]=y;
        ans=tmp*(stp^(n-1));
        printf("%lld\n",ans.h[1][1]);
    }
}