hdu 3306 Another kind of Fibonacci 矩阵快速幂

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)².

s( n ) = s(n-1) + An^2 

An^2 = X^2*An-1 ^2 + Y^2*An-2^2 + 2*X*Y*An-1*An-2

An*An-1  = X*An-1^2 + Y*An-1 *An-2 

矩阵  【Sn-1 ， An-1*An-2  An-1^2    An-2 ^2】 

矩阵2

{ 1     0   0   0 }

{2XY  Y 2XY 0}

 {X^2  X  X^2 0 }

 {Y^2  0  Y^2  0}

通过矩阵就可以递推。

这样就可以 矩阵快速幂。

#include <vector>
#include <list>
#include <map>
#include <set>
#include <deque>
#include <stack>
#include <cstring>
#include <bitset>
#include <algorithm>
#include <functional>
#include <numeric>
#include <utility>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <assert.h>
#include <queue>
#define REP(i,n) for(int i=0;i<n;i++)
#define TR(i,x) for(typeof(x.begin()) i=x.begin();i!=x.end();i++)
#define ALLL(x) x.begin(),x.end()
#define SORT(x) sort(ALLL(x))
#define CLEAR(x) memset(x,0,sizeof(x))
#define FILLL(x,c) memset(x,c,sizeof(x))
using namespace std;
const double eps = 1e-9;
#define LL long long 
#define pb push_back
const int  maxn = 4;
const int size = 4;
const LL MOD =  10007 ;
LL X, Y;
LL n;
struct Mat{
	long long mat[maxn][maxn];
	void init(){
		 memset(mat,0,sizeof(mat));
		 mat[0][0] = 1;
		 mat[1][0] = 2*X%MOD*Y%MOD;  mat[1][1] = Y ; mat[1][2] = 2*X%MOD*Y%MOD;
		 mat[2][0] = X*X%MOD;        mat[2][1] = X%MOD;  mat[2][2] = X*X%MOD; mat[2][3] =1;
		 mat[3][0] = Y*Y%MOD;        mat[3][2] =  Y*Y%MOD;       
	}
	void show(){
		for(int i=0;i<=4;i++){
			for(int j=0;j<=4;j++){
				 cout << mat[i][j]<< " ";
			}
			cout <<endl;
		}
	}
	LL get_ans(){
		return (mat[1][0] + mat[2][0]+ mat[3][0])%MOD;
	}
}E,g;
void init_E(){
   memset(E.mat,0,sizeof(E.mat));
   for(int i=0;i<size ;i++){
      E.mat[i][i] =1 ;
   }
}
Mat operator*(const Mat &a,const Mat &b){
	Mat c;
	for(int i=0;i<size;i++)
		for(int j=0;j<size;j++){
			c.mat[i][j] = 0;
			for(int k=0;k<size;k++){
				if(a.mat[i][k] && b.mat[k][j])
					c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j]%MOD)%MOD;
			}
			
		}
	return c;
}

Mat operator^(Mat a,LL x){
	Mat c=E;
	for(;x;x>>=1){
		if(x&1)
			c=c*a;
		a=a*a;
	}
	return c;
}
Mat pow2(Mat a,LL x){
	Mat c = E;
	while(x){
		if(x&1){
			c = c*a;
		}
		a= a*a;
		x= x>>1;
	}
	return c;
}
void solve(){
	g.init();
	g = g^(n-1);
	LL ans = g.get_ans();
	ans +=2;
	ans %= MOD;
	cout << ans<<endl;
}
int main(){
	init_E();
    while(~scanf("%lld",&n)){
    	scanf("%lld%lld",&X,&Y);
    	
    	solve();
    }
    return 0;
}