HDU 3117 Fibonacci Numbers（矩阵快速幂求斐波那契）

思路：

对于小于1e8的斐波那契数，直接打表输出即可；
对于大于1e8的数，我们需要分别计算它的前四位和后四位；
计算前四位：
1.斐波那契数的通项公式为：
$$f(n)=\frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2}{)}^n-(\frac{1-\sqrt{5}}{2}{)}^n]$$
记我们要算s的前四位d，s的长度为len,则
$$s=d.xxx\ast 10^{len-4}$$
两边取lg，得
$$\lg s=\lg d.xxx+len-4$$
从而得到
$$d.xxx=10^{\lg s+4-len}$$
根据斐波那契数的通项，我们可以求出当n足够大时，s的近似值
$$s=\frac{1}{\sqrt{5}}(\frac{1+\sqrt{5}}{2}{)}^n$$
而s的长度len可以表示为
$$len=(int)\lg s+1$$
因此我们可以得到最终d的求解公式为
$$\{\begin{array}{r}\lg s=n\ast \lg \frac{1+\sqrt{5}}{2}-lg\sqrt{5}\\ \\ d=(int)10^{\lg s-(int)\lg s+3}\end{array}$$
计算后四位：
由斐波那契数列的性质我们可以得到以下递推矩阵：
$$\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\ast \left[\begin{array}{c}f(n)\\ f(n-1)\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]}^{n-1}\ast \left[\begin{array}{c}1\\ 0\end{array}\right]=\left[\begin{array}{c}f(n)\\ f(n-1)\end{array}\right]$$
由于右边的行列式是[1,0]因此我们只需计算出n-1次方后取(0,0)位置的值即是答案（当然要在过程中取余）；
计算方法就用矩阵的快速幂，来降低时间复杂度~

代码：

#include<iostream>
#include<cmath>
#include<vector>
using namespace std;
typedef long long LL;
const int MAX_F=1e8;
vector<int> f(2);
double sq5=sqrt(5);
int first_4(int n){//返回前四位数字 
	double log_s=n*log10((1+sq5)/2)-log10(sq5);
	return (int)pow(10,log_s-(int)log_s+3);
}
typedef vector<int> vec;
typedef vector<vec> mat;
const int MOD=1e4;
mat mul(mat &A,mat &B){ 
	mat C(A.size(),vec(B[0].size()));
	for(int i=0;i<A.size();i++){
		for(int k=0;k<B.size();k++){
			for(int j=0;j<B[0].size();j++){
				C[i][j]=(C[i][j]+A[i][k]*B[k][j])%MOD;
			}
		}
	}
	return C;
}
mat pow(mat A,LL n){
	mat B(A.size(),vec(A.size()));
	for(int i=0;i<A.size();i++) B[i][i]=1;
	while(n>0){
		if(n&1) B=mul(B,A);
		A=mul(A,A);
		n>>=1;
	}
	return B;
}
int main(){
//	freopen("in.txt","r",stdin);
	f[0]=0; f[1]=1;	
	for(int i=2;f[i-1]<MAX_F;i++) f.push_back(f[i-1]+f[i-2]); 
	int n,sz=f.size()-1;
	while(~scanf("%d",&n)){
		if(n<sz) printf("%d\n",f[n]);
		else{
			printf("%d...",first_4(n));
			mat A(2,vec(2));
			A[0][0]=A[0][1]=A[1][0]=1; A[1][1]=0;
			mat B=pow(A,n-1);
			printf("%04d\n",B[0][0]);
		}
	}
	return 0;
}