矩阵在计算斐波那契数列的运用

斐波那契数列的定义为：f(n) = f(n-1) + f(n-2)，其中f(1)=0,f(2)=1.
按照一般的迭代算法,时间复杂度是O(n)。
还有一种算法叫矩阵算法。
如果用矩阵，理论上会快点，但是实际上会很慢，除非n很大。原因主要在于其系数很大，所以终究是属于学院派算法。
矩阵算式：A(n)=[f(n),f(n-1)]=[1,1,1,0]*[f(n-1),f(n-2)];
通过递归，可得A(n) = [1,1,1,0]^(n-2) * [1,0];
现在，重点就是计算矩阵M(n) = [1,1,1,0]^n的值。
这个算式时间复杂度为O(logn);
所以总的时间复杂度为O(logn);
空间复杂度为O(1).
目前没看到有比这个更加快的计算方式，用公式的不算。

其实通过一般化，对类型f(n)=af(n-1)+bf(n-2)的方程，都可以用矩阵算法，
A(n)=[f(n),f(n-1)]=[a,b,1,0]*[f(n-1),f(n-2)];
从而将时间复杂度降低到O(logn);

先给出代码实现，在进行算法分析。

#ifndef _HPP_MATRIX_
#define _HPP_MATRIX_
#include <memory>

int g_arr[64] = { 0 };
// 矩阵乘法
template<typename T1, typename T2, typename T3 >
bool Matrixmulti(
    T1  *Lhs,
    int nlCol,
    int nlRow,
    T2  *Rhs,
    int nrCol,
    int nrRow,
    T3  *Result )
{
    if ( nlRow != nrCol
         && Lhs == nullptr
         && Rhs == nullptr
         && Result == nullptr )
    {
        return false;
    }

    for ( int i = 0; i < nlCol; i++ )
        for ( int j = 0; j < nrRow; j++ )
            for ( int k = 0; k < nrCol; k++ )
            {
                Result[i * nrRow + j] 
                    += Lhs[i * nlRow + k] * Rhs[k * nrRow + j];
            }

    return true;
}

// 矩阵平方
template< typename T>
bool MatrixSquare( T *Src, int nOrder, T *Result )
{
    return Matrixmulti(
        Src,
        nOrder,
        nOrder,
        Src,
        nOrder,
        nOrder,
        Result
        );
}

// 矩阵幂，采用递归法
template< typename T >
bool Matrices( T *Src, int nOrder, int nPower, T *Result )
{
    if ( nPower == 2 )
    {
        return MatrixSquare( Src, nOrder, Result );
    }
    if ( nPower == 1 )
    {
        memcpy( Result, Src, nOrder * nOrder );
        return true;
    }

    if ( nPower <= 0 )
    {
        return false;
    }

    T* pTemp = new T[nOrder * nOrder];
    
    int nLen = sizeof( T ) * nOrder * nOrder;
    memset( pTemp, 0, nLen );
    if ( 0 == nPower % 2 )
    {
        if ( g_arr[nPower/2] )
        {
            memcpy( pTemp, (void*)g_arr[nPower / 2], nLen );
        }
        else
        {
            bool bRet = Matrices( Src, nOrder, nPower / 2, pTemp );
            if ( !bRet )
            {
                return false;
            }
            g_arr[nPower / 2] = (int)pTemp;
        }     
        return MatrixSquare( pTemp, nOrder, Result );
    }
    else
    {
        if ( g_arr[nPower - 1] )
        {
            memcpy( pTemp, (void*)g_arr[nPower - 1], nLen );
        }
        else
        {
            bool bRet = Matrices( Src, nOrder, nPower - 1, pTemp );
            if ( !bRet )
            {
                return false;
            }
            g_arr[nPower - 1] = (int)pTemp;
        }
        return Matrixmulti(
            pTemp,
            nOrder,
            nOrder,
            Src,
            nOrder,
            nOrder,
            Result );
    }
    return true;
}

// 计算斐波那契数的第n项大小。
uint64_t Fibonacci( int n )
{
    typedef uint64_t Value_type;
    Value_type m[2][2] = { 1, 1, 1, 0 };
    Value_type init[2][1] = { 1, 0 };
    Value_type *result = new Value_type[2];
    Value_type temp[2][2] = {0};
    if ( n < 0 || n > 93 )
    {
        return -1;
    }
    if ( n < 2 )
    {
        return init[1-n][0];
    }
    memset( result, 0, sizeof(Value_type)* 2 );
#ifdef Recursion // 采用递归法
    bool bRet = Matrices( (Value_type*)m, 2, n - 1, (Value_type*)temp );
    if ( !bRet )
    {
        return -1;
    }
    bRet = Matrixmulti( (Value_type*)temp, 2, 2, (Value_type*)init, 2, 1, (Value_type*)result );
    if ( !bRet )
    {
        return -1;
    }

    return *result;

#else // 采用迭代法

    g_arr[0] = (int)m;
    n--;

    for ( int i = 1; (1 << i) <= n; i++ )
    {
        if ( !g_arr[i] )
        {
            Value_type* pTemp = new Value_type[4];
            memset( pTemp, 0, sizeof(Value_type)* 4 );

            MatrixSquare( (Value_type*)g_arr[i - 1], 2, pTemp );
            g_arr[i] = (int)pTemp;
        }
    }
    
    bool begin = false;
    for ( int i = 0; (1<<i) <= n; i++ )
    {
        if ( (1 << i) & n ) //第i位是否为0
        {
            if ( begin == false )
            {
                begin = true;
                memcpy( temp, (void*)g_arr[i], sizeof(Value_type) * 4 );
            }
            else
            {
                Value_type temp1[2][2] = { 0 }; 
                Matrixmulti( (Value_type*)temp, 2, 2, (Value_type*)g_arr[i], 2, 2, (Value_type*)temp1 );
                memcpy( temp, temp1, sizeof(Value_type) * 4 );
            }
        }
    }

    int bRet = Matrixmulti( (Value_type*)temp, 2, 2, (Value_type*)init, 2, 1, (Value_type*)result );
    if ( !bRet )
    {
        return -1;
    }

    return *result;
#endif
}
#endif

算法分析：
Fibonacci函数内含两种实现，一种是递归法，另一种是迭代法。实际测试迭代法比较高效。这里有一个缓存矩阵2的n次幂的计算结果，记在全局变量g_arr。首先对参数n进行预处理，将n的二进制形式的最高bit上的1所在位置i算出来，计算矩阵[1,i]次幂，结果保存到g_arr。比如给出33=b100001 那么i=5(从0计数)。对于int类型的n来说，共有32位，那么最多要计算32次幂，计算量为32*Matrixmulti,其中Matrixmulti计算量固定。所以预处理的时间复杂度为O(logn).
后面一个循环主要是用n中每一bit为1的位置i从g_arr取出，并相乘。这个的时间复杂度不必预处理高，因为其最坏的情况就是每一bit都为1，此时跟预处理复杂度一样。所以总的时间复杂度为O(logn);

这意味着，如果实现好大数相乘的情况下，当n很大的时候，可以达到很高的效率。

测试代码：

DWORD nt = GetTickCount();
    for ( int i = 0; i < 94; i++ )
    {
        cout << Fibonacci( i ) << '\t';
    }
    cout << endl << "运行时间  " << GetTickCount() - nt << endl;

计算结果

0 1 1 2 3 5 8 13 21 34

55 89 144 233 377 610 987 1597 2584 4181

6765 10946 17711 28657 46368 75025 121393 196418 317811 514229

832040 1346269 2178309 3524578 5702887 9227465 14930352 24157817

39088169 63245986 102334155 165580141 267914296

433494437 701408733 1134903170 1836311903 2971215073

4807526976 7778742049 12586269025 20365011074 32951280099

53316291173 86267571272 139583862445 225851433717 365435296162

591286729879 956722026041 1548008755920 2504730781961 4052739537881

6557470319842 10610209857723 17167680177565 27777890035288 44945570212853

72723460248141 117669030460994 190392490709135 308061521170129 498454011879264

806515533049393 1304969544928657 2111485077978050 3416454622906707

5527939700884757 8944394323791464 14472334024676221

23416728348467685 37889062373143906 61305790721611591 99194853094755497
160500643816367088 259695496911122585 420196140727489673
679891637638612258 1100087778366101931 1779979416004714189

2880067194370816120 4660046610375530309 7540113804746346429 12200160415121876738

运行时间 140

结果说明，在64位下，最高能计算前93位斐波数。这个数其实很小，不足以体现出优势。