快速幂算法

在四则运算的基础上,朴素算法计算 x^a 是将x连乘a次, 算法时间复杂度为O(n), 而快速幂算法可以优化到O(logn), 效率得到极大提升。

例如我们想计算5^{10}的值。 首先有一个显然的事:任何整数都可以用二进制表示, 也就是说,任何整数都可以按 2 的升幂展开,我们把指数 10 展开, 10 = (1010)_210 = 0\times2^0 + 1\times2^1+ 0\times2^2 + 1\times2^3.

那么5^{10} = 5^{0\times2^0 + 1\times2^1+ 0\times2^2 + 1\times2^3} 考虑到x^{ab} = ({x^b})^a 于是:

5^{10} \\= 5^{(0\times2^0 + 1\times2^1+ 0\times2^2 + 1\times2^3)} \\={(5^{2^0})}^0 \times {(5^{2^1})}^1 \times {(5^{2^2})}^0 \times {(5^{2^3})}^1 \\=5^0 \times {25}^1 \times {625}^0 \times {390625}^1

巧的是式子中 5 的平方是 25,25 的平方是 625,625 的平方是 390625, 那么整个式子的计算只需要 1 次乘法和 3 次平方 (5^0{625}^0都是 1 不用参与计算), 原来 10 次计算就被我们简化到 4 次运算。

不失一般性,我们把指数 a 按 2 的升幂展开:

a = a_0\times2^0 + a_1\times2^1 + a_2\times2^2 + a_3\times2^3 + \dots + a_n\times2^n.

指数 a 用这个展开式代替,即

x^a = x^ { a_0\times2^0 + a_1\times2^1 + a_2\times2^2 + a_3\times2^3 + \dots + a_n\times2^n }.

那么,

x^a \\= x^ { a_0\times2^0 + a_1\times2^1 + a_2\times2^2 + a_3\times2^3 + \dots + a_n\times2^n } \\= ({x^{2^0}})^{a_0} \times ({x^{2^1}})^{a_1} \times ({x^{2^2}})^{a_2} \times \dots \times ({x^{2^n}})^{a_n} \\=(x)^{a_0} \times ({x^2})^{a_1} \times ({x^4})^{a_2} \times \dots \times ({x^{2^n}})^{a_n}

不要感觉最后这个式子很复杂,最后一个式子中 a_n 要么是 0 要么是 1。 如果a_n为 0,那么这一项就是1,不用参与计算。 如果a_n为 1,那么这一项就是{x^{2^n}}。 而且,括号里面第一项是 x, 第二项是x^2, 第三项是x^4, 即后一项都是一项的平方,程序运行需要乘法和平方就行了, 这就是快速幂算法基本思想。

有了基本原理,代码实现就很简单了。 式子 (x)^{a_0} \times ({x^2})^{a_1} \times ({x^4})^{a_2} \times \dots \times ({x^{2^n}})^{a_n} a_0 , a_1 , a_2 , a_3 , \dots , a_n就是二进制每一位, C语言取一个数二进制每一位可以使用位运算迭代。 而括号里面底数 x 不断对自身平方更新,即 x = x*x(或 x *= x)。

			int x;//底数
			int a;//指数
			int res = 1;//结果
			scanf("%d%d",&x,&a);
			
			while(a)
			{
				if(a&1)//取二进制最后一位,二进制最后一位为 1 参与运算才有意义
				{
					res*=x;
				}
				else//二进制最后一位为 0,不用计算
				{
					
				}
				
				x *= x;
				//x 自身平方
				a>>=1;
				//指数每一位位数右移一次,这样才能拿到 a 每一个二进制位,相当于 a/=2;
			}
			
			printf("%d",res);
		

为了方便使用,把快速幂封装成一个函数, 而且程序也很容易扩展到底数 x 为实数、a 为任意整数情况。

			double fast_pow(double x,int a)
			{
				if(a<0)
				{
					return 1/fast_pow(x,-a);
				}
				
				double res=1;
				while(a)
				{
					if(a&1)
					{
						res*=x;
					}
			
					x*=x;
					a>>=1;
				}
				return res;
			}
		

快速幂的思想很容易延伸到矩阵幂运算, 只不过我们需要自己定义矩阵的乘法函数 matrix_multiplication()。

			int ** matrix_multiplication(int ** A,int ARow,int AColumn,int **B,int BRow,int BColumn)
			{
				//ARow 是矩阵 A 的行数
				//AColumn 是 A 的列数
			
				//BRow 是矩阵 B 的行数
				//BColumn 是矩阵 B 的列数
			
				//R 是计算结果矩阵
				int RRow = ARow;
				int RColumn = BColumn;
				int **R = new int *[RRow];
				for(int i=0;i<RRow;i++)
				{
					R[i] = new int[RColumn];
				}
			
				for(int i=0;i<RRow;i++)
				{
					for(int j=0;j<RColumn;j++)
					{
						int sum=0;
						for(int k=0;k<AColumn;k++)
						{
							sum+=A[i][k]*B[k][j];
						}
						R[i][j]=sum;
					}
				}
				
				return R;
			}
			
			int ** matrix_fast_pow(int **A, int n, int k)
			{
				//矩阵 A 必为方阵
				//n 是阶数
				//k 是指数
				int **R = new int *[n];
				for (int i = 0; i<n; i++)
				{
					R[i] = new int[n];
				}
				//初始化单位矩阵
				for (int i = 0; i < n; i++)
				{
					for(int j=0;j<n;j++)
					{
						if(i==j)//主对角线上元素
						{
							R[i][j] = 1;
						}
						else
						{
							R[i][j] = 0;
						}
					}
				}
				
				
				while (k)
				{
					if (k & 1)
					{
						R = matrix_multiplication(R, n, n, A, n, n);
					}
			
					A = matrix_multiplication(A, n, n, A, n, n);
					k >>= 1;
				}
				return R;
			}
		

矩阵快速幂还可以计算斐波那契数列。

我们来看看一个特殊的矩阵 A = \begin{bmatrix} 1 & 1\\\\ 1 & 0 \end{bmatrix} , 我们对 A 进行幂运算:

A^0 = \begin{bmatrix} 1 & 0\\\\ 0 & 1 \end{bmatrix} A^1 = \begin{bmatrix} 1 & 1\\\\ 1 & 0 \end{bmatrix} A^2 = \begin{bmatrix} 2 & 1\\\\ 1 & 1 \end{bmatrix} A^3 = \begin{bmatrix} 3 & 2\\\\ 2 & 1 \end{bmatrix} A^4 = \begin{bmatrix} 5 & 3\\\\ 3 & 2 \end{bmatrix} A^5 = \begin{bmatrix} 8 & 5\\\\ 5 & 3 \end{bmatrix} A^6 = \begin{bmatrix} 13 & 8\\\\ 8 & 5 \end{bmatrix} A^7 = \begin{bmatrix} 21 & 13\\\\ 13 & 8 \end{bmatrix} ...

太巧了吧,矩阵 \begin{bmatrix} 1 & 1\\\\ 1 & 0 \end{bmatrix} 的 k 次幂运算得出的新矩阵的第一个元素就是斐波那契数列的第 k项。 这并不是巧合,网上有证明材料。 利用矩阵快速幂是高效计算斐波那契数列的常用方式, 尤其是在项数很大的时候,这种方式比迭代或递归高效得多。

还有一个问题,为什么我们把指数按 2 的升幂展开?为什么不按其他数展开?难道仅仅为了利用位运算?

假如我们把指数 a 按 10 的升幂展开:

x^a \\= x^ {( a_0\times{10}^0 + a_1\times{10}^1 + a_2\times{10}^2 + a_3\times{10}^3 + \dots + a_n\times{10}^n )} \\= {(x^{{10}^0})}^{a_0} \times {(x^{{10}^1})}^{a_1} \times {(x^{{10}^2})}^{a_2} \times {(x^{{10}^3})}^{a_3} \times \dots \times {(x^{{10}^n})}^{a_n} \\={(x)}^{a_0} \times {(x^{10})}^{a_1} \times {(x^{100})}^{a_2} \times {(x^{1000})}^{a_2} \times \dots \times {(x^{10^n})}^{a_n}

那么问题来了,这里a_n的取值为 0、1、2、3、4、5、6、7、8、9, 导致({x^{10^n}})^{a_n}不好计算; 而且底数由 x 变到x^{10}变到x^{100} 再变到x^{1000} ... 实现起来有点麻烦。 这样做不如经典方法简单,虽然这两个问题能借助循环克服:

			double fast_pow_10(double x, int a)
			{
				if (a<0)
				{
					return 1/fast_pow_10(x, -a);
				}
				double res = 1;
				while (a)
				{
					int lastbit = a % 10;
					if (lastbit != 0)
					{
						for (int i = 1; i <= lastbit; i++)
						{
							res *= x;
						}
					}
					
					x = x * x  x * x * x * x * x * x * x * x;
					a /=10;
				}
				return res;
			}
		

最后,我们可以看到快速幂算法执行进行了 n 次循环,即指数的二进制的位数,即指数按 2 的升幂展开式的项数。 假如指数为 a,可以算出循环执行了 \log_2a 次左右。 而朴素算法循环次数就是 a 次,这就是为什么快速幂算法那么高效,尤其在指数很大的时候。