hdu3306矩阵快速幂讲解

来源：互联网发布：mac book air能做什么编辑：程序博客网时间：2024/05/17 15:37

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

比如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

很显然我们先推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);    }}

0 0