HDU 3117 Fibonacci Numbers(矩阵快速幂+公式)

题目地址：HDU 3117

对于后四位可以用矩阵快速幂快速求出来，但前四位就没办法了。要知道斐波那契数列是有通项公式的，所以只能通过通项公式来求前四位，但公式不能求后四位，因为公式使用浮点数求的，精度显然不够，求前四位要用到对数。

通项公式为：

f(n)=1/sqrt(5)(((1+sqrt(5))/2)^n+((1-sqrt(5))/2)^n)
假设F
可以表示成 t * 10^k（t是一个小数），那么对于F
取对数log10，答案就为log10 t + K，此时很明显log10 t<1，于是我们去除整数部分，就得到了log10 t ，

再用pow（10，log10 t)我们就还原回了t。将t×1000就得到了F
的前四位。 具体实现的时候Log10 F
约等于((1+sqrt(5))/2)^n/sqrt(5)，这里我们把((1-sqrt(5))/2)^n这一项忽略了，

因为当N>=40时，这个数已经小的可以忽略。于是log10 F
就可以化简成log10 1/sqrt(5) + n*log10 (1+sqrt(5))/2。
于是就可以用矩阵快速幂求后四位，用公式求前四位。
代码如下：
#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <queue>
#include <map>
#include <set>
#include <algorithm>

using namespace std;
#define LL __int64
const int mod=1e4;
struct matrix
{
int ma[3][3];
} init, res;
matrix Mult(matrix x, matrix y)
{
int i, j, k;
matrix tmp;
memset(tmp.ma,0,sizeof(tmp.ma));
for(i=0; i<2; i++)
{
for(k=0; k<2; k++)
{
for(j=0; j<2; j++)
{
tmp.ma[i][j]=(tmp.ma[i][j]+x.ma[i][k]*y.ma[k][j])%mod;
}
}
}
return tmp;
}
matrix Pow(matrix x, int k)
{
matrix tmp;
int i, j;
for(i=0; i<2; i++) for(j=0; j<2; j++) tmp.ma[i][j]=(i==j);
while(k)
{
if(k&1) tmp=Mult(tmp,x);
x=Mult(x,x);
k>>=1;
}
return tmp;
}
int main()
{
int k, fib[40], i;
fib[0]=0;
fib[1]=1;
for(i=2; i<40; i++)
{
fib[i]=fib[i-1]+fib[i-2];
}
while(scanf("%d",&k)!=EOF)
{
if(k<40)
{
printf("%d\n",fib[k]);
continue ;
}
init.ma[0][0]=init.ma[0][1]=init.ma[1][0]=1;
init.ma[1][1]=0;
res=Pow(init,k);
double ans;
ans=-0.5*(log10(5.0))+k*log10((sqrt(5.0)+1.0)/2);
ans-=(int)ans;
ans=pow(10.0,ans);
ans*=1000;
printf("%d...",(int)ans);
printf("%4.4d\n",res.ma[0][1]);
}
return 0;
}