快速幂算法

递归快速幂

template<typename T>
T qpow(T a, T n)
{
    if (n == 0)
        return T(1);
    else if (n % 2 == 1)
        return qpow(a, n - 1) * a;
    else
    {
        T temp = qpow(a, n / 2);
        return temp * temp;
    }
}

非递归快速幂

$$a^n=a^{b_i{2}^i+...+b_1{2}^1+b_0}=a^{b_i{2}^i}\times ...\times a^{b_1{2}^1}\times a^{b_0}$$
因此用各项$a^{b_i{2}^i}=\{\begin{array}{cc}{a}^{2^i}& ,如果b_i=1\\ 1& ,如果b_i=0\end{array}$的积来计算$a^n$，也就是连续项$a^{2^i}$的乘积，其中跳过了那些二进制位$b_i$为0的项。因为$a^{2^i}=(a^{2^{i-1}}{)}^2$。所以我们可以利用二进制位串从右到左来计算$a^n$。

template<typename T>
T RightLeftBinaryExponentiation(T a, vector<T> b)
// 输入：一个数字a和存放二进制位b_i,...,b_0的数组b(n)
// 		这些位来自于非负整数n的二进制展开式，但是在存储时b_0存放在b[0],b_i存放在b[i]
// 输出：a^n的值
{
    T term = a;
    T product;
    if(b[0] == 1) {
    	product = a;
    }
    else {
    	product = 1;
    }
    for(int i=1; i< b.size(); ++i) {
		term = term * term;
		if(b[i] == 1) {
			product = product * term;
		}
	}
	return product;
}

如果我们现在不输入位串数组，直接输入n，那么就可以得到非递归快速幂的写法。

template<typename T>
T qpow(T a, T n)
{
	T term = a;
    T product = 1;
    while(n) {
	    if(n&1) { // 按位与运算，对应位置上有一个为0结果为0，两个都是1结果为1，n&1==1，说明n的二进制位的最低位为1
	    	product = product * term;
    	}
    	term = term * term;
    	n >>= 1;
	}
	return product;
}

进一步因为我们只想获得最后的值，所以我们可以直接对a进行操作。

template<typename T>
T qpow(T a, T n)
{
    T product = 1;
    while(n) {
	    if(n&1) { // 按位与运算，对应位置上有一个为0结果为0，两个都是1结果为1，n&1==1，说明二进制位的最低位为1
	    	product = product * a;
    	}
    	a = a * a;
    	n >>= 1;
	}
	return product;
}

应用：
用快速幂计算斐波那契数
斐波那契数列：
$$F_n=\{\begin{array}{cc}{F}_{n-1}+F_{n-2}& ,n\ge 2,n\in N^{\ast }\\ 0& ,n=0\\ 1& ,n=1\end{array}$$
首先我们可以构建这样一个递推关系
$$\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right)=\left(\begin{array}{c}{F}_n+F_{n-1}\\ F_n\end{array}\right)=\left(\begin{array}{c}{F}_{n+1}\\ F_n\end{array}\right)$$
所以：
$$\left(\begin{array}{c}{F}_{n+1}\\ F_n\end{array}\right)={\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)}^n\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)$$
令$M=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)$，则$\left(\begin{array}{c}{F}_{n+1}\\ F_n\end{array}\right)=M^n\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)$，那么$\left(\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right)=M^{n-1}\left(\begin{array}{c}{F}_1\\ F_0\end{array}\right)$
我们可以利用快速幂来计算$M^{n-1}$.
我们这里就处理n比较小的情形，如果涉及n较大的情形，要涉及到高精度的问题。

#include<vector>
#include<iostream>
using namespace std;

vector<vector<int>> matrix_multiply(vector<vector<int>>& a, vector<vector<int>>& b) {
    vector<vector<int>> c{{0, 0}, {0, 0}};
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j];
        }
    }
    return c;
}

vector<vector<int>> matrix_pow(vector<vector<int>>& a, int n) {
    vector<vector<int>> ret{{1, 0}, {0, 1}};
    while (n) {
        if (n & 1) {
            ret = matrix_multiply(ret, a);
        }
        n >>= 1;
        a = matrix_multiply(a, a);
    }
    return ret;
}

int fib(int n) {
    if (n < 2) {
        return n;
    }
    vector<vector<int>> q{{1, 1}, {1, 0}};
    vector<vector<int>> res = matrix_pow(q, n - 1);
    return res[0][0];
}

int main() {
    for (int i = 0; i < 10;++i) {
        cout << fib(i) << endl;
    }
    system("pause");
    return 0;
}