編程之美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]
#include <iostream> 
 
using namespace std; 
 
#define MAXSIDE 10 
struct Matrix{int side; long long a[MAXSIDE][MAXSIDE];}; 
 
// 方陣的乘法 
void MatrixMulti(const Matrix x, const Matrix y, Matrix& z) 
{ 
    if (x.side != y.side) return; 
    z.side = x.side; 
    for (int i=0; i<z.side; i++) 
        for (int j=0; j<z.side; j++) 
        { 
            z.a[i][j] = 0; 
            for (int k=0; k<z.side; k++) 
                z.a[i][j] += x.a[i][k]*y.a[k][j]; 
        } 
} 
 
// 方陣的冪次運算 
void MatrixPower(const Matrix x, Matrix &y, int n) 
{ 
    y.side = x.side; 
    memset(y.a, 0, y.side*y.side*sizeof(long long)); 
    for (int i=0; i<y.side; i++) 
        y.a[i][i] = 1; 
    Matrix tmp = x; 
    // 只需要O(logn)的復雜度就能算出x的n次冪 
    while (n) 
    { 
        // 如果n的最低二進制位為1，則乘上對應的冪次tmp 
        if ((n&1) == 1) 
            MatrixMulti(y, tmp, y); 
        MatrixMulti(tmp, tmp, tmp); 
        // n移位 
        n >>= 1; 
    } 
} 
 
int main() 
{ 
    int n; 
    Matrix tmp; 
    tmp.side = 2; 
    tmp.a[0][0] = 1; tmp.a[0][1] = 1; 
    tmp.a[1][0] = 1; tmp.a[1][1] = 0; 
    while (cin >> n) 
    { 
        Matrix x = tmp, y; 
        MatrixPower(x, y, n-1); 
        long long fn = y.a[0][0]; 
        cout << fn << endl; 
    } 
} 

如果想輸出更高項的值，就要用大整數來表示矩陣的元素。這裡我們可以表示1000位的整數，表達能力是上面算法的50倍，能計算出第4000多個元素，而且計算速度非常快。
[cpp]
#include <iostream> 
using namespace std; 
 
// 大整數類型 
#define MAXLEN 1000 
struct HP {int len, s[MAXLEN];}; 
 
void PrintHP(const HP &x)  
{ 
    for (int i=x.len; i>=1; i--) 
        cout << x.s[i]; 
} 
 
// 字符串轉大整數 
void Str2HP(const char *s, HP &x) 
{ 
    x.len = strlen(s); 
    for (int i=1; i<=x.len; i++) 
        x.s[i] = s[x.len-i] - '0'; 
    if (x.len == 0) 
    { 
        x.len = 1; 
        x.s[1] = 0; 
    } 
} 
 
// 大整數的加法 
void Plus(const HP a, const HP b, HP &c) 
{ 
    int i; c.s[1] = 0; 
    // 大整數a,b的加法操作和結果c的進位操作 
    for (i=1; i<=a.len || i<=b.len || c.s[i]; i++) 
    { 
        if (i <= a.len) c.s[i] += a.s[i]; 
        if (i <= b.len) c.s[i] += b.s[i]; 
        c.s[i+1] = c.s[i]/10; c.s[i] %= 10; 
    } 
    // 退出循環到原因是c.s[i]==0，所以取前一位 
    c.len = i-1;  
    if (c.len == 0) c.len = 1; 
} 
 
// 大整數的減法 
void Subtract(const HP a, const HP b, HP &c) 
{ 
    int i, j; 
    for (i=1,j=0; i<=a.len; i++) 
    { 
        // j表示是否要對高位進行借位 
        c.s[i] = a.s[i] - j; 
        if (i <= b.len) c.s[i] -= b.s[i]; 
        if (c.s[i] < 0)  
        { 
            // 向高位借位，補10 
            j = 1; 
            c.s[i] += 10; 
        } 
        else j = 0; 
    } 
    c.len = a.len; 
    while (c.len > 1 && !c.s[c.len]) c.len--; 
} 
 
// 大整數的比較 
int HPCompare(const HP &x, const HP &y) 
{ 
    if (x.len > y.len) return 1; 
    if (x.len < y.len) return -1; 
    int i = x.len; 
    while (i>1 && (x.s[i]==y.s[i])) i--; 
    return x.s[i] - y.s[i]; 
} 
 
// 大整數的乘法 
void Multi(const HP a, const HP b, HP &c) 
{ 
    int i, j; 
    // 對乘法結果賦初值，以方便之後的+=運算 
    c.len = a.len + b.len; 
    for (i=1; i<=c.len; i++) c.s[i] = 0; 
    for (i=1; i<=a.len; i++) 
        for (j=1; j<=b.len; j++) 
            c.s[i+j-1] += a.s[i]*b.s[j]; 
    // 運算結果進位 
    for (i=1; i<c.len; i++) {c.s[i+1] += c.s[i]/10; c.s[i] %= 10;} 
    // 最高位繼續進位 
    while (c.s[i]) {c.s[i+1] = c.s[i]/10; c.s[i] %= 10; i++;} 
    // 確保最高位不為0 
    while (i>1 && !c.s[i]) i--; 
    c.len = i; 
} 
 
// 大整數的除法 
void Divide(const HP a, const HP b, HP &c, HP &d) 
{ 
    int i, j; 
    // 用余數d存被除數a的前i位數據，用來多次減去除數b，以得到商c 
    d.len = 1; d.s[1] = 0; 
    for (i=a.len; i>0; i--) 
    { 
        if (!(d.len == 1 && d.s[1] == 0)) 
        { 
            // i沒移一位，余數d也移位 
            for (j=d.len; j>0; j--) 
                d.s[j+1] = d.s[j]; 
            d.len++; 
        } 
        d.s[1] = a.s[i]; 
        c.s[i] = 0; 
        // 余數d大於除數b時，才可以進行減操作 
        while ((j=HPCompare(d,b)) >= 0) 
        { 
            Subtract(d, b, d); 
            c.s[i]++; 
            if (j == 0) break; 
        } 
    } 
    c.len = a.len; 
    while (c.len > 1 && c.s[c.len] == 0) 
        c.len--; 
}  
 
 
#define MAXSIDE 3 
struct Matrix{int side; HP a[MAXSIDE][MAXSIDE];}; 
 
// 方陣的乘法 
void MatrixMulti(const Matrix x, const Matrix y, Matrix& z) 
{ 
    if (x.side != y.side) return; 
    HP tmp; 
    z.side = x.side; 
    for (int i=0; i<z.side; i++) 
        for (int j=0; j<z.side; j++) 
        { 
            z.a[i][j].len=1; z.a[i][j].s[1]=0; 
            for (int k=0; k<z.side; k++) 
            { 
                Multi(x.a[i][k], y.a[k][j], tmp); 
                Plus(z.a[i][j], tmp, z.a[i][j]); 
            } 
        } 
} 
 
// 方陣的冪次運算 
void MatrixPower(const Matrix x, Matrix &y, int n) 
{ 
    y.side = x.side; 
    int i, j; 
    for (i=0; i<y.side; i++) 
        for (j=0; j<y.side; j++) 
        { 
            y.a[i][j].len = 1; 
            y.a[i][j].s[1] = 0; 
        } 
    for (i=0; i<y.side; i++) 
    { 
        y.a[i][i].len = 1; 
        y.a[i][i].s[1] = 1; 
    } 
    Matrix tmp = x; 
    // 只需要O(logn)的復雜度就能算出x的n次冪 
    while (n) 
    { 
        // 如果n的最低二進制位為1，則乘上對應的冪次tmp 
        if ((n&1) == 1) 
            MatrixMulti(y, tmp, y); 
        MatrixMulti(tmp, tmp, tmp); 
        // n移位 
        n >>= 1; 
    } 
} 
 
int main() 
{ 
    int n; 
    Matrix tmp; 
    tmp.side = 2; 
    tmp.a[0][0].len=tmp.a[0][1].len=tmp.a[1][0].len=tmp.a[1][1].len=1; 
    tmp.a[0][0].s[1] = 1; tmp.a[0][1].s[1] = 1; 
    tmp.a[1][0].s[1] = 1; tmp.a[1][1].s[1] = 0; 
    while (cin >> n) 
    { 
        Matrix x = tmp, y; 
        MatrixPower(x, y, n-1); 
        HP fn = y.a[0][0]; 
        PrintHP(fn); 
        cout << endl; 
    } 
} 
作者：linyunzju