「アルゴリズム導論」学習ノート-最大サブ配列(分治戦略、動的計画)
21345 ワード
一、分治策略
分治法の思想
原問題をいくつかの規模が小さいが原問題に類似するサブ問題に分解し,これらのサブ問題を再帰的に解き,その後,これらのサブ問題の解を統合して原問題の解を確立する.
再帰式
再帰式は、再帰式を用いることで、再帰式アルゴリズムの実行時間を自然に描画することができるため、分治方法と密接に関連している.
分治戦略では,各層の再帰に次の3つのステップを適用する問題を再帰的に解いた.
分解:問題をいくつかのサブ問題に分け、サブ問題の形式は元の問題と同じで、規模が小さいだけです.
解決:サブ問題を再帰的に解く.サブ問題の規模が十分に小さい場合は,再帰を停止し,直接解く.
≪マージ|Merge|emdw≫:サブ問題の解を元の問題の解に組み合わせます.
学習の心得.
再帰形式のアルゴリズムは典型的に分治思想に従っている.しかし、再帰アルゴリズムを用いて再帰を行う具体的な計算方法の違い(再帰方程式の違い)は、アルゴリズムの時間的複雑さの違いをもたらすため、アルゴリズムが再帰的な方法を使用しているとしても、アルゴリズムの時間的複雑さを減少させるより良いテクニックがあるとは限らない.再帰アルゴリズムは、より小さな規模のサブ問題にしか分けられない可能性がある.(分治思想はここでは明らかに表現されていない.サブ問題は複数ではなく、1つしかないからだが、問題の規模は確かに小さくなっている)、複数の規模がより小さいサブ問題に分解される可能性もある(分治思想の表現が明らかである).
二、最大サブ配列問題
問題の説明:与えられた配列Aに対して、Aのおよび最大の非空連続サブ配列を探す.
例:配列{2,4,−7,5,2,−1,2,−4,3}の最大連続サブ配列は{5,2,−1,2}であり、最大連続サブ配列の和は5+2−1+2=8である.
バージョン1:暴力の解法
原理:Aの各シード配列の和を計算し、最大のサブ配列を見つけます.
時間複雑度:ネストされたループ反復(各シード配列の開始位置について、各シード配列の終了位置の和を解く)が必要であるため、時間複雑度T(n)=O(n^2)となる.
体得:暴力求解法の時間複雑度を計算する際に,配列と組合せ式を用いて計算を補助することができる.また,アルゴリズムは実際には動的計画の考え方も用いられている:アルゴリズムの和を求める際に累積を採用し,前のステップの結果を格納し,後続のステップで直接使用し,重複労働を減らす.累積を用いなければ,時間複雑度はT(n)=O(n^3)となる.
C言語実現コアコード:
バージョン2:分治法
原理:分治の思想を採用し、手順は以下の通りである.
分解:配列をより小さい規模でできるだけ等しいサブ配列に分割し、ここで配列を中間位置からサブ配列に分割します.
解決:これらのサブ配列の最大サブ配列を分治法を用いて再帰的に解く.最小サブ配列の場合、配列の長さが1の場合、その要素の値が返されます.
マージ:2つのサブ配列にまたがる最も大きなサブ配列、すなわち、この隣接するサブ配列にまたがって合成された配列の中で点の最も大きなサブ配列を解き、2つのサブ配列のそれぞれの最も大きなサブ配列と比較し、選択し、最大者を求める.
時間複雑度:再帰アルゴリズムを用いて,T(n)=2 T(n/2)+O(n)(n>1)を実現し,T(n)=O(nlgn)を導いた.
体得:分治の思想は最も主要な考慮の2つの方面です:1.どのように問題をサブ問題に分解するか、ここには実際に多くの考えがあり、例えばデータが半分になっている.例えば、1つのデータを減らすだけで、合理的な区分が分治思想のメリットを体現することができる.2.2つのサブ問題をどのようにマージするか、ここでサブ問題が何であるかを明確にする必要があります.マージ後に解決される問題は現在の問題であり、上位レベルにとっても実際にはサブ問題です.
C言語実現コアコード:
バージョン3:ダイナミックプランニングアルゴリズム(線形時間)
原理:動的計画の思想を利用して、しかし前の状態を保存して、重なり合うサブ問題を解決して、状態の移行方程式と結びつけて解く.後進解析:n番目の要素a[n]を加えた配列a[0--n]の最大子配列と、前の配列a[0--n-1]を加えた最大子配列との関係を考慮して、状態遷移の方程式を見つける.それらの関係は次の2つの場合があります.
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]);
時間複雑度:1回の遍歴のみを設計する:T(n)=O(n).
体得:動的計画の方法は主にサブ問題と問題の間の転態転移関係を考慮し、状態転移方程式と初期条件を利用して解き、肝心なのは状態がどのように転移したのかを明確にし、このステップを分析すれば、アルゴリズムは比較的に実現しやすい.
学習:このアルゴリズムの1つに付属しているのは、データを1回だけスキャンし、要素が読み込まれて処理されると、記憶される必要がなくなります.したがって、配列がディスクまたはテープ上にある場合、配列の任意の部分をプライマリ・ストレージに格納する必要はありません.それだけでなく、任意の時点で、このアルゴリズムはすでに読み込まれたデータに対して最も大きなサブ配列(前の2つのアルゴリズムはこのような特性を持たない)を与えることができる.このような特性を持つアルゴリズムをオンラインアルゴリズムと呼ぶ.
C言語実現コアコード:
バージョン4:ダイナミックプランニングアルゴリズム最適化(線形時間)
バージョン3で得られた状態遷移方程式は以下の通りである.
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]の2つの状態を別々に記録することなく,判断文を1つにまとめることができる.
このバージョンは、累積を行う過程で、現在得られている和が負数であれば、これと次の累積で捨てて再クリアすべきであり、そうでなければ、この負数は次の和を減少させることを別の考え方で理解することもできます.
時間複雑度:1回の遍歴のみを設計する:T(n)=O(n).
体得:ループに前回の結果と現在の結果を格納する際に、必ずしも別々に格納してからループが終了した後に値を割り当てる必要はありません.彼らの関係をさらに簡略化し、前回の結果を利用して直接現在の結果として格納することができます.
C言語実現コアコード:
この問題の詳細ソースコードの表示とダウンロード(UBUNTUシステム+GCCコンパイルとテストに合格):https://github.com/chenxilinsidney/funnycprogram/blob/master/introduction_to_algorithms/maximum_subarray
三、行列乗算問題
質問説明:入力Aはm×n行列とBはn×p行列,それらの積ABはmである×1≦i≦m,1≦j≦pの元素の1つのp行列は,積ABを求める.
例1:
入力:
マトリックスA:1 0 3 5マトリックスB:3 1 2出力(C=A*B):
行列C:31 1 19 13
例2:
マトリックスA:1 0 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
バージョン1:暴力の解法
原理:マトリクス乗算定義を用いて直接解く.
時間複雑度:関数は3層のネストループを含み,T(m,n,p)=O(m*n*p),行列A,Bがn次方程式である場合,T(n)=O(n^3)である.
体得:ループネストの状況から見ると、時間の複雑さが大きいまで見ることができます.C言語実現コアコード:
バージョン2:分治法
原理:ここでは,A,Bはn次方程式であり,nは2のべき乗であると仮定する.これにより,n*n行列を4つの同じ大きさの(n/2)*(n/2)行列に分けることができ,問題を4つの規模の小さいサブ問題に分解し,行列乗算式を再利用して解くことができる.
時間複雑度:再帰過程における行列加算を考慮し,再帰方程式:T(n)=8 T(n/2)+O(n^2),n>1;T(n)=1,n=1,計算時間複雑度T(n)=O(n^3)
体得:プログラムは1次元配列格納マトリクスを用いるため,小さい4つのマトリクスを構築する際,この4つのマトリクスは元の配列では不連続であるため,新しい配列格納マトリクスを構築し,4つのサブマトリクスに基づいてこのマトリクスを解く方法を採用した.もう1つの方法は、ポインタインデックスを修正することによって対応するサブマトリクスを得ることであり、各マトリクス行の長さ値を再帰入力しながらヘッダポインタ位置を修正すればよい.C言語実現コアコード:
バージョン3:Strassenアルゴリズム
原理:学習補充を待つ.
時間の複雑さ:学習の補充を待つ.
体得:勉強しなければならない.C言語実現コアコード:
この問題の詳細ソースコードの表示とダウンロード(UBUNTUシステム+GCCコンパイルとテストに合格):https://github.com/chenxilinsidney/funnycprogram/tree/master/introduction_to_algorithms/square_matrix_multiply
四、読書を広げ、深く勉強する
1.問題2配列の先頭と末尾がつながっている場合を考慮して、最大のサブ配列を解く.他の類似の問題も含まれています.
参考資料:http://www.ahathinking.com/archives/120.html
2.タイトル2はすべて負の値の配列の戻り値を考慮することができ、上のアルゴリズムはすべて最大の負の値を返し、0の値を返すには、個別変数の初期化時に1または最後に判断を加えるだけでよい.
分治法の思想
原問題をいくつかの規模が小さいが原問題に類似するサブ問題に分解し,これらのサブ問題を再帰的に解き,その後,これらのサブ問題の解を統合して原問題の解を確立する.
再帰式
再帰式は、再帰式を用いることで、再帰式アルゴリズムの実行時間を自然に描画することができるため、分治方法と密接に関連している.
分治戦略では,各層の再帰に次の3つのステップを適用する問題を再帰的に解いた.
分解:問題をいくつかのサブ問題に分け、サブ問題の形式は元の問題と同じで、規模が小さいだけです.
解決:サブ問題を再帰的に解く.サブ問題の規模が十分に小さい場合は,再帰を停止し,直接解く.
≪マージ|Merge|emdw≫:サブ問題の解を元の問題の解に組み合わせます.
学習の心得.
再帰形式のアルゴリズムは典型的に分治思想に従っている.しかし、再帰アルゴリズムを用いて再帰を行う具体的な計算方法の違い(再帰方程式の違い)は、アルゴリズムの時間的複雑さの違いをもたらすため、アルゴリズムが再帰的な方法を使用しているとしても、アルゴリズムの時間的複雑さを減少させるより良いテクニックがあるとは限らない.再帰アルゴリズムは、より小さな規模のサブ問題にしか分けられない可能性がある.(分治思想はここでは明らかに表現されていない.サブ問題は複数ではなく、1つしかないからだが、問題の規模は確かに小さくなっている)、複数の規模がより小さいサブ問題に分解される可能性もある(分治思想の表現が明らかである).
二、最大サブ配列問題
問題の説明:与えられた配列Aに対して、Aのおよび最大の非空連続サブ配列を探す.
例:配列{2,4,−7,5,2,−1,2,−4,3}の最大連続サブ配列は{5,2,−1,2}であり、最大連続サブ配列の和は5+2−1+2=8である.
バージョン1:暴力の解法
原理: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;
}
}
}バージョン2:分治法
原理:分治の思想を採用し、手順は以下の通りである.
分解:配列をより小さい規模でできるだけ等しいサブ配列に分割し、ここで配列を中間位置からサブ配列に分割します.
解決:これらのサブ配列の最大サブ配列を分治法を用いて再帰的に解く.最小サブ配列の場合、配列の長さが1の場合、その要素の値が返されます.
マージ:2つのサブ配列にまたがる最も大きなサブ配列、すなわち、この隣接するサブ配列にまたがって合成された配列の中で点の最も大きなサブ配列を解き、2つのサブ配列のそれぞれの最も大きなサブ配列と比較し、選択し、最大者を求める.
時間複雑度:再帰アルゴリズムを用いて,T(n)=2 T(n/2)+O(n)(n>1)を実現し,T(n)=O(nlgn)を導いた.
体得:分治の思想は最も主要な考慮の2つの方面です:1.どのように問題をサブ問題に分解するか、ここには実際に多くの考えがあり、例えばデータが半分になっている.例えば、1つのデータを減らすだけで、合理的な区分が分治思想のメリットを体現することができる.2.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.
");
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.
");
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.
");
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.
");
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.
");
return data_cross;
}
}
}
バージョン3:ダイナミックプランニングアルゴリズム(線形時間)
原理:動的計画の思想を利用して、しかし前の状態を保存して、重なり合うサブ問題を解決して、状態の移行方程式と結びつけて解く.後進解析:n番目の要素a[n]を加えた配列a[0--n]の最大子配列と、前の配列a[0--n-1]を加えた最大子配列との関係を考慮して、状態遷移の方程式を見つける.それらの関係は次の2つの場合があります.
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]);
時間複雑度:1回の遍歴のみを設計する:T(n)=O(n).
体得:動的計画の方法は主にサブ問題と問題の間の転態転移関係を考慮し、状態転移方程式と初期条件を利用して解き、肝心なのは状態がどのように転移したのかを明確にし、このステップを分析すれば、アルゴリズムは比較的に実現しやすい.
学習:このアルゴリズムの1つに付属しているのは、データを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 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;
}バージョン4:ダイナミックプランニングアルゴリズム最適化(線形時間)
バージョン3で得られた状態遷移方程式は以下の通りである.
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]の2つの状態を別々に記録することなく,判断文を1つにまとめることができる.
このバージョンは、累積を行う過程で、現在得られている和が負数であれば、これと次の累積で捨てて再クリアすべきであり、そうでなければ、この負数は次の和を減少させることを別の考え方で理解することもできます.
時間複雑度:1回の遍歴のみを設計する: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コンパイルとテストに合格):https://github.com/chenxilinsidney/funnycprogram/blob/master/introduction_to_algorithms/maximum_subarray
三、行列乗算問題
質問説明:入力Aはm×n行列とBはn×p行列,それらの積ABはmである×1≦i≦m,1≦j≦pの元素の1つのp行列は,積ABを求める.
例1:
入力:
マトリックスA:1 0 3 5マトリックスB:3 1 2出力(C=A*B):
行列C:31 1 19 13
例2:
マトリックスA:1 0 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
バージョン1:暴力の解法
原理:マトリクス乗算定義を用いて直接解く.
時間複雑度:関数は3層のネストループを含み,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
#include
// #define NDEBUG
#include
#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 =
");
for (i = 0; i < M; i++) {
for (j = 0; j < P; j++) {
printf("%3d ", matrix_c[i * P + j]);
}
printf("
");
}
return EXIT_SUCCESS;
} バージョン2:分治法
原理:ここでは,A,Bはn次方程式であり,nは2のべき乗であると仮定する.これにより,n*n行列を4つの同じ大きさの(n/2)*(n/2)行列に分けることができ,問題を4つの規模の小さいサブ問題に分解し,行列乗算式を再利用して解くことができる.
時間複雑度:再帰過程における行列加算を考慮し,再帰方程式:T(n)=8 T(n/2)+O(n^2),n>1;T(n)=1,n=1,計算時間複雑度T(n)=O(n^3)
体得:プログラムは1次元配列格納マトリクスを用いるため,小さい4つのマトリクスを構築する際,この4つのマトリクスは元の配列では不連続であるため,新しい配列格納マトリクスを構築し,4つのサブマトリクスに基づいてこのマトリクスを解く方法を採用した.もう1つの方法は、ポインタインデックスを修正することによって対応するサブマトリクスを得ることであり、各マトリクス行の長さ値を再帰入力しながらヘッダポインタ位置を修正すればよい.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
#include
#include
#include
#include
#define NDEBUG
#include
#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.
");
/// 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.
");
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 =
");
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
printf("%3d ", matrix_c[i * N + j]);
}
printf("
");
}
return EXIT_SUCCESS;
} バージョン3:Strassenアルゴリズム
原理:学習補充を待つ.
時間の複雑さ:学習の補充を待つ.
体得:勉強しなければならない.C言語実現コアコード:
この問題の詳細ソースコードの表示とダウンロード(UBUNTUシステム+GCCコンパイルとテストに合格):https://github.com/chenxilinsidney/funnycprogram/tree/master/introduction_to_algorithms/square_matrix_multiply
四、読書を広げ、深く勉強する
1.問題2配列の先頭と末尾がつながっている場合を考慮して、最大のサブ配列を解く.他の類似の問題も含まれています.
参考資料:http://www.ahathinking.com/archives/120.html
2.タイトル2はすべて負の値の配列の戻り値を考慮することができ、上のアルゴリズムはすべて最大の負の値を返し、0の値を返すには、個別変数の初期化時に1または最後に判断を加えるだけでよい.