経路数を計算する
それではおためし。単位行列Aとfig-2に示した隣接行列Bを用意し、A = A * Bを三回繰り返してB^3を求めます。
list-4 01_adjacency\01_adjacency.cpp
#include "sparse/matrix.h"
#include "sparse/multiplies.h"
#include "sparse/fill_diagonal.h"
#include <iostream> // cout, endl
int main() {
using namespace sparse;
matrix<float> A;
fill_diagonal(A, 4, 1.0f); // 4x4の単位行列
matrix<float> B; // 隣接行列
B.insert(0, 0, 1.0f);
B.insert(0, 1, 1.0f);
B.insert(0, 2, 1.0f);
B.insert(1, 2, 1.0f);
B.insert(2, 3, 1.0f);
B.insert(3, 1, 1.0f);
B.finalize();
// A = A * B を3回行う
A = multiplies(A, B, 4, 4, 4);
A = multiplies(A, B, 4, 4, 4);
A = multiplies(A, B, 4, 4, 4);
A.enumerate([](int row, int col, float val) {
std::cout << '[' << row << ',' << col << "] "
<< val << std::endl;});
std::cout << std::endl;
}
隣接行列を3回かけたときの成分（i,j）が何を意味するかというと、「頂点iから頂点jへ3つの辺をたどって到達する経路の数」なんです。例えば「[0,3] 2」とありますから、島0から3本の橋を渡って島3にたどり着くルートは2通りです。確かに0→0→2→3と0→1→2→3の2通りですね。
経路の最短／最長距離を求める
さて、ここから面白くなりますよ。
K行L列の行列AとL行M列の行列Bの掛け算関数multipliesは、その成分（row,col）を以下のアルゴリズムで求めています。
val = 0;
for ( i = 0; i < L; ++i ) {
val = val + A(row, i) * B(i, col);
for-loopの中で乗算と加算が1つずつ現れてますね。この乗算と加算を引数で与えるmultipliesを用意します。
list-5 02_adjacency\02_adjacency.cpp
template<typename T, typename U, typename MulOp, typename AddOp>
matrix<T> multiplies(const matrix<T>& A, const matrix<U>& B,
int k, int l, int m,
MulOp nmul, AddOp nadd) {
matrix<T> out;
// out[row, col] = Σ(A[row,i] * B[i,col]
for ( int row = 0; row < k; ++row ) {
for ( int col = 0; col < m; ++col ) {
T s; T* sp = nullptr;
T v; T* vp;
T a; T* ap;
U b; U* bp;
for ( int i = 0; i < l; ++i ) {
ap = A.at(row, i, &a); // a = A[row,i]
bp = B.at(i, col, &b); // b = B[i,col]
vp = nmul(ap, bp, &v); // v = a * b
sp = nadd(sp, vp, &s); // s = s + v;
}
if ( sp != nullptr ) {
out.insert(row, col, s);
}
}
}
out.finalize();
return out;
}
ここで与える乗算／加算関数MulOo／AddOpは、以下のプロトタイプに従うこととします。
T* operation(const T* a, const T* b, T* c)
*aと*bを演算子の両辺とし、演算結果を*cにセットしてcを返します。ただし左辺／右辺が無効値であったときはa,bはnullptrですし、演算結果が無効値ならばnullptrを返します。
ここで乗算と加算を以下のように定義します:
- 乗算：両辺の和を返す。ただしどちらかが無効値であるなら無効値を返す。
- 加算：両辺の「大きくない方」を返す。ただし一方が無効値なら他方を、両辺とも無効値なら無効値を返す。
隣接行列にも手を加え、各成分を（1ではなく）橋の長さに置き換えます。
list-6 03_shortestpath\03_shortestpath.cpp
#include "sparse/matrix.h"
#include "sparse/multiplies.h"
#include "sparse/fill_diagonal.h"
// 乗算 <*>
// 0 <*> 0 = 0
// x <*> 0 = 0
// 0 <*> x = 0
// x <*> y = x+y
template<typename T>
struct nullable_multiplies {
T* operator()(const T* ap, const T* bp, T* result) const {
if ( ap == nullptr ) return nullptr;
if ( bp == nullptr ) return nullptr;
*result = (*ap) + (*bp);
return result;
}
};
// 加算 <+>
// 0 <+> 0 = 0
// 0 <+> x = x
// x <+> 0 = x
// x <+> y = min(x,y)
template<typename T>
struct nullable_min {
T* operator()(const T* ap, const T* bp, T* result) const {
if ( ap == nullptr && bp == nullptr ) return nullptr;
if ( ap == nullptr ) *result = *bp;
else *result = ( bp == nullptr ) ? *ap : ((*ap) < (*bp) ? (*ap) : (*bp));
return result;
}
};
// 加算 <+>
// 0 <+> 0 = 0
// 0 <+> x = x
// x <+> 0 = x
// x <+> y = max(x,y)
template<typename T>
struct nullable_max {
T* operator()(const T* ap, const T* bp, T* result) const {
if ( ap == nullptr && bp == nullptr ) return nullptr;
if ( ap == nullptr ) *result = *bp;
else *result = ( bp == nullptr ) ? *ap : ((*ap) > (*bp) ? (*ap) : (*bp));
return result;
}
};
#include <iostream> // cout, endl
int main() {
using namespace sparse;
matrix<float> A;
fill_diagonal(A, 4, 0.0f);
matrix<float> B;
B.insert(0, 0, 0.8f);
B.insert(0, 1, 1.5f);
B.insert(0, 2, 1.8f);
B.insert(1, 2, 1.1f);
B.insert(2, 3, 2.1f);
B.insert(3, 1, 2.2f);
B.finalize();
nullable_multiplies<float> mulop;
nullable_min<float> addop;
A = multiplies(A, B, 4, 4, 4, mulop, addop);
A = multiplies(A, B, 4, 4, 4, mulop, addop);
A = multiplies(A, B, 4, 4, 4, mulop, addop);
A.enumerate([](int row, int col, float val) {
std::cout << '[' << row << ',' << col << "] "
<< val << std::endl;});
std::cout << std::endl;
}
この演算で得られる成分（i,j）は、「3本の橋をたどって島iから島jに到達するときの最短路長」なんです。加算の定義を「小さくない方」に差し替えれば最長路長が求まります。
