[POJ3070]Fibonacci sequence——矩阵+快速幂

#斐波那契数列#
题目描述
斐波那契数列是由如下递推式定义的数列
$F_0=0$
$F_1=1$
$F_{n+1}=F_{n+1}+F_n$
求这个数列第n项的值对$10^4$取余后的结果。
限制条件
$\cdot$ $0⩽n⩽10^{16}$

    我也是刚刚才搞懂了矩阵乘法(如果你不知道什么是矩阵乘法的话,右转百度百科)，于是来应用一下新知识，如有表述不到位的地方请见谅。
下面进入正文

首先，我们先介绍一下对于斐波那契数列如何求解。把斐波那契数列的递推式表示成矩阵就得到下面的式子
$\left(\begin{array}{c}{F}_{n+2}\\ F_{n+1}\end{array}\right)=$ $\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)$ $\left(\begin{array}{c}{F}_{n+1}\\ F_n\end{array}\right)$

我们发现式子里有个固定的矩阵 $\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)$

记这个矩阵为A，则有
$\left(\begin{array}{c}{F}_{n+1}\\ F_n\end{array}\right)=$$A^n$$\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)=$$A^n$$\left(\begin{array}{c}1\\ 0\end{array}\right)$

因此只要求出$A^n$就可以求出$F_n$了。关于$A^n$的计算我们可以采用类似快速幂的算法,在$O(logn)$时间里求出第n项的值。

用结构体定义矩阵$\left(\begin{array}{cc}a0& a1\\ a2& a3\end{array}\right)$，并重载矩阵乘法。

struct Matrix{
	int a[4];
	Matrix operator *(const Matrix &tmp)const{
		Matrix res;
		res.a[0]=(a[0]*tmp.a[0]+a[1]*tmp.a[2])%P;
		res.a[1]=(a[0]*tmp.a[1]+a[1]*tmp.a[3])%P;
		res.a[2]=(a[2]*tmp.a[0]+a[3]*tmp.a[2])%P;
		res.a[3]=(a[2]*tmp.a[1]+a[3]*tmp.a[3])%P;
		return res;	
	}
	Matrix(){memset(a,0,sizeof(a));}
}A[Max];

计算前要先初始化

void Init(){
	A[0].a[0]=1;
	A[0].a[1]=1;
	A[0].a[2]=1;
	for(int i=1;i<Max;i++)
	A[i]=A[i-1]*A[i-1];
}

求答案时要注意初始的res为$\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$

int main(){
	Init();
	Rd(n);
	Matrix res;
	res.a[0]=1;
	res.a[3]=1;
	for(int i=Max-1;i>=0;i--)
	if((n>>i)&1)res=res*A[i];
	printf("%d\n",res.a[2]);
	return 0;
}

完整的代码如下：

#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<iostream>
#define P 10000
#define M 1000005
#define Max 25
using namespace std;
void Rd(int &res){
	char c;res=0;int k=1;
	while(c=getchar(),c<48&&c!='-');
	if(c=='-'){k=-1;c='0';}
	do{
		res=(res<<3)+(res<<1)+(c^48);
	}while(c=getchar(),c>=48);
	res*=k;
}
struct Matrix{
	int a[4];
	Matrix operator *(const Matrix &tmp)const{
		Matrix res;
		res.a[0]=(a[0]*tmp.a[0]+a[1]*tmp.a[2])%P;
		res.a[1]=(a[0]*tmp.a[1]+a[1]*tmp.a[3])%P;
		res.a[2]=(a[2]*tmp.a[0]+a[3]*tmp.a[2])%P;
		res.a[3]=(a[2]*tmp.a[1]+a[3]*tmp.a[3])%P;
		return res;	
	}
	Matrix(){memset(a,0,sizeof(a));}
}A[Max];
void Init(){
	A[0].a[0]=1;
	A[0].a[1]=1;
	A[0].a[2]=1;
	for(int i=1;i<Max;i++)
	A[i]=A[i-1]*A[i-1];
}
int n;
int main(){
	Init();
	Rd(n);
	Matrix res;
	res.a[0]=1;
	res.a[3]=1;
	for(int i=Max-1;i>=0;i--)
	if((n>>i)&1)res=res*A[i];
	printf("%d\n",res.a[2]);
	return 0;
}

第一次写博客，谢谢各位支持。O(∩_∩)O