今回は反復深化で 8 パズルを解いてみましょう。経路の探索 で説明したように、反復深化は最短手数を求めることができるアルゴリズムです。幅優先探索と違って局面を保存する必要が無いため、必要となるメモリは深さ優先探索と同程度で済みます。また、プログラムも深さ優先探索と同じくらい簡単に作成することができます。
ただし、同じ探索を何度も繰り返すため実行時間が増大する、という欠点があります。ようするに、使用するメモリは少ないが実行時間が長くなるアルゴリズムなのです。実行時間が長くなるといっても、枝刈りを工夫することでパズルを高速に解くことができます。メモリ不足になる場合には、積極的に使ってみたいアルゴリズムといえるでしょう。
幅優先探索では全ての局面を保存しましたが、反復深化ではその必要はありません。盤面は配列 board で表します。駒の移動は board を書き換えて、バックトラックする時は元に戻すことにします。動かした駒はリスト move_list に格納します。動かした駒がわかれば局面を再現できるので、それで移動手順を表すことにします。
それでは、解を求める関数 solver を作りましょう。次のリストを見てください。
リスト 1 : 単純な反復深化による解法 let solver board goal = let count = ref 0 in let rec ids n limit space move_list = if n = limit then if board = goal then begin count := !count + 1; print_answer (List.tl (List.rev move_list)) end else () else List.iter (fun x -> let p = board.(x) in if p <> List.hd move_list then begin (* 駒を動かす *) board.(space) <- p; board.(x) <- 0; (* 再帰呼び出し *) ids (n+1) limit x (p::move_list); (* 元に戻す *) board.(space) <- 0; board.(x) <- p end else ()) adjacent.(space) in let i = ref 1 in while (!i <= 31 && !count = 0) done ids 0 !i (position 0 board) [-1]; i := !i + 1 done
探索処理は局所関数 ids で行います。ids の引数 n が手数、limit が反復深化の上限値、space が空き場所の位置、move_list が移動手順を表します。n が limit に達したら、パズルが解けたかチェックします。solver の引数 goal が完成形を表す配列です。完成形に到達したら、変数 count の値を +1 してから関数 print_answer で手順を表示します。上限値に達していない場合は、駒を移動して新しい局面を作ります。
8 パズルのように、元の局面に戻すことが可能(可逆的)なパズルの場合、単純な深さ優先探索では同じ移動手順を何度も繰り返すことがあります。そうなると、とんでもない解を出力するだけではなく、再帰呼び出しが深くなるとスタックがオーバーフローしてプログラムの実行が停止してしまいます。
このような場合、局面の履歴を保存しておいて同じ局面がないかチェックすることで、解を求めることができるようになります。ただし、同一局面をチェックする分だけ時間が余分にかかりますし、最初に見つかる解が最短手数とは限りません。
反復深化では深さが制限されているため、同一局面のチェックを行わなくてもスタックオーバーフローが発生することはありません。そのかわり、無駄な探索はどうしても避けることができません。8 パズルの場合、1 手前に動かした駒を再度動かすと 2 手前の局面に戻ってしまいます。完全ではありませんが、このチェックを入れるだけでもかなりの無駄を省くことができます。
プログラムでは、リスト move_list に移動した駒を格納しているので、1 手前と同じ駒は動かさないようにチェックしています。なお、move_list の初期値はダミーデータを入れて [-1] としています。
最後に、関数 ids を呼び出します。変数 i が上限値を表します。i を 1 手ずつ増やして関数 ids を呼び出します。変数 count が 0 より大きい場合、解が見つかったので while ループを終了します。プログラムはこれで完成です。
それでは実行してみましょう。
let () = let a = Sys.time () in solver [|8;6;7;2;5;4;3;0;1|] [|1;2;3;4;5;6;7;8;0|]; print_float (Sys.time () -. a);;
当然ですが最短手数は 31 手で 40 通りの手順が表示されました。実行時間は 53.8 秒 (ocamlc version 4.05.0, Ubunts 18.04 (Windows Subsystem for Linux), Intel Core i5-6200U 2.30GHz) かかりました。1 分ちかくかかるのですから、やっぱり遅いですね。反復深化の場合、枝刈りを工夫しないと高速に解くことはできません。そこで、反復深化の常套手段である「下限値枝刈り法」を使うことにしましょう。
下限値枝刈り法は難しいアルゴリズムではありません。たとえば、5 手進めた局面を考えてみます。探索の上限値が 10 手とすると、あと 5 手だけ動かすことができますね。この時、パズルを解くのに 6 手以上かかることがわかれば、ここで探索を打ち切ることができます。
このように、必要となる最低限の手数が明確にわかる場合、この値を「下限値 (Lower Bound)」と呼びます。この下限値を求めることができれば、「今の移動手数+下限値」が探索手数を超えた時点で、枝刈りすることが可能になります。これが下限値枝刈り法の基本的な考え方です。
さて、下限値を求める方法ですが、これにはいろいろな方法が考えられます。今回は、各駒が正しい位置へ移動するまでの手数 (移動距離) [*1] を下限値として利用することにしましょう。次の図を見てください。
図 1 : 下限値の求め方
たとえば、右下にある 1 の駒を左上の正しい位置に移動するには、最低でも 4 手必要です。もちろん、ほかの駒との関連で、それ以上の手数が必要になる場合もあるでしょうが、4 手より少なくなることは絶対にありません。同じように、各駒について最低限必要な手数を求めることができます。そして、その合計値はパズルを解くのに最低限必要な手数となります。これを下限値として利用することができます。ちなみに、図 1 (2) の初期状態の下限値は 21 手になります。
下限値枝刈り法を使う場合、下限値の計算を間違えると正しい解を求めることができなくなります。たとえば、10 手で解ける問題の下限値を 11 手と計算すれば、最短手数を求めることができなくなります。それどころか、10 手の解しかない場合は、答えを求めることすらできなくなります。下限値の計算には十分に注意してください。
それでは、プログラムを作りましょう。下限値の求め方ですが、駒を動かすたびに各駒の移動距離を計算していたのでは時間がかかります。8 パズルの場合、1 回に一つの駒しか移動しないので、初期状態の下限値を求めておいて、動かした駒の差分だけ計算すればいいでしょう。
また、駒の移動距離はいちいち計算するのではなく、あらかじめ計算した結果を配列に格納しておきます。この配列を distance とすると、盤面から移動距離を求めるプログラムは次のようになります。
リスト 2 : 移動距離を求める let distance = [| [|0; 0; 0; 0; 0; 0; 0; 0; 0|]; [|0; 1; 2; 1; 2; 3; 2; 3; 4|]; [|1; 0; 1; 2; 1; 2; 3; 2; 3|]; [|2; 1; 0; 3; 2; 1; 4; 3; 2|]; [|1; 2; 3; 0; 1; 2; 1; 2; 3|]; [|2; 1; 2; 1; 0; 1; 2; 1; 2|]; [|3; 2; 1; 2; 1; 0; 3; 2; 1|]; [|2; 3; 4; 1; 2; 3; 0; 1; 2|]; [|3; 2; 3; 2; 1; 2; 1; 0; 1|] |] (* 移動距離を求める *) let get_distance board = Array.fold_left (fun a n -> a + distance.(board.(n)).(n)) 0 board
distance は 2 次元配列で「駒の種類×駒の位置」を表しています。空き場所は関係ないので、distance.(0) はダミー (要素が 0 の配列) となります。関数 get_distance は盤面 board にある駒と位置から移動距離を求めます。この処理は distance.(board.(n)).(n) の合計値を Array.fold_left で求めるだけです。
次は、下限値枝刈り法による反復深化を行う関数 solver を作ります。次のリストを見てください。
リスト 3 : 下限値枝刈り法 let solver board goal = let count = ref 0 in let rec ids n limit space move_list lower = if n = limit then if board = goal then begin count := !count + 1; print_answer (List.tl (List.rev move_list)) end else () else List.iter (fun x -> let p = board.(x) in if p <> List.hd move_list then let new_lower = lower - distance.(p).(x) + distance.(p).(space) in if new_lower + n <= limit then begin (* 駒を動かす *) board.(space) <- p; board.(x) <- 0; (* 再帰呼び出し *) ids (n+1) limit x (p::move_list) new_lower; (* 元に戻す *) board.(space) <- 0; board.(x) <- p end else () else ()) adjacent.(space) in let lower = get_distance board in let i = ref lower in while (!i <= 31 && !count = 0) do ids 0 !i (position 0 board) [-1] lower; i := !i + 1 done
局所関数 ids の引数 lower は現在の盤面 board の下限値を表しています。動かす駒の差分を計算して、新しい下限値 new_lower を求めます。そして、new_lower + n が上限値 limit を越えたら枝刈りを行います。limit 以下であれば ids を再帰呼び出しします。
最後に ids を呼び出す処理を修正します。関数 get_distance で初期状態の下限値 lower を求めます。下限値がわかるのですから、上限値 i は 1 手からではなく下限値 lower からスタートします。あとは ids に下限値 lower を渡して呼び出すだけです。
プログラムの主な修正はこれだけです。実際に実行してみると、実行時間は 0.063 秒でした。約 850 倍という高速化に驚いてしまいました。下限値枝刈り法の効果は極めて高いですね。
(* * eight3.ml : 8 Puzzle (反復深化) * * Copyright (C) 2008 Makoto Hiroi *) (* 隣接リスト *) let adjacent = [| [1; 3]; [0; 2; 4]; [1; 5]; [0; 4; 6]; [1; 3; 5; 7]; [2; 4; 8]; [3; 7]; [4; 6; 8]; [5; 7] |] (* 駒の位置を返す *) let position x ary = let rec iter n = if n = Array.length ary then raise Not_found else if x = ary.(n) then n else iter (n + 1) in iter 0 (* 手順の表示 *) let print_answer ls = List.iter (fun x -> print_int x; print_string " ") ls; print_newline () (* 反復深化 *) let solver board goal = let count = ref 0 in let rec ids n limit space move_list = if n = limit then if board = goal then print_answer (List.tl (List.rev move_list)) else () else List.iter (fun x -> let p = board.(x) in if p <> List.hd move_list then begin (* 駒を動かす *) board.(space) <- p; board.(x) <- 0; (* 再帰呼び出し *) ids (n+1) limit x (p::move_list); (* 元に戻す *) board.(space) <- 0; board.(x) <- p end else ()) adjacent.(space) in let i = ref 1 in while (!i <= 31 && !count = 0) do ids 0 !i (position 0 board) [-1]; i := !i + 1 done (* 実行 *) let () = let a = Sys.time () in solver [|8;6;7;2;5;4;3;0;1|] [|1;2;3;4;5;6;7;8;0|]; print_float (Sys.time () -. a)
(* * eight4.ml : 8 Puzzle (反復深化 + 下限値枝刈り法) * * Copyright (C) 2008 Makoto Hiroi *) (* 隣接リスト *) let adjacent = [| [1; 3]; [0; 2; 4]; [1; 5]; [0; 4; 6]; [1; 3; 5; 7]; [2; 4; 8]; [3; 7]; [4; 6; 8]; [5; 7] |] (* 移動距離 *) let distance = [| [|0; 0; 0; 0; 0; 0; 0; 0; 0|]; [|0; 1; 2; 1; 2; 3; 2; 3; 4|]; [|1; 0; 1; 2; 1; 2; 3; 2; 3|]; [|2; 1; 0; 3; 2; 1; 4; 3; 2|]; [|1; 2; 3; 0; 1; 2; 1; 2; 3|]; [|2; 1; 2; 1; 0; 1; 2; 1; 2|]; [|3; 2; 1; 2; 1; 0; 3; 2; 1|]; [|2; 3; 4; 1; 2; 3; 0; 1; 2|]; [|3; 2; 3; 2; 1; 2; 1; 0; 1|] |] (* 駒の位置を返す *) let position x ary = let rec iter n = if n = Array.length ary then raise Not_found else if x = ary.(n) then n else iter (n + 1) in iter 0 (* 手順の表示 *) let print_answer ls = List.iter (fun x -> print_int x; print_string " ") ls; print_newline () (* 移動距離を求める *) let get_distance board = Array.fold_left (fun a n -> a + distance.(board.(n)).(n)) 0 board (* 反復深化 *) let solver board goal = let count = ref 0 in let rec ids n limit space move_list lower = if n = limit then if board = goal then begin count := !count + 1; print_answer (List.tl (List.rev move_list)) end else () else List.iter (fun x -> let p = board.(x) in if p <> List.hd move_list then let new_lower = lower - distance.(p).(x) + distance.(p).(space) in if new_lower + n <= limit then begin (* 駒を動かす *) board.(space) <- p; board.(x) <- 0; (* 再帰呼び出し *) ids (n+1) limit x (p::move_list) new_lower; (* 元に戻す *) board.(space) <- 0; board.(x) <- p end else () else ()) adjacent.(space) in let lower = get_distance board in let i = ref lower in while (!i <= 31 && !count = 0) do ids 0 !i (position 0 board) [-1] lower; i := !i + 1 done let () = let a = Sys.time () in solver [|8;6;7;2;5;4;3;0;1|] [|1;2;3;4;5;6;7;8;0|]; print_float (Sys.time () -. a)
次のパズルを反復深化を使って解いてください。
ペグ・ソリテアは、盤上に配置されたペグ (駒) を、最後にはひとつ残るように取り除いていく、古典的なパズルです。ルールの説明は拙作のページ Puzzle DE Programming ペグ・ソリティア をお読みください。Hoppers は芦ヶ原伸之氏が考案されたペグ・ソリテアです。次の図を見てください。
図 : Hoppers
Hoppers は穴を 13 個に減らしていて、遊ぶのに手頃な大きさになっています。上図に示したように、最初に中央のペグを取り除きます。この状態から始めて、最後のペグが中央の位置に残る跳び方の最小手数を求めてください。
今回はペグの有無を bool (true, false) で表します。盤面のデータ型は bool array になります。盤面と配列の対応は下図を見てください。
図 : Hoppers の盤面
ペグの移動は跳び先表 (配列 jump_table) を用意すると簡単です。配列の要素はリストであることに注意してください。リストの要素は、跳び越されるペグの位置と跳び先の位置を格納したタプルです。たとえば、0 番の位置にあるペグは、1 番を跳び越して 2 番へ移動する場合と、3 番を跳び越して 6 番へ移動する場合と、5 番を飛び越して 10 番へ移動する場合の 3 通りがあります。これをリスト [(1, 2); (3, 6); (5, 10)] で表しています。
あとは単純な反復深化で最短手順を求めるだけです。詳細はプログラムリストをお読みくださいませ。
リスト : Hoppers の解法 (* 跳び先表 *) let jump_table = [| [(1, 2); (3, 6); (5, 10)]; [(3, 5); (6, 11); (4, 7)]; [(1, 0); (4, 6); (7, 12)]; [(6, 9)]; [(6, 8)]; [(3, 1); (6, 7); (8, 11)]; [(3, 0); (4, 2); (8, 10); (9, 12)]; [(4, 1); (6, 5); (9, 11)]; [(6, 4)]; [(6, 3)]; [(5, 0); (8, 6); (11, 12)]; [(8, 5); (6, 1); (9, 7)]; [(11, 10); (9, 6); (7, 2)] |] let size = 13 (* 盤面の大きさ *) let max_jump = 11 (* ペグがジャンプする回数 *) let hole = 6 (* 中央の穴の位置 *) (* 手順の表示 *) (* 連続跳びは [i,j,k,...] で表す *) let print_move1 move = let rec _print_move prev = function [] -> print_endline "]" | (from, to_)::xs -> if prev = from then Printf.printf ",%d" to_ else Printf.printf "][%d,%d" from to_; _print_move to_ xs in match move with [] -> () | (from, to_)::xs -> Printf.printf "[%d,%d" from to_; _print_move to_ xs (* 反復深化 *) let solver1 () = let count = ref 0 in let board = Array.make size true in let rec dfs jc limit move = if jc <= limit then if List.length move = max_jump then if board.(hole) then (count := !count + 1; print_move1 (List.rev move)) else () else for from = 0 to size - 1 do if not board.(from) then () else List.iter (fun (del, to_) -> if board.(del) && not board.(to_) then (board.(from) <- false; board.(del) <- false; board.(to_) <- true; dfs (if snd (List.hd move) = from then jc else jc + 1) limit ((from, to_)::move); board.(to_) <- false; board.(del) <- true; board.(from) <- true) else ()) jump_table.(from) done in (* 初手を 0 -> (3) -> 6: hole に限定 *) board.(0) <- false; board.(3) <- false; board.(hole) <- true; let i = ref 2 in while !i <= max_jump && !count = 0 do Printf.printf "----- %d -----\n" !i; dfs 1 !i [(0, hole)]; i := !i + 1 done
# solver1 ();; ----- 2 ----- ----- 3 ----- ----- 4 ----- ----- 5 ----- ----- 6 ----- ----- 7 ----- [0,6][9,3][2,0,6][11,1][10,0,2,6][8,4][12,2,6] [0,6][9,3][2,0,6][11,1][10,6][4,8][12,2,0,10,6] [0,6][9,3][2,0,6][11,1][12,2,6][8,4][10,0,2,6] [0,6][9,3][2,6][8,4][10,0,2,6][7,5][12,10,0,6] [0,6][9,3][2,6][8,4][10,0,2,6][11,1][12,2,0,6] [0,6][9,3][2,6][8,4][10,0,6][7,5][12,10,0,2,6] [0,6][9,3][2,6][8,4][12,2,0,6][5,7][10,12,2,6] [0,6][9,3][2,6][8,4][12,2,0,6][11,1][10,0,2,6] [0,6][9,3][2,6][8,4][12,2,6][5,7][10,12,2,0,6] [0,6][9,3][10,0,6][7,5][2,0,10,6][4,8][12,10,6] [0,6][9,3][10,0,6][7,5][2,6][8,4][12,10,0,2,6] [0,6][9,3][10,0,6][7,5][12,10,6][4,8][2,0,10,6] [0,6][9,3][10,6][4,8][2,0,6][11,1][12,2,0,10,6] [0,6][9,3][10,6][4,8][2,0,10,6][7,5][12,10,0,6] [0,6][9,3][10,6][4,8][2,0,10,6][11,1][12,2,0,6] [0,6][9,3][10,6][4,8][12,10,0,6][1,11][2,12,10,6] [0,6][9,3][10,6][4,8][12,10,0,6][7,5][2,0,10,6] [0,6][9,3][10,6][4,8][12,10,6][1,11][2,12,10,0,6] - : unit = ()
7 手で解くことができました。解は全部で 18 通りになりました。最近のパソコンは高性能なので、穴の数が少ない盤面であれば、単純な反復深化でも高速に解くことができるようです。
組み合わせの生成は拙作のページ 順列と組み合わせ で説明しました。このほかに、n 個の中から m 個を選ぶ組み合わせは、ビットのオンオフで表すことができます。たとえば、5 個の数字 (0 - 4) から 3 個を選ぶ場合、数字を 0 bit から 4 bit に対応させます。すると、1, 3, 4 という組み合わせは 11010 と表すことができます。今回はビットを使って組み合わせを求めてみましょう。最初に OCaml のビット演算について説明します。
OCaml にはビット演算を行う演算子が用意されています。主な演算子を表に示します。
演算子 | 型 | 機能 |
---|---|---|
land | int -> int -> int | ビットごとの論理積を返す |
lor | int -> int -> int | ビットごとの論理和を返す |
lxor | int -> int -> int | ビットごとの排他的論理和を返す |
lnot | int -> int | ビットごとの論理的な否定を返す |
lsl | int -> int -> int | m lsl n は m を n ビットだけ左シフトする |
lsr | int -> int -> int | m lsr n は m を n ビットだけ右シフトする |
asr | int -> int -> int | m lsr n は m を n ビットだけ算術右シフトする |
land はビットごとの論理積を返します。
# 5 land 3;; - : int = 1
0101 and 0011 --------- 0001
lor はビットごとの論理和を返します。
# 5 lor 3;; - : int = 7
0101 or 0011 -------- 0111
lxor はビットごとの排他的論理和を返します。
# 5 lxor 3;; - : int = 6
0101 xor 0011 --------- 0110
lnot はビットごとの論理的な否定を返します。
# lnot 0;; - : int = -1 # lnot 1;; - : int = -2
m lsl n は m を n ビット左シフトします。m lsr n は m を n ビット右シフトします。
# 1 lsl 8;; - : int = 256 # 256 lsr 4;; - : int = 16
このほかに、算術右シフトを行う演算子 asr もあります。
組み合わせを求めるプログラムは次のようになります。
リスト 4 : 組み合わせの生成 let rec combination n m a = if m = 0 then Printf.printf "%x\n" a else if m = n then Printf.printf "%x\n" (a lor ((1 lsl m) - 1)) else begin combination (n - 1) m a; combination (n - 1) (m - 1) (a lor (1 lsl (n - 1))) end
関数 combination は n 個の中から m 個を選ぶ組み合わせを生成して出力します。組み合わせは引数 a にセットします。m が 0 になったら、組み合わせがひとつできたので a を出力します。n が m と等しくなったならば、残り m 個を全て選びます。(1 lsl m) - 1 で m 個のビットをオンにして出力します。
あとは combination を再帰呼び出しします。最初の呼び出しは n 番目の数字を選ばない場合です。n - 1 個の中から m 個を選びます。次の呼び出しが n 番目の数字を選ぶ場合で、a の n - 1 ビットをオンにします。そして、n - 1 個の中から m - 1 個を選びます。
それでは 5 個の中から 3 個を選ぶ combination 5 3 0 の実行例を示します。
7 (00111) b (01011) d (01101) e (01110) 13 (10011) 15 (10101) 16 (10110) 19 (11001) 1a (11010) 1c (11100)
この場合、最小値は 0x07 (00111) で最大値は 0x1c (11100) になります。このように、combination は組み合わせを表す数を昇順で出力します。
次は、N 通りある組み合わせに 0 から N - 1 までの番号を付ける方法を紹介しましょう。たとえば、6 個の中から 3 個を選ぶ組み合わせは 20 通りありますが、この組み合わせに 0 から 19 までの番号を付けることができます。1 1 1 0 0 0 を例題に考えてみましょう。次の図を見てください。
図 : 6C3 の組み合わせ
最初に 5 をチェックします。5 を選ばない場合は 5C3 = 10 通りありますね。この組み合わせに 0 から 9 までの番号を割り当てることにすると、5 を選ぶ組み合わせの番号は 10 から 19 までとなります。
次に、4 をチェックします。4 を選ばない場合は、4C2 = 6 通りあります。したがって、5 を選んで 4 を選ばない組み合わせに 10 から 15 までの番号を割り当てることにすると、5 と 4 を選ぶ組み合わせには 16 から 19 までの番号となります。
最後に、3 をチェックします。同様に 3 を選ばない場合は 3 通りあるので、これに 16 から 18 までの番号を割り当て、5, 4, 3 を選ぶ組み合わせには 19 を割り当てます。これで組み合わせ 1 1 1 0 0 0 の番号を求めることができました。
では、0 0 0 1 1 1 はどうなるのでしょうか。左から順番にチェックしていくと、最初の 1 が見つかった時点で、その数字を選ばない組み合わせは存在しません。つまり、残りの数字をすべて選ぶしかないわけです。したがって、これが 0 番目となります。
このように、数字を選ぶときに、数字を選ばない場合の組み合わせの数を足し算していけば、その組み合わせの番号を求めることができるのです。
組み合わせを番号に変換するプログラムは次のようになります。
リスト 5 : 組み合わせを番号に変換 (* 組み合わせの数を求める *) let rec comb n r = if n = r || r = 0 then 1 else (comb n (r - 1)) * (n - r + 1) / r (* 組み合わせを番号に変換 *) let rec comb_to_num c n r v = if r = 0 || n = r then v else if c land (1 lsl (n - 1)) <> 0 then comb_to_num c (n - 1) (r - 1) (v + comb (n - 1) r) else comb_to_num c (n - 1) r v
関数 comb_to_num の引数 c はビットのオンオフで表した組み合わせ、引数 n と r は nCr の n と r を表しています。引数 v は求める番号を表します。n と r の値が同じになるか、もしくは r が 0 になれば、組み合わせの番号を計算できたので value を返します。
そうでない場合、c の n - 1 ビットの値を調べます。ビットがオンであれば、v に comb (n - 1) r の値を足し算し、r を -1 して comb_to_num を再帰呼び出しします。そうでなければ、v と r の値はそのままで comb_to_num を再帰呼び出しします。
逆に、番号から組み合わせを求めるプログラムも簡単に作ることができます。次のリストを見てください。
リスト 6 : 番号を組み合わせに変換 let rec num_to_comb v n r c = if r = 0 then c else if n = r then c lor ((1 lsl n) - 1) else let k = comb (n - 1) r in if v >= k then num_to_comb (v - k) (n - 1) (r - 1) (c lor (1 lsl (n - 1))) else num_to_comb v (n - 1) r c
引数 v が番号で、引数 n と r は nCr の n と r を表しています。引数 c が求める組み合わせです。たとえば、n = 5, r = 3 の場合、ビットが 1 になるのは 4C2 = 6 通りあり、0 になるのは 4C3 = 4 通りあります。したがって、数値が 0 - 3 の場合はビットを 0 にし、4 - 9 の場合はビットを 1 にすればいいわけです。
ビットを 0 にした場合、残りは 4C3 = 4 通りになるので、同様に次のビットを決定します。ビット 1 にした場合、残りは 4C2 = 6 通りになるので、v から 4 を引いて num_to_comb を再帰呼び出しして次のビットを決定します。
r が 0 の場合は、組み合わせが完成したので c を返します。n と r が等しい場合は、残りのビットをすべて 1 にセットしてから c を返します。それ以外の場合は、n-1Cr の値を comb (n - 1) r で求めて変数 k にセットします。v が k 以上であれば変数 c のビットを 1 にセットし、v から k を引き算して comb_to_num を再帰呼び出しします。そうでなければ、num_to_comb を再帰呼び出しするだけです。
それでは、n = 5, r = 3 の場合の実行例を示します。
# for i = 0 to 9 do let v = num_to_comb i 5 3 0 in let n = comb_to_num v 5 3 0 in Printf.printf "%d => %x => %d\n" i v n done;; 0 => 7 => 0 1 => b => 1 2 => d => 2 3 => e => 3 4 => 13 => 4 5 => 15 => 5 6 => 16 => 6 7 => 19 => 7 8 => 1a => 8 9 => 1c => 9
正常に動作していますね。この方法を使うと、n 個ある組み合わせの中の i 番目 (0 <= i < n) の組み合わせを簡単に求めることができます。
最も右側 (LSB 側) にある 1 を 0 にクリアする、逆に最も右側にある 0 を 1 にセットすることは簡単にできます。
(1) 右側にある 1 をクリア => x AND (- x) x : 1 1 1 1 x - 1 : 1 1 1 0 ---------------- AND : 1 1 1 0 x : 1 0 0 0 x - 1 : 0 1 1 1 ---------------- AND : 0 0 0 0 (2) 右側にある 0 を 1 にセット => x OR (x + 1) x : 0 0 0 0 x + 1 : 0 0 0 1 ---------------- OR : 0 0 0 1 x : 0 1 1 1 x + 1 : 1 0 0 0 ---------------- OR : 1 1 1 1
上図 (1) を見てください。x から 1 を引くと、右側から連続している 0 は桁借りにより 1 になり、最初に出現する 1 が 0 になります。したがって、x AND (x - 1) を計算すると、最も右側にある 1 を 0 にクリアすることができます。(2) の場合、x に 1 を足すと、右側から連続している 1 は桁上がりにより 0 になり、最初に出現する 0 が 1 になります。x OR (x + 1) を計算すれば、最も右側にある 0 を 1 にセットすることができます。
また、最も右側にある 1 を取り出すことも簡単にできます。簡単な例として 4 ビットの整数値を考えてみます。負の整数を 2 の補数で表した場合、4 ビットで表される整数は -8 から 7 になります。次の図を見てください。
0 : 0000 1 : 0001 -1 : 1111 1 AND (-1) => 0001 2 : 0010 -2 : 1110 2 AND (-2) => 0010 3 : 0011 -3 : 1101 3 AND (-3) => 0001 4 : 0100 -4 : 1100 4 AND (-4) => 0100 5 : 0101 -5 : 1011 5 AND (-5) => 0001 6 : 0110 -6 : 1010 6 AND (-6) => 0010 7 : 0111 -7 : 1001 7 AND (-7) => 0001 -8 : 1000 図 : 最も右側にある 1 を取り出す方法
2 の補数はビットを反転した値 (1 の補数) に 1 を加算することで求めることができます。したがって、x と -x の論理積 x AND (-x) は、最も右側にある 1 だけが残り、あとのビットはすべて 0 になります。
次は、ビットが 1 の個数を数える処理を作ってみましょう。次のリストを見てください。
リスト : ビットカウント let bit_count n = let rec _bit_count n c = if n = 0 then c else _bit_count (n land (n - 1)) (c + 1) in _bit_count n 0
整数 n の右側から順番に 1 をクリアしていき、0 になるまでの回数を求めます。とても簡単ですね。
整数を 32 bit に限定する場合、次の方法で 1 の個数をもっと高速に求めることができます。
リスト : ビットカウント (2) let bit_count32 n = let a = (n land 0x55555555) + ((n lsr 1) land 0x55555555) in let b = (a land 0x33333333) + ((a lsr 2) land 0x33333333) in let c = (b land 0x0f0f0f0f) + ((b lsr 4) land 0x0f0f0f0f) in let d = (c land 0x00ff00ff) + ((c lsr 8) land 0x00ff00ff) in (d land 0xffff) + (d lsr 16)
最初に、整数を 2 bit ずつに分割して 1 の個数を求めます。たとえば、整数 n を 4 bit で考えてみましょう。5 を 2 進数で表すと 0101 になります。n と論理積を計算すると 0, 2 番目のビットが 1 であれば、結果の 0, 2 番目のビットは 1 になります。同様に n を 1 ビット右シフトして論理積を計算すると、1, 3 番目のビットが 1 であれば、結果の 0, 2 番目のビットは 1 になります。あとは、それを足し算すれば 2 bit の中にある 1 の個数を求めることができます。
変数 a には 2 ビットの中の 1 の個数が格納されています。左隣の 2 ビットの値を足し算すれば、4 ビットの中の 1 の個数を求めることができます。次に、左隣の 4 ビットの値を足し算して 8 ビットの中の 1 の個数を求め、左隣の 8 ビットの値を足し算して、というように順番に値を加算していくと 32 ビットの中にある 1 の個数を求めることができます。
bit_count は 1 の個数が多くなると遅くなりますが、bit_count32 は 1 の個数に関係なく高速に動作します。興味のある方は試してみてください。