HDU3306-Another kind of Fibonacci（矩阵构造）

Another kind of Fibonacci

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

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

6 196

 

Author

wyb

 

虽然以前接触过矩阵，但是拿到这题后还是无从下手，后来上网看了好多有关矩阵的东西，总算明白如何用矩阵来求解问题了，这道题算是自己的一个矩阵入门吧。

这题与Fibonacci求f(n)很像，但是这题求得是前n项的平方和，已知,A(n)=x*A(n-1)+y*A(n-2)--(1),S(n)=A(0)^2+A(1)^2+……+A(n)^2=S(n-1)+A(n)^2--(2)

将(1)式带入(2)得，S(n)=S(n-1)+x*x*A(n-1)^2+y*y*A(n-2)^2+2*x*y*A(n-1)*A(n-2)--(3),根据(3)式可得矩阵B={S(n-1),A(n-1)^2,A(n-2)^2,A(n-1)*A(n-2)},

现在需要一个四阶矩阵A，使得B乘以A后得到{S(n),A(n)^2,A(n-1)^2,A(n)*A(n-1)},接下来就是构造这个函数A，不难得出

        |1       0      0     0|

        |x*x  x*x    1     x|

  A=  |y*y  y*y    0     0|

        |2xy  2xy    0     y|

最后{S(1) , A(1)^2 , A(0)^2 , A(1)*A(2)}*A^(n-1)={S(n) , A(n)^2 , A(n-1)^2 , A(n)*A(n-1)},该矩阵第一个元素即为所求答案。其中A^(n-1)用快速幂求解。

#include<stdio.h>
#include<string.h>
void fac1(__int64 a[4][4],__int64 b[4][4])
{
 __int64 i,j,k,c[4][4];
 memset(c,0,sizeof(c));
 for(i=0;i<4;i++)
 {
  for(j=0;j<4;j++)
  {
   for(k=0;k<4;k++)
   {
    c[i][j]+=(a[i][k]*b[k][j])%10007;
   }
  }
 }
 for(i=0;i<4;i++)
  for(j=0;j<4;j++)
   a[i][j]=c[i][j]%10007;
}
void fac2(__int64 a[4][4],__int64 b[4])
{
 __int64 i,j,c[4];
 memset(c,0,sizeof(c));
 for(i=0;i<4;i++)
 {
  for(j=0;j<4;j++)
   c[i]+=(b[j]*a[j][i])%10007;
 }
 for(i=0;i<4;i++)
  b[i]=c[i]%10007;
}
void pow(__int64 a[4][4],__int64 n)
{
 __int64 i,j,ans[4][4];
 memset(ans,0,sizeof(ans));
 for(i=0;i<4;i++)
  ans[i][i]=1;
 while(n>=1)
 {
  if(n%2)fac1(ans,a);
  n/=2;
  fac1(a,a);
 }
 for(i=0;i<4;i++)
  for(j=0;j<4;j++)
   a[i][j]=ans[i][j]%10007;
}
int main()
{
 __int64 n,x,y;
 while(scanf("%I64d%I64d%I64d",&n,&x,&y)!=EOF)
 {
  __int64 a[4][4],b[4];
  a[0][0]=a[1][2]=1;
  a[0][1]=a[0][2]=a[0][3]=a[2][2]=a[2][3]=a[3][2]=0;
  a[1][0]=a[1][1]=(x*x)%10007;
  a[1][3]=x%10007;
  a[2][0]=a[2][1]=(y*y)%10007;
  a[3][0]=a[3][1]=(2*x*y)%10007;
  a[3][3]=y%10007;
  b[0]=2;b[1]=1;b[2]=1;b[3]=1;
  pow(a,n-1);
  fac2(a,b);
  printf("%I64d ",b[0]%10007);
 }
 return 0;
}