HDU 3306 矩阵构造

 

Another kind of Fibonacci

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

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

6196

 

Author

wyb

 

Source

HDOJ Monthly Contest – 2010.02.06

 

以前使用矩阵二分幂来处理那种线性组合的式子是在matrix67的博客中看到的。

里面讲了那种类似于斐波那契的第n项的求法

 

但是那里面讲解的方法只是说了等式两边都只有关于f的式子

这个题等式两边又有S又有A

但是求解方法还是一样的

 

这里就把矩阵哥的方法推广了一下

只要两边的未知量有递推的关系，那么就可以用矩阵来构造

 

这个题

A[N]=X*A[N-1]+Y*A[N-2]

S[N]=S[N-1]+A[N]^2

两个式子联立一下

变成了：

S[N]=S[N-1]+X^2*A[N-1]^2+Y^2*A[N-2]^2+2XYA[N-1]*A[N-2]

这里可以构造一个矩阵

分别对应着S[N]  A[N]^2   A[N-1]^2   A[N]*A[N-1]

构造的矩阵为：

1          0         0     0

 x^2      x^2     1     x

 y^2      y^2     0     0

2*x*y   2*x*y  0     y

 

注意到矩阵:S[N-1],A[N-1]^2,A[N-2]^2,A[N-1]*A[N-2]

乘以这个矩阵之后就可以得到目标矩阵了

所以转化为了

类似于斐波那契数列的那种求第n项的问题！

 

我的代码先开始一直超时，后来我把很多取模的冗余操作去掉，然后再把某一些变量由原来的int变成了__int64

就可以过了，大概是500++MS

 

我的代码：

#include<stdio.h>#include<algorithm>#include<string.h>using namespace std;struct mart{__int64 mat[5][5];};__int64 mod=10007;mart kk;mart multi(mart a,mart b){mart c;int i,j,k;for(i=1;i<=4;i++)for(j=1;j<=4;j++){c.mat[i][j]=0;for(k=1;k<=4;k++)c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;}return c;}mart power(int k){mart p,q;int i,j;for(i=1;i<=4;i++)for(j=1;j<=4;j++){p.mat[i][j]=kk.mat[i][j];if(i==j)q.mat[i][j]=1;elseq.mat[i][j]=0;}if(k==0)return q;while(k!=1){if(k&1){k--;q=multi(p,q);}else{k=k>>1;p=multi(p,p);}}p=multi(p,q);return p;}void init(__int64 x,__int64 y){kk.mat[1][1]=1;kk.mat[1][2]=(x%mod)*(x%mod)%mod;kk.mat[1][3]=(y%mod)*(y%mod)%mod;kk.mat[1][4]=(2*x%mod)*(y%mod)%mod;kk.mat[2][1]=0;kk.mat[2][2]=(x%mod)*(x%mod)%mod;kk.mat[2][3]=(y%mod)*(y%mod)%mod;kk.mat[2][4]=(2*x%mod)*(y%mod)%mod;kk.mat[3][1]=0;kk.mat[3][2]=1;kk.mat[3][3]=0;kk.mat[3][4]=0;kk.mat[4][1]=0;kk.mat[4][2]=x%mod;kk.mat[4][3]=0;kk.mat[4][4]=y%mod;}int main(){int n,i;__int64 a[5],x,y;mart xx;while(scanf("%d%I64d%I64d",&n,&x,&y)!=EOF){a[1]=2,a[2]=1,a[3]=1,a[4]=1;init(x,y);xx=power(n-1);__int64 ans=0;for(i=1;i<=4;i++)ans=(ans+a[i]*xx.mat[1][i])%mod;printf("%I64d\n",ans);}return 0;}