编程之美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