strassen算法-现实2次幂,偶数奇数

合集下载

strassen算法实现矩阵相乘,当输入为偶数,奇数,任意数的实现

strassen算法实现矩阵相乘,当输入为偶数,奇数,任意数的实现

strassen算法实现矩阵相乘,当输⼊为偶数,奇数,任意数的实现1.使⽤随机数⽣成两个矩阵,并实现输出⽅法。

public void makeMatrix(int[][] matrixA, int[][] matrixB, int length){ //⽣成矩阵Random random=new Random();for(int i=0;i<length;i++){for (int j=0;j<length;j++){matrixA[i][j]=random.nextInt(5);matrixB[i][j]=random.nextInt(5);}}}public void printMatrix(int[][] matrixA,int length){ //输出for(int i=0;i<length;i++){for (int j=0;j<length;j++){System.out.print(matrixA[i][j]+" ");if((j+1)%length==0)System.out.println();}}}2.使⽤Strassen算法需要涉及到矩阵的加减。

所以先准备好⽅法。

public void add(int[][] matrixA,int[][] matrixB,int[][] matrixC,int length){for(int i=0;i<length;i++) {for (int j = 0; j < length; j++) {matrixC[i][j]= matrixA[i][j]+ matrixB[i][j];}}}public void jian(int[][] matrixA,int[][] matrixB,int[][] matrixC,int length){for(int i=0;i<length;i++) {for (int j = 0; j < length; j++) {matrixC[i][j]= matrixA[i][j] - matrixB[i][j];}}}public void cheng(int[][] matrixA,int[][] matrixB,int[][] matrixC,int length){for(int i=0;i<length;i++) {for (int j = 0; j < length; j++) {matrixC[i][j]=0;for(int k=0;k<length;k++){matrixC[i][j] = matrixC[i][j]+ matrixA[i][k] * matrixB[k][j];}}}}3.当矩阵的阶数为 2的k次⽅//阶数为 2 的 K 次⽅的时候 Strassen算法public void strassen(int[][] matrixA,int[][] matrixB,int[][] matrixC,int N){int newsize=N/2;if(N==2){cheng(matrixA,matrixB,matrixC,N);return;}int[][] A11=new int[newsize][newsize];int[][] A12=new int[newsize][newsize];int[][] A21=new int[newsize][newsize];int[][] A22=new int[newsize][newsize];int[][] B11=new int[newsize][newsize];int[][] B12=new int[newsize][newsize];int[][] B21=new int[newsize][newsize];int[][] B22=new int[newsize][newsize];int[][] C11=new int[newsize][newsize];int[][] C12=new int[newsize][newsize];int[][] C21=new int[newsize][newsize];int[][] C22=new int[newsize][newsize];int[][] M1=new int[newsize][newsize];int[][] M2=new int[newsize][newsize];int[][] M3=new int[newsize][newsize];int[][] M4=new int[newsize][newsize];int[][] M5=new int[newsize][newsize];int[][] M6=new int[newsize][newsize];int[][] M7=new int[newsize][newsize];int[][] Aresult=new int[newsize][newsize];int[][] Bresult=new int[newsize][newsize];//分别给 A11 A12 A21 A22赋值for(int i=0;i<N/2;i++){for(int j=0;j<N/2;j++){A11[i][j]=matrixA[i][j];A12[i][j]=matrixA[i][j+N/2];A21[i][j]=matrixA[i+N/2][j];A22[i][j]=matrixA[i+N/2][j+N/2];B11[i][j]=matrixB[i][j];B12[i][j]=matrixB[i][j+N/2];B21[i][j]=matrixB[i+N/2][j];B22[i][j]=matrixB[i+N/2][j+N/2];}}//计算M1 到M7add(A11,A22,Aresult,newsize);add(B11,B22,Bresult,newsize);strassen(Aresult,Bresult,M1,newsize);//M2add(A21,A22,Aresult,newsize);strassen(Aresult,B11,M2,newsize);//M3jian(B12,B22,Bresult,newsize);strassen(A11,Bresult,M3,newsize);//M4jian(B21,B11,Bresult,newsize);strassen(A22,Bresult,M4,newsize);//M5add(A11,A12,Aresult,newsize);strassen(Aresult,B22,M5,newsize);//M6jian(A21,A11,Aresult,newsize);add(B11,B12,Bresult,newsize);strassen(Aresult,Bresult,M6,newsize);//M7jian(A12,A22,Aresult,newsize);add(B21,B22,Bresult,newsize);strassen(Aresult,Bresult,M7,newsize);//C11add(M1,M4,Aresult,newsize);jian(M5,M7,Bresult,newsize);jian(Aresult,Bresult,C11,newsize);//C12add(M3,M5,C12,newsize);//C21add(M2,M4,C21,newsize);//C22add(M1,M3,Aresult,newsize);jian(M2,M6,Bresult,newsize);jian(Aresult,Bresult,C22,newsize);//把C的值填充for(int i=0;i<N/2;i++){for(int j=0;j<N/2;j++){matrixC[i][j]=C11[i][j];matrixC[i][j+N/2]=C12[i][j];matrixC[i+N/2][j]=C21[i][j];matrixC[i+N/2][j+N/2]=C22[i][j];}}}4.当阶数为偶数的时候假设阶数为 n ,所以 n=m*2^k 。

常用算法大全-分而治之算法

常用算法大全-分而治之算法

君主和殖民者们所成功运用的分而治之策略也可以运用到高效率的计算机算法的设计过程中。

本章将首先介绍怎样在算法设计领域应用这一古老的策略,然后将利用这一策略解决如下问题:最小最大问题、矩阵乘法、残缺棋盘、排序、选择和计算一个几何问题——找出二维空间中距离最近的两个点。

本章给出了用来分析分而治之算法复杂性的数学方法,并通过推导最小最大问题和排序问题的复杂性下限来证明分而治之算法对于求解这两种问题是最优的(因为算法的复杂性与下限一致)。

2.1 算法思想分而治之方法与软件设计的模块化方法非常相似。

为了解决一个大的问题,可以:1) 把它分成两个或多个更小的问题; 2) 分别解决每个小问题; 3) 把各小问题的解答组合起来,即可得到原问题的解答。

小问题通常与原问题相似,可以递归地使用分而治之策略来解决。

例2-1 [找出伪币] 给你一个装有1 6个硬币的袋子。

1 6个硬币中有一个是伪造的,并且那个伪造的硬币比真的硬币要轻一些。

你的任务是找出这个伪造的硬币。

为了帮助你完成这一任务,将提供一台可用来比较两组硬币重量的仪器,利用这台仪器,可以知道两组硬币的重量是否相同。

比较硬币1与硬币2的重量。

假如硬币1比硬币2轻,则硬币1是伪造的;假如硬币2比硬币1轻,则硬币2是伪造的。

这样就完成了任务。

假如两硬币重量相等,则比较硬币3和硬币4。

同样,假如有一个硬币轻一些,则寻找伪币的任务完成。

假如两硬币重量相等,则继续比较硬币5和硬币6。

按照这种方式,可以最多通过8次比较来判断伪币的存在并找出这一伪币。

另外一种方法就是利用分而治之方法。

假如把1 6硬币的例子看成一个大的问题。

第一步,把这一问题分成两个小问题。

随机选择8个硬币作为第一组称为A组,剩下的8个硬币作为第二组称为B组。

这样,就把1 6个硬币的问题分成两个8硬币的问题来解决。

第二步,判断A和B组中是否有伪币。

可以利用仪器来比较A组硬币和B组硬币的重量。

假如两组硬币重量相等,则可以判断伪币不存在。

算法设计与分析课后习题

算法设计与分析课后习题

第一章1. 算法分析题算法分析题1-1 求下列函数的渐进表达式(1). 3n^2 + 10n < 3n^2 + 10n^2 = 13n^2 = O(n^2)(2). n^2 / 10 + 2^n当n>5是,n^2 〈2 ^n所以,当n >= 1时,n^2/10 〈2 ^n故:n^2/10 + 2^n < 2 ^n + 2^n = 2*2^n = O(2^n)(3). 21 + 1/n < 21 + 1 = 22 = O(1)(4). log(n^3)=3log(n)=O(log(n))(5). 10log(3^n)= (10log3)n = O(n)算法分析题1—6(1)因为:f(n)=log(n^2) = 2log(n); g(n) = log(n) + 5所以:f(n)=Θ(log(n)+5)=Θ(g(n))(2)因为:log(n) 〈√n ; f(n)= 2log(n);g(n)=√n所以:f(n)= O(g(n))(3)因为:log(n)< n;f(n) = n;g(n) = log(n^2) = 2log(n)所以;f(n)= Ω(g(n))(4)因为:f(n) = nlogn +n; g(n) = logn所以:f(n) =Ω(g(n))(5)因为: f(n)= 10;g(n) = log(10)所以:f(n)=Θ(g(n))(6)因为: f(n)=log^2(n);g(n)= log(n)所以:f(n) ==Ω(g(n))(7)因为: f(n) = 2^n < 100*2^n; g(n)=100n^2; 2^n > n ^2所以:f(n)= Ω(g(n))(8)因为:f(n)= 2^n; g(n)= 3 ^n;2 ^n 〈3 ^n所以:f(n)= O(g(n))习题1-9 证明:如果一个算法在平均情况下的计算时间复杂性为Θ(f(n)),该算法在最坏情况下所需的计算时间为Ω(f(n)).分析与解答:因此,Tmax(N) = Ω(Tavg(N)) = Ω(Θ(f(n)))=Ω(f(n))。

【论文】Solovay-Strassen素数判定算法

【论文】Solovay-Strassen素数判定算法

1前言网络信息安全的研究主要有两个研究方向,一个是网络协议的分析,通过分析网络协议检测出网络协议中存在的漏洞,然后有效的修改这些网络协议,弥补这些漏洞。

另一个就是通过对信息进行加密,把需要通过网络传输的信息加密后再进行传输和把机密的信息加密后再进行存储。

近年来,随着计算机运算速度的不断提升利黑客攻击手段的不断改进,保密通信和信息安全越来越受到人们的重视,直接促进了密码学的发展和应用。

在现代密码系统中,大素数的判别和寻找对一些加密系统起着关键作用。

很多公钥密码算法都用到素数,特别是160位(二进制)以上的大素数。

熟知,RSA 的公共模数就是两个512(近年来又有了对1024,乃至2048)位以上的大素数的乘F上离散对数的公钥密码体制,其中素数p要求在512积;又如,基于有限域p位以上,且p-1有一个大素因子q(160比特以上);再如,基于椭圆曲线离散对数的公钥密码体制(ECC)中,一般要求基点的阶是一个160位以上的大整数,且是一个素数。

由此可见大数的素性判断的重要性。

但当前既没有一个有效的确定性素数检测算法,也没有一个高效确定性素数产生算法。

当前使用的高效的素数判定算法都是概率算法,速度虽快,但可能会错判误判,即把合数判定成素数,这将会严重影响系统的可靠性以及安全性。

而RSA 算法是目前最优秀的公钥解决方案之一, 其安全性建立在大整数分解为两个素数之积的困难性假设基础之上。

因此, RSA 算法的核心问题是要解决通过何种方式能快速的找到两个大的随机素数。

这样既有利于提高RSA 加密的安全性, 又有利于提高加密效率。

因此,开展在多项式时间内判断一个正整数是否是素数的确定性算法研究,以及快速素数产生算法的研究是非常必要和急需的,这对于促进公钥密码体制的应用以及增强保密通信的可靠性以及安全性都有重要意义。

2 算法概述素数就是一个除了1和它自身以外不能被其它数整除的数。

关于素数的一个基本问题是如何有效地确定一个数是否是一个素数,即素性测试问题。

矩阵乘法之 strassen算法

矩阵乘法之 strassen算法

矩阵乘法之strassen 算法一般情况下矩阵乘法需要三个for循环,时间复杂度为O(n^3),现在我们将矩阵分块如图:( 来自MIT算法导论)一般算法需要八次乘法r = a * e + b * g ;s = a * f + b * h ;t = c * e + d * g;u = c * f + d * h;strassen将其变成7次乘法,因为大家都知道乘法比加减法消耗更多,所有时间复杂更高!strassen的处理是:令:p1 = a * ( f - h )p2 = ( a + b ) * hp3 = ( c +d ) * ep4 = d * ( g - e )p5 = ( a + d ) * ( e + h )p6 = ( b - d ) * ( g + h )p7 = ( a - c ) * ( e + f )那么我们可以知道:r = p5 + p4 + p6 - p2s = p1 + p2t = p3 + p4u = p5 + p1 - p3 - p7我们可以看到上面只有7次乘法和多次加减法,最终达到降低复杂度为O( n^lg7 ) ~= O( n^2.81 );代码实现如下:[cpp] view plaincopyprint?// strassen 算法:将矩阵相乘的复杂度降到O(n^lg7) ~= O(n^2.81)// 原理是将8次乘法减少到7次的处理// 现在理论上的最好的算法是O(n^2,367),仅仅是理论上的而已////// 下面的代码仅仅是简单的实例而已,不必较真哦,呵呵~// 下面的空间可以优化的,此处就不麻烦了~#include <stdio.h>#define N 10//matrix + matrixvoid plus( int t[N/2][N/2], int r[N/2][N/2], int s[N/2][N/2] ) {int i, j;for( i = 0; i < N / 2; i++ ){for( j = 0; j < N / 2; j++ ){t[i][j] = r[i][j] + s[i][j];}}}//matrix - matrixvoid minus( int t[N/2][N/2], int r[N/2][N/2], int s[N/2][N/2] ) {int i, j;for( i = 0; i < N / 2; i++ ){for( j = 0; j < N / 2; j++ ){t[i][j] = r[i][j] - s[i][j];}}}//matrix * matrixvoid mul( int t[N/2][N/2], int r[N/2][N/2], int s[N/2][N/2] ) {int i, j, k;for( i = 0; i < N / 2; i++ ){for( j = 0; j < N / 2; j++ ){t[i][j] = 0;for( k = 0; k < N / 2; k++ ){t[i][j] += r[i][k] * s[k][j];}}}}int main(){int i, j, k;int mat[N][N];int m1[N][N];int m2[N][N];int a[N/2][N/2],b[N/2][N/2],c[N/2][N/2],d[N/2][N/2];int e[N/2][N/2],f[N/2][N/2],g[N/2][N/2],h[N/2][N/2];int p1[N/2][N/2],p2[N/2][N/2],p3[N/2][N/2],p4[N/2][N/2];int p5[N/2][N/2],p6[N/2][N/2],p7[N/2][N/2];int r[N/2][N/2], s[N/2][N/2], t[N/2][N/2], u[N/2][N/2], t1[N/2][N/2], t2[N/2][N/2];printf("\nInput the first matrix...:\n");for( i = 0; i < N; i++ ){for( j = 0; j < N; j++ ){scanf("%d", &m1[i][j]);}}printf("\nInput the second matrix...:\n");for( i = 0; i < N; i++ ){for( j = 0; j < N; j++ ){scanf("%d", &m2[i][j]);}}// a b c d e f g hfor( i = 0; i < N / 2; i++ ){for( j = 0; j < N / 2; j++ ){a[i][j] = m1[i][j];b[i][j] = m1[i][j + N / 2];c[i][j] = m1[i + N / 2][j];d[i][j] = m1[i + N / 2][j + N / 2];e[i][j] = m2[i][j];f[i][j] = m2[i][j + N / 2];g[i][j] = m2[i + N / 2][j];h[i][j] = m2[i + N / 2][j + N / 2];}}//p1minus( r, f, h );mul( p1, a, r );//p2plus( r, a, b );mul( p2, r, h );//p3plus( r, c, d );mul( p3, r, e );//p4minus( r, g, e );mul( p4, d, r );//p5plus( r, a, d );plus( s, e, f );mul( p5, r, s );//p6minus( r, b, d );plus( s, g, h );mul( p6, r, s );//p7minus( r, a, c );plus( s, e, f );mul( p7, r, s );//r = p5 + p4 - p2 + p6plus( t1, p5, p4 );minus( t2, t1, p2 );plus( r, t2, p6 );//s = p1 + p2plus( s, p1, p2 );//t = p3 + p4plus( t, p3, p4 );//u = p5 + p1 - p3 - p7 = p5 + p1 - ( p3 + p7 )plus( t1, p5, p1 );plus( t2, p3, p7 );minus( u, t1, t2 );for( i = 0; i < N / 2; i++ ){for( j = 0; j < N / 2; j++ ){mat[i][j] = r[i][j];mat[i][j + N / 2] = s[i][j];mat[i + N / 2][j] = t[i][j];mat[i + N / 2][j + N / 2] = u[i][j];}}printf("\n下面是strassen算法处理结果:\n"); for( i = 0; i < N; i++ ){for( j = 0; j < N; j++ ){printf("%d ", mat[i][j]);}printf("\n");}//下面是朴素算法处理printf("\n下面是朴素算法处理结果:\n");for( i = 0; i < N; i++ ){for( j = 0; j < N; j++ ){mat[i][j] = 0;for( k = 0; k < N; k++ ){mat[i][j] += m1[i][j] * m2[i][j];}}}for( i = 0; i < N; i++ ){for( j = 0; j < N; j++ ){printf("%d ", mat[i][j]);}printf("\n");}return 0;}现在最好的计算矩阵乘法的复杂度是O( n^2.376 ),不过只是理论上的结果。

Strassen算法及其python实现

Strassen算法及其python实现

Strassen 算法及其python 实现题⽬描述请编程实现矩阵乘法,并考虑当矩阵规模较⼤时的优化⽅法。

思路分析根据wikipedia上的介绍:两个矩阵的乘法仅当第⼀个矩阵B的列数和另⼀个矩阵A的⾏数相等时才能定义。

如A是m×n矩阵和B是n×p矩阵,它们的乘积AB是⼀个m×p矩阵,它的⼀个元素其中 1 ≤ i ≤ m, 1 ≤ j ≤ p。

值得⼀提的是,矩阵乘法满⾜结合律和分配率,但并不满⾜交换律,如下图所⽰的这个例⼦,两个矩阵交换相乘后,结果变了:下⾯咱们来具体解决这个矩阵相乘的问题。

解法⼀、暴⼒解法其实,通过前⾯的分析,我们已经很明显的看出,两个具有相同维数的矩阵相乘,其复杂度为O(n^3),参考代码如下:解法⼆、Strassen 算法在解法⼀中,我们⽤了3个for循环搞定矩阵乘法,但当两个矩阵的维度变得很⼤时,O(n^3)的时间复杂度将会变得很⼤,于是,我们需要找到⼀种更优的解法。

⼀般说来,当数据量⼀⼤时,我们往往会把⼤的数据分割成⼩的数据,各个分别处理。

遵此思路,如果丢给我们⼀个很⼤的两个矩阵呢,是否可以考虑分治的⽅法循序渐进处理各个⼩矩阵的相乘,因为我们知道⼀个矩阵是可以分成更多⼩的矩阵的。

如下图,当给定⼀个两个⼆维矩阵A B时:1 //矩阵乘法,3个for 循环搞定2 void Mul(int ** matrixA, int ** matrixB, int ** matrixC)3 {4 for (int i = 0; i < 2; ++i)5 {6 for (int j = 0; j < 2; ++j)7 {8 matrixC[i][j] = 0;9 for (int k = 0; k < 2; ++k)10 {11 matrixC[i][j] += matrixA[i][k] * matrixB[k][j];12 }13 }14 }15 }这两个矩阵A B相乘时,我们发现在相乘的过程中,有8次乘法运算,4次加法运算:矩阵乘法的复杂度主要就是体现在相乘上,⽽多⼀两次的加法并不会让复杂度上升太多。

斯特拉森算法

斯特拉森算法

斯特拉森算法斯特拉森算法(Strassen algorithm)是一种矩阵乘法算法,由德国数学家弗莱戈·斯特拉森(Frigyes Riesz)和德国计算机科学家赫尔曼·斯特拉森(Hermann Strassen)在1969年提出。

传统的矩阵乘法的时间复杂度为 $O(n^3)$,但是斯特拉森算法只需要 $O(n^{log_2 7})$ 的时间复杂度,其中 $n$ 表示矩阵的尺寸。

斯特拉森算法的优点是,对于大矩阵的乘法,它的计算速度比传统的矩阵乘法要快很多。

斯特拉森算法的核心思想是将两个 $n\times n$ 的矩阵 $A$ 和 $B$ 分别划分为四个$\frac{n}{2}\times\frac{n}{2}$ 的子矩阵,即:$$A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}, \quad B = \begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix}$$然后,按照下面的公式计算出 $C = AB$:$$\begin{aligned} P_1 &= A_{11}(B_{12}-B_{22}) \\ P_2 &= (A_{11}+A_{12})B_{22} \\ P_3 &= (A_{21}+A_{22})B_{11} \\ P_4 &= A_{22}(B_{21}-B_{11}) \\ P_5 &=(A_{11}+A_{22})(B_{11}+B_{22}) \\ P_6 &= (A_{12}-A_{22})(B_{21}+B_{22}) \\ P_7 &= (A_{11}-A_{21})(B_{11}+B_{12}) \end{aligned}$$这个公式看起来很神秘,但是在具体实现时,可以递归应用斯特拉森算法来计算$P_1$ 到 $P_7$,直到遇到 $1\times 1$ 的矩阵,或者达到某个设定的阈值。

分治法-大整数乘法和Strassen矩阵乘法

分治法-大整数乘法和Strassen矩阵乘法

分治法-⼤整数乘法和Strassen矩阵乘法4.5.1 ⼤整数乘法对于100位甚⾄更多位的⼗进制之间的乘法运算还是⽐较复杂的。

我们使⽤经典的笔算算法来对两个n位整数相乘,第⼀个数中的n个数字都要被第⼆个数中的n个数字相乘,这样就需要做n2次相乘,⽽使⽤分治技术,我们就能设计出乘法次数少于n2次的算法。

先来看下这个简单公式:令,则我们实际上要处理的就是中间的这⼀部分,就是将这两次乘法转为⼀次乘法,具体实现可由下⾯这个公式得到:我们令,所以,原式为:额,这个算法还是有点复杂,代码不知道该怎么写。

4.5.2 S t rassen矩阵乘法V.Strassen在1969年发表了这个算法,它的成功依赖于这个发现:计算两个2阶⽅阵A和B的积C只需要进⾏7次乘法运算,⽽不是蛮⼒算法所需要的8次。

公式参照如下:其中,因此,对于两个2阶⽅阵相乘时,Strassen算法执⾏了7次乘法和18次加减法,⽽蛮⼒法需要执⾏8次乘法和4次加法。

虽然只是减少了⼀次乘法,但当矩阵的阶趋于⽆穷⼤时,算法卓越的效率就渐渐表现出来了。

代码实现这个算法对我来说感觉还是有点复杂:-),毕竟考虑的因素有很多,因为进⾏乘法运算的矩阵并不都是2n阶的,⽽且矩阵之间是⽆法进⾏乘法运算的,总之,思路感觉有点多啊。

以下代码是我排除了各种不定因素,且进⾏乘法运算的矩阵都是2n阶的⽅阵(好像是有点low哦,不过不管啦)。

代码实现:/*** Strassen算法进⾏矩阵相乘* @author xiaofeig* @since 2015.9.19* @param marix1 要进⾏相乘的矩阵1* @param marix2 要进⾏相乘的矩阵2* @return返回相乘的结果* */public static int[][] strassenMultiplyMatrix(int[][] marix1, int[][] marix2){if(marix1.length==1){return new int[][]{{marix1[0][0]*marix2[0][0]}};}int xLen=marix1[0].length;int yLen=marix1.length;int[][] a00=copyArrayOfRange(marix1, 0, 0, yLen/2, xLen/2);int[][] a01=copyArrayOfRange(marix1, 0, xLen/2, yLen/2, xLen);int[][] a10=copyArrayOfRange(marix1, yLen/2, 0, yLen, xLen/2);int[][] a11=copyArrayOfRange(marix1, yLen/2, xLen/2, yLen, xLen);xLen=marix2[0].length;yLen=marix2.length;int[][] b00=copyArrayOfRange(marix2, 0, 0, yLen/2, xLen/2);int[][] b01=copyArrayOfRange(marix2, 0, xLen/2, yLen/2, xLen);int[][] b10=copyArrayOfRange(marix2, yLen/2, 0, yLen, xLen/2);int[][] b11=copyArrayOfRange(marix2, yLen/2, xLen/2, yLen, xLen);int[][] m1=strassenMultiplyMatrix(plusMarix(a00, a11), plusMarix(b00, b11));int[][] m2=strassenMultiplyMatrix(plusMarix(a10, a11), b00);int[][] m3=strassenMultiplyMatrix(a00, minusMarix(b01, b11));int[][] m4=strassenMultiplyMatrix(a11, minusMarix(b10, b00));int[][] m5=strassenMultiplyMatrix(plusMarix(a00, a01), b11);int[][] m6=strassenMultiplyMatrix(minusMarix(a10, a00), plusMarix(b00, b01));int[][] m7=strassenMultiplyMatrix(minusMarix(a01, a11), plusMarix(b10, b11));int[][] newMarix1=plusMarix(minusMarix(plusMarix(m1, m4), m5), m7);int[][] newMarix2=plusMarix(m3, m5);int[][] newMarix3=plusMarix(m2, m4);int[][] newMarix4=plusMarix(minusMarix(plusMarix(m1, m3), m2), m6);return mergeMarix(newMarix1, newMarix2, newMarix3, newMarix4);}/*** 复制指定矩阵的某范围内的数据到以新的数组* @author xiaofeig* @since 2015.9.19* @param array ⽬标数组* @param i,j 左上⾓元素下标(包含)* @param m,n 右下⾓元素下标(不包含)* @return返回指定数组某范围的新数组* */public static int[][] copyArrayOfRange(int[][] array,int i,int j,int m,int n){int[][] result=new int[m-i][n-j];int index=0;while(i<m){result[index]=Arrays.copyOfRange(array[i], j, n);index++;i++;}return result;}/*** 进⾏矩阵之间的加法运算* @author xiaofeig* @since 2015.9.19* @param marix1 加数矩阵1* @param marix2 加数矩阵2* @return返回结果矩阵* */public static int[][] plusMarix(int[][] marix1,int[][] marix2){int[][] result=new int[marix1.length][marix1[0].length];for(int i=0;i<marix1.length;i++){for(int j=0;j<marix1[0].length;j++){result[i][j]=marix1[i][j]+marix2[i][j];}}return result;}/*** 进⾏矩阵之间的减法运算* @author xiaofeig* @since 2015.9.19* @param marix1 减数矩阵* @param marix2 被减数矩阵* @return返回结果矩阵* */public static int[][] minusMarix(int[][] marix1,int[][] marix2){int[][] result=new int[marix1.length][marix1[0].length];for(int i=0;i<marix1.length;i++){for(int j=0;j<marix1[0].length;j++){result[i][j]=marix1[i][j]-marix2[i][j];}}return result;}/*** 将四个矩阵合并为⼀个矩阵* @param marix1 数组1* @param marix2 数组2* @param marix3 数组3* @param marix4 数组4* @return返回合并之后的新矩阵* */public static int[][] mergeMarix(int[][] marix1,int[][] marix2,int[][] marix3,int[][] marix4){ int m=marix1.length,n=marix1[0].length;int[][] marix=new int[m*2][n*2];for(int i=0;i<marix.length;i++){for(int j=0;j<marix[i].length;j++){if(i<m){if(j<n){marix[i][j]=marix1[i][j];}else{marix[i][j]=marix2[i][j-n];}}else{if(j<n){marix[i][j]=marix3[i-m][j];}else{marix[i][j]=marix4[i-m][j-n];}}}}return marix;}算法分析:上⾯的代码我⽤了两个23阶的矩阵测试过,结果是正确的,其它阶数的矩阵我没测试,估计会有很多错误。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

7
实现功能:
• 1、实现Strassen算法。
• 2、Strassen要求阶n是2的幂,但这样的情况很少。对任 一偶数n,求n阶矩阵想乘。 • 3、对任意的正整数n,求n阶矩阵相乘。
8
当矩阵的维数2的幂,例如4阶矩阵相乘:C=A*B
a11 a 21 C a31 a41
a12 a22 a32 a42
B12 B22
当矩阵的维数为任一偶数n,总有n=m*2k,将矩阵分为m*m个 的2k阶矩阵。例如n=6,把矩阵分为3*3个21小矩阵。小矩阵用 Strassen相乘,大矩阵用传统算法。
*
=
当矩阵的维数n为奇数时,例如n=5时,分别给矩阵加一行一列 0。
算法代码
指向指针的指针: int **C; C=(int**)malloc(n*sizeof(int)); for(int i=0;i<n;i++) C[i]=(int *)malloc(n*sizeof(int)); 之后,就可以用到C[i][j]取到元素了。
运行结果
20
Thank you
23
#include <iostream> #include <time.h> using namespace std; class Matrix{ public: int **m; int n; Matrix(){ m=NULL; n=0; } Matrix(int n){ this->n=n; m=new int*[n]; for(int i=0;i<n;i++){ m[i]=new int[n]; for(int j=0;j<n;j++) m[i][j]=0; } } void InitMatrix(){ for(int i=0;i<n;i++) for(int j=0;j<n;j++) m[i][j]=rand()%10; }
Strassen算法
姓名: 王军袖 学号: 31609059
1
Hale Waihona Puke 目录1算法思想
2
3
算法代码
运行结果
2
算法思想
分治算法
其核心思想是将大的问题分解成若干小的子问题进行 处理,从而使复杂的问题变得简单。表现在算法中就 是递归算法。
3
一般的矩阵乘法:
c[i, j ] k 1 a[i, k ] * b[k , j ]
k n
C11=A11B11+A12B21 C12=A11B12+A12B22 C21=A21B11+A22B21 C22=A21B12+A22B22
算法分析:
T(n)=n3
传统方法2个2阶方阵相乘需要8次乘法。按照分治的思想可以看 出,要想减少乘法运算次数,关键在于计算2个2阶方阵的乘积时,能 否用少于8次的乘法运算。Strassen提出了一种新的算法来计算2个2阶 方阵的乘积。此算法只用了7次乘法运算,但增加了加、减法的运算 次数。 先计算下面7个量:
15
//阶为2的n次幂阶矩阵相乘
Matrix Muti2n(Matrix A,Matrix B){ if(A.n==2)//如果阶是2,直接相乘 return MutiFor2(A,B); //将矩阵分成四块 //用Strassen相乘 Matrix M1=Muti2n(Plus(A11,A22),Plus(B11,B22)); Matrix M2=Muti2n(Plus(A21,A22),B11); Matrix M3=Muti2n(A11,Minus(B12,B22)); Matrix M4=Muti2n(A22,Minus(B21,B11)); Matrix M5=Muti2n(Plus(A11,A12),B22); Matrix M6=Muti2n(Minus(A21,A11),Plus(B11,B12)); Matrix M7=Muti2n(Minus(A12,A22),Plus(B21,B22)); Matrix C(A.n); for (int i=0;i<C.n/2;i++) for (int j=0;j<C.n/2;j++){ C.m[i][j]=M1.m[i][j]+M4.m[i][j]-M5.m[i][j]+M7.m[i][j]; C.m[i][j+C.n/2]=M3.m[i][j]+M5.m[i][j]; C.m[i+C.n/2][j]=M2.m[i][j]+M4.m[i][j]; C.m[i+C.n/2][j+C.n/2]=M1.m[i][j]+M3.m[i][j]-M2.m[i][j]+M6.m[i][j]; } return C;}
6
递归中最了7次(n/2)阶矩阵相乘,18次加法,认为 是一个常数c。 设n阶矩阵相乘需T(n),则 T(1)=1 T(n)=7T(n/2)+c 解之得:T(n)=O(nlog7)≈O(n2.81)
一般矩阵乘法的时间复杂度为:T(n)=n3。由此可见,Strassen矩阵乘 法的计算时间复杂性比普通矩阵乘法有改进。
13
//分块函数,分割矩阵A,将第x行第y列的元素作为新矩阵的元素, 分割出一个n阶矩阵。 Matrix Split(Matrix A,int x,int y,int n){ Matrix C(n); for (int i=0;i<n;i++) for (int j=0;j<n;j++) C.m[i][j]=A.m[i+x][j+y]; return C; }
14
//阶为2的方阵相乘: Matrix MutiFor2(Matrix A,Matrix B){//二阶矩阵相乘,返回积 Matrix C(A.n); int M1=(A.m[0][0]+A.m[1][1])*(B.m[0][0]+B.m[1][1]); int M2=(A.m[1][0]+A.m[1][1])*B.m[0][0]; int M3=A.m[0][0]*(B.m[0][1]-B.m[1][1]); int M4=A.m[1][1]*(B.m[1][0]-B.m[0][0]); int M5=(A.m[0][0]+A.m[0][1])*B.m[1][1]; int M6=(A.m[1][0]-A.m[0][0])*(B.m[0][0]+B.m[0][1]); int M7=(A.m[0][1]-A.m[1][1])*(B.m[1][0]+B.m[1][1]); C.m[0][0]=M1+M4-M5+M7; C.m[0][1]=M3+M5; C.m[1][0]=M2+M4; C.m[1][1]=M1+M3-M2+M6; return C;}
//给矩阵加一行和一列0 void ReinitMatris(){ int i,j; n=n+1; for( i=0;i<n-1;i++){ delete m[i]; } m=new int*[n]; for( i=0;i<n;i++){ m[i]=new int[n]; for( j=0;j<n;j++) m[i][j]=0; } for( i=0;i<n-1;i++) for( j=0;j<n-1;j++) m[i][j]=rand()%10;}
16
Matrix **cr=new Matrix*[n];//初始化cr,储存ar*br结果 //阶为偶数矩阵想乘 for ( i=0;i<n;i++){ Matrix MutiEven(Matrix A,Matrix B){ cr[i]=new Matrix[n]; int n=A.n,m=1; for ( j=0;j<n;j++) int i,j; cr[i][j]=Matrix(m); //将矩阵阶拆分成m*n,n为奇数 } m为2的幂次 for( i=0;i<n;i++)//用传统方法将ar和br相乘,储存在cr中 while(!(n&1)){ for( j=0;j<n;j++) n>>=1; for(int k=0;k<n;k++) m<<=1; cr[i][j]=Plus(cr[i][j],Muti2n(ar[i][k],br[k][j])); Matrix C(m*n);//声明C } for( i=0;i<n;i++)//将cr合并成一个矩阵C //将A拆分为n^2个m阶矩阵,在ar中储存。 for( j=0;j<n;j++) Matrix **ar=new Matrix*[n]; for(int x=0;x<m;x++) for ( i=0;i<n;i++){ for(int y=0;y<m;y++) ar[i]=new Matrix[n]; for ( j=0;j<n;j++) C.m[i*m+x][j*m+y]=cr[i][j].m[x][y]; ar[i][j]=Split(A,i*m,j*m,m); return C;} } Matrix **br=new Matrix*[n];//同A for (i=0;i<n;i++){ br[i]=new Matrix[n]; for ( j=0;j<n;j++) 17 br[i][j]=Split(B,i*m,j*m,m);}
5
相关文档
最新文档