矩阵快速幂

来源:互联网 发布:tensorflow 中文 编辑:程序博客网 时间:2024/06/11 03:58

快速幂:

求某数的n次方,如A^9,直接A*A*A*A*A*A*A*A*A会很慢,如果不觉得慢,试试求A^999,A^9999吧!

这样考虑:
A^2 = A*A
A^4 = (A^2)*(A^2)
A^8 = (A^4)*(A^4)
A^9 = (A^8)*A
要简单很多,因为A^2,A^4,A^8都在重复利用,只需要算4次乘法(以A^3考虑也需要算4次乘法)。
我们以二分法考虑,任意n次幂有:
若n为偶数, A^n = A^(n/2) * A^(n/2); 
若n为奇数, A^n = A * A^((n-1)/2) * A^((n-1)/2);

这样我们能递归地从A^2计算到A^n,而乘法次数从n降到了log(n)。

十进制数除以2,以2为底的对数,二分等运算在计算机上都能和二进制运算挂钩,

我们将9写成二进制为1001,也就是说A^9=A^8*A^1。

我们从二进制表示的右边(低位)向左计算,每计算一次就把A平方一次后存起来以便重复利用,在当前位数是1时再将这个存起来的数累乘到答案中,这样计算下去就得到我们要的答案。

接下来是代码实现A^n:

while(n>0) {
if(n%2==1) {
ans*=A;
}
A*=A;
n/=2;//相当于n>>1
}


矩阵快速幂:

当我们对许一大堆数字进行一系列运算时,我们通常采用矩阵运算来一次性计算这大堆数字,若把这堆数字看成对事物的属性描述,那么这堆数字组成的矩阵代表这个事物,而事物产生的一系列变化就看做对代表事物矩阵进行的一系列矩阵运算(例如图像处理中对图像平移、旋转、拉伸、放大等操作都是图像的矩阵运算)。

我们都知道斐波拉契数列,并知道相关计算方法,比如迭代着求、递归着求,其实求斐波拉契数列采用矩阵乘法求解更为高效。若把数列的第n项和第n-1项组成矩阵(F(n),F(n-1))T,这里T是转置符号,由斐波拉契数列计算公式:F(n)=F(n-1)+F(n-2),则有:

(F(n),F(n-1))T=((1,1),(1,0))T * (F(n-1),F(n-2))T,

当然递推下去有:

(F(n),F(n-1))T=[((1,1),(1,0))T]^(n-2) * (F(2),F(1))T,

则工作重点就转移到了求矩阵((1,1),(1,0))T的n-2次幂了,而求高次幂就可以用上面的快速幂运算,只是将运算数字变为矩阵,数字乘法变为矩阵乘法了。


矩阵乘法:

矩阵乘法的算法实现有很多种,比如朴素算法(就是我们平常的计算方法,时间复杂度为o(n^3))、strassen算法(采用分治思想,挺复杂的,时间复杂度o(n^(log7)),即o(n^2.81))等。

我们讨论下一个简单的小优化,进行幂运算的矩阵为稀疏矩阵或0的数量很多的矩阵时(其实在很多运用中,用来给原矩阵做变换的矩阵有不少是0很多的矩阵)有不少的0存在,而对应位置的0的乘法运算会反复做多次(跟行数列数相关),其实只要判断到对应位置的0,就将之后的乘法给剪掉就行了。而方便计算,我们将正常的“i行乘j列得出答案ans[ i ] [ j ]”的朴素方法改造为:i行分别乘被乘矩阵的每一行,再将结果累加到答案矩阵的i行各数字对应的位置,这样,当i行中某一位是0就可剪掉后续相关乘法运算了。

实现代码如下:

private static int[][] matrix_mult(int a[][], int b[][]) {

if(a[0].length!=b.length)
return null;

int c[][]=new int[a.length][b[0].length];
for(int i=0;i<a.length;i++) {
for(int j=0;j<a[i].length;j++) {
if(a[i][j]==0)
continue;
for(int k=0;k<b[j].length;k++) {
c[i][k]+=a[i][j]*b[j][k];
}
}
}

return c;
}


接下来,我们来个实际例子:

网易2017年春招笔试题(来自牛客网):

小易拥有一个拥有魔力的手环上面有n个数字(构成一个环),当这个魔力手环每次使用魔力的时候就会发生一种奇特的变化:每个数字会变成自己跟后面一个数字的和(最后一个数字的后面一个数字是第一个),一旦某个位置的数字大于等于100就马上对100取模(比如某个位置变为103,就会自动变为3).现在给出这个魔力手环的构成,请你计算出使用k次魔力之后魔力手环的状态。 

输入描述:
输入数据包括两行:第一行为两个整数n(2 ≤ n ≤ 50)和k(1 ≤ k ≤ 2000000000),以空格分隔第二行为魔力手环初始的n个数,以空格分隔。范围都在0至99.


输出描述:
输出魔力手环使用k次之后的状态,以空格分隔,行末无空格。

输入例子1:
3 21 2 3

输出例子1:
8 9 7

Java代码:

import java.util.*;


public class Main {


public static void main(String[] args) {
// TODO Auto-generated method stub
Scanner sysin = new Scanner(System.in);
int n = sysin.nextInt(), k = sysin.nextInt();
/*手环数字符文列表*/
int[][] list = new int[n][1];
for(int i=0;i<n;i++)
list[i][0]=sysin.nextInt();
sysin.close();

/*使用一次的变换矩阵*/
int[][] use = new int[n][n];
for(int i=0;i<n-1;i++) {
use[i][i]=1;
use[i][i+1]=1;
}
use[n-1][n-1]=1;
use[n-1][0]=1;

/*初始化单位矩阵*/
int[][] tmp = new int[n][n];
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
tmp[i][j]= (i==j?1:0);

/*快速幂运算*/
while(k>0) {
if(k%2==1) {
tmp=matrix_mult(tmp, use);
}
use=matrix_mult(use,use);
k/=2;
}

list=matrix_mult(tmp,list);
for(int i=0;i<list.length-1;i++) {
System.out.print(list[i][0]%100+" ");
}
System.out.println(list[list.length-1][0]%100);
}


/*矩阵相乘(基于快速矩阵幂运算的剪枝)*/
private static int[][] matrix_mult(int a[][], int b[][]) {

if(a[0].length!=b.length)
return null;

int c[][]=new int[a.length][b[0].length];
for(int i=0;i<a.length;i++) {
for(int j=0;j<a[i].length;j++) {
if(a[i][j]==0)
continue;
for(int k=0;k<b[j].length;k++) {
c[i][k]+=a[i][j]*b[j][k]%100;
}
}
}

return c;
}

}


对于100求模应该还有更好的处理方式。