题解 HDU3306 Another kind of Fibonacci

Posted 2020-12-18 point-king

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

题目大意

众所周知，斐波那契数列为 (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)) 的递推式：

[egin{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[egin{matrix} S(n)&A^2(n)&A^2(n-1)&A(n)A(n-1)\end{matrix} ight] ]

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

[egin{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[egin{matrix} 1&0&0&0x^2&x^2&1&xy^2&y^2&0&02xy&2xy&0&y\end{matrix} ight] ]

即：

[egin{array}{} left[egin{matrix} S(n-1)&A^2(n-1)&A^2(n-2)&A(n-1)A(n-2)\end{matrix} ight] * left[egin{matrix} 1&0&0&0x^2&x^2&1&xy^2&y^2&0&02xy&2xy&0&y\end{matrix} ight] = left[egin{matrix} S(n)&A^2(n)&A^2(n-1)&A(n)A(n-1)\end{matrix} ight]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
",ans.h[1][1]);
    }
}