矩阵二分乘

题目描述：
给定a0,a1,以及an=p*a(n-1) + q*a(n-2)中的p,q。这里n >= 2。 求第k个数对10000的模。

输入：
输入包括5个整数：a0、a1、p、q、k。

输出：
第k个数a(k)对10000的模。

样例输入：

20 1 1 14 5

样例输出：

8359

来源：2009年清华大学计算机研究生机试真题看到这个题目，如果按照最普通的方法（根据递推公式，带入每组输入的数，最后对a(k)求余），这样多半会超时。因为结果要对10000取模，所以我想到找规律，测试了几组数据发现很难找到什么规律。哎~~输入的变量太多了，不好控制。没办法，上网看到有人说用矩阵二分乘来做。What？矩阵二分乘是什么东东？弱爆了~~~下面引用了jzlikewei的文章里的内容，从快速取幂到矩阵二分乘。
快速取幂：

利用分治思想进行:

A^p=(A^(p/2))*(A^(p/2)) //p为偶数

A^p=(A^(p/2))*(A^(p/2)) *A //p为奇数

那么如何求A^(p/2),依旧用上面的两个式子，直至p=1。

一个递归程序就得到了：

longlong Qexp (int A, int p )

{

if (p==1)

return(A );

// if(0==p)

//　　　　return 1;

longlong tmp=Qexp(A,p/2);

if ( p %2 ==0)

return(tmp *tmp );

return(tmp*tmp*A);

}

矩阵二分乘

很多递推关系通常并不是一个前项对应一个后项的，二对一甚至多对一都是可能，一个臭名昭著著名的例子就是斐波那契数列(F[i]=f[i-1]+f[i-2],f[0]=0;f[1]=1)，这些式子很难写出通项公式。

以斐波那契数列数列举例,如何快速求出第N项？

考虑矩阵乘法:

#include <iostream>
#include <cstdio>
using namespace std;
struct Matrix_2_2
{
int a[2][2];
Matrix_2_2 operator*(Matrix_2_2 &m)
{
Matrix_2_2 new_m;
new_m.a[0][0]=(a[0][0]*m.a[0][0]+a[0][1]*m.a[1][0])%10000;
new_m.a[0][1]=(a[0][0]*m.a[0][1]+a[0][1]*m.a[1][1])%10000;
new_m.a[1][0]=(a[1][0]*m.a[0][0]+a[1][1]*m.a[1][0])%10000;
new_m.a[1][1]=(a[1][0]*m.a[0][1]+a[1][1]*m.a[1][1])%10000;
return new_m;
}
};

Matrix_2_2 MatrixPower(Matrix_2_2 &m,int n)
{
if(0==n)
{
Matrix_2_2 m0;
m0.a[0][0]=1;
m0.a[0][1]=0;
m0.a[1][0]=0;
m0.a[1][1]=1;
return m0;
}
Matrix_2_2 temp;
temp=MatrixPower(m,n/2);
if(n&1)
return temp*temp*m;
else
return temp*temp;
}
int main()
{
int a0,a1,p,q,k;
Matrix_2_2 m,m_last;
while(scanf("%d%d%d%d%d",&a0,&a1,&p,&q,&k)!=EOF)
{
if(0==k)
{
printf("%d\n",a0);
continue;
}
if(1==k)
{
printf("%d\n",a1);
continue;
}
m.a[0][0]=p;
m.a[0][1]=q;
m.a[1][0]=1;
m.a[1][1]=0;
m_last=MatrixPower(m,k-1);
int res=((m_last.a[0][0]*a1)%10000+(m_last.a[0][1]*a0)%10000)%10000;
printf("%d\n",res);
}
return 0;
}

View Code
2013,Goodbye.