编程之美2.9----斐波那契数列

问题：

斐波那契数列由如下递推关系式定义：F(0) = 0，F(1)=1，F(n)=F(n-1)+F(n-2) if n>1。

解法：

斐波那契数列是二阶递推数列，所以存在一个2*2的矩阵A，使得：

(Fn, Fn-1) = (Fn-1, Fn-2)*A

求得A=(1   1)

            (1   0)

那么求数列的第n项就是等于求矩阵A的第n-1次幂，计算的速度非常快，时间复杂度为O(logn)。

首先我们用long long 型表示数列中的元素，它只能表示20位的整数，能表示的范围太小，最多第92个元素。

[cpp] view
plaincopy

#include <iostream>  

#include <cassert>  

using namespace std;  

  

const int MAXLENGTH = 10;  

  

struct Matrix   

{  

    unsigned side;  

    long long dat[MAXLENGTH*MAXLENGTH];  

};  

  

// 方阵的乘法   

void MatrixMult(const Matrix a, const Matrix b, Matrix &m)  

{  

    assert(a.side == b.side);  

    m.side = a.side;  

    for (int i=0; i<m.side; ++i)  

        for (int j=0; j<m.side; ++j)  

        {  

            m.dat[i*m.side+j] = 0;  

            for (int k=0; k<m.side; ++k)  

                m.dat[i*m.side+j] += a.dat[i*a.side+k]*b.dat[k*b.side+j];  

        }  

}  

  

long long Fibonaci(unsigned n)  

{  

    if (n==0) return 0;  

    --n;    // 计算矩阵prod的n-1次幂  

    Matrix res; res.side = 2;  

    res.dat[0] = 1; res.dat[1] = 0;   

    res.dat[2] = 0; res.dat[3] = 1;  

    Matrix prod;        prod.side = 2;  

    prod.dat[0] = 1;  prod.dat[1] = 1;  

    prod.dat[2] = 1;  prod.dat[3] = 0;  

    // 只需要O(logn)的复杂度就能算出x的n次幂    

    while (n)  

    {  

        // 如果n的最低二进制位为1，则乘上对应的幂次prod  

        if (n&1) MatrixMult(res, prod, res);  

        MatrixMult(prod, prod, prod);  

        n >>= 1;  

    }  

    return res.dat[0];  

}  

  

int main()  

{  

    long long res = Fibonaci(93);  

    cout << res << endl;  

}  

 

如果想输出更高项的值，就要用大整数来表示矩阵的元素。这里我们可以表示1000位的整数，表达能力是上面算法的50倍，能计算出第4000多个元素，而且计算速度非常快。

[cpp] view
plaincopy

#include <iostream>  

#include <cstring>  

#include <cassert>  

using namespace std;  

  

// 大整数类型  

const int MAXLENGTH = 1000;  

struct HP{int len, s[MAXLENGTH];};  

  

// 字符串转大整数  

void Str2HP(const char* s, HP& x)  

{  

    x.len = strlen(s);  

    for (int i=0; i<x.len; ++i)  

    {  

        assert(s[i]>='0' && s[i]<='9');  

        x.s[x.len-1-i] = s[i]-'0';  

    }  

    if (x.len == 0)  

    {  

        x.len = 1;  

        x.s[0] = 0;  

    }  

}  

  

void PrintHP(const HP& x)  

{  

    for (int i=x.len-1; i>=0; --i)  

        cout << x.s[i];  

}  

  

// 大整数的加法  

void PlusHP(const HP x, const HP y, HP &z)  

{  

    int i;  z.s[0] = 0;  

    // 大整数x,y的加法操作和输出大整数z的进位操作  

    for (i=0; i<x.len || i<y.len; ++i)  

    {  

        if (i<x.len) z.s[i] += x.s[i];  

        if (i<y.len) z.s[i] += y.s[i];  

        z.s[i+1] = z.s[i]/10;  

        z.s[i] %= 10;  

    }  

    // 第i位不会再进位，这里可省去  

    while (z.s[i])  

    {  

        z.s[i+1] = z.s[i]/10;  

        z.s[i] %= 10;  

        ++i;  

    }  

    z.len = i;  

}  

  

// 大整数的减法   

void SubtractHP(const HP x, const HP y, HP &z)  

{  

    int i, j;  

    // j表示是否要对高位进行借位，1表示借1位，0表示不借位  

    for (i=0, j=0; i<x.len; ++i)  

    {  

        z.s[i] = x.s[i] - j;  

        if (i<y.len) z.s[i] -= y.s[i];  

        if (z.s[i]<0)   

        {  

            // 向高位借位，该位补10  

            j=1;   

            z.s[i] += 10;  

        }  

        else  

            j=0;  

    }  

    do --i;  

    while (i>0 && !z.s[i]);  

    z.len = i+1;  

}  

  

// 大整数的比较  

int CompareHP(const HP &x, const HP &y)  

{  

    if (x.len > y.len) return 1;  

    if (x.len < y.len) return -1;  

    int i = x.len-1;  

    while (i>0 && x.s[i]==y.s[i]) --i;  

    return x.s[i]-y.s[i];  

}  

  

// 大整数的乘法  

void MultiHP(const HP x, const HP y, HP &z)  

{  

    int i, j;  

    // 对乘法结果大整数z初值化为0，以方便之后的+=运算  

    z.len = x.len + y.len;  

    for (i=0; i<z.len; ++i) z.s[i] = 0;  

    for (i=0; i<x.len; ++i)  

        for (j=0; j<y.len; ++j)  

            z.s[i+j] += x.s[i]*y.s[j];  

    // 大整数z进位  

    for (i=0; i<z.len-1; ++i) {z.s[i+1] += z.s[i]/10; z.s[i] %= 10;}  

    // 最高位继续进位，这步不会执行可省去  

    while (z.s[i]) {z.s[i+1] = z.s[i]/10; z.s[i] %= 10; i++;}  

    // 直到最高位不为0 ,以确定大整数的长度  

    while (i>0 && !z.s[i]) --i;  

    z.len = i+1;  

}  

  

  

void DivideHP(const HP x, const HP y, HP &u, HP &v)  

{  

    int i, j;  

    v.len = 1; v.s[0] = 0;  

    // u表示x的前i位除y的商  

    // v表示x的前i位除y的余数  

    for (i=x.len-1; i>=0; --i)  

    {  

        // 余数v先向左移一位，再加上x.s[i]  

        if (!(v.len==1 && v.s[0]==0))  

        {  

            for (j=v.len-1; j>=0; --j)  

                v.s[j+1] = v.s[j];  

            ++v.len;  

        }  

        v.s[0] = x.s[i];  

        // 每次循环都计算出商的第i位u.s[i]  

        u.s[i] = 0;  

        while ((j=CompareHP(v, y))>=0)  

        {  

            // 余数v大于等于除数y时，就进行减操作  

            SubtractHP(v, y, v);  

            ++u.s[i];  

            if (j==0) break;  

        }  

    }  

    i = x.len - 1;  

    while (i>0 && !u.s[i]) --i;  

    u.len = i+1;  

}  

  

// 大整数置0  

inline void ZeroHP(HP& x)  

{  

    x.len = 1; x.s[0] = 0;  

}  

  

// 大整数置1  

inline void OneHP(HP& x)  

{  

    x.len = 1; x.s[0] = 1;  

}  

  

// 二维数组（方阵）  

const int MAXSIDE = 2;  

struct Matrix{int side; HP a[MAXSIDE*MAXSIDE];};  

  

// 方阵的乘法  

void MultiMatrix(const Matrix x, const Matrix y, Matrix &z)  

{  

    assert(x.side == y.side);  

    z.side = x.side;  

    HP tmp;  

    for (int i=0; i<z.side; ++i)  

        for (int j=0; j<z.side; ++j)  

        {  

            ZeroHP(z.a[i*z.side+j]);  

            for (int k=0; k<z.side; ++k)  

            {  

                MultiHP(x.a[i*x.side+k], y.a[k*y.side+j], tmp);  

                PlusHP(z.a[i*z.side+j], tmp, z.a[i*z.side+j]);  

            }  

        }  

}  

  

const HP& Fibonaci(int n)  

{  

    HP zero;  

    ZeroHP(zero);  

    if (n==0) return zero;  

    --n;  

    Matrix tmp; tmp.side = 2;  

    OneHP(tmp.a[0]);    OneHP(tmp.a[1]);  

    OneHP(tmp.a[2]);    ZeroHP(tmp.a[3]);  

    Matrix res; res.side = 2;  

    OneHP(res.a[0]);    ZeroHP(res.a[1]);  

    ZeroHP(res.a[2]);   OneHP(res.a[3]);  

    while (n)  

    {  

        if (n&1) MultiMatrix(res, tmp, res);  

        MultiMatrix(tmp, tmp, tmp);  

        n >>= 1;  

    }  

    return res.a[0];  

}  

  

int main()  

{  

    const HP& res = Fibonaci(4000);  

    PrintHP(res);  

    cout << endl;  

}  

转载自：http://blog.csdn.net/linyunzju/article/details/7706896