今回は n 個の中から r 個を選ぶ組み合わせをビットのオンオフで表してみましょう。たとえば、5 個の数字 (0 - 4) から 3 個を選ぶ場合、数字を 0 bit から 4 bit に対応させます。すると、1, 3, 4 という組み合わせは 11010 と表すことができます。これを ISLisp でプログラムすると次のようになります。
リスト : 組み合わせの生成
(import "bit")
;;; 2 進数表示
(defun display2 (n)
(format (standard-output) "~B~%" n))
(defun comb (fn n r a)
(cond
((= r 0) (funcall fn a))
((= n 0) nil)
(t
(comb fn (1- n) r a)
(comb fn (1- n) (1- r) (logior (ash 1 (1- n)) a)))))
(defun combination (fn n r) (comb fn n r 0))
関数 combination は n 個の中から r 個を選ぶ組み合わせを生成して出力します。実際の処理は関数 comb で行います。組み合わせは引数 a にセットします。r が 0 になったら、組み合わせがひとつできたので関数 fn を呼び出します。n が 0 になったら選ぶ数字がなくなったので nil を返します。
あとは comb を再帰呼び出しします。最初の呼び出しは n 番目の要素を選ばない場合です。n - 1 個の中から r 個を選びます。次の呼び出しが n 番目の要素を選ぶ場合です。(ash 1 (- n 1)) で 1 を n - 1 ビット左シフトして、a との論理和を計算します。これで、n - 1 番目のビットをオンにすることができます。そして、n - 1 個の中から r - 1 個を選びます。
それでは 5 個の中から 3 個を選ぶ combination の実行例を示します。
$ eisl Easy-ISLisp Ver2.65 > (load "bitcomb.lsp") T > (combination #'display2 5 3) 111 1011 1101 1110 10011 10101 10110 11001 11010 11100 NIL
この場合、最小値は 111 で最大値は 11100 になります。このように、combination は組み合わせを表す数を昇順で出力します。ところで、参考文献『Ｃ言語による最新アルゴリズム事典』の「組み合わせの生成」には、再帰呼び出しを使わずに同じ結果を得る方法が解説されてます。とても巧妙な方法なので、興味のある方は読んでみてください。
次は、N 通りある組み合わせに 0 から N - 1 までの番号を付ける方法を紹介しましょう。たとえば、6 個の中から 3 個を選ぶ組み合わせは 20 通りありますが、この組み合わせに 0 から 19 までの番号を付けることができます。1 1 1 0 0 0 を例題に考えてみましょう。次の図を見てください。
5 4 3 2 1 0
─────────
0 0 0 1 1 1 ↑
0 0 1 0 1 1 │
0 0 1 1 0 1 │
0 0 1 1 1 0 │
0 1 0 0 1 1 │
0 1 0 1 0 1 5Ｃ3 = 10 通り
0 1 0 1 1 0 │
0 1 1 0 0 1 │
0 1 1 0 1 0 │
0 1 1 1 0 0 ↓
─────────
1 0 0 0 1 1 ↑
1 0 0 1 0 1 │
1 0 0 1 1 0 │
1 0 1 0 0 1 4Ｃ2 = 6 通り
1 0 1 0 1 0 │
1 0 1 1 0 0 ↓
────────
1 1 0 0 0 1 ↑
1 1 0 0 1 0 3Ｃ1 = 3 通り
1 1 0 1 0 0 ↓
───────
1 1 1 0 0 0 19 番目
─────────
図 : 6Ｃ3 の組み合わせ
最初に 5 をチェックします。5 を選ばない場合は \({}_5 \mathrm{C}_3 = 10\) 通りありますね。この組み合わせに 0 から 9 までの番号を割り当てることにすると、5 を選ぶ組み合わせの番号は 10 から 19 までとなります。
次に、4 をチェックします。4 を選ばない場合は、\({}_4 \mathrm{C}_2 = 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 番目となります。
このように、数字を選ぶときに、数字を選ばない場合の組み合わせの数を足し算していけば、その組み合わせの番号を求めることができるのです。プログラムは次のようになります。
リスト : 組み合わせに番号を付ける
;;; 組み合わせの数
(defun combination-number (n r)
(if (or (= n r) (= r 0))
1
(div (* (combination-number n (- r 1)) (+ (- n r) 1)) r)))
;;; 組み合わせに番号を付ける
(defun comb-to-num (n r c v)
(cond
((or (= r 0) (= n r)) v)
((logbitp (- n 1) c)
(comb-to-num (- n 1) (- r 1) c (+ (combination-number (- n 1) r) v)))
(t
(comb-to-num (- n 1) r c v))))
(defun to-number (n r c)
(comb-to-num n r c 0))
関数 to-number は n 個の中から r 個を選ぶ組み合わせ c を番号に変換します。実際の処理は関数 comb-to-num で行います。この関数は c の上位ビットから順番にチェックしていきます。
n 個の中から r 個を選ぶ場合、n - 1 番目のビットが 1 であれば \({}_{n-1} \mathrm{C}_r\) の値を v に加算して comb-to-num を再帰呼び出しします。このとき、n 個の中から一つ選んだので、r の値を -1 することをお忘れなく。ビットが 0 であれば、v の値はそのままで comb-to-num を再帰呼び出しします。n = r または r = 0 になったら v を返します。これが再帰呼び出しの停止条件になります。
次は、番号から組み合わせを求める関数 from-number を作ります。次のリストを見てください。
リスト : 番号から組み合わせを求める
(defun num-to-comb (n r v c)
(if (= n r)
(logior (- (expt 2 n) 1) c)
(let ((k (combination-number (- n 1) r)))
(if (>= v k)
(num-to-comb (- n 1) (- r 1) (- v k) (logior (ash 1 (- n 1)) c))
(num-to-comb (- n 1) r v c)))))
(defun from-number (n r v)
(num-to-comb n r v 0))
関数 from-number は、組み合わせ \({}_n \mathrm{C}_r\) の番号 v を組み合わせ c に変換します。実際の処理は関数 num-to-comb で行います。この関数は組み合わせ c の上位ビットから決定していきます。
たとえば、n = 6, r = 3 の場合、ビットが 1 になるのは \({}_5 \mathrm{C}_2 = 10\) 通りあり、0 になるのは \({}_5 \mathrm{C}_3 = 10\) 通りあります。したがって、数値が 0 - 9 の場合はビットを 0 にし、10 - 19 の場合はビットを 1 にすればいいわけです。
ビットを 0 にした場合、残りは \({}_5 \mathrm{C}_3 = 10\) 通りになるので、同様に次のビットを決定します。ビットを 1 にした場合、残りは \({}_5 \mathrm{C}_3 = 10\) 通りになります。数値から \({}_5 \mathrm{C}_3 = 10\) を引いて次のビットを決定します。
プログラムでは、\({}_{n-1} \mathrm{C}_r\) の値を関数 combination-number で求めて変数 k にセットします。v が k 以上であれば c の n - 1 番目のビットを 1 にセットし、v から k を引き算します。そして、次のビットを決めればいいわけです。r = 0 になったら c を返します。また、r が n と等しくなったら、c の残りのビットを全て 1 にセットした値を返します。
それでは実際に試してみましょう。
リスト : 簡単なテスト
(defun test (n r)
(combination (lambda (c)
(let ((m (to-number n r c)))
(format (standard-output) "~B => ~D => ~B~%" c m (from-number n r m))))
n r))
> (test 5 3) 111 => 0 => 111 1011 => 1 => 1011 1101 => 2 => 1101 1110 => 3 => 1110 10011 => 4 => 10011 10101 => 5 => 10101 10110 => 6 => 10110 11001 => 7 => 11001 11010 => 8 => 11010 11100 => 9 => 11100 NIL > (test 6 3) 111 => 0 => 111 1011 => 1 => 1011 1101 => 2 => 1101 1110 => 3 => 1110 10011 => 4 => 10011 10101 => 5 => 10101 10110 => 6 => 10110 11001 => 7 => 11001 11010 => 8 => 11010 11100 => 9 => 11100 100011 => 10 => 100011 100101 => 11 => 100101 100110 => 12 => 100110 101001 => 13 => 101001 101010 => 14 => 101010 101100 => 15 => 101100 110001 => 16 => 110001 110010 => 17 => 110010 110100 => 18 => 110100 111000 => 19 => 111000 NIL
正常に動作していますね。この方法を使うと、n 個ある組み合わせの中の i 番目 (0 <= i < n) の組み合わせを簡単に求めることができます。
それでは簡単な例題として、Puzzel DE Programming で取り上げた「ライツアウト」というパズルを ISLisp で解いてみましょう。ライツアウトは光っているボタンをすべて消すことが目的のパズルです。ルールはとても簡単です。あるボタンを押すと、そのボタンと上下左右のボタンの状態が反転します。つまり、光っているボタンは消灯し消えていたボタンは点灯します。次の図を見てください。
□□□□□ □□□□□ ０ １ ２ ３ ４
□□□□□ □□■□□ ５ ６ ７ ８ ９
□□□□□ ─→ □■■■□ 10 11 12 13 14
□□□□□ □□■□□ 15 16 17 18 19
□□□□□ □□□□□ 20 21 22 23 24
（Ａ）中央のボタンを押した場合 （Ｂ）座標
図 : ライツアウトの点灯パターン
ボタンは 5 行 5 列に配置されています。上図に示したように、中央のボタン 12 を押すとそのボタンと上下左右のボタンの状態が反転します。
ライツアウトはライトオン・オフの 2 種類の状態しかないので、盤面はリストよりもビットを使って表した方が簡単です。ライトオン・オフの状態を 1 と 0 で表し、各ビットとボタンの座標を対応させると、盤面は 0 から 33554431 の整数値で表すことができます。
ボタンを押してライトの状態を反転する処理も簡単です。たとえば、中央のボタン 12 を押した場合、7, 11, 12, 13, 17 のライトを反転させます。この場合、5 つのボタンのビットをオンにした値 #x23880 と、盤面を表す整数値の排他的論理和 (xor) を求めれば、5 つのライトの状態を反転することができます。次の例を見てください。
0 xor #x23880 => #x23880 % 消灯の状態でボタン 12 を押す（点灯する） #x23880 xor #x23880 => 0 % もう一度同じボタンを押す（消灯する）
このように、ライツアウトは同じボタンを二度押すと元の状態に戻ります。したがって、同じボタンは二度押さなくてよいことがわかります。また、実際にボタンを押してみるとわかりますが、ボタンを押す順番は関係がないことがわかります。たとえば、ボタン 0 と 1 を押す場合、0 -> 1 と押すのも 1 -> 0 と押すのも同じ結果になります。この 2 つの法則から、ボタンを押す組み合わせは全部で 2 ^ 25 通りになります。
ライツアウトを解くいちばん単純な方法は、ボタンを押す組み合わせを生成して、実際にライトが全部消えるかチェックすることです。プログラムは次のようになります。
リスト : ライツアウトの解法
(import "macro")
;;; ボタンを押したときのパターン
(defglobal pattern
#(#x0000023 #x0000047 #x000008e #x000011c #x0000218
#x0000461 #x00008e2 #x00011c4 #x0002388 #x0004310
#x0008c20 #x0011c40 #x0023880 #x0047100 #x0086200
#x0118400 #x0238800 #x0471000 #x08e2000 #x10c4000
#x0308000 #x0710000 #x0e20000 #x1c40000 #x1880000))
;;; ボタンを押す
(defun push-buttons (xs board n)
(cond
((= xs 0) board)
((logtest xs 1)
(push-buttons (ash xs -1) (logxor (aref pattern n) board) (+ n 1)))
(t
(push-buttons (ash xs -1) board (+ n 1)))))
;;; 解の表示
(defun print-answer (xs)
(dotimes (n 25)
(format (standard-output)
"~A "
(if (logbitp n xs) "O" "X"))
(if (= (mod (+ n 1) 5) 0)
(format (standard-output) "~%"))))
;;; 解法
(defun solver (board)
(block
exit
(for ((n 1 (+ n 1)))
((> n 25))
(format (standard-output) "----- ~D -----~%" n)
(combination
(lambda (xs)
(cond
((= (push-buttons xs board 0) 0)
(print-answer xs)
(return-from exit))))
25 n))))
最初に、ボタンを押したときにライトの状態を反転させるための値を配列 pattern に定義します。そして、関数 solver で盤面 board の解を求めます。combination で 25 個のボタン (0 - 24) から n 個を選ぶ組み合わせを生成し、関数 push-buttons で選んだボタンを押します。その結果が 0 であれば、関数 print-answer で解を出力して、return-from で処理を終了します。
関数 push-buttons は引数 xs の n 番目のビットがオンであれば、n 番目のボタンを押して新しい盤面を生成します。これを再帰定義で行っています。関数 print-answer は押すボタンを O で、押さないボタンを X で表示します。5 行 5 列に出力するため、(mod x 5) の値が 0 であれば format で改行を出力します。
それでは実行してみましょう。ライトが全部点灯している状態 (#x1ffffff) を解いてみます。
Easy-ISLisp Ver2.65 > (compile-file "bitcomb.lsp") ... ワーニング (省略) ... initialize pass1 pass2 compiling DISPLAY2 compiling COMBINATION-NUMBER compiling COMB compiling COMBINATION compiling COMB-TO-NUM compiling TO-NUMBER compiling NUM-TO-COMB compiling FROM-NUMBER compiling TEST compiling PUSH-BUTTONS compiling PRINT-ANSWER compiling SOLVER finalize invoke CC T > (load "bitcomb.o") T > (time (solver #x1ffffff)) ----- 1 ----- ----- 2 ----- ----- 3 ----- ----- 4 ----- ----- 5 ----- ----- 6 ----- ----- 7 ----- ----- 8 ----- ----- 9 ----- ----- 10 ----- ----- 11 ----- ----- 12 ----- ----- 13 ----- ----- 14 ----- ----- 15 ----- X O O X O X O O O X X X O O O O O X O O O O X X X Elapsed Time(second)=268.465406 <undef> 実行環境 : Ubunts 20.04 (WSL1), Intel Core i5-6200U 2.30GHz
15 手で解くことができましたが、時間がかかりますね。実は、もっと高速に解く方法があるのです。
下図を見てください。
ＡＢＣＤＥ ＡＢＣＤＥ ＡＢＣＤＥ ＡＢＣＤＥ
１□■□■□ １□□□□□ １□□□□□ １□□□□□
２□□□□□ ２■■□■■ ２□□□□□ ２□□□□□
３□□□□□ ３□■□■□ ３□■□■□ ３□□□□□
４□□□□□ ４□□□□□ ４■■□■■ ４□□□□□
５□□□□□ ５□□□□□ ５□□□□□ ５□■□■□
（１） （２） （３） （４）
図 : 1 行ずつボタンを消灯していく方法
(1) では、1 行目のボタンが 2 つ点灯しています。このボタンを消すには、真下にある 2 行目の B と D のボタンを押せばいいですね。すると (2) の状態になります。次に、2 行目のボタンを消します。3 行目の A, B, D, E のボタンを押して (3) の状態になります。
あとはこれを繰り返して 4 行目までのボタンを消したときに、5 行目のボタンも全部消えていれば成功となります。(4) のように、5 行目のボタンが消えない場合は失敗です。この場合は、1 行目のボタンを押して、点灯パターンを変更します。
2 - 5 行目のボタンの押し方は、1 行目の点灯パターンにより決定されるので、けっきょく 1 行目のボタンの押し方により、解けるか否かが決まります。この場合、ボタンの押し方は、2 ^ 5 = 32 通りしかありせん。つまり、たった 32 通り調べるだけでライツアウトの解を求めることができるのです。
それではプログラムを作りましょう。次のリストをみてください。
リスト : ライツアウトの解法 (高速版)
(defun solver-fast (board)
(dotimes (i 32)
(let ((xs i) ; 押したボタン (ビット)
(new-board (push-buttons i board 0))) ; 1 行目を押す (32 通り)
;; 1 行ずつライトを消していく
(dotimes (j 20)
(when (logbitp j new-board)
(setq new-board (logxor (aref pattern (+ j 5)) new-board))
(setq xs (logior xs (ash 1 (+ j 5))))))
(when (= new-board 0)
(print-answer xs)
(format (standard-output) "~%")))))
1 行目のボタンの押し方は 32 通りあるので、ボタンの押し方を 0 から 31 までの数値で表すことにします。これらの値は 5 ビットで表すことができるので、ビットがオンであればそのボタンを押すことにします。押したボタンは変数 xs にセットします。最初に、1 行目のボタンを押して新しい盤面を生成し、それを変数 new-board にセットします。
次は、1 行ずつ new-board のライトを消していきます。変数 j がチェックするボタンの位置を表します。logbitp で j の位置のビットを調べます。それがオンであればライトが点灯しいるので、1 行下のボタンを押します。押すボタンの位置は (+ j 5) で求めることができます。そして最後に new-board の値をチェックします。new-board が 0 であればライトが全部消えているので、print-answer で解を出力します。
これでプログラムは完成です。それでは実行してみましょう。
> (load "bitcomb.lsp") T > (time (solver-fast #x1ffffff)) O O X X X O O X O O X X O O O X O O O X X O O X O O X O O X X O O O X O O O X X O O X O O X X X O O X O O X O X O O O X X X O O O O O X O O O O X X X X X X O O O O X O O O O O X X X O O O X O X O O X Elapsed Time(second)=0.041613 <undef>
4 通りの解が出力されました。インタプリタでも高速に解くことができます。
ボタンを押した回数はどの解も 15 回になります。実は、これがライツアウトの最長手数なのです。ライツアウトの場合、ライトの点灯パターンは 2 ^ 25 = 33554432 通りありますが、実際に解が存在するパターンは、その 1 / 4 の 8388608 通りしかありません。その中で最短回数が 15 回で解けるパターンは 7350 通りあり、そのうちのひとつがライトが全部点灯しているパターンなのです。
ライツアウトの最長手数に興味のある方は、「Puzzle DE Programming: ライツアウト」の「最長手数を求める」をお読みくださいませ。
高橋謙一郎さんの「コンピュータ＆パズル」では、細江万太郎さんが考案されたライツアウトを連立方程式で解く方法が紹介されています。また、拙作のページ「お気楽 Numpy プログラミング超入門」や「Julia Programming: Puzzle DE Julia!!」でも連立方程式によるライツアウトの解法を取り上げています。よろしければお読みくださいませ。
;;;
;;; bitcomb.lsp : ビット操作による組み合わせの生成
;;;
;;; Copyright (C) 2022 Makoto Hiroi
;;;
(import "bit")
(import "macro")
(defun display2 (n)
(format (standard-output) "~B~%" n))
;;; 組み合わせの数
(defun combination-number (n r)
(if (or (= n r) (= r 0))
1
(div (* (combination-number n (- r 1)) (+ (- n r) 1)) r)))
;;; 組み合わせの生成
(defun comb (fn n r a)
(cond
((= r 0) (funcall fn a))
((= n 0) nil)
(t
(comb fn (- n 1) r a)
(comb fn (- n 1) (- r 1) (logior (ash 1 (- n 1)) a)))))
(defun combination (fn n r) (comb fn n r 0))
;;; 組み合わせに番号を付ける
(defun comb-to-num (n r c v)
(cond
((or (= r 0) (= n r)) v)
((logbitp (- n 1) c)
(comb-to-num (- n 1) (- r 1) c (+ (combination-number (- n 1) r) v)))
(t
(comb-to-num (- n 1) r c v))))
(defun to-number (n r c)
(comb-to-num n r c 0))
;;; 番号から組み合わせを求める
(defun num-to-comb (n r v c)
(if (= n r)
(logior (- (expt 2 n) 1) c)
(let ((k (combination-number (- n 1) r)))
(if (>= v k)
(num-to-comb (- n 1) (- r 1) (- v k) (logior (ash 1 (- n 1)) c))
(num-to-comb (- n 1) r v c)))))
(defun from-number (n r v)
(num-to-comb n r v 0))
(defun test (n r)
(combination (lambda (c)
(let ((m (to-number n r c)))
(format (standard-output) "~B => ~D => ~B~%" c m (from-number n r m))))
n r))
;;;
;;; ライツアウトの解法
;;;
;;; ボタンを押したときのパターン
(defglobal pattern
#(#x0000023 #x0000047 #x000008e #x000011c #x0000218
#x0000461 #x00008e2 #x00011c4 #x0002388 #x0004310
#x0008c20 #x0011c40 #x0023880 #x0047100 #x0086200
#x0118400 #x0238800 #x0471000 #x08e2000 #x10c4000
#x0308000 #x0710000 #x0e20000 #x1c40000 #x1880000))
;;; ボタンを押す
(defun push-buttons (xs board n)
(cond
((= xs 0) board)
((logtest xs 1)
(push-buttons (ash xs -1) (logxor (aref pattern n) board) (+ n 1)))
(t
(push-buttons (ash xs -1) board (+ n 1)))))
;;; 解の表示
(defun print-answer (xs)
(dotimes (n 25)
(format (standard-output)
"~A "
(if (logbitp n xs) "O" "X"))
(if (= (mod (+ n 1) 5) 0)
(format (standard-output) "~%"))))
;;; 解法
(defun solver (board)
(block
exit
(for ((n 1 (+ n 1)))
((> n 25))
(format (standard-output) "----- ~D -----~%" n)
(combination
(lambda (xs)
(cond
((= (push-buttons xs board 0) 0)
(print-answer xs)
(return-from exit))))
25 n))))
(defun solver-fast (board)
(dotimes (i 32)
(let ((xs i) ; 押したボタン (ビット)
(new-board (push-buttons i board 0))) ; 1 行目を押す (32 通り)
;; 1 行ずつライトを消していく
(dotimes (j 20)
(when (logbitp j new-board)
(setq new-board (logxor (aref pattern (+ j 5)) new-board))
(setq xs (logior xs (ash 1 (+ j 5))))))
(when (= new-board 0)
(print-answer xs)
(format (standard-output) "~%")))))
n 個の整数 0, 1, ..., n - 1 の順列を考えます。このとき、i 番目の要素が整数 i ではない順列を「完全順列 (derangement)」といいます。たとえば、0, 1, 2 の完全順列は次の 2 通りになります。
1, 2, 0 2, 0, 1
整数 0 は 0 番目に、1 は 1 番目に、 2 は 2 番目に現れていません。これが完全順列になります。そして、完全順列の総数を「モンモール数 (Montmort number)」といいます。モンモール数は次の漸化式で求めることができます。
今回は以下の関数を ISLisp で作ってみましょう。
リスト : 完全順列
(defun derangement-sub (fn m n a)
(if (= n m)
(funcall fn (reverse a))
(for ((x 0 (+ x 1)))
((>= x m))
(if (and (not (= x n))
(not (member x a)))
(derangement-sub fn m (+ n 1) (cons x a))))))
(defun derangement (func m)
(derangement-sub func m 0 nil))
> (derangement #'print 3) (1 2 0) (2 0 1) NIL > (derangement #'print 4) (1 0 3 2) (1 2 3 0) (1 3 0 2) (2 0 3 1) (2 3 0 1) (2 3 1 0) (3 0 1 2) (3 2 0 1) (3 2 1 0) NIL
実際の処理は関数 derangement-sub で行います。derangement-sub は順列を生成するプログラムに条件を追加するだけです。数字 x を選択するとき、その値が位置 n と等しくないことを確認するだけです。
リスト : モンモール数
(defun montmort-number (n)
(cond
((= n 1) 0)
((= n 2) 1)
(t
(* (- n 1) (+ (montmort-number (- n 1))
(montmort-number (- n 2)))))))
;;; 別解
(defun montmort-number1 (n)
(for ((i 1 (+ i 1))
(a 0)
(b 1))
((>= i n) a)
(let ((c a))
(setq a b)
(setq b (* (+ i 1) (+ c b))))))
> (import "macro") T > (dotimes (n 10) (format (standard-output) "~D, ~D~%" (montmort-number (+ n 1)) (montmort-number1 (+ n 1)))) 0, 0 1, 1 2, 2 9, 9 44, 44 265, 265 1854, 1854 14833, 14833 133496, 133496 1334961, 1334961 NIL > (montmort-number1 20) 895014631192902121 > (montmort-number1 40) 300158458444475693321518926221316715906770469041 > (montmort-number1 60) 3061120089030075932414107959955432155359786926462466522532734900460828338404172761
関数 montmort は公式をそのままプログラムしただけです。二重再帰になっているので、実行速度はとても遅くなります。これを繰り返しに変換すると別解のようになります。考え方はフィボナッチ数列と同じです。累積変数 a に i 番目の値を、b に i + 1 番目の値を保存しておきます。すると、i + 2 番目の値は (i + 1) * (a + b)) で計算することができます。あとは、b の値を a に、新しい値を b にセットして処理を繰り返すだけです。
リストや配列で表した集合 xs を分割することを考えます。たとえば、集合 [1, 2, 3] は次のように分割することができます。
1 分割 : [[1, 2, 3]] 2 分割 : [[1, 2], [3]], [[1, 3], [2]], [[1], [2, 3]] 3 分割 : [[1], [2], [3]]
このように、分割した集合 ys は元の集合 xs の部分集合になります。分割した部分集合の積は空集合になり、分割した部分集合のすべての和を求めると元の集合になります。そして、集合を分割する方法の総数を「ベル数 (Bell Number)」といいます。ベル数は次の漸化式で求めることができます。
今回は以下に示す関数を ISLisp で作ってみましょう。
集合を分割するアルゴリズムは簡単です。たとえば、n -1 個の要素 x1, ..., xn-1 を持つ集合を分割したところ、i 個の部分集合 S1, ..., Si が生成されたとしましょう。ここに、n 番目の要素 xn を追加すると、要素が n 個の集合を分割することができます。
新しい要素を追加する場合は次に示す手順で行います。
簡単な例を示しましょう。次の図を見てください。
[] ─ [[1]] ─┬─ [[1, 2]] ─┬─ [[1, 2, 3]]
│ │
│ └─ [[1, 2], [3]]
│
└─ [[1], [2]] ─┬─ [[1, 3] [2]]
│
├─ [[1,], [2, 3]]
│
└─ [[1], [2], [3]]
図 : 集合 [1, 2, 3] を分割する
部分集合を格納する配列を用意します。最初、部分集合は空集合なので空の配列に初期化します。次に、要素 1 を追加します。部分集合は空なので、手順 1 は適用できません。手順 2 を適用して新しい部分集合 [1] を追加します。
次に要素 2 を追加します。[[1]] に 手順 1 を適用すると、部分集合 [1] に要素を追加して [[1, 2]] になります。手順 2 を適用すると、新しい部分集合 [2] を追加して [[1], [2]] になります。最後に 3 を追加します。[[1, 2]] に手順 1 を適用すると [[1, 2, 3]] に、手順 2 を適用すると [[1, 2], [3]] になります。[[1], [2]] に手順 1 を適用すると [[1, 3] [2]] と [[1], [2, 3]] になり、手順 2 を適用すると [[1], [2], [3]] になります。
このように、簡単な方法で集合を分割することができます。実際にプログラムを作る場合、上図を木と考えて、深さ優先で木をたどると簡単です。次のリストを見てください。
リスト : 集合の分割
(defun part-set-sub (fn xs a)
(cond
((null xs)
(funcall fn a))
(t
(mapc (lambda (ys)
(part-set-sub
fn
(cdr xs)
(cons (cons (car xs) (car ys)) (cdr ys))))
(selects a))
(part-set-sub fn (cdr xs) (cons (list (car xs)) a)))))
(defun partition-of-set (fn xs)
(part-set-sub fn (cdr xs) (list (list (car xs)))))
関数 partition-of-set の引数 xs が集合、生成した部分集合は累積変数 a に格納します。xs が空リストの場合、追加する要素がなくなったので fn を評価します。要素がある場合、各々の部分集合に先頭要素 (car xs) を追加して、partition-of-set を再帰呼び出しします。すべての部分集合に要素を追加したら、(car xs) を要素として持つ部分集合を生成して累積変数 a に追加します。
簡単な実行例を示します。
> (partition-of-set #'print (reverse '(a b c))) ((A B C)) ((A) (B C)) ((A B) (C)) ((A C) (B)) ((A) (B) (C)) NIL > (partition-of-set #'print (reverse '(a b c d))) ((A B C D)) ((A) (B C D)) ((A B) (C D)) ((A C D) (B)) ((A) (B) (C D)) ((A B C) (D)) ((A D) (B C)) ((A) (B C) (D)) ((A B D) (C)) ((A C) (B D)) ((A) (B D) (C)) ((A B) (C) (D)) ((A C) (C) (D)) ((A D) (C) (C)) ((A) (B) (C) (D)) NIL
リスト : ベル数
(import "combination")
(import "list")
;;; 添え字付き畳み込み
(defun fold-left-with-index (fn a xs)
(let ((idx (iota 0 (- (length xs) 1))))
(fold-left fn a idx xs)))
;;; ベル数
(defun bell-number (n)
(for ((i 0 (+ i 1))
(bs (list 1)))
((= i n) (car bs))
(setq bs (cons (fold-left-with-index
(lambda (a k x) (+ (* (combination-number i k) x) a))
0
bs)
bs))))
bell-number は公式をそのままプログラムするだけです。累積変数 bs にベル数を逆順で格納します。nＣk はライブラリ combination の関数 combination-number で求めます。
nＣk * B(k) の総和は関数 fold-left-with-index で計算します。fold-left-with-index は添字を関数に渡して畳み込みを行います。ラムダ式の引数 x がリストの要素、k が添字、a が累積変数です。bs は逆順になっていますが、二項係数は nＣi と nＣn - i の値が同じになるので、そのまま計算しても大丈夫です。もちろん、reverse で bs を逆順にしてから計算してもかまいません。
簡単な実行例を示します。
> (import "macro") T > (dotimes (n 11) (print (bell-number n))) 1 1 2 5 15 52 203 877 4140 21147 115975 NIL > (bell-number 20) 51724158235372 > (bell-number 30) 846749014511809332450147 > (bell-number 40) 157450588391204931289324344702531067 > (bell-number 50) 185724268771078270438257767181908917499221852770