一、分治策略
分治法的思想
将原问题分解为几个规模较小但类似于原问题的子问题,递归地求解这些子问题,然后再合并这些子问题的解来建立原问题的解。
递归式
递归式与分治方法是紧密相关的,因为使用递归式可以很自然地刻画分治算法的运行时间。
在分治策略中,我们递归地求解一个问题,在每层递归中应用如下的三个步骤:
分解:将问题划分为一些子问题,子问题的形式与原问题一样,只是规模更小。
解决:递归地求解出子问题。如果子问题的规模足够小,则停止递归,直接求解。
合并:将子问题的解组合成原问题的解。
学习心得
递归形式的算法典型的遵循了分治思想,然而采用递归算法进行递归采用具体的计算方式的不同(递归方程不同)会导致算法时间复杂度的不同,因而即使算法都是使用了递归的方法,仍说不定有些更好的技巧使算法时间复杂度减小。一个递归算法,它可能只分成一个规模更小的子问题(分治思想在这里表现的不明显,因为子问题并不是多个,只有一个,但是问题规模确实是变小了),也可能分解成多个规模更小的子问题(分治思想表现明显)。
二、最大子数组问题
问题描述:对给定数组A,寻找A的和最大的非空连续子数组。
例子:数组{2,4,-7,5,2,-1,2,-4,3}的最大连续子数组为{5,2,-1,2},最大连续子数组的和为5+2-1+2=8。
版本一:暴力求解法
原理:计算A中每种子数组的和,然后找到最大的子数组。
时间复杂度:需要一次嵌套循环迭代进行求解(对于每种子数组的起始位置,再求解每种子数组的结束位置下的和),故时间复杂度T(n) = O(n^2)。
体会:在计算暴力求解法的时间复杂度时,可以利用排列和组合公式辅助计算。另外,算法实际上也用到了动态规划的思想:算法求和时采用累加,它将前一步的结果存储起来并在后续的步骤中直接使用,减少重复性劳动。若不采用累加,时间复杂度将会变为T(n) = O(n^3)。
C语言实现核心代码:
// directly cacluating TYPE j; TYPE index_left, index_right; TYPE diff_max = INT_MIN; // 起始累加位置循环 for(i = 0; i < count - 1; i++) { temp = array[i]; // 结束累加位置循环 for(j = i + 1; j < count; j++) { // 存储最大子数组 if((temp += array[j]) > diff_max) { diff_max = temp; index_left = i; index_right = j; } } }
版本二:分治法
原理:采用分治的思想,步骤如下
分解:将数组划分为规模更小且规模尽量相等的子数组,这里把数组从中间位置划分成子数组。
解决:使用分治法递归地求解这些子数组的最大子数组。最小子数组时当数组长度为一的情况,此时返回该元素的值。
合并:求解跨越两个子数组的最大子数组,即跨越这个相邻子数组合成的数组中点的最大子数组,并与两个子数组各自的最大子数组进行比较,选取和最大者。
时间复杂度:采用递归算法实现:T(n) = 2T(n/2) + O(n) (n > 1),得出 T(n)=O(nlgn)。
体会:分治的思想最主要考虑的是两个方面:1.如何把问题分解成子问题,这里实际上有很多种想法,比如数据对半分,比如只减少一个数据,合理的划分才能体现出分治思想的好处。2.如何合并两个子问题,这里需要明确子问题是什么,合并后解决的问题是当前问题,对于上一级来说实际上也是一个子问题。
C语言实现核心代码:
/// used for save data struct subarray { TYPE index_begin; ///< index begin from an array when get max sum TYPE index_end; ///< index end from an array when get max sum TYPE sum_max; ///< max sum of the subarray }; /** * @brief get the max subarray when cross the middle point of the array * * @param[in] array * @param[in] index_begin * @param[in] index_end * * @return subarray */ struct subarray find_max_crossing_subarray(TYPE* array, TYPE index_begin, TYPE index_end) { // should have special input value assert(index_begin < index_end); // get middle index TYPE index_middle = (unsigned)(index_begin + index_end) >> 1; struct subarray return_subarray; TYPE i; TYPE sum_temp = INT_MIN; TYPE sum = 0; for(i = index_middle; i >= index_begin; i--) { sum += array[i]; if(sum > sum_temp) { return_subarray.index_begin = i; sum_temp = sum; } } return_subarray.sum_max = sum_temp; sum = 0; sum_temp = INT_MIN; for(i = index_middle + 1; i <= index_end; i++) { sum += array[i]; if(sum > sum_temp) { return_subarray.index_end = i; sum_temp = sum; } } return_subarray.sum_max += sum_temp; return return_subarray; } /** * @brief find maximum subarray by divide and conquer method. * * @param TYPE * @param index_begin * @param index_end */ struct subarray find_maximum_subarray(TYPE* array, TYPE index_begin, TYPE index_end) { DEBUG_PRINT_STRING("In recursion now.\n"); DEBUG_PRINT_VALUE("%d", index_begin); DEBUG_PRINT_VALUE("%d", index_end); assert(index_begin <= index_end); if(index_begin == index_end) { struct subarray return_subarray = {index_begin, index_end, array[index_begin]}; DEBUG_PRINT_VALUE("%d", return_subarray.sum_max); DEBUG_PRINT_VALUE("%d", return_subarray.index_begin); DEBUG_PRINT_VALUE("%d", return_subarray.index_end); DEBUG_PRINT_VALUE("%d", index_begin); DEBUG_PRINT_VALUE("%d", index_end); DEBUG_PRINT_STRING("Out recursion now.\n"); return return_subarray; } else { TYPE index_middle = (unsigned)(index_begin + index_end) >> 1; struct subarray data_left = find_maximum_subarray(array, index_begin, index_middle); struct subarray data_right = find_maximum_subarray(array, index_middle + 1, index_end); struct subarray data_cross = find_max_crossing_subarray(array, index_begin, index_end); if(data_left.sum_max >= data_right.sum_max && data_left.sum_max >= data_cross.sum_max) { DEBUG_PRINT_VALUE("%d", data_left.sum_max); DEBUG_PRINT_VALUE("%d", data_left.index_begin); DEBUG_PRINT_VALUE("%d", data_left.index_end); DEBUG_PRINT_VALUE("%d", index_begin); DEBUG_PRINT_VALUE("%d", index_end); DEBUG_PRINT_STRING("Out recursion now.\n"); return data_left; } else if(data_right.sum_max >= data_left.sum_max && data_right.sum_max >= data_cross.sum_max) { DEBUG_PRINT_VALUE("%d", data_right.sum_max); DEBUG_PRINT_VALUE("%d", data_right.index_begin); DEBUG_PRINT_VALUE("%d", data_right.index_end); DEBUG_PRINT_VALUE("%d", index_begin); DEBUG_PRINT_VALUE("%d", index_end); DEBUG_PRINT_STRING("Out recursion now.\n"); return data_right; } else { DEBUG_PRINT_VALUE("%d", data_cross.sum_max); DEBUG_PRINT_VALUE("%d", data_cross.index_begin); DEBUG_PRINT_VALUE("%d", data_cross.index_end); DEBUG_PRINT_VALUE("%d", index_begin); DEBUG_PRINT_VALUE("%d", index_end); DEBUG_PRINT_STRING("Out recursion now.\n"); return data_cross; } } }
版本三:动态规划算法(线性时间)
原理:利用动态规划的思想,保存但前的状态,解决重叠子问题,结合状态转移方程求解。从后往前分析:考虑加入第n的元素a[n]后的数组a[0--n]的最大子数组与加入前的数组a[0--n-1]的最大子数组之间的关系,找到状态转移的方程。它们的关系可能如下两种情况,
1.数组a[0--n]的最大子数组以a[n]元素结尾,可能由单个a[n]元素组成;
2.数组a[0--n]的最大子数组与a[n]无关,与数组a[0--n-1]的最大子数组相同。
若我们存储以a[n]元素结尾的最大子数组状态为end[n],数组a[0--n]的最大子数组的状态为all[n],
则得出状态转移方程如下:
end[n] = max(end[n-1] + array[n], array[n]);
all[n] = max(all[n-1], end[n]);
时间复杂度:只设计一次遍历:T(n)=O(n)。
体会:动态规划的方法主要考虑子问题和问题之间的转态转移关系,利用状态转移方程和初始条件进行求解,关键是要明确状态是如何转移的,分析好这一步,算法就比较好实现了。
学习:该算法的一个附带的有点是:它只对数据进行一次扫描,一旦元素被读入并被处理,它就不再需要被记忆。因此,如果数组在磁盘或磁带上,他就可以被顺序读入,在主存中不必存储数组的任何部分。不仅如此,在任意时刻,该算法都能对它已经读入的数据给出最大子数组(前面两种算法不具有这种特性)。具有这种特性的算法叫做联机算法。
C语言实现核心代码:
/// used for save data struct subarray { TYPE index_begin; ///< index begin from an array when get max sum TYPE index_end; ///< index end from an array when get max sum TYPE sum_max; ///< max sum of the subarray }; /** * @brief find maximum subarray in linear time. * * @param TYPE * @param index_begin * @param index_end */ struct subarray find_maximum_subarray(TYPE* array, TYPE index_begin, TYPE index_end) { assert(index_begin <= index_end); struct subarray end_current; struct subarray all_current; struct subarray end_previous = {index_begin, index_begin, array[index_begin]}; struct subarray all_previous = {index_begin, index_begin, array[index_begin]}; TYPE i; for(i = index_begin + 1; i <= index_end; i++) { if(end_previous.sum_max < 0) { end_current.sum_max = array[i]; end_current.index_begin = end_current.index_end = i; } else { end_current.sum_max = end_previous.sum_max + array[i]; end_current.index_end++; } if(end_current.sum_max >= all_previous.sum_max) { all_current = end_current; } else { all_current = all_previous; } all_previous = all_current; end_previous = end_current; } return all_current; }
版本四:动态规划算法优化(线性时间)
版本三得出的状态转移方程如下:
end[n] = max(end[n-1] + array[n], array[n]);
all[n] = max(all[n-1], end[n]);
分析第一个方程,可以发现只要end[n-1]小于0,那么end[n]将重新开始加入array[n]元素再进行累加,那么我们可以通过判断end[n]在累加过程中是否小于零,来觉得继续累加还是重新开始,这样我们就不用分开记录前一次end[n-1]和当前end[n]两个状态而通过判断语句合二为一。
这个版本也可以用另外一种思路去理解:在进行累加过程中,如果当前得到的和是个负数,那么这个和在接下来的累加中应该抛弃并重新清零,不然的话这个负数将会减少接下来的和。
时间复杂度:只设计一次遍历:T(n)=O(n)。
体会:循环中存储上一次结果和当前结果时不一定非要分开存储然后一次循环结束后赋值,可以通过他们之间的关系进行进一步简化,利用上一次结果直接存储为当前结果。
C语言实现核心代码:
/// used for save data struct subarray { TYPE index_begin; ///< index begin from an array when get max sum TYPE index_end; ///< index end from an array when get max sum TYPE sum_max; ///< max sum of the subarray }; /** * @brief find maximum subarray in linear time. * * @param TYPE * @param index_begin * @param index_end */ struct subarray find_maximum_subarray(TYPE* array, TYPE index_begin, TYPE index_end) { assert(index_begin <= index_end); struct subarray return_subarray = {index_begin, index_begin, array[index_begin]}; TYPE sum_end = array[index_begin]; TYPE sum_all = array[index_begin]; TYPE index_begin_current = index_begin; TYPE i; for(i = index_begin + 1; i <= index_end; i++) { if(sum_end < 0) { sum_end = array[i]; index_begin_current = i; } else { sum_end += array[i]; } if(sum_end >= sum_all) { sum_all = sum_end; return_subarray.index_begin = index_begin_current; return_subarray.index_end = i; } } return_subarray.sum_max = sum_all; return return_subarray; }
此问题详细源代码查看及下载(UBUNTU系统+ GCC编译并测试通过):
三、矩阵乘法问题
问题描述:输入A是m×n矩阵和B是n×p矩阵,它们的乘积AB是一个m×p矩阵,它的一个元素其中 1 ≤ i ≤ m, 1 ≤ j ≤ p,求乘积AB。
例子一:
输入:
矩阵A:1 0
3 5
矩阵B:3 1
2 2
输出(C=A*B):
矩阵C:3 1
19 13
例子二:
矩阵A:1 0 3 3
3 5 3 4
2 4 1 5
7 5 2 0
矩阵B:3 1 3 4
2 2 1 5
3 4 4 2
6 4 9 4
输出(C=A*B):
矩阵C:30 25 42 22
52 41 62 59
47 34 59 50
37 25 34 57
版本一:暴力求解法
原理:利用矩阵乘法定义直接求解。
时间复杂度:函数包含三层嵌套循环,T(m,n,p) = O(m*n*p),若矩阵A,B为n阶方阵,则T(n) = O(n^3)。
体会:从循环嵌套情况来看可以直到时间复杂度很大。
C语言实现核心代码:
/** * @file square_matrix_multiply_by_directly_caculation.c * @brief caculate matrix C = A*B;A:M*N B:N*P C:M*P * with directly caculation method * @author chenxilinsidney * @version 1.0 * @date 2015-01-05 */ #include <stdio.h> #include <stdlib.h> // #define NDEBUG #include <assert.h> #include "memory.h" // #define NDBG_PRINT #include "debug_print.h" typedef int TYPE; /// matrix size #define M 2 #define N 2 #define P 2 /// matrix A, B, C TYPE matrix_a[M*N] = {1, 0, 3, 5}; TYPE matrix_b[N*P] = {3, 1, 2, 2}; TYPE matrix_c[M*P] = {0}; /** * @brief get matrix C from C = A*B * * @param matrix_A matrix input A * @param matrix_B matrix input B * @param matrix_C matrix output C * @param size_M matrix A size:M * N * @param size_N matrix B size:N * P * @param size_P matrix C size:M * P */ void square_matrix_multiply(TYPE* matrix_A, TYPE* matrix_B, TYPE* matrix_C, TYPE size_M, TYPE size_N, TYPE size_P) { TYPE i, j, k; for (i = 0; i < size_M; i++) { for (j = 0; j < size_P; j++) { matrix_C[i * size_P + j] = 0; for (k = 0; k < size_N; k++) { matrix_C[i * size_P + j] += matrix_A[i * size_N + k] * matrix_B[k * size_P + j]; } } } return; } int main(void) { /// square matrix multiply square_matrix_multiply(matrix_a, matrix_b, matrix_c, M, N, P); /// output result TYPE i, j; printf("matrix c = \n"); for (i = 0; i < M; i++) { for (j = 0; j < P; j++) { printf("%3d ", matrix_c[i * P + j]); } printf("\n"); } return EXIT_SUCCESS; }
版本二:分治法
原理:这里假定A,B是n阶方阵,并且n为2的幂。这样就可以把n*n矩阵分为4个相同大小的(n/2)*(n/2)矩阵,从而把问题分解为4个规模更小的子问题,再利用矩阵乘法公式求解。
时间复杂度:考虑递归过程中涉及到矩阵相加,递归方程:T(n) = 8T(n/2) + O(n^2), n > 1;T(n) = 1,n=1,计算得时间复杂度T(n)=O(n^3)
体会:由於程序使用一维数组存储矩阵,所以在构造小的4个矩阵时,由于这4个矩阵在原来的数组中是不连续的,采用了方法是构造新的数组存储这些矩阵,再根据4个子矩阵求解该矩阵。另外一个方法是通过修改指针索引来获得相应的子矩阵,在递归输入的再传入各矩阵行长度值同时修改头指针位置即可。
C语言实现核心代码:
/** * @file square_matrix_multiply_by_divide_conquer.c * @brief caculate matrix C = A*B;A:M*N B:N*P C:M*P * with divide and conquer method. * @author chenxilinsidney * @version 1.0 * @date 2015-01-05 */ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <float.h> #include <string.h> #define NDEBUG #include <assert.h> #include "memory.h" #define NDBG_PRINT #include "debug_print.h" typedef int TYPE; /// matrix size #define N 4 /// matrix A, B, C TYPE matrix_a[N*N] = {1, 0, 3, 3, 3, 5, 3, 4, 2, 4, 1, 5, 7, 5, 2, 0}; TYPE matrix_b[N*N] = {3, 1, 3, 4, 2, 2, 1, 5, 3, 4, 4, 2, 6, 4, 9, 4}; TYPE matrix_c[N*N] = {0}; /** * @brief get matrix C from C = A*B * * @param matrix_A matrix input A * @param matrix_B matrix input B * @param matrix_C matrix output C * @param size_N matrix C size:N * N * matrix B size:N * N * matrix C size:N * N */ void square_matrix_multiply(TYPE* matrix_A, TYPE* matrix_B, TYPE* matrix_C, TYPE size_N) { DEBUG_PRINT_STRING("In recursion now.\n"); /// demand size_N = 2^N assert((log2((double)size_N) - (int)log2((double)size_N)) <= FLT_EPSILON); DEBUG_PRINT_VALUE("%d", size_N); if(size_N == 1) { *matrix_C = *matrix_A * *matrix_B; DEBUG_PRINT_VALUE("%d", *matrix_C); } else { /// get four smaller size matrix for each A,B matrix TYPE size_half_N = (unsigned)size_N >> 1; DEBUG_PRINT_VALUE("%d", size_half_N); TYPE i, j; /// malloc memory TYPE* matrix_A11 = SMALLOC(size_half_N * size_half_N, TYPE); TYPE* matrix_A12 = SMALLOC(size_half_N * size_half_N, TYPE); TYPE* matrix_A21 = SMALLOC(size_half_N * size_half_N, TYPE); TYPE* matrix_A22 = SMALLOC(size_half_N * size_half_N, TYPE); TYPE* matrix_B11 = SMALLOC(size_half_N * size_half_N, TYPE); TYPE* matrix_B12 = SMALLOC(size_half_N * size_half_N, TYPE); TYPE* matrix_B21 = SMALLOC(size_half_N * size_half_N, TYPE); TYPE* matrix_B22 = SMALLOC(size_half_N * size_half_N, TYPE); /// copy data for(i = 0; i < size_half_N; i++) { memcpy(matrix_A11 + i * size_half_N, matrix_A + i * size_N, size_half_N * sizeof(TYPE)); memcpy(matrix_A12 + i * size_half_N, matrix_A + size_half_N + i * size_N, size_half_N * sizeof(TYPE)); memcpy(matrix_A21 + i * size_half_N, matrix_A + (size_half_N + i) * size_N, size_half_N * sizeof(TYPE)); memcpy(matrix_A22 + i * size_half_N, matrix_A + size_half_N +(size_half_N + i) * size_N, size_half_N * sizeof(TYPE)); memcpy(matrix_B11 + i * size_half_N, matrix_B + i * size_N, size_half_N * sizeof(TYPE)); memcpy(matrix_B12 + i * size_half_N, matrix_B + size_half_N + i * size_N, size_half_N * sizeof(TYPE)); memcpy(matrix_B21 + i * size_half_N, matrix_B + (size_half_N + i) * size_N, size_half_N * sizeof(TYPE)); memcpy(matrix_B22 + i * size_half_N, matrix_B + size_half_N +(size_half_N + i) * size_N, size_half_N * sizeof(TYPE)); DEBUG_PRINT_VALUE("%d", *(matrix_A11 + i * size_half_N)); DEBUG_PRINT_VALUE("%d", *(matrix_A12 + i * size_half_N)); DEBUG_PRINT_VALUE("%d", *(matrix_A21 + i * size_half_N)); DEBUG_PRINT_VALUE("%d", *(matrix_A22 + i * size_half_N)); DEBUG_PRINT_VALUE("%d", *(matrix_B11 + i * size_half_N)); DEBUG_PRINT_VALUE("%d", *(matrix_B12 + i * size_half_N)); DEBUG_PRINT_VALUE("%d", *(matrix_B21 + i * size_half_N)); DEBUG_PRINT_VALUE("%d", *(matrix_B22 + i * size_half_N)); } /// get result for each smaller matrix multiply TYPE* matrix_C1 = SMALLOC(size_half_N * size_half_N, TYPE); TYPE* matrix_C2 = SMALLOC(size_half_N * size_half_N, TYPE); square_matrix_multiply(matrix_A11, matrix_B11, matrix_C1, size_half_N); square_matrix_multiply(matrix_A12, matrix_B21, matrix_C2, size_half_N); for(i = 0; i < size_half_N; i++) { for(j = 0; j < size_half_N; j++) { matrix_C[i * size_N + j] = matrix_C1[i * size_half_N + j] + matrix_C2[i * size_half_N + j]; } } square_matrix_multiply(matrix_A11, matrix_B12, matrix_C1, size_half_N); square_matrix_multiply(matrix_A12, matrix_B22, matrix_C2, size_half_N); for(i = 0; i < size_half_N; i++) { for(j = 0; j < size_half_N; j++) { matrix_C[i * size_N + size_half_N + j] = matrix_C1[i * size_half_N + j] + matrix_C2[i * size_half_N + j]; } } square_matrix_multiply(matrix_A21, matrix_B11, matrix_C1, size_half_N); square_matrix_multiply(matrix_A22, matrix_B21, matrix_C2, size_half_N); for(i = 0; i < size_half_N; i++) { for(j = 0; j < size_half_N; j++) { matrix_C[(size_half_N + i) * size_N + j] = matrix_C1[i * size_half_N + j] + matrix_C2[i * size_half_N + j]; } } square_matrix_multiply(matrix_A21, matrix_B12, matrix_C1, size_half_N); square_matrix_multiply(matrix_A22, matrix_B22, matrix_C2, size_half_N); for(i = 0; i < size_half_N; i++) { for(j = 0; j < size_half_N; j++) { matrix_C[(size_half_N + i) * size_N + size_half_N + j] = matrix_C1[i * size_half_N + j] + matrix_C2[i * size_half_N + j]; } } DEBUG_PRINT_VALUE("%d", *matrix_C); /// free memory SFREE(&matrix_A11); SFREE(&matrix_A12); SFREE(&matrix_A21); SFREE(&matrix_A22); SFREE(&matrix_B11); SFREE(&matrix_B12); SFREE(&matrix_B21); SFREE(&matrix_B22); SFREE(&matrix_C1); SFREE(&matrix_C2); } DEBUG_PRINT_STRING("Out recursion now.\n"); return; } int main(void) { /// square matrix multiply square_matrix_multiply(matrix_a, matrix_b, matrix_c, N); /// output result TYPE i, j; printf("matrix c = \n"); for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { printf("%3d ", matrix_c[i * N + j]); } printf("\n"); } return EXIT_SUCCESS; }
版本三:Strassen算法
原理:待学习补充。
时间复杂度:待学习补充。
体会:待学习补充。
C语言实现核心代码:
此问题详细源代码查看及下载(UBUNTU系统+ GCC编译并测试通过):
https://github.com/chenxilinsidney/funnycprogram/tree/master/introduction_to_algorithms/square_matrix_multiply
四、拓展阅读及深入学习
1.题目二考虑数组首尾相连的情况来求解最大子数组。包括其他类似问题。
参考资料:http://www.ahathinking.com/archives/120.html
2.题目二可以考虑全是负值的数组的返回值情况,上面算法都是返回最大的一个负值,若要返回0值,只需要个别变量初始化时置为1或最后加入判断即可。