程式師世界 >> 編程語言 >> C語言 >> C++ >> C++入門知識 >> HDU 3117 Fibonacci Numbers(矩陣快速冪+公式)

HDU 3117 Fibonacci Numbers(矩陣快速冪+公式)

編輯：C++入門知識

HDU 3117 Fibonacci Numbers(矩陣快速冪+公式)

題目地址：HDU 3117

對於後四位可以用矩陣快速冪快速求出來，但前四位就沒辦法了。要知道斐波那契數列是有通項公式的，所以只能通過通項公式來求前四位，但公式不能求後四位，因為公式使用浮點數求的，精度顯然不夠，求前四位要用到對數。

通項公式為：

f(n)=1/sqrt(5)(((1+sqrt(5))/2)^n+((1-sqrt(5))/2)^n)

假設F[n]可以表示成 t * 10^k（t是一個小數），那麼對於F[n]取對數log10，答案就為log10 t + K，此時很明顯log10 t<1，於是我們去除整數部分，就得到了log10 t ，
再用pow（10，log10 t)我們就還原回了t。將t×1000就得到了F[n]的前四位。具體實現的時候Log10 F[n]約等於((1+sqrt(5))/2)^n/sqrt(5)，這裡我們把((1-sqrt(5))/2)^n這一項忽略了，
因為當N>=40時，這個數已經小的可以忽略。於是log10 F[n]就可以化簡成log10 1/sqrt(5) + n*log10 (1+sqrt(5))/2。

於是就可以用矩陣快速冪求後四位，用公式求前四位。

代碼如下：

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

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;
}