02
分 治 法
divide-and-conquer algorithm
无论人们在祈祷什么，他们总是在祈祷一个奇迹。每一个祈祷都可以简化为：伟大的上帝呀，请让两个二相加不等于四吧。
——伊万·屠格涅夫（1818 － 1883），俄国作家和短篇小说家
关于分治法，本公司的另一位老板在一年前也写过啦（老板真强，老板真强。。。），大家可以看看那篇，很详细（不是说讲解详细）：
经典优化算法之分治法（Divide-and-Conque Algorithm）
转载 |【算法】分治法（Divide-and-Conquer Algorithm）经典例子分析
从字面上分析就可以看到有哪些步骤：
分-分解-将问题分解为规模更小的子问题，子问题最好相同或相似；
治-求解-将这些规模更小的子问题逐个击破；
合-合并-将已解决的子问题合并，最终得出原问题的解；
（有时原始问题的解只存在于分解出的某一个（或某几个）子问题中，则只需要在这一（或这几个）子问题中求解即可。例如在图书馆里找书。）
从上述步骤中我们可以看出，分治算法一般适用满足以下条件的场景：
1）问题规模缩小到一定的程度就可以很容易解决；
2）问题可以分解为若干个规模较小的相同问题；
3）问题分解出的若干子问题的解可以合并为该问题的解；
4）每个子问题都是独立的，相互之间没有交集。（这是区别分治法与减）
在“分”的过程中，我们尽可能让分解出的子问题与原始问题相似，而规模更小。这刚好符合递归的特性。因此，分治法往往与递归联系在一起。
说到这里，是不是有小伙伴觉得分治法与减治法很相似，傻傻分不清？我们这里再举f(n)=a^n的栗子。
n>1时：f(n)=a^(n/2)*a^(n/2)
n=1时：f(n)=a
比较一下减常数因子：
n是偶数时：f(n)= f(n/2)^2
n是大于１的奇数时：f(n)= f((n-1)/2)^2* a
n = 1时：f(n)=a
也就是说，分治法对分治出的部分需要分别处理，进行分开的单独计算，而减治法则利用了"一个问题给定实例的解和同样问题较小实例的解之间的关系"，只针对部分子问题求解，减治掉的那部分就不需要了。
其实，减常因子的减治法也可以看做是分治的变种。
需要注意的是，不是所有的分治算法都一定比简单蛮干更有效，前面的减治法也是，就比方说这里的栗子，时间复杂度仍为o（n）。不过，“通常我们向算法女神所做的祈祷都被应允了”，使用分治法往往比使用其他方法效率更高。
分治法的时间效率T(n)一般满足方程T(n)=aT(n/b)+f(n)。因此，分治法对于并行算法是非常理想的。
我们介绍一个和数学有关的算法：
Strassen矩阵乘法。
（虽然在骁老板的文章里有提到，但我还是不辞辛苦地再写一遍八）
考虑到大家都是小白，先说明矩阵乘法的定义。（其实是我自己不懂才先去查的）
这张图很清楚的说明了矩阵乘法的计算公式。为了方便讲解，我们先以n*n的偶数阶方阵为例，之后再拓展到一般的矩阵乘法。
我们从数学中回到算法来。这个问题如果直接暴力计算，需要循环三次：关于i,j,k分别循环。时间复杂度为o（n^3）。
我们采用分治的思想，把偶数阶方阵如图划分为四份（这里的ABC可以是矩阵）：
按照矩阵中的知识，可以得到：
C11 = A11 • B11 + A12 • B21
C12 = A11 • B12 + A12 • B22
C21 = A21 • B11 + A22 • B21
C22 = A21 • B12 + A22 • B22
我们估算一下时间复杂度：T(n) = 8T(n/2) + o(n^2)（前面是不是见过这个式子？）递归的结果是时间复杂度并没有降低。算法女神去死八。
聪明的Strassen不甘心。他发明了一种新的方法，通过降低递归式中T(n/2)的系数来降低时间复杂度。
仍然把每个矩阵分割为4份，然后创建如下10个中间矩阵：
S1 = B12 - B22
S2 = A11 + A12
S3 = A21 + A22
S4 = B21 - B11
S5 = A11 + A22
S6 = B11 + B22
S7 = A12 - A22
S8 = B21 + B22
S9 = A11 - A21
S10 = B11 + B12
接着，计算7次矩阵乘法，在建立7个中间的中间矩阵（中中间矩阵）：
M1 = A11 • S1
M2 = S2 • B22
M3 = S3 • B11
M4 = A22 • S4
M5 = S5 • S6
M6 = S7 • S8
M7 = S9 • S10
最后，根据这7个中中间矩阵就可以计算出C矩阵：
C11 = M5 + M4 – M2 + M6
C12 = M1 + M2
C21 = M3 + M4
C22 = M5 + M1 - M3 - M7
证明就留给线性代数的同学们吧，应该不难。当然能不能想到就另说咯。
但不管怎么说，我们利用这个结论，成功地将系数8变成了7，算法的时间复杂度也就降低到了o（n^log7）。
我们再回到一般性矩阵。对于一般性矩阵，我们还是认真看定义吧。。。
对于非偶数阶方阵，我们可以用0将其填充为偶数阶方阵：
如果是奇数阶方阵，我们也可以在找到最近的偶数阶方阵，其余部分直接暴力计算。
（代码中我只实现第一种方式）
很明显这种优化是为了处理大型问题的，既然如此我们在写代码时也需要开辟足够大的规模。一种方法是用指针创建它（这就是代码冗长的主要原因）。因为是二维数组，所以需要用到指向指针的指针，再用数组表示指针，然后就可以用熟悉的处理数组的方式处理数据啦。
代码小长，实际上也没什么内容：
//strassen算法（在此感谢互联网，减少工作量） #include <iostream> #include <iomanip> #include <cmath> using namespace std; void Strassen(int N, int** MatrixA, int ** MatrixB, int ** MatrixC); //进行矩阵乘法 void ADD(int** MatrixA, int** MatrixB, int** MatrixResult, int length ); //矩阵加法 void SUB(int** MatrixA, int** MatrixB, int** MatrixResult, int length ); //矩阵减法 void MUL(int** MatrixA, int** MatrixB, int** MatrixResult ); //计算规模最小的二阶矩阵（常规方法） void FillMatrix( int** matrix, int, int,int); //输入原始矩阵 void PrintMatrix( int **Matrix, int length,int width ); //输出矩阵乘法结果 int findN(int l1,int l2,int l3,int l4); //非偶数阶方阵，寻找扩充边界 int main() { int length1,length2,width1,width2; int** MatrixA; //指向指针的指针，存放矩阵数据 int** MatrixB; int** MatrixC; cout<<"请输入两个矩阵的规模(先长后宽，第一个组的宽与第二组的长相等): "<<endl; cin>>length1>>width1; cin>>length2>>width2; int N =findN(length1,length2,width1,width2); MatrixA = new int *[N]; //开辟一维指针数组，每个数组容量为N MatrixB = new int *[N]; MatrixC = new int *[N]; for (int i = 0; i < N; i++) { MatrixA[i] = new int [N]; //开辟二维数组 ，每个数组容量为N*N MatrixB[i] = new int [N]; MatrixC[i] = new int [N]; } cout<<"请分别输入两个矩阵的数据："<<endl; FillMatrix(MatrixA,length1,width1,N); FillMatrix(MatrixB,length2,width2,N); Strassen( N, MatrixA, MatrixB, MatrixC ); cout<<"运算结果为："; int length3=length1>length2?length1:length2; PrintMatrix(MatrixC,length1,width2); return 0; } //寻找拓展规模N int findN(int l1,int l2,int l3,int l4) { int max,k=1,max1,max2; max1=l1>l2?l1:l2; max2=l3>l4?l3:l4; max=max1>max2?max1:max2; while (true) { if (max<=pow(2,k)) return pow(2,k); k++; } } //计算矩阵乘法 void Strassen(int N, int **MatrixA, int **MatrixB, int **MatrixC) { int HalfSize = N/2; int newSize = N/2; if ( N == 2 ) //边界最小规模：二阶矩阵 { MUL(MatrixA,MatrixB,MatrixC); } else { //和主函数里一样，定义指针，开辟空间 int** A11; int** A12; int** A21; int** A22; int** B11; int** B12; int** B21; int** B22; int** C11; int** C12; int** C21; int** C22; int** M1; int** M2; int** M3; int** M4; int** M5; int** M6; int** M7; int** AResult; //result返回结果 int** BResult; A11 = new int *[newSize]; A12 = new int *[newSize]; A21 = new int *[newSize]; A22 = new int *[newSize]; B11 = new int *[newSize]; B12 = new int *[newSize]; B21 = new int *[newSize]; B22 = new int *[newSize]; C11 = new int *[newSize]; C12 = new int *[newSize]; C21 = new int *[newSize]; C22 = new int *[newSize]; M1 = new int *[newSize]; M2 = new int *[newSize]; M3 = new int *[newSize]; M4 = new int *[newSize]; M5 = new int *[newSize]; M6 = new int *[newSize]; M7 = new int *[newSize]; AResult = new int *[newSize]; BResult = new int *[newSize]; int newLength = newSize; //同上 for ( int i = 0; i < newSize; i++) { A11[i] = new int[newLength]; A12[i] = new int[newLength]; A21[i] = new int[newLength]; A22[i] = new int[newLength]; B11[i] = new int[newLength]; B12[i] = new int[newLength]; B21[i] = new int[newLength]; B22[i] = new int[newLength]; C11[i] = new int[newLength]; C12[i] = new int[newLength]; C21[i] = new int[newLength]; C22[i] = new int[newLength]; M1[i] = new int[newLength]; M2[i] = new int[newLength]; M3[i] = new int[newLength]; M4[i] = new int[newLength]; M5[i] = new int[newLength]; M6[i] = new int[newLength]; M7[i] = new int[newLength]; AResult[i] = new int[newLength]; BResult[i] = new int[newLength]; } //拆分A、B矩阵为四个小矩阵 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]; } } //计算7个中中间矩阵 M //M1=(A11+A22)*(B11+B22) ADD( A11,A22,AResult, HalfSize); ADD( B11,B22,BResult, HalfSize); Strassen( HalfSize, AResult, BResult, M1 ); //M2=(A21+A22)*B11 ADD( A21,A22,AResult, HalfSize); Strassen(HalfSize, AResult, B11, M2); //M3=A11*(B12-B22) SUB( B12,B22,BResult, HalfSize); Strassen(HalfSize, A11, BResult, M3); //M4=A22(B21-B11) SUB( B21, B11, BResult, HalfSize); Strassen(HalfSize, A22, BResult, M4); //M5=(A11+A12)B22 ADD( A11, A12, AResult, HalfSize); Strassen(HalfSize, AResult, B22, M5); //M6=(A21-A11)(B11+B12) SUB( A21, A11, AResult, HalfSize); ADD( B11, B12, BResult, HalfSize); Strassen( HalfSize, AResult, BResult, M6); //M7=(A12-A22)(B21+B22) SUB(A12, A22, AResult, HalfSize); ADD(B21, B22, BResult, HalfSize); Strassen(HalfSize, AResult, BResult, M7); //计算出新矩阵的四个部分C //C11 = M1 + M4 - M5 + M7; ADD( M1, M4, AResult, HalfSize); SUB( M7, M5, BResult, HalfSize); ADD( AResult, BResult, C11, HalfSize); //C12 = M3 + M5; ADD( M3, M5, C12, HalfSize); //C21 = M2 + M4; ADD( M2, M4, C21, HalfSize); //C22 = M1 + M3 - M2 + M6; ADD( M1, M3, AResult, HalfSize); SUB( M6, M2, BResult, HalfSize); ADD( AResult, BResult, C22, HalfSize); //合并为新矩阵 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]; } } // 释放空间 for (int i = 0; i < newLength; i++) { delete[] A11[i];delete[] A12[i];delete[] A21[i];delete[] A22[i]; delete[] B11[i];delete[] B12[i];delete[] B21[i];delete[] B22[i]; delete[] C11[i];delete[] C12[i];delete[] C21[i];delete[] C22[i]; delete[] M1[i];delete[] M2[i];delete[] M3[i];delete[] M4[i];delete[] M5[i];delete[] M6[i];delete[] M7[i]; delete[] AResult[i];delete[] BResult[i] ; } delete[] A11;delete[] A12;delete[] A21;delete[] A22; delete[] B11;delete[] B12;delete[] B21;delete[] B22; delete[] C11;delete[] C12;delete[] C21;delete[] C22; delete[] M1;delete[] M2;delete[] M3;delete[] M4;delete[] M5; delete[] M6;delete[] M7; delete[] AResult; delete[] BResult ; } } //计算矩阵加法 ，记过输入result void ADD(int** MatrixA, int** MatrixB, int** MatrixResult, int MatrixSize ) { for ( int i = 0; i < MatrixSize; i++) { for ( int j = 0; j < MatrixSize; j++) { MatrixResult[i][j] = MatrixA[i][j] + MatrixB[i][j]; } } } //矩阵减法，同上 void SUB(int** MatrixA, int** MatrixB, int** MatrixResult, int MatrixSize ) { for ( int i = 0; i < MatrixSize; i++) { for ( int j = 0; j < MatrixSize; j++) { MatrixResult[i][j] = MatrixA[i][j] - MatrixB[i][j]; } } } //暴力法求解，计算二阶矩阵 void MUL( int** MatrixA, int** MatrixB, int** MatrixResult ) { for (int i=0;i<2 ;i++) { for (int j=0;j<2 ;j++) { MatrixResult[i][j]=0; for (int k=0;k<2 ;k++) { MatrixResult[i][j]=MatrixResult[i][j]+MatrixA[i][k]*MatrixB[k][j]; } } } } //输入数据 ，填充为偶数阶方阵 void FillMatrix( int** Matrix, int length,int width,int N) { for(int row = 0; row<N; row++) { if (row<length) { for(int column = 0; column<N; column++) { if (column<width) cin>>Matrix[row][column]; else Matrix[row][column]=0; } } else { for(int column = 0; column<N; column++) Matrix[row][column]=0; } } } //输出正确规模的新矩阵 void PrintMatrix(int **MatrixA,int length,int width) { cout<<endl; for(int row = 0; row<length; row++) { for(int column = 0; column<width; column++) { cout<<MatrixA[row][column]<<"\t"; if ((column+1)%((width)) == 0) cout<<endl; } } cout<<endl; }
（检验用了小一点的数组，别介意。。。）
