hdu3306矩阵快速幂讲解

矩阵快速幂可以说是数论题中较为常见也是较为简单的一部分。原理很直接，就是根据题目所给的公式推出递推公式，用矩阵乘法表示，然后套用快速幂模板得出结果。

比如hdu3306--Another kind of Fibonacci

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
    
    

 

很显然我们先推s[n]=s[n-1]+A[n]^2;A[n]^2=x*x*A[n-1]^2+y*y*A[n-2]^2+2*x*y*A[n-1][n-2];于是我们发现又涉及到了A[n-1]*A[n-2]。A[n]*A[n-1]=x*A[n-1]^2+y*A[n-1]*A[n-2]；

这样矩阵就好构造了：分列四项s[n-1]、A[n-1]^2、A[n-2]^2、A[n-1]*A[n-2]；另一边的：s[n]、A[n]^2、A[n-1]^2、A[n]*A[n-1]。

而矩阵内容就是各项系数了。代码如下：

#include "cstdio"
#include "cstring"
#include "algorithm"
using namespace std;
typedef struct
{
  int matrix[4][4];
}Matrix;

Matrix multi(Matrix x,Matrix y)
{
  Matrix res;
  int i,j,k;
  int sum;

  for(i = 0;i<4;i++)
    for(j = 0;j<4;j++)
    {
      sum = 0;
      for(k = 0;k<4;k++)
        sum+=(x.matrix[i][k]*y.matrix[k][j])%10007;
      res.matrix[i][j] = sum%10007;
    }
    return res;
}

Matrix powermod(Matrix x,int n)
{
  Matrix res;
  int i,j;

  for(i = 0;i<4;i++)
    for(j = 0;j<4;j++)
    {
       if(i==j)
        res.matrix[i][j] = 1;
       else
        res.matrix[i][j] = 0;
    }
    for(;n;n>>=1)
    {
      if(n&1)
        res = multi(res,x);
       x = multi(x,x);
    }
    return res;
}

int main()
{
    int N,X,Y;
    while(~scanf("%d %d %d",&N,&X,&Y))
    {
        if(N <= 1)
        {
            printf("%d\n",N+1);
            continue;
        }
        X%=10007,Y%=10007;
        Matrix res;
        res.matrix[0][0]=1;
        res.matrix[0][1]=(X*X)%10007;
        res.matrix[0][2]=(Y*Y)%10007;
        res.matrix[0][3]=(2*X*Y)%10007;
        res.matrix[1][0]=0;
        res.matrix[1][1]=(X*X)%10007;
        res.matrix[1][2]=(Y*Y)%10007;
        res.matrix[1][3]=(2*X*Y)%10007;
        res.matrix[2][0]=0;
        res.matrix[2][1]=1;
        res.matrix[2][2]=0;
        res.matrix[2][3]=0;
        res.matrix[3][0]=0;
        res.matrix[3][1]=X;
        res.matrix[3][2]=0;
        res.matrix[3][3]=Y;
        res=powermod(res,N-1);
        printf("%d\n",(2*res.matrix[0][0]+res.matrix[0][1]+res.matrix[0][2]+res.matrix[0][3])%10007);
    }
}