Knuth Morris Pratt (KMP) algorithm

Knuth, Morris, Pratt によって考案された KMP アルゴリズムは、T[n], P[m] に対して最悪の場合でも n + m 回の比較で文字列照合を行うことができます（O(n + m)）。

Brute force アルゴリズムでは、不一致 T[i] ≠ P[j] が検出されると、i = i - j + 1 として i のポインタをリセットしていました。KMP アルゴリズムは、不一致 T[i] ≠ P[j] を検出した時点で、テキストとパタンがどこまで一致しているかが分かっていることを有効に利用します。つまり不一致 T[i] ≠ P[j] を検出した時点で 1 ≦ k ≦ j について T[i - k] = P[j - k] ということが既に分かっているのです（下図参照）。

具体的な例をいくつか見てみましょう。

Case 1:

右図の例において、不一致 T[i] ≠ P[j] が検出された時点で、T[i - j] = b, T[i - j + 1] = a, ... T[i - 1] = a ということが既に分かっています。
つまり、P[0] = b ということから ( T[1] = T[2] = T[3] ) ≠ P[0] ということが既に分かっているので、Brute force のように i のポインタを i - j + 1 に後戻りさせるのは明らかに無駄な処理を発生させます。
この例の場合、i はそのまま、j = 0 としてパタンのチェックを再開することができます。

Case 2:

右図の例において、不一致 T[i] ≠ P[j] が検出された時点で、T[0] = T[2] = a, T[1] = T[3] = b ということが既に分かっています。
P[0] ≠ P[1] = T[1] ということが既に分かっているので、Brute force のように i のポインタを i - j + 1 に後戻りさせるのは無駄になります。
では、テキトウに i を１つだけ後戻りさせてみましょう（i = i - 1）。しかしこれは、P[0] = T[2], P[1] = T[3] という事実、つまりパタンを見逃してしまいます。これは、少なくとも i を１つしか後戻りさせなかったら、パタンを見逃してしまうことを意味します。
この例の場合、i を２つ後戻りさせ（i = i - 2）、チェックを再開することが正しいのです。

KMP アルゴリズムでは、この「P[j] の位置で不一致が検出された時、i をどこまで後戻りさせるか」をテーブル N として記録しておき有効利用します。N の値は P にのみ依存することから、N は照合の処理をする前にあらかじめ計算しておくことができます。

j > 0 に対する N[j] の値は、k < j かつ パタンの先頭から j 文字の部分列において、最初の k 文字と最後の k 文字が一致するような最大の k です。

このパタンをあるテキストと照合してみましょう。

P[8] で不一致が検出されたとすると、k = N[8] = 5 なので、i を 5 つ後戻りさせ照合を再開すればよいことになります。

しかし、上で述べたように、 1 ≦ k ≦ j について T[i - k] = P[j - k] ということが既に分かっているので、元の i の位置から再開することができます。従って、N[j] は「i をどこまで後戻りさせるか」と表現するよりは「j をどこから再開するか」と表現した方が直感的に分かりやすいと思います。

ただし、P[0] で不一致が検出された場合は、無条件に i = i + 1, j = 0 となります。この例外処理を行うために、以下のプログラムのように N[0] = -1 としておくと便利です。

void initNext( string P, int m, int N[]){
    int i, j;
    N[0] = -1;
    for ( i = 0, j = -1; i < m; i++, j++, N[i] = j )
    while( ( j >= 0 ) && ( P[i] != P[j] ) ) j = N[j];
}

int KMPSearch( string T, int n, string P, int m ){
    int N[MAX];
    initNext( P, m, N );
    
    int i, j;
    i = j = 0;
    while ( i < n ){
        while( i < n && j < m ){
            while ( j >= 0 && T[i] != P[j] ) j = N[j];
            i++; j++;
        }
        if ( j == m ) {
            print( i - j );
            j = N[j];
        }
    }
}

上記のプログラムは、若干巧妙なので、P[0] で不一致が検出されたときの処理を明示的に書いたプログラム例を以下に示します。

void initNext( string P, int m, int N[]){
    int i, j;
    N[0] = -1;
    for ( i = 0, j = -1; i < m; i++, j++, N[i] = j )
    while( ( j >= 0 ) && ( P[i] != P[j] ) ) j = N[j];
}

int KMPSearch( string T, int n, string P, int m ){
    int N[MAX];
    initNext( P, m, N );
    
    int i, j;
    i = j = 0;
    
    while ( i < n ){
        while( i < n && j < m ){
            if ( j == 0 ){
                if ( T[i] != P[j] ){ j = 0; i++; }
                else { i++; j++; }
            } else {
                if ( T[i] != P[j] ){ j = N[j]; }
                else { i++; j++; }
            }
        }
        if ( j == m ) {
            printf("P appers at %d\n", i - j); j = N[j];
        }
    }
}

Brute force algorithm 力任せな文字列照合アルゴリズム

文字列照合を行うための最も自明（素朴）な方法は、テキスト内のパタンが存在する可能性のある全ての位置を順番にチェックするものです。テキストの位置 i が 0 から n - m について、T[i + j] = P[j] ( 0 ≦ j ≦ m - 1 ) を満たすかをチェックします。

T[i] ≠ P[j] が検出されたときに、i = i - j + 1 (次の i の位置) j = 0 (パタンの先頭)にリセットします。（下図参照） 文字列照合アルゴリズムでは、文字を比較する回数により計算効率を評価します。Brute force algorithm は最悪の場合 m( n - m + 1) 回の比較が必要なので、m が n に比べて十分小さいと考えるとオーダーは O(nm) となります。例えば、テキスト "000000001" に対して、パタン "001" の照合がこの最悪のケースを引き起こします（試してみましょう）。 しかし、実際のテキストエディタなどのアプリケーションで、この最悪のケースが発生することはないと考えて良いでしょう。なぜなら、T[i + j] と P[j] をチェックする内側のループは高確率で j = 0 のときに終了することが考えられるからです。従って、文字コードが多い実際の文字列に対しては、Brute force algorithm で十分実用的なのです。 しかし一方で、先ほどの例のように、バイナリ（S = {0, 1}）のテキストなど、ある特定の文字列に対しては非常に効率が悪くなります。

Brute force algorithm の流れ

int bruteForceSearch( string T, int n, string P, int m ){
    int i, j;
    i = j = 0;
    while ( i <= n - m ){
        while( i < n && j < m ){
            if ( T[i] != P[j] ){
                i = i - j + 1; j = 0;
            } else { i++; j++; }
        }
        if ( j == m ) {
            match( i - j );
            i = i - j + 1; j = 0;
        }
    }
}

Warshall Floyd's algorithm ワーシャル フロイド

All Pairs Shortest Path (APSP) problem (全対最短経路問題) は、グラフ G(V, E) に対して、G に含まれるノードの全てのペアの最短経路（距離）を求める問題です。G に負の重みのあるエッジがなければ、この問題は各ノードを起点として Dijkstra's Algorithm を |V| 回行うことによって解くことができます。この方法のオーダーは O(|V|³) または O(|V|(|E| + |V|)log|V|) となります。

Floyd's algorithm は APSP 問題を O(|V|³) で解くことができます。しかも Floyd's algorithm は、G に負の閉路がない限り、G に負の重みを持つエッジが存在しても正しく動作します。負の閉路とは、その閉路を成すエッジの重みの合計が負となっている閉路です（負の閉路をもつとき G に最短経路は存在しません）。さらに、Floyd's algorithm は、G に負の閉路が存在するか否かを判定することができます。アルゴリズムが終了した時点で、G のあるノード v からノード v (それ自身) への最短コストが負になっていれば、G に負の閉路が存在します。

Floyd's algorithm は、ノードの全てのペア [i, j] に対する最短コストを与える配列 A[i, j] を以下の方法で求めます。まず、G(V, E) のノードを {1, 2, 3, ... |V|} とします。

A^k[i, j] をノード V^k = {1, 2, 3, ... k} のみを経由するノード i とノード j の最短コスト、その時の経路の１つを P^k[i, j] とします。Floyd's algorithm は A^k-1 を利用して A^k を求めるということを k = {1, 2, 3, ... |V|} に対して行い、A^|V| つまり A[i, j] (および P[i, j])を求めます（Dynamic Programming (動的計画法)）。Floyd's algorithm は以下のように数学的帰納法によって証明・説明することができます。

1.	A0 ( k = 0 ) の時の証明 A0[i, j] は、ノード i とノード j 以外に経由ノードが存在しないので、単にノード i とノード j を繋ぐエッジの重みになります。つまり、入力されたグラフ G に対応させます。G のノード i からノード j に向かって重み d のエッジがあれば A[i, j] = d とし、エッジがなければ A[i, j] = ∞ とます。これを全対の [i, j] に対して施します。ただし A[i, i] = 0 とします。
	この時点で、A0[i, j] がノード i からノード j への最短コストであることは明白なので、帰納法の基礎が確立できました。
2.	Ak ( k > 0 ) の時の証明 k = {1, 2, 3, ... |V|} について、Ak-1 に基づいて Ak が正しく求められるか？を示します。
	Pk[i, j] がノード k を経由しない場合とする場合を考えます。 (i) Pk[i, j] がノード k を経由しない場合、Pk[i, j] は端点であるノード i とノード j 以外では Vk-1 = {1, 2, 3, ... k-1} に属するノードしか通らないので、Pk-1[i, j] と等しいと考えることができます。よって、この場合は Ak[i, j] = Ak-1[i, j] となります。 (ii) Pk[i, j] がノード k を経由する場合、Pk[i, j] はノード k によって　i → k と k → j の２つの部分路に分けられます。この２つの部分路はどちらも Vk-1 = {1, 2, 3, ... k-1} に属するノードしか経由しません（下図参照）。従ってノード k を経由する最短経路の部分路は Pk-1[i, k] と Pk-1[k, j] となります。つまり、Ak[i, j] = Ak-1[i, k] + Ak-1[k, j] となります。
	(i)(ii) より、i, j の全ての対に対して、      Ak[i, j] = min( Ak-1[i, j], Ak-1[i, k] + Ak-1[k, j]) を行い Ak を求めます。

このアルゴリズムを実装する時に、A^k[i, j] を求めるために、３次元のメモリ空間を必要と思うかもしれません。つまり A^k[i, j] の計算中に A^k-1[i, j] の値が変化してしまうのでは？と考えてしまいます。しかし、これは以下の理由により問題はありません。

A^k[k, k] = 0 より、

     A^k[i, k] = min( A^k-1[i, k], A^k-1[i, k] + A^k-1[k, k] ) = A^k-1[i, k]
     A^k[k, j] = min( A^k-1[k, j], A^k-1[k, j] + A^k-1[k, k] ) = A^k-1[k, j]

よって、ノードの全てのペア [i, j] について A^k-1[i, j] を A^k[i, j] に上書きしても不都合はありません。

int N;
int A[MAX][MAX];

void floyd(){
    for ( int k = 0; k < N; k++ ){
        for ( int i = 0; i < N; i++ ){
            for ( int j = 0; j < N; j++ ){
                A[i][j] = min( A[i][j], A[i][k] + A[k][j] );
            }
        }
    }
}

Priority First Search 順位優先探索

Dijkstra's Algorithm のオーダーは、O(|V|^2) (ノードの数の２乗)です。これは、u を S に追加するためのループ（|V| 回）の中で、u を選択するためのループ（|V| 回）を行っていることから、簡単に導くことができます。このオーダーは、隣接行列で実装しても、隣接リストで実装しても同じですが、疎なグラフに対しては隣接リストで実装すべきです。（以下のテーブルを参照して下さい）。

タイプ	演算回数	オーダー
隣接行列	|V|^2 + |V|^2	|V|^2
隣接リスト	|V|^2 + |E|	|V|^2

Dijkstra's Algorithm は、隣接リストによる実装と、Priority Queue (Heap) を応用することによって、効率化することができます。このアルゴリズムは、Priority First Search (順位優先探索) とも呼ばれ、以下のように実装します。

	グラフ G(V, E) のノード全体の集合を V、 source をスタート地点、 S を「 source から S 内の各ノード v への最短経路が必ず S 内にあるような、最短経路の部分的な解の集合」とします（Dijkstra's Algorithm 参照）。
	各計算ステップにおいて、V - S の各ノード v について、S 内のノードのみを経由した source から v への最短経路のコスト(距離)を d[v] とします。また、各ノード間のコスト（距離） (x, y) は、行列 M[x][y] に記録されているものとします。
1.	初期状態で、S は空です。 source に対して    d[s] = 0 source 以外の V に属する全てのノード v に対して    d[v] = ∞ と初期化します。
	d[v] をキーとして、V - S に属するノードを minimum heap として構築します。これを heap とします。
2.	heap から d[u] が最小であるノード u を取り出します。   u = heap.extractMin()
	u を S に追加すると同時に、u に隣接しかつ V - S に属する全てのノード v に対する値を以下のように更新します。
	d[v] = min( d[v], d[u] + M[u][v] )     d[v] が更新された場合は、heap を更新します。         heap.update( v )　(v から upheap を施します)
	この処理を S = V となるまで繰り返します。

Priority First Search のオーダーは、O((|V| + |E|)log|V|) となります。u を heap から取り出すために |V|log|V|、d[v] を更新するために |E|log|V| の演算回数が必要です。

Priority First Search のプログラム例を以下に示します。

class Graph{
    public:
    vector<vector<int> > adjList;
    vector<vector<int>::iterator> pos;
    
    Graph(){}
    Graph( int size ){ resize( size );}
    
    void resize( int size ){
        adjList.resize( size ); pos.resize( size );
        for ( int i = 0; i < size; i++ ) adjList[i].clear();
    }
    
    int size(){ return adjList.size(); }
    
    void insert( int i, int j ){
        adjList[i].push_back( j );
    }
    
    int next( int i ){
        if ( pos[i] == adjList[i].end() ) return -1;
        return *pos[i]++;
    }
    
    void reset( int i ){ pos[i] = adjList[i].begin(); }
    void reset(){ for( int i = 0; i < size(); i++ ) reset(i); }
};

class Heap{
    public:
    int size;
    int H[MAX + 1]; // heap buffer
    int index[MAX + 1]; // index of v in the heap buffer
    int *key; // keys for heap condition
    
    Heap(){}
    Heap( int size, int d[] ): size(size){
        key = d;
    }
    
    void construct(){
        for ( int i = 0; i < size; i++ ){
            H[i + 1] = i; // heap is 1 base not 0
            index[i] = i + 1;
        }
        for ( int i = size / 2; i >= 1; i-- ) downheap( i );
    }
    
    int extractMin(){
        int v = H[1];
        H[1] = H[size];
        index[ H[size] ] = 1;
        size--;
        downheap(1);
        return v;
    }
    
    int update( int v ){
        upheap( index[v] );
    }
    
    private:
    void downheap( int k ){
        int j, v;
        v = H[k];
        while( k <= size/2 ){
            j = k + k;
            if ( j < size && key[ H[j] ] > key[ H[j+1] ] ) j++;
            if ( key[v] <= key[ H[j] ] ) break;
            H[k] = H[j];
            index[ H[j] ] = k;
            k = j;
        }
        H[k] = v;
        index[ v ] = k;
    }
    
    void upheap( int k){
        int v;
        v = H[k];
        while ( k > 1 && key[ H[k/2] ] >= key[v] ){
            H[k] = H[k/2];
            index[ H[k/2] ] = k;
            k = k/2;
        }
        H[k] = v;
        index[v] = k;
    }
};

int size;
int M[MAX][MAX];
Graph graph;

void dijkstra( int s, int d[] ){
    bool visited[MAX]; // S に属するノードは true
    for ( int i = 0; i < size; i++ ){
        d[i] = INT_MAX;
        visited[i] = false;
    }
    
    d[s] = 0; // 最初に s が u として選ばれる
    
    Heap heap = Heap( size, d );
    heap.construct();
    
    graph.reset();
    
    while( heap.size >= 1 ){
        int u = heap.extractMin();
        
        visited[u] = true; // u を S に追加
        
        int v;
        while ( ( v = graph.next( u ) ) != -1 ){
            if ( visited[v] ) continue;
            if ( d[u] + M[u][v] < d[v] ){
                d[v] =  d[u] + M[u][v];
                heap.update( v );
            }
        }
    }
}

このプログラム例ですと、自作の Heap を使っているので、プログラムが分かりづらくなってしまいました。私は、普段は自作の Heap を使うことはなく、以下のように既存の Priority queue に候補となるノードを突っ込んでいくアルゴリズムを実装しています。

class Graph{
    public:
    vector<vector<int> > adjList;
    vector<vector<int>::iterator> pos;
    
    Graph(){}
    Graph( int size ){ resize( size );}
    
    void resize( int size ){
        adjList.resize( size ); pos.resize( size );
        for ( int i = 0; i < size; i++ ) adjList[i].clear();
    }
    
    int size(){ return adjList.size(); }
    
    void insert( int i, int j ){
        adjList[i].push_back( j );
    }
    
    int next( int i ){
        if ( pos[i] == adjList[i].end() ) return -1;
        return *pos[i]++;
    }
    
    void reset( int i ){ pos[i] = adjList[i].begin(); }
    void reset(){ for( int i = 0; i < size(); i++ ) reset(i); }
};

class Node{
    public:
    int id, cost;
    Node(){}
    Node( int id, int cost ): id(id), cost(cost){}
    
    bool operator < ( const Node &n ) const{
        return cost > n.cost;
    }
};

int size;
int M[MAX][MAX];
Graph graph;

void dijkstra( int s, int d[] ){
    bool visited[MAX]; // S に属するノードは true
    for ( int i = 0; i < size; i++ ){
        d[i] = INT_MAX;
        visited[i] = false;
    }
    
    d[s] = 0;
    
    priority_queue<Node> PQ;
    PQ.push( Node( s, 0 ) );// 最初に s が u として選ばれる
    
    graph.reset();
    
    while( !PQ.empty() ){
        int u = PQ.top().id; PQ.pop();
        
        visited[u] = true; // u を S に追加
        
        int v;
        while ( ( v = graph.next( u ) ) != -1 ){
            if ( visited[v] ) continue;
            if ( d[u] + M[u][v] < d[v] ){
                d[v] =  d[u] + M[u][v];
                PQ.push( Node( v, d[v] ) );
            }
        }
    }
}

純粋に Heap を用いるのと比べて効率はやや劣りますが（なぜなら priority queue の中に余分なゴミが溜まっていくから）、オーダーがたかだか O((|V| + |E|)log|E|) になるだけなので、私の経験上、実用には問題ないようです。

String matching problem 文字列照合

ある文字列（テキスト）に対して、指定された照合文字列（パタン）の存在を発見するための機能は、多くのテキストエディタに実装されています。「テキストが変更中のドキュメントで、パタンがユーザーが探したい特定のキーワード」というのが典型的な例でしょう。

この文字列照合問題のための効率の良いアルゴリズムは、テキストエディタのプログラムには不可欠です。また、文字列照合アルゴリズムは、ＤＮＡ配列中の特定のパタンを発見したり、インターネットにおけるキーワード検索など、様々な場面で応用されます。

文字列照合問題は以下のように形式化されます。
文字列 T[n] を長さ n のテキスト、文字列 P[m] を長さ m のパタンとし、T と P の各要素はある有限の集合 S に属するとします。例えば、S = {0, 1}、S = {0, 1, 2, ..., a, b, ... z} などが考えられます。

「 0 ≦ i ≦ n - m 」 かつ 「 0 ≦ j ≦ m - 1 について P[j] = T[i + j] 」を満たすとき、「テキスト T の i の位置にパタン P が存在する」と言います（下図参照）。

例: i = 3 P[0] = T[3 + 0] = T[3] P[1] = T[3 + 1] = T[4] P[2] = T[3 + 2] = T[5]

文字列照合問題は、テキストの中に存在する指定されたパタンの全ての位置を発見する問題です。

これから紹介していこうと思うアルゴリズムです↓

• Brute-force algorithm 力任せの文字列照合アルゴリズム
• Knuth-Morris-Pratt (KMP) algorithm
• Boyer-Moore algorithm

バブルソート Bubble Sort

選択ソートはある決められた範囲から、最小値を見つけ要素の交換を行いました。バブルソートは、隣り合う要素の比較のみを繰り返すことによって並び替えを行います。

バブルソートのアルゴリズムを以下に示します。

各計算ステップにおいて、配列は「ソート済みの部分」と「未ソートの部分」とに分けられます。 未ソートの部分がなくなるまで、以下の処理を繰り返します。 データの先頭から隣どうしを比較して、大小関係が逆ならばそれらを交換していきます。この交換処理を未ソートの部分の末尾まで順番に行います。

バブルソートはその名前からイメージできるように、泡（Bubble）が水面に上がっていくようにデータが動いていきます。上記のアルゴリズムでは、ソートされたデータが配列の後ろから順番に決定されていきます。

具体的な例で、バブルソートの流れを初期状態から見てみましょう。上図に示す処理が終了すると、データの中で一番大きい要素が配列の末尾に移動します。続きの処理を下図に示します。



上図のステップ1.～5.が終了すると、データの中で2番目に大きい要素が配列の後ろから2番目に移動します。以下同様に、ステップ6.～9.、ステップ10.～12.、ステップ13.～14.でそれぞれソート済みの部分に追加されるデータが決まっていきます。

以下のプログラムはは、データの数 size　を読み込み、size個の整数のデータを配列 data に読み込んで、配列の要素をバブルソートによって昇順に整列し、それを出力します。プログラムの本体である bubbleSort を見てみましょう。下図に示すように、07行目の最初のループ変数であるendは未ソートの部分の末尾を示します。08行目のforループによってデータの先頭からendまでの範囲で、隣どうしを比較して逆順ならば交換する処理を繰り返します。この２番目のループが終わると、その範囲のなかで一番大きい値がソートされていない部分の末尾に移動します。

#include<iostream>
using namespace std;

#define MAX 100

void bubbleSort( int size, int data[] ){
    for ( int end = size - 1; end > 0; end-- ){
        for ( int j = 0; j < end; j++ ){
            // 隣どうしを比較して、逆順ならば交換
            if ( data[j+1] < data[j] ){
                int tmp = data[j];
                data[j] = data[j+1];
                data[j+1] = tmp;
            }
        }
    }
}

バブルソートの計算量を考えてみましょう。データの数をnとすると、バブルソートは未ソートの部分における隣どうしの比較演算をn-1回、n-2回、n-3回、・・・、1回行います。よって合計の比較演算の回数は(n²-n)/2となり、バブルソートはO(n²)のアルゴリズムとなります。

練習問題

エラトステネスの篩（ふるい） The sieve of Eratosthenes

素数の項目で紹介した素数判定プログラムは、与えられた１つの数 x が素数であるか否かを判定するだけでしたが、多くの問題を解く場合、あらかじめ作成された素数表を用いるのが効率的です。

与えられた範囲内の全ての素数を列挙するための高速なアルゴリズムが、エラトステネスの篩（ふるい） The sieve of Eratosthenes です。エラトステネスの篩は、以下の手順で素数表を生成します。ここで、素数表を prime[x] が真ならば x は素数、偽ならば x は合成数（素数ではない整数）であるようなブール型の１次元配列とします。

	2 以上の整数を列挙しておく
	最小である 2 を残して、その倍数をすべて消す
	残った最小の 3 を残して、その倍数をすべて消す
	残った最小の 5 を残して、その倍数をすべて消す
	以下同様に、まだ消えていない最小の数を残してその倍数を消すことを繰り返せば、素数だけが残る

エラトステネスの篩を実装したプログラム例を以下に示します。

001 void eratos ( int n, bool prime[]){
002     for ( int i = 0; i <= n; i++ ) prime[i] = false;
003     // 奇数を素数の候補とする
004     for ( int i = 3; i <= n; i += 2 ) prime[i] = true;
005     prime[2] = true;
006     int limit = (int)sqrt((double)n) + 1;
007     // 合成数でない奇数だけを残す
008     for ( int i = 3; i <= limit; i += 2 ){
009         if ( !prime[i] ) continue;
010         // i は素数　i の倍数をすべて消す
011         for ( int j = i + i; j <= n; j += i ) prime[j] = false;
012     }
013 }

このプログラムは、以下のように多少改良することができます。

void eratos ( int n, bool prime[]){
    for ( int i = 0; i <= n; i++ ) prime[i] = false;
    // 奇数を素数の候補とする
    for ( int i = 3; i <= n; i += 2 ) prime[i] = true;
    prime[2] = true;
    int limit = (int)sqrt((double)n) + 1;
    // 合成数でない奇数だけを残す
    for ( int i = 3; i <= limit; i += 2 ){
        if ( !prime[i] ) continue;
        for ( int j = i * i, k = i * 2; j <= n; j += k )
        prime[j] = false;
    }
}

素数

約数が「1 とその数自身だけ」であるような自然数を素数と言います。素数を小さい順に列挙すると、

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, , ,

となります。1 は素数ではありません。素数は数論において重要な数ですが、コンピュータにおいても、暗号化のアルゴリズムに応用されるなど、重要な役割を果たします。

まずは、与えられた数が素数であるか否かを判定する素数判定のアルゴリズムを考えてみましょう。

誰もが思いつく素朴なアルゴリズムの１つを実装したプログラムを以下に示します。

001 bool isPrime( int x ){
002     if ( x <= 1 ) return false;
003     // 2 から x - 1 までの数で x を割ってみて
004     // 割り切れたら素数ではないので false を返す
005     for ( int i = 2; i < x - 1; i++ ){
006         if ( x % i == 0 ) return false;
007     }
008     return true;
009 }

この素朴な素数判定のアルゴリズムは、以下の２つのポイントに着目することで高速化が図れます。

	2 以外の偶数は素数ではないこと
	2 から x - 1 までではなく、2 から x の平方根までについて割り切れるかどうかをチェックすれば十分である

1 つめのポイントは明らかですが、2 つめのポイントは直感的に分かりにくいかもしれません。例を考えてみましょう。与えられた x が 31 とします。31 の平方根は 5.56... なので、31 が素数であるかどうかは、31 を 2 から 6 までの数で割ってみれば十分ということになります。ではなぜ 6 までで十分なのでしょうか。それは 7 から 30 までの数でもし 31 を割り切れる数が存在するならば、すでにチェックした 2 から 6 までの数に 31 を割り切れる数が必ず存在するからです。この改良を施したアルゴリズムを実装したプログラムを以下に示します。

001 bool isPrime( int x ){
002     if ( x < 2 ) return false;
003     else if ( x == 2 ) return true;
004     if ( x % 2 == 0 ) return false; // 偶数は素数ではない
005     // 3 から開始して 2 づつ増やす→奇数のみチェックする
006     int limit = (int)sqrt((double)x) + 1;
007     for ( int i = 3; i <= limit; i += 2 ){
008         if ( x % i == 0 ) return false;
009     }
010     
011     return true;
012 }

Dijkstra's Algorithm ダイクストラのアルゴリズム

グラフ G(V, E) における Single Source Shortest Path を求めるためのアルゴリズムの１つが、 Dijkstra によって発見された Dijkstra's Algorithm (ダイクストラのアルゴリズム) です。

Dijkstra's Algorithm は Greedy Algorithm で、以下のように実装します。

	グラフ G(V, E) のノード全体の集合を V、 s をスタート地点 (source)、 S を「 s から S 内の各ノード v への最短経路が必ず S 内にあるような、最短経路の部分的な解の集合」とします（下図を参照して下さい）。
	各計算ステップにおいて、V - S の各ノード v について、S 内のノードのみを経由した s から v への最短経路のコスト(距離)を d[v] とします。また、各ノード間のコスト（距離） (x, y) は、行列 M[x][y] に記録されているものとします。
1.	初期状態で、S は空です。 s に対して     d[s] = 0 s 以外の V に属する全てのノード v に対して     d[v] = ∞ と初期化します。
2.	ノードの集合 V - S の中から、d[u] が最小であるノード u を選択します。
	u を S に追加すると同時に、u に隣接しかつ V - S に属する全てのノード v に対する値を以下のように更新します。
	d[v] = min( d[v], d[u] + M[u][v] )
	この処理を S = V となるまで繰り返します。

以下の図は、u が選ばれた直後、v を更新する様子を示しています。

2. の計算ステップ直後（つまり次に u を選ぶ直前）において、 d[v] には、s から S 内のノードのみを経由した v までの最短コストが記録されています。

すなわち、処理が終了した時点で、V に属する全てのノード d[v] には、s から v までの最短コスト（距離）が記録されています。

コスト（距離）だけでなく、経路も求めたい場合は、Shortest spanning tree における v の親である pi[v] を用意し、pi[s] = -1と初期化し、更新部分を以下のように書き換えます。

もし、d[u] + M[u][v] ＜ d[v] ならば、     d[v] = d[u] + M[u][v]     pi[v] = u;

ダイクストラのアルゴリズムのプログラム例です。

int N;
int M[MAX][MAX];

void dijkstra( int s, int d[] ){
    bool visited[MAX]; // S に属するノードは true
    for ( int i = 0; i < N; i++ ){
        d[i] = INT_MAX;
        visited[i] = false;
    }
    
    d[s] = 0; // 最初に s が u として選ばれる
    
    while( 1 ){
        int u; // 最適なノード u を選ぶ
        int mincost = INT_MAX;
        for ( int i = 0; i < N; i++ ){
            if ( !visited[i] && d[i] < mincost ){
                mincost = d[i]; u = i;
            }
        }
        
        // u が存在しなかったとき、つまり S がこれ以上増えないとき、終了
        if ( mincost == INT_MAX ) break;
        
        visited[u] = true; // u を S に追加
        
        for ( int v = 0; v < N; v++ ){
            // v が S に含まれる場合　または　エッジ(u, v)がない場合は　無視
            if ( visited[v] || M[u][v] == INT_MAX ) continue;
            
            d[v] = min( d[v], d[u] + M[u][v] );
        }
    }
}

経路を生成するために pi を加えたプログラムです。

001 int N;
002 int M[MAX][MAX];
003 
004 void dijkstra( int source, int d[], int pi[] ){
005     bool visited[MAX]; // S に属するノードは true
006     for ( int i = 0; i < N; i++ ){
007         d[i] = INT_MAX;
008         visited[i] = false;
009     }
010     
011     d[source] = 0; // 最初に s が u として選ばれる
012     pi[source] = -1;
013     
014     while( 1 ){
015         int u; // 最適なノード u を選ぶ
016         int mincost = INT_MAX;
017         for ( int i = 0; i < N; i++ ){
018             if ( !visited[i] && d[i] < mincost ){
019                 mincost = d[i]; u = i;
020             }
021         }
022         
023         // u が存在しなかったとき、つまり S がこれ以上増えないとき、終了
024         if ( mincost == INT_MAX ) break;
025         
026         visited[u] = true; // u を S に追加
027         
028         for ( int v = 0; v < N; v++ ){
029             // v が S に含まれる場合　または　エッジ(u, v)がない場合は　無視
030             if ( visited[v] || M[u][v] == INT_MAX ) continue;
031             if ( d[u] + M[u][v] < d[v] ){
032                 d[v] = d[u] + M[u][v];
033                 pi[v] = u;
034             }
035         }
036     }
037 }

Disjoint Sets, Union Find

Disjoint Sets, Union Find とは？
Disjoint Sets は、多くのデータを互いに素な集合（１つの要素が複数の集合に属することがない集合）に分類して管理するためのデータ構造です。このデータ構造は、動的に以下の操作を高速に処理します。

	makeSet( x ): 要素が x ただ１つである新しい集合を作る
	findSet( x ): 要素 x が属する集合の"代表"の要素を求める つまり、要素 x がどの集合に属するかを調べることができる
	union( x, y ): 指定された２つの要素 x, y を unite（合併）する

Disjoint Sets に対して、「指定された２つの要素 x, y が、同じ集合に属するか？」どうかを調べる操作、つまり findSet(x) == findSet(y) は Union Find と呼ばれます。

Disjoint Sets の表現
Disjoint Sets を実装する方法はいくつかありますが、ここでは Disjoint Sets forests と呼ばれる"森"の構造を用います。森を構成する tree (木)が各集合を表し、tree の各ノードが集合内の各要素を表します。

findSet( x )
集合を区別するために、各 tree の根を代表 (representative) とし、集合を識別するために用います。従って、findSet( x ) は、要素 x が属する tree （集合）の根を値として返します。これを実現するために、各ノードはそれ自身から代表まで到達できるように、親へのポインタを持ちます（ただし、代表 はそれ自身へのポインタを持ちます）。

例えば、下図の Disjoint Sets で findSet( 5 ) の結果は 1 、findSet( 0 ) の結果も 1 なので、5 と 0 が同じ集合に属することが分かります。

findSet( x ) の計算量は、各ノードから代表までの間にたどるポインタの数 = 木の高さです。ここで、findSet( x ) は、単に代表を求めるだけではなく、次に実行される findSet( x ) の効率を高めるために、ある工夫をしています。それは、ある要素の代表を求めるとき、その要素から代表に至る経路の全てのノードについて、ポインタが直接代表を指すように変更します。例えば、下図はある集合 S に対して findSet( 0 ) を実行した結果を示しています。

S		S'
	→  findSet( 0 )  →

S において、各ノードは親へのポインタをもち、要素 0 の代表は、0 → 1 → 2 → 3 → 4 とたどることによって 4 と求められます。その経路上のノードのポインタを直接 4 を指すように処理し S' を得ます。

この工夫によって、極めて高さの低い木が生成されることがわかります。つまり、次回実行されるであろう S' に属する要素 x における findSet( x ) の操作は極めて少ない計算量で行うことができます。これを path compression (経路圧縮)といいます。

再帰を用いて path compression を行う findSet の実装は、以下のようになります。

    int findSet( int x ){
        if ( x != p[x] ) p[x] = findSet( p[x] );
        return p[x];
    }

union( x, y )
指定された２つの要素 x, y を合併する操作は、x の代表と y の代表を求め、どちらか一方を新しい代表として選びます。つまり、代表にならなかったほうのポインタを新しい代表を指すように更新します。例えば、下図はある Disjoint Sets S に対して union( 2, 4 ) を実行した結果を示しています。

S		S'
	→  union( 2, 4 )  →

ここで重要なことが、"どちらの代表を新しい代表として選ぶか"です。ポイントは、合併する集合を表す木の高さです。低い方の木を高い方の木に合併すれば、新しい木の高さが高くなることはありません。従って、高い木の代表を新しい代表にします。各ノード x を根としたときの木の高さに関する情報を rank[x] という変数に記録します。１つの要素が１つの集合を成している初期状態では、rank[x] は全て 0 としておきます。そして同じ rank の木を合併するときに rank を１つ増やします。以下に例を示します。

S		S'
	→  merge( 0, 5 )  →

S		S'
	→  merge( 0, 5 )  →

merge や経路圧縮を伴う findSet によって、rank の値は変化していきます。従って、rank は木の正確な高さを与えるものではなく、木の高さの上限を与えます。

以下に、DisjointSets のデータ構造を実装したプログラムを示します。

class DisjointSet{
    public:
    DisjointSet(){}
    DisjointSet( int size ){
        rank.resize( size, 0 );
        p.resize( size, 0 );
    }
    
    void makeSet( int x ){
        p[x] = x;
        rank[x] = 0;
    }
    
    void union( int x, int y ){
        link( findSet(x), findSet(y) );
    }
    
    int findSet( int x ){
        if ( x != p[x] ) p[x] = findSet( p[x] );
        return p[x];
    }
    
    bool isSameSet( int x, int y ){
        return ( findSet(x) == findSet(y) );
    }
    
    private:
    vector<int> rank, p;
    
    void link ( int x, int y ){
        if ( rank[x] > rank[y] ){
            p[y] = x;
        } else {
            p[x] = y;
            if ( rank[x] == rank[y] ) rank[y]++;
        }
    }
};

計算量の解析
ここで紹介した DisjointSets における Union Find は、経路圧縮と rank を用いることによって極めて高速になります。（解析は省略します）