今回は Prolog で「ナンバープレース (数独)」の解法プログラムを作ってみましょう。最近のパソコンはハイスペックなので、9 行 9 列盤のナンバープレースであれば、特に工夫しなくても単純な深さ優先探索で解くことができます。ただし、普通の Prolog は配列をサポートしていないので、他のプログラミング言語よりもちょっとだけ難しくなります。
SWI-Prolog の場合、一番簡単な方法は制約論理プログラミング用のライブラリ clpfd を使うことです。clpfd のマニュアルにはナンバープレースの解法プログラムが掲載されています。プログラムはとても簡単、高速にパズルを解くことができます。実際に動かしてみると、制約プログラミングのパワーを実感することができるでしょう。制約論理プログラミングについては拙作のページ 制約論理プログラミング超入門 をお読みくださいませ。
clpfd を使わない場合、いろいろな方法が考えられると思いますが、盤面をリストで表すことにすると、空き場所を数字 0 で表すよりも自由変数で表したほうが Prolog らしいプログラムになります。そして、空き場所に置くことができる数字のリストと、空き場所が属する縦横枠の自由変数のリストをあらかじめ求めておいて、深さ優先探索で矛盾しないように自由変数の値を決定していくことにします。
数字のリストを求める場合、集合演算を使うと簡単ですが、自由変数のリストを求めるときは注意が必要です。Prolog の自由変数はどんな値にでもマッチングするので、今回の処理に SWI-Prolog の集合演算を使うことはできません。簡単な例を示しましょう。
?- union([1,2,3,4], [3,4,5,6], X). X = [1, 2, 3, 4, 5, 6]. ?- union([A,B,C,D], [C,D,E,F], X). A = B, B = C, C = D, X = [D, D, E, F].
和集合を求める述語 union/3 は、要素が数字だと正常に動作しますが、自由変数だと正常に動作しません。この場合、自由変数のリスト [A,B,C,D,E,F] がほしいわけです。
そこで、等値の判定に述語 == を使う集合演算を定義することにします。次のリストを見てください。
リスト : 集合演算
% 述語 == による member
memq(X, [Y | _]) :- X == Y.
memq(X, [_ | Ys]) :- memq(X, Ys).
% memq による集合述語
% 和集合
unionq([], Ys, Ys).
unionq([X | Xs], Ys, Zs) :- memq(X, Ys), !, unionq(Xs, Ys, Zs).
unionq([X | Xs], Ys, [X | Zs]) :- unionq(Xs, Ys, Zs).
% 積集合
intersectq([], _, []).
intersectq([X | Xs], Ys, [X | Zs]) :-
memq(X, Ys), !, intersectq(Xs, Ys, Zs).
intersectq([_ | Xs], Ys, Zs) :- intersectq(Xs, Ys, Zs).
% 差集合
differenceq([], _, []).
differenceq([X | Xs], Ys, Zs) :-
memq(X, Ys), !, differenceq(Xs, Ys, Zs).
differenceq([X | Xs], Ys, [X | Zs]) :- differenceq(Xs, Ys, Zs).
述語 memq は演算子 == で等値を判定します。名前は Scheme から拝借しました。member と違って第 1 引数を自由変数にしてもマッチングは行われません。あとは、拙作のページ Yet Another Prolog Problems (1) の問題 18, 19, 20 のプログラムで、member を memq に変更するだけです。
簡単な実行例を示しましょう。
?- memq(1, [1,2,3,4,5]). true ; false. ?- memq(A, [1,2,3,4,5]). false. ?- unionq([1,2,3,4],[3,4,5,6], X). X = [1, 2, 3, 4, 5, 6]. ?- unionq([A,B,C,D], [C,D,E,F], X). X = [A, B, C, D, E, F]. ?- intersectq([1,2,3,4],[3,4,5,6], X). X = [3, 4]. ?- intersectq([A,B,C,D], [C,D,E,F], X). X = [C, D]. ?- differenceq([1,2,3,4],[3,4,5,6], X). X = [1, 2]. ?- differenceq([A,B,C,D], [C,D,E,F], X). X = [A, B].
正常に動作していますね。
次は盤面を表すデータ構造を定義しましょう。最初にナンバープレースの盤面 (9 行 9 列) を下図に示します。
列 0 1 2 3 4 5 6 7 8 行 +-------+-------+-------+ 0 | | | | 1 | 枠 0 | 1 | 2 | 2 | | | | +-------+-------+-------+ 3 | | | | 4 | 3 | 4 | 5 | 5 | | | | +-------+-------+-------+ 6 | | | | 7 | 6 | 7 | 8 | 8 | | | | +-------+-------+-------+ 図 : 数独 (9 * 9) の盤面
この盤面を次のようにリストのリストで表すことにすると、横 (行) の関係はすぐに求めることができます。
リスト : 問題 (出典: 数独 - Wikipedia の問題例)
problem(0, [[5,3,_, _,7,_, _,_,_],
[6,_,_, 1,9,5, _,_,_],
[_,9,8, _,_,_, _,6,_],
[8,_,_, _,6,_, _,_,3],
[4,_,_, 8,_,3, _,_,1],
[7,_,_, _,2,_, _,_,6],
[_,6,_, _,_,_, 2,8,_],
[_,_,_, 4,1,9, _,_,5],
[_,_,_, _,8,_, _,7,9]]).
要素のリストの中で 1 から 9 の数字が重複せずに一つずつ入ればいいわけです。
次は、縦 (列) の関係を求める述語を作ります。盤面を行列と考えると、列の関係は転置行列を求めることと同じになります。次の図を見てください。
[[1,2,3], [[1,4,7], [4,5,6], = 転置行列 => [2,5,8], [7,8,9]] [3,6,9]]
このように、行列の行と列を入れ替えた行列を「転置行列 (transposed matrix)」といいます。SWI-Prolog のライブラリ clpfd には転置行列を求める述語 transpose が用意されていますが、maplist を使うと私たちでも簡単に定義することができます。次のリストを見てください。
リスト : 転置行列
head([X | _], X).
tail([_ | Xs], Xs).
transpose(Xs, []) :- member([], Xs), !.
transpose(Xs, [Y | Ys]) :-
maplist(head, Xs, Y),
maplist(tail, Xs, Xs1),
transpose(Xs1, Ys).
述語 head はリストの先頭要素を取り出します。tail はリストの先頭要素を取り除きます。これらの述語は Lisp / Scheme の car, cdr と同じです。
transpose の最初の節が再帰呼び出しの停止条件で、要素のリストが空リストになったか member でチェックします。次の節で maplist に head を渡して、各リストの先頭要素を格納したリスト、つまり列を表すリスト Y を作ります。そして、maplist に tail を渡して、先頭要素を取り除いたリストを格納したリスト Xs1 を作ります。この Xs1 に transpose を適用すれば、次の列の要素を格納したリストを作ることができます。
簡単な実行例を示します。
?- transpose([[1,2],[3,4]], Xs). Xs = [[1, 3], [2, 4]]. ?- transpose([[1,2],[3,4],[5,6]], Xs). Xs = [[1, 3, 5], [2, 4, 6]]. ?- transpose([[1,2,3],[4,5,6],[7,8,9]], Xs). Xs = [[1, 4, 7], [2, 5, 8], [3, 6, 9]].
次は枠 (ブロック) の関係を求めるプログラムを作りましょう。次のリストを見てください。
リスト : 枠 (ブロック) の関係を求める
make_block_sub([],[],[],[]).
make_block_sub([X1, X2, X3 | Xs],
[Y1, Y2, Y3 | Ys],
[Z1, Z2, Z3 | Zs],
[[X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3] | Bs]) :-
make_block_sub(Xs, Ys, Zs, Bs).
make_block([], []).
make_block([X, Y, Z | Ls], Gs) :-
make_block_sub(X, Y, Z, Gs1),
make_block(Ls, Gs2),
append(Gs1, Gs2, Gs).
述語 make_block はリストの先頭から 3 行 (X, Y, Z) ずつ取り出して述語 make_block_sub に渡します。make_block_sub では、各リストの先頭から 3 つの要素を取り出して、それらを一つのリストに格納します。これで同じ枠内にある要素を一つのリストにまとめることができます。
それでは簡単なテストを行ってみましょう。
リスト : テスト
test :-
Ls = [
[11, 12, 13, 14, 15, 16, 17, 18, 19],
[21, 22, 23, 24, 25, 26, 27, 28, 29],
[31, 32, 33, 34, 35, 36, 37, 38, 39],
[41, 42, 43, 44, 45, 46, 47, 48, 49],
[51, 52, 53, 54, 55, 56, 57, 58, 59],
[61, 62, 63, 64, 65, 66, 67, 68, 69],
[71, 72, 73, 74, 75, 76, 77, 78, 79],
[81, 82, 83, 84, 85, 86, 87, 88, 89],
[91, 92, 93, 94, 95, 96, 97, 98, 99]],
maplist(writeln, Ls),
nl,
transpose(Ls, Cs),
maplist(writeln, Cs),
nl,
make_block(Ls, Gs),
maplist(writeln, Gs).
?- test. [11,12,13,14,15,16,17,18,19] [21,22,23,24,25,26,27,28,29] [31,32,33,34,35,36,37,38,39] [41,42,43,44,45,46,47,48,49] [51,52,53,54,55,56,57,58,59] [61,62,63,64,65,66,67,68,69] [71,72,73,74,75,76,77,78,79] [81,82,83,84,85,86,87,88,89] [91,92,93,94,95,96,97,98,99] [11,21,31,41,51,61,71,81,91] [12,22,32,42,52,62,72,82,92] [13,23,33,43,53,63,73,83,93] [14,24,34,44,54,64,74,84,94] [15,25,35,45,55,65,75,85,95] [16,26,36,46,56,66,76,86,96] [17,27,37,47,57,67,77,87,97] [18,28,38,48,58,68,78,88,98] [19,29,39,49,59,69,79,89,99] [11,12,13,21,22,23,31,32,33] [14,15,16,24,25,26,34,35,36] [17,18,19,27,28,29,37,38,39] [41,42,43,51,52,53,61,62,63] [44,45,46,54,55,56,64,65,66] [47,48,49,57,58,59,67,68,69] [71,72,73,81,82,83,91,92,93] [74,75,76,84,85,86,94,95,96] [77,78,79,87,88,89,97,98,99] true.
正常に動作していますね。
次は空き場所に置くことができる数字と、空き場所が属する縦横枠にある自由変数 (空き場所) を求める述語 analysis を作りましょう。次のリストを見てください。
リスト : 盤面の解析
% 数字を取り出す
get_number(X, Y, Ls, N) :- nth0(Y, Ls, L), nth0(X, L, N).
% 解析
analysis(_, 9, _, _, _, []).
analysis(9, Y, Ls, Cs, Gs, Zs) :-
Y1 is Y + 1, analysis(0, Y1, Ls, Cs, Gs, Zs).
analysis(X, Y, Ls, Cs, Gs, Zs) :-
get_number(X, Y, Ls, N),
nonvar(N),
X1 is X + 1,
analysis(X1, Y, Ls, Cs, Gs, Zs).
analysis(X, Y, Ls, Cs, Gs, [Z | Zs]) :-
get_number(X, Y, Ls, N),
var(N),
analysis_sub(X, Y, N, Ls, Cs, Gs, Z),
X1 is X + 1,
analysis(X1, Y, Ls, Cs, Gs, Zs).
analysys(X, Y, Ls, Cs, Gs, Zs) の引数 X, Y は盤面の列と行、Ls が盤面 (行)、Cs が盤面を転置したもの (列)、Gs が枠を表します。Zs の要素はリストで要素は次のようになります。
[空き場所 (自由変数), 置くことができる数字のリスト, 空き場所が属する縦横枠の自由変数のリスト]
このリストは述語 analysys_sub で作ります。最初に、述語 get_number で (X, Y) にある数字 N を求めます。N が自由変数でなければ、次の場所を調べます。自由変数の場合、analysis_sub で数字と自由変数のリストを求めます。
次は述語 analysis_sub を作ります。プログラムは次のようになります。
リスト : 盤面の解析 (2)
analysis_sub(X, Y, N, Ls, Cs, Gs, [N, As, Bs]) :-
nth0(Y, Ls, Ys),
nth0(X, Cs, Xs),
G is (Y // 3) * 3 + X // 3,
nth0(G, Gs, Zs),
% 数字と変数に分ける
partition(integer, Ys, Ys1, Ys2),
partition(integer, Xs, Xs1, Xs2),
partition(integer, Zs, Zs1, Zs2),
unionq(Ys1, Xs1, As1),
unionq(Zs1, As1, As2),
% 未確定の数字
differenceq([1,2,3,4,5,6,7,8,9], As2, As),
% 変数のリスト
unionq(Ys2, Xs2, Bs1),
unionq(Zs2, Bs1, Bs).
最初に nth0 で Y 行のリストを Ys に、X 列のリストを Xs に、G 番目の枠のリストを Zs にもと目増す。次に、それぞれのリストの要素を述語 partition/4 で数字と自由変数に分けます。partition の説明は拙作のページ 高階プログラミング をお読みくださいませ。
リスト Xs, Ys, Zs を partition で Xs1, Xs2, Ys1, Ys2, Zs1, Zs2 に分離します。ヒントの数字は Xs1, Ys1, Zs1 の和集合を unionq で求めるだけです。この場合、SWI-Prolog の述語 union を使っても問題ありません。この値を As2 とすると、置くことができる数字は、1 から 9 までの数字の集合から As2 を引き算する、つまり differenceq で差集合を求めるだけです。自由変数のリストは Xs2, Ys2, Zs2 の和集合を unionq で求めるだけです。この場合、SWI-Prolog の union を使ってはいけません。ご注意ください。
ここまで準備が整ったら、あとは簡単です。ナンバープレースの解法プログラムは次のようになります。
リスト : ナンバープレースの解法
% 深さ優先探索
dfs([]).
dfs([[N, As, Bs] | Vs]) :-
member(I, As),
maplist(\==(I), Bs),
N = I,
dfs(Vs).
% 解法
solver(N) :-
problem(N, Ls),
transpose(Ls, Cs),
make_block(Ls, Gs),
analysis(0, 0, Ls, Cs, Gs, Vs),
dfs(Vs),
maplist(writeln, Ls).
述語 solver の引数 N は問題番号を表します。problem から問題 (盤面) を取り出して、transpose で列の関係を、make_block で枠の関係を求めます。次に、述語 analysis で盤面を解析して、述語 dfs で深さ優先探索します。解が見つかったら maplist で盤面を表示します。
述語 dfs も簡単です。引数 N が空き場所を表す自由変数、As が置くことができる数字のリスト、Bs が自由変数のリストです。最初に、member で As から数字 I を選びます。次に、Bs の中で I と同じ値がないか maplist でチェックします。最初、Bs の要素は自由変数しかありませんが、探索を進めていくと数字と自由変数が混在するようになります。述語には \== を使うことに注意してください。
同じ値がない場合は N = I で N の値を I に決定し、次の空き場所の数字を決めるため dfs を再帰呼び出しします。引数 Vs が空リストになったならば解を求めることができました。これが再帰呼び出しの停止条件になります。
それでは実行してみましょう。
?- time(solver(0)). [5,3,4,6,7,8,9,1,2] [6,7,2,1,9,5,3,4,8] [1,9,8,3,4,2,5,6,7] [8,5,9,7,6,1,4,2,3] [4,2,6,8,5,3,7,9,1] [7,1,3,9,2,4,8,5,6] [9,6,1,5,3,7,2,8,4] [2,8,7,4,1,9,6,3,5] [3,4,5,2,8,6,1,7,9] % 257,012 inferences, 0.034 CPU in 0.035 seconds (98% CPU, 7575796 Lips) true ; % 25,008 inferences, 0.013 CPU in 0.013 seconds (93% CPU, 1992928 Lips) false. 実行環境 : Lubuntu 16.10 on VirtualBox, Core i7-2670QM 2.20GHz, SWI-Prolog Version 7.2.3
0.1 秒もかからずに解くことができました。ただし、問題によっては時間がかかることがあります。Puzzle Generater Japan にある Java版標準問題集 の問題 8-a, 8-b, 9-a, 9-b, 10-a, 10-b を解いてみました。解を見つけるまでの時間を示します。セミコロン ( ; ) を入力してバックトラックすると、もっと時間がかかることに注意してください。
表 : 実行時間 (秒)
問題 : Prolog : clpfd
------+--------+-------
8-a : 0.32 : 0.11
8-b : 0.88 : 0.13
9-a : 2.37 : 0.11
9-b : 0.65 : 0.11
10-a : 0.12 : 0.12
10-b : 0.34 : 0.10
実行環境 : Lubuntu 16.10 on VirtualBox, Core i7-2670QM 2.20GHz, SWI-Prolog Version 7.2.3
今回のプログラムは、問題 9-a を解くのに 3 秒ちかくかかりますが、clpfd はどの問題でも 0.1 秒程度で解くことができました。 制約プログラミングは凄いですね。この結果には M.Hiroi も大変驚きました。興味のある方はいろいろ試してみてください。
%
% numplace.swi : ナンバープレースの解法
%
% Copyright (C) 2016 Makoto Hiroi
%
% 述語 == による member
memq(X, [Y | _]) :- X == Y.
memq(X, [_ | Ys]) :- memq(X, Ys).
% memq による集合述語
unionq([], Ys, Ys).
unionq([X | Xs], Ys, Zs) :- memq(X, Ys), !, unionq(Xs, Ys, Zs).
unionq([X | Xs], Ys, [X | Zs]) :- unionq(Xs, Ys, Zs).
intersectq([], _, []).
intersectq([X | Xs], Ys, [X | Zs]) :-
memq(X, Ys), !, intersectq(Xs, Ys, Zs).
intersectq([_ | Xs], Ys, Zs) :- intersectq(Xs, Ys, Zs).
differenceq([], _, []).
differenceq([X | Xs], Ys, Zs) :-
memq(X, Ys), !, differenceq(Xs, Ys, Zs).
differenceq([X | Xs], Ys, [X | Zs]) :- differenceq(Xs, Ys, Zs).
% 行列の転置
head([X | _], X).
tail([_ | Xs], Xs).
transpose(Xs, []) :- member([], Xs), !.
transpose(Xs, [Y | Ys]) :-
maplist(head, Xs, Y),
maplist(tail, Xs, Xs1),
transpose(Xs1, Ys).
% ブロックの生成
make_block_sub([],[],[],[]).
make_block_sub([X1, X2, X3 | Xs],
[Y1, Y2, Y3 | Ys],
[Z1, Z2, Z3 | Zs],
[[X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3] | Bs]) :-
make_block_sub(Xs, Ys, Zs, Bs).
%
make_block([], []).
make_block([X, Y, Z | Ls], Gs) :-
make_block_sub(X, Y, Z, Gs1),
make_block(Ls, Gs2),
append(Gs1, Gs2, Gs).
% 数字を取り出す
get_number(X, Y, Ls, N) :- nth0(Y, Ls, L), nth0(X, L, N).
% 解析
analysis_sub(X, Y, N, Ls, Cs, Gs, [N, As, Bs]) :-
nth0(Y, Ls, Ys),
nth0(X, Cs, Xs),
G is (Y // 3) * 3 + X // 3,
nth0(G, Gs, Zs),
% 数字と変数に分ける
partition(integer, Ys, Ys1, Ys2),
partition(integer, Xs, Xs1, Xs2),
partition(integer, Zs, Zs1, Zs2),
unionq(Ys1, Xs1, As1),
unionq(Zs1, As1, As2),
% 未確定の数字
differenceq([1,2,3,4,5,6,7,8,9], As2, As),
% 変数のリスト
unionq(Ys2, Xs2, Bs1),
unionq(Zs2, Bs1, Bs).
analysis(_, 9, _, _, _, []).
analysis(9, Y, Ls, Cs, Gs, Zs) :-
Y1 is Y + 1, analysis(0, Y1, Ls, Cs, Gs, Zs).
analysis(X, Y, Ls, Cs, Gs, Zs) :-
get_number(X, Y, Ls, N),
nonvar(N),
X1 is X + 1,
analysis(X1, Y, Ls, Cs, Gs, Zs).
analysis(X, Y, Ls, Cs, Gs, [Z | Zs]) :-
get_number(X, Y, Ls, N),
var(N),
analysis_sub(X, Y, N, Ls, Cs, Gs, Z),
X1 is X + 1,
analysis(X1, Y, Ls, Cs, Gs, Zs).
% 深さ優先探索
dfs([]).
dfs([[N, As, Bs] | Vs]) :-
member(I, As),
maplist(\==(I), Bs),
N = I,
dfs(Vs).
% 解法
solver(N) :-
problem(N, Ls),
transpose(Ls, Cs),
make_block(Ls, Gs),
analysis(0, 0, Ls, Cs, Gs, Vs),
dfs(Vs),
maplist(writeln, Ls),
fail.
% 問題 (出典: 数独 - Wikipedia の問題例)
problem(0, [[5,3,_, _,7,_, _,_,_],
[6,_,_, 1,9,5, _,_,_],
[_,9,8, _,_,_, _,6,_],
[8,_,_, _,6,_, _,_,3],
[4,_,_, 8,_,3, _,_,1],
[7,_,_, _,2,_, _,_,6],
[_,6,_, _,_,_, 2,8,_],
[_,_,_, 4,1,9, _,_,5],
[_,_,_, _,8,_, _,7,9]]).