hdu 3306 Another kind of Fibonacci（矩阵快速幂）

来源：互联网发布：人工智能与医疗编辑：程序博客网时间：2024/04/29 00:40

Another kind of Fibonacci

Time Limit: 3000/1000 MS (Java/Others) Memory Limit: 65536/65536 K (Java/Others)
Total Submission(s): 1691 Accepted Submission(s): 660

Problem Description

As we all known , the Fibonacci series : F(0) = 1, F(1) = 1, F(N) = F(N - 1) + F(N - 2) (N >= 2).Now we define another kind of Fibonacci : A(0) = 1 , A(1) = 1 , A(N) = X * A(N - 1) + Y * A(N - 2) (N >= 2).And we want to Calculate S(N) , S(N) = A(0)² +A(1)²+……+A(n)².

Input

There are several test cases.
Each test case will contain three integers , N, X , Y .
N : 2<= N <= 2³¹ – 1
X : 2<= X <= 2³¹– 1
Y : 2<= Y <= 2³¹ – 1

Output

For each test case , output the answer of S(n).If the answer is too big , divide it by 10007 and give me the reminder.

Sample Input

2 1 1 3 2 3

Sample Output

主要是递推矩阵，然后矩阵快速幂：

根据S(n)=S(n-1)+A(n)^2，A(n)^2=X^2*A(n-1)^2+Y^2*A(n-2)^2+2*X*Y*A(n-1)*A(n-2)，A(n)*A(n-1)=Y*A(n-1)*A(n-2)+X*A(n-1)^2，可以构造如下的矩阵，然后用二分矩阵的方法求解即可。

递推矩阵

以上摘自：http://www.cnblogs.com/staginner/archive/2012/04/26/2471507.html

代码：

#include <iostream>#include <cstdio>#include <cstring>#include <algorithm>using namespace std;const int mod=10007;struct matrix{    long long ma[5][5];};matrix multi(matrix x,matrix y)//矩阵相乘{    matrix ans;    memset(ans.ma,0,sizeof(ans.ma));    for(int i=1;i<=4;i++)    {        for(int j=1;j<=4;j++)        {            if(x.ma[i][j])//稀疏矩阵优化            for(int k=1;k<=4;k++)            {                ans.ma[i][k]=(ans.ma[i][k]+(x.ma[i][j]*y.ma[j][k])%mod)%mod;            }        }    }    return ans;}matrix pow(matrix a,long long m){    matrix ans;    for(int i=1;i<=4;i++)    {        for(int j=1;j<=4;j++)        {            if(i==j)            ans.ma[i][j]=1;            else            ans.ma[i][j]=0;        }    }    while(m)    {        if(m&1)        ans=multi(ans,a);        a=multi(a,a);        m=m>>1;    }    return ans;}int main(){    long long x,y,n;    while(~scanf("%I64d%I64d%I64d",&n,&x,&y))    {        matrix a,b;        memset(a.ma,0,sizeof(a.ma));        memset(b.ma,0,sizeof(b.ma));        a.ma[1][1]=1;        a.ma[1][2]=1;        a.ma[2][2]=(x*x)%mod;        a.ma[2][3]=(y*y)%mod;        a.ma[2][4]=(2*x*y)%mod;        a.ma[3][2]=1;        a.ma[4][2]=x;        a.ma[4][4]=y;        b.ma[1][1]=1;        b.ma[2][1]=1;        b.ma[3][1]=1;        b.ma[4][1]=1;        a=pow(a,n);        a=multi(a,b);        printf("%I64d\n",a.ma[1][1]);    }    return 0;}

0 0