斐波那契数列的分治法计算

           《编程之美》一书中讲述斐波那契数列的问题，之前大学本科的时候就接触这个问题，那时候开始就知道使用递归来计算，可是一直没有考虑过改进下该算法。。。囧~~菜   直到看到这本书，发现原来之前好多问题都可以优化，斐波那契就是其中之一，其中书本中讲述了三种方法、：

1. 第一种就是对平时的递归算法进行优化，增加了数组专门记录每个子问题的解，实际上是动态规划的思想；
2. 第二种利用数学中通项公式，利用F0=0，F1=1，Fn=F(n-1)+F(n-2)（n>=2，n∈N*）递推公式得到F（n）的特征方程，然后直接利用公式计算。。。好bug哇~ 不过存在无理数，数据的精度就无法保证了~~
3. 利用分治策略，将问题转换为矩阵乘幂的问题。

  其余两种大家可以去网上找答案脑补下，我这里关键阐述下第三种哈~

基本的原理是运用矩阵：

(Fn Fn-1) = (Fn-1 Fn-2)*A，求解得到

A = { {1,1}, {1,0}}

然后通过递推式

(Fn Fn-1) = (Fn-1 Fn-2)*A = (Fn-2 Fn-3)*A^2 = .... = (F1 F0)*A^(n-1)

然后只要计算A^(n-1)，再与矩阵(F1 F0)相乘，就可以的得到 Fn的值，我非常佩服作者的是，能把重复计算转化为对A的n-1乘法计算，因为计算A的n-1乘法有快速相乘的方法，比如计算m的10000次方，其实最少的计算次数，等于10000的最高比特位的位置与零的个数，即14~28次乘法运算，并不需要10000次乘法。

对于n有 n = ak*2^k + ak-1*2^k-1 + ... + a1*2 + a0，其中ai = 0 或1 ，i = 0,1,2... k ，将幂数变为二进制形式表示，此时我们只需要进行logn次的运算(以2为底)。 也即 快速指数相乘法。

下面是自己实现的代码，可能写的比较挫~ 关键是matrix矩阵的写法哈~

#include <iostream>

using namespace std;

//
struct Matrix 
{
    Matrix(int n)
    {
        if (n >= 1)
        {
            nLength = n;
            data = new long long*[nLength];
            for (int i = 0; i!= n; i++)
            {
                data[i] = new long long[nLength];
            }
        }
        else
        {
            data = NULL;
        }
    }
    Matrix()
    {
        data = NULL;
    }
    // 拷贝构造函数
    Matrix(const Matrix &m)
    {
        int nLen = m.GetMatrixLength();
        if (nLen >= 1)
        {
            nLength = nLen;
            data = new long long*[nLength];
            for (int i = 0; i!= nLen; i++)
            {
                data[i] = new long long[nLength];
            }
            //
            for (int j = 0; j!=nLen; j++)
            {
                for (int k = 0; k != nLen; k++)
                {
                    data[j][k] = m.GetMatrixValue(j,k);
                }
            }
        }
        else
        {
            data = NULL;
        }
    }

    // 重载赋值操作符
    // 错误都在这里。。。注意因为涉及到指针，一定要深拷贝
    Matrix& operator = (const Matrix& m)
    {
        // 如果是自己 直接返回即可
        if (this == &m)
        {
            return *this;
        }
        // 如果当前对象已有空间也需要释放掉 重新分配
        if (data)
        {
            for (int i = 0; i!=nLength ; i++)
            {
                if (data[i])
                {
                    delete[] data[i];
                }
            }
            delete[] data;
            data = NULL;
        }
        //
        if (data == NULL)
        {
            nLength = m.GetMatrixLength();
            data = new long long*[nLength];
            for (int i = 0; i!= nLength; i++)
            {
                data[i] = new long long[nLength];
            }
        }
        for (int i=0; i!= nLength; i++)
        {
            for (int j = 0; j != nLength; j++)
            {
                this->data[i][j] = m.GetMatrixValue(i,j);
            }
        }

        // 方法返回时 调用析构 函数  释放空间
        return *this;
    }


    // 释放所有分配的空间
    ~Matrix()
    {
        if (data)
        {
            for (int i = 0; i!=nLength ; i++)
            {
                if (data[i])
                {
                    delete[] data[i];
                }
            }
            delete[] data;
            data = NULL;
        }
    }

    void setIdentity()
    {
        for (int i = 0; i != nLength; i++)
        {
            for (int j = 0; j!= nLength; j++)
            {
                data[i][j] =( (i ==j) ? 1 : 0 );
            }
        }
    }

    // 重载乘等操作符 一定要记得深拷贝哦 亲
    Matrix& operator *= (const Matrix& m)
    {
        if (this->nLength == m.nLength)
        {
            Matrix  tmp(m.nLength);
            for (int i = 0; i!= nLength; i++)
            {
                for (int j = 0; j!= nLength; j++)
                {
                    // 计算每个位置的值
                    long long nValue = 0;
                    for (int p = 0; p != nLength; p++)
                    {
                        nValue += data[i][p] * (m.GetMatrixValue(p,j));
                    }
                    tmp.SetMatrixValue(i,j,nValue);
                }
            }
            //
            *this = tmp;
            // 返回时会调用析构 释放tmp对象
        }
        return *this;
    }

    //
    void SetMatrixValue(const int i,const int j, const long long nValue)
    {
        data[i][j] = nValue;
    }    long long  GetMatrixValue(const int i, const int j) const 
    {
        long long n = data[i][j];
        return n;
    }
    unsigned int GetMatrixLength() const    {
        return nLength;
    }

private:
    unsigned int nLength;
    long long **data;
};




void GetMatrixPow(const Matrix &m, int index,Matrix &result)
{
    // 先单位化结果矩阵
    result.setIdentity();
//    Matrix  tmp(m.GetMatrixLength());        //
//    tmp = m;
    Matrix tmp = m;    // C语言风格的赋值 会调用拷贝构造函数

    for (;index;index >>=1)
    {
        if (index & 1)
        {
            result *= tmp;
        }
        tmp *= tmp;
    }

    // 该函数结束 调用matrix析构函数 释放tmp所占空间
}

long long Fibonacci(int n)
{
    if (n == 0)
    {
        return 0;
    }
    else if (n == 1)
    {
        return 1;
    }

    long long nValue = 0;
    Matrix A(2);
    A.SetMatrixValue(0,0,1);
    A.SetMatrixValue(0,1,1);
    A.SetMatrixValue(1,0,1);
    A.SetMatrixValue(1,1,0);
    // 结果矩阵一定还是2*2的
    Matrix result(2);
    GetMatrixPow(A , n-1,result);

    nValue = result.GetMatrixValue(1,0);//F0* result.GetMatrixValue(0,0) + F1* result.GetMatrixValue(1,0);
    
    return nValue;
}


int main()
{
    int  N;
    cout<<"输入你要计算的斐波那契数列的项(目前可以计算到93~)：";
    cin>>N;

    for (int i = 0; i <= N; i++)
    {
        long long nResult = Fibonacci(i);
        cout<<nResult<<endl;
    }

    return 0;
}

程序的运行结果如下所示：

当然这里矩阵每个元素我用了long long 64位，经测试可以计算到93~ 需要改进的话，需要选择合适的变量类型。