簡単なプログラム
●ラテン方陣
「ラテン方陣」は数独の枠の条件を無くした方陣です。ラテン方陣の定義を 参考文献 より引用します。
『ラテン方陣を一般的にいうなら、n 行 n 列の正方形の枡に n 種類の記号を n 個ずつ配列して、各行各列に記号の重複のないものを n 次のラテン方陣というのです。』
このラテン方陣をパズルに応用したものが「ナンバープレース (数独)」というわけです。
簡単な例を示しましょう。3 次のラテン方陣は次に示す 12 通りになります。
0 1 2 0 1 2 0 2 1 0 2 1 1 0 2 1 0 2
1 2 0 2 0 1 1 0 2 2 1 0 0 2 1 2 1 0
2 0 1 1 2 0 2 1 0 1 0 2 2 1 0 0 2 1
標準形
1 2 0 1 2 0 2 0 1 2 0 1 2 1 0 2 1 0
0 1 2 2 0 1 0 1 2 1 2 0 0 2 1 1 0 2
2 0 1 0 1 2 1 2 0 0 1 2 1 0 2 0 2 1
図 : 3 次のラテン方陣
この中で、最初の行と列の要素を昇順に並べたものを「標準形」といいます。3 次のラテン方陣の場合、標準形は 1 種類しかありません。ラテン方陣は任意の行を交換する、または任意の列を交換してもラテン方陣になります。3 次のラテン方陣の場合、標準形から行または列を交換することで、残りの 11 種類のラテン方陣を生成することができます。
今回は標準形ラテン方陣の総数を求めるプログラムを作ります。
リスト : ラテン方陣
function check(mat, x, y, n, z)
for i = 1 : n
if mat[x, i] == z || mat[i, y] == z
return false
end
end
true
end
function solver_sub(mat, x, y, n)
global cnt
if y > n
cnt += 1
elseif x > n
solver_sub(mat, 2, y + 1, n)
else
for z = 1 : n
if check(mat, x, y, n, z)
mat[x, y] = z
solver_sub(mat, x + 1, y, n)
mat[x, y] = 0
end
end
end
end
function solver(n)
global cnt = 0
mat = zeros(Int, n, n)
for i = 1 : n
mat[1, i] = i
mat[i, 1] = i
end
solver_sub(mat, 2, 2, n)
println(cnt)
end
for n = 3 : 7
@time solver(n)
end
プログラムは簡単なので説明は割愛します。実行結果は次のようになりました。
$ julia latin.jl 1 0.037558 seconds (22.72 k allocations: 1.398 MiB, 98.58% compilation time) 4 0.000097 seconds (8 allocations: 432 bytes) 56 0.000204 seconds (8 allocations: 512 bytes) 9408 0.012512 seconds (8.91 k allocations: 139.609 KiB) 16942080 30.332830 seconds (16.94 M allocations: 258.509 MiB, 0.07% gc time) 実行環境 : Julia ver 1.6.4, Windows 10, Intel Core i5-6200U 2.30GHz
単純なプログラムなので実行時間は遅いですね。高次の標準形ラテン方陣の総数は、簡単に求めることができない非常にハードな問題だといわれています。興味のある方は挑戦してみてください。
大村平, 『数理パズルのはなし』, 日科技連出版社, 1998
●三目並べ
皆さんお馴染みのゲーム「三目並べ」で、両者が最善を尽くすと引き分けになることを示すプログラムです。詳しい説明は拙作のページ Puzzle DE Programming 三目並べ をお読みください。
リスト : 三目並べ
# 盤面
# 1 2 3
# 4 5 6
# 7 8 9
# 直線
const line_table = [
[1, 2, 3], [4, 5, 6], [7, 8, 9], [1, 4, 7],
[2, 5, 8], [3, 6, 9], [1, 5, 9], [3, 5, 7]
]
# 盤面
# 0 : 空き場所
# 1 : ○
# -1 : ×
const MARU = 1
const BATU = -1
const board = [0,0,0, 0,0,0, 0,0,0]
const SIZE = 9
const MIN_VALUE = -2
const MAX_VALUE = 2
# ３つ並んだか？
function check_line(x, y, z)
p = board[x]
if p == board[y] == board[z]
p
else
0
end
end
# 勝敗の判定
function check_winner()
for line = line_table
r = check_line(line...)
if r != 0
return r
end
end
0
end
# 先手
function think_maru(n)
value = MIN_VALUE
for i = 1 : SIZE
if board[i] != 0; continue; end
# MARU を書く
board[i] = MARU
# 勝敗の判定
v = check_winner();
# 決着していなければ相手の手番へ
if v == 0 && n < SIZE
v = think_batu(n + 1)
end
# ミニマックス
if v > value; value = v; end
# 元に戻す
board[i] = 0
end
value
end
# 後手
function think_batu(n)
value = MAX_VALUE
for i = 1 : SIZE
if board[i] != 0; continue; end
# BATU を書く
board[i] = BATU;
# 勝敗の判定
v = check_winner()
# 決着していなければ相手の手番へ */
if v == 0 && n < SIZE
v = think_maru(n + 1)
end
# ミニマックス
if v < value; value = v; end
# 元に戻す
board[i] = 0
end
value
end
function solver()
for i in 1 : SIZE
# 初手
board[i] = MARU
# 相手の手番
v = think_batu(2)
# 結果
println("pos = $i, value = $v")
# 元に戻す
board[i] = 0
end
end
@time solver()
$ julia tictactoe.jl pos = 1, value = 0 pos = 2, value = 0 pos = 3, value = 0 pos = 4, value = 0 pos = 5, value = 0 pos = 6, value = 0 pos = 7, value = 0 pos = 8, value = 0 pos = 9, value = 0 0.574310 seconds (174 allocations: 7.922 KiB, 1.22% compilation time) 実行環境 : Julia ver 1.6.4, Windows 10, Intel Core i5-6200U 2.30GHz
ところで、ルールを 3 つ並んだほうが負けに変更すると、初手が中央以外は後手必勝になります。プログラムの修正は関数 check_line() の返り値 r を -r に変更するだけです。
$ julia tictactoe.jl pos = 1, value = -1 pos = 2, value = -1 pos = 3, value = -1 pos = 4, value = -1 pos = 5, value = 0 pos = 6, value = -1 pos = 7, value = -1 pos = 8, value = -1 pos = 9, value = -1 0.562648 seconds (174 allocations: 7.922 KiB, 1.18% compilation time)
●分割数
整数 n を 1 以上の自然数の和で表すことを考えます。これを「整数の分割」といいます。整数を分割するとき、同じ自然数を何回使ってもかまいませんが、並べる順序が違うだけのものは同じ分割とします。そして、分割の仕方の総数を「分割数」といいます。詳しい説明は拙作のページ Puzzle DE Programming 分割数 をお読みください。
●分割数の求め方
整数 n を k 以下の整数で分割する総数を求める関数を \(p(n, k)\) とすると、\(p(n, k)\) は次のように定義することができます。
たとえば、p(6, 6) は次のように計算することができます。
p(6, 6) => p(0, 6) + p(6, 5)
=> 1 + p(1, 5) + p(6, 4)
=> 1 + 1 + p(2, 4) + p(6, 3)
=> 1 + 1 + 2 + 7
=> 11
p(2, 4) => p(-2, 4) + p(2, 3)
=> 0 + p(-1, 3) + p(2, 2)
=> 0 + 0 + p(0, 2) + p(2, 1)
=> 0 + 0 + 1 + 1
=> 2
p(6, 3) => p(3, 3) + p(6, 2)
=> p(0, 3) + p(3, 2) + p(4, 2) + p(6, 1)
=> 1 + p(1, 2) + p(3, 1) + p(2, 2) + p(4, 1) + 1
=> 1 + 1 + 1 + p(0, 2) + p(2, 1) + 1 + 1
=> 1 + 1 + 1 + 1 + 1 + 1 + 1
=> 7
分割数を求める関数を partition_number() とすると、関数 p(n, k) を使って次のようにプログラムすることができます。
リスト : 分割数
function part_num(n, k)
if n == 1 || k == 1
1
elseif n < 0 || k < 1
0
else
part_num(n - k, k) + part_num(n, k - 1)
end
end
partition_number(n) = part_num(n, n)
@time println(partition_number(40))
@time println(partition_number(60))
@time println(partition_number(80))
37338 0.032928 seconds (13.52 k allocations: 839.243 KiB, 97.76% compilation time) 966467 0.006672 seconds (8 allocations: 240 bytes) 15796476 0.098773 seconds (8 allocations: 240 bytes) 実行環境 : Julia ver 1.6.4, Windows 10, Intel Core i5-6200U 2.30GHz
関数 part_num は p(n, k) の定義をそのままプログラムしただけです。ただし、このプログラムは二重再帰で何度も同じ値を求めているため実行速度はとても遅くなります。関数 part_num() をメモ化すると高速化することができますが、大きな値を計算すると Julia のスタックがオーバーフローしてしまいます。
●動的計画法による高速化
動的計画法を使うと、大きな値でも高速に計算することができます。次の図を見てください。
k 1 : [1, 1, 1, 1, 1, 1, 1] 2 : [1, 1, 1+1=2, 1+1=2, 2+1=3, 2+1=3, 3+1=4] => [1, 1, 2, 2, 3, 3, 4] 3: [1, 1, 2, 1+2=3, 1+3=4, 2+3=5, 3+4=7] => [1, 1, 2, 3, 4, 5, 7] 4: [1, 1, 2, 3, 1+4=4, 1+5=6, 2+7=9] => [1, 1, 2, 3, 5, 6, 9 5: [1, 1, 2, 3, 5, 1+6=7, 1+9=10] => [1, 1, 2, 3, 5, 7, 10] 6: [1, 1, 2, 3, 5, 7, 10+1=11] => [1, 1, 2, 3, 5, 7, 11]
最初に大きさ n + 1 のベクタを用意します。上図ではベクタの添字が 0 から始まることに注意してください。ベクタの添字が n を表していて、p(n, 1) から順番に値を求めていきます。p(n, 1) の値は 1 ですから、ベクタの要素は 1 に初期化します。次に、p(n, 2) の値を求めます。定義により p(n, 2) = p(n - 2, 2) + p(n, 1) なので、2 番目以降の要素に n - 2 番目の要素を加算すれば求めることができます。あとは、k の値をひとつずつ増やして同様の計算を行えば p(n, n) の値を求めることができます。
プログラムは次のようになります。
リスト : 分割数 (動的計画法)
function partition_number2(n)
table = fill(big(1), n + 1)
for k = 2 : n, m = k : n
table[m + 1] += table[m - k + 1]
end
table[n + 1]
end
@time println(partition_number2(1000))
@time println(partition_number2(2000))
@time println(partition_number2(4000))
Julia の配列は添字が 1 から始まるので、配列にアクセスするときは +1 していることに注意してください。あとは説明をそのままプログラムしただけなので、とくに難しいところはないと思います。
それでは実際に試してみましょう。1000, 2000, 4000 の分割数を求めてみました。
24061467864032622473692149727991 0.121129 seconds (1.01 M allocations: 26.751 MiB, 19.59% compilation time) 4720819175619413888601432406799959512200344166 0.365793 seconds (4.00 M allocations: 112.827 MiB, 25.82% gc time) 1024150064776551375119256307915896842122498030313150910234889093895 1.338935 seconds (16.00 M allocations: 501.862 MiB, 29.99% gc time)
動的計画法の効果はとても高いですね。ですが、Puzzle DE Programming で使ったプログラミング言語 PyPy より遅くなりました。Julia は多倍長整数 (BigInt) の計算がちょっと苦手なのかもしれません。
●さらなる高速化
ところで、数がもっと大きくなると動的計画法を使ったプログラムでも遅くなります。参考 URL 1 によると、次の漸化式を使うと分割数を高速に求めることができるそうです。
漸化式の説明を 参考 URL 1 より引用します。
『ここで p(0) = 1 および負の整数 k に対して p(k) = 0 とし、和は (1/2)n(3n - 1) の形（ただし n は正または負の整数全体を走る）の一般五角数全体にわたってとるものとする（順に n = 1, -1, 2, -2, 3, -3, 4, -4 ..., とすると、値として 1, 2, 5, 7, 12, 15, 22, 26, 35, 40, 51, ... が得られる）。和における符号は交互に +, +, -, -, +, +, ... と続く。』
分割数 p(k) は k - 1 以下の分割数がわかれば求めることができます。この漸化式も動的計画法を使えば簡単にプログラムできます。次のリストを見てください。
リスト : 分割数 (オイラーの五角数定理)
# 五角数
pentagon(n) = div(n * (3 * n - 1), 2)
function partition_number3(n)
p = zeros(BigInt, n + 1)
p[1] = big(1)
for i = 1 : n
j = 1
s = 1
while true
k = pentagon(j)
if i < k break end
p[i + 1] += p[i - k + 1] * s
k = pentagon(-j)
if i < k break end
p[i + 1] += p[i - k + 1] * s
j += 1
s *= -1
end
end
p[n + 1]
end
@time println(partition_number3(5000))
@time println(partition_number3(10000))
@time println(partition_number3(20000))
配列 p は分割数 p(k) を記憶するために使います。配列の添字は 1 から始まるので、p をアクセスするときは添字を +1 していることに注意してください。p[1] を 1 に初期化したあと、for ループで 1 から n までの分割数を順番に求めていきます。あとは、漸化式をそのままプログラムするだけです。変数 s は符号 (+. -) を表していて、j が奇数のとき s は 1 になり、j が偶数のときは -1 になります。
それでは実際に 5000, 10000, 20000 の分割数を求めてみましょう。
169820168825442121851975101689306431361757683049829233322203824652329144349 0.148793 seconds (1.90 M allocations: 46.527 MiB, 22.44% gc time) 36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144 0.485117 seconds (5.39 M allocations: 150.787 MiB, 21.62% gc time) 25211481381252969791661953323047045228132894960181159343685031410803428442380156495662397073168982436919232 4789351994903016411826230578166735959242113097 1.694703 seconds (15.30 M allocations: 500.695 MiB, 19.64% gc time)
このように、20000 の分割数でも 2 秒未満で求めることができます。
●分割の仕方の列挙
最後に整数の分割の仕方を列挙するプログラムを示します。
リスト : 整数の分割
function part_int_sub(f, n, k, a)
if n == 0
f(a)
elseif n == 1
push!(a, 1)
f(a)
pop!(a)
elseif k == 1
for _ in 1 : n; push!(a, 1); end
f(a)
for _ in 1 : n; pop!(a); end
else
if n >= k
push!(a, k)
part_int_sub(f, n - k, k, a)
pop!(a)
end
part_int_sub(f, n, k - 1, a)
end
end
partition_of_int(f, n) = part_int_sub(f, n, n, Int[])
partition_of_int(println, 5)
partition_of_int(println, 6)
partition_of_int(println, 7)
基本的な考え方は partition_number() と同じです。関数 part_int_sub() は選んだ数値を累積変数 a の配列に格納していくだけです。n が 0 の場合は f(a) を評価し、n が 1 の場合は a に [1] を追加してから f(a) を評価します。k が 1 の場合は n 個の 1 を a に追加してから f(a) を評価します。
5, 6, 7 の分割の仕方は次のようになります。
[5] [4, 1] [3, 2] [3, 1, 1] [2, 2, 1] [2, 1, 1, 1] [1, 1, 1, 1, 1] [6] [5, 1] [4, 2] [4, 1, 1] [3, 3] [3, 2, 1] [3, 1, 1, 1] [2, 2, 2] [2, 2, 1, 1] [2, 1, 1, 1, 1] [1, 1, 1, 1, 1, 1] [7] [6, 1] [5, 2] [5, 1, 1] [4, 3] [4, 2, 1] [4, 1, 1, 1] [3, 3, 1] [3, 2, 2] [3, 2, 1, 1] [3, 1, 1, 1, 1] [2, 2, 2, 1] [2, 2, 1, 1, 1] [2, 1, 1, 1, 1, 1] [1, 1, 1, 1, 1, 1, 1]
●参考 URL
- 分割数 - Wikipedia
- 自然数の分割 - Wikipedia
- 分割数 - 桃の天然水, (inamori さん)
●集合の分割 (1)
配列で表した集合 ls を分割することを考えます。たとえば、集合 [1, 2, 3] は次のように分割することができます。
1 分割 : [[1, 2, 3]] 2 分割 : [[1, 2], [3]], [[1, 3], [2]], [[1], [2, 3]] 3 分割 ; [[1], [2], [3]]
このように、分割した集合 xs は元の集合 ls の部分集合になります。分割した部分集合の積は空集合になり、分割した部分集合のすべての和を求めると元の集合になります。今回は分割の仕方をすべて求める関数 parititon_of_set() を作ります。
集合を分割するアルゴリズムは簡単です。たとえば、n -1 個の要素 x1, ..., xn-1 を持つ集合を分割したところ、i 個の部分集合 S1, ..., Si が生成されたとしましょう。ここに、n 番目の要素 xn を追加すると、要素が n 個の集合を分割することができます。
新しい要素を追加する場合は次に示す手順で行います。
- 部分集合 Sk (k = 1 から i まで) に要素 xn を追加する
- 新しい部分集合 Si+1 (要素が xn だけの集合) を生成する
簡単な例を示しましょう。次の図を見てください。
部分集合を格納する配列を用意します。最初、部分集合は空集合なので空の配列に初期化します。次に、要素 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]] になります。
このように、簡単な方法で集合を分割することができます。実際にプログラムを作る場合、上図を木と考えて、深さ優先で木をたどると簡単です。次のリストを見てください。
リスト : 集合の分割
function part_sub(f, ls, a)
if length(ls) == 0
f(a)
else
for i = 1 : length(a)
push!(a[i], ls[1])
part_sub(f, ls[2:end], a)
pop!(a[i])
end
push!(a, [ls[1]])
part_sub(f, ls[2:end], a)
pop!(a)
end
end
partition_of_set(f, ls) = part_sub(f, ls[2:end], [[ls[1]]])
partition_of_set(println, [1, 2, 3])
println("")
partition_of_set(println, [1, 2, 3, 4])
実際の処理は関数 part_sub() で行います。生成した部分集合は累積変数 a に格納します。ls が空の配列の場合、追加する要素がなくなったので f(a) を評価します。要素がある場合、i 番目の部分集合に要素 ls[1] を追加して、part_sub() を再帰呼び出しします。すべての部分集合に要素を追加したら、ls[1] を要素として持つ部分集合を生成して累積変数 a に追加します。
[[1, 2, 3]] [[1, 2], [3]] [[1, 3], [2]] [[1], [2, 3]] [[1], [2], [3]] [[1, 2, 3, 4]] [[1, 2, 3], [4]] [[1, 2, 4], [3]] [[1, 2], [3, 4]] [[1, 2], [3], [4]] [[1, 3, 4], [2]] [[1, 3], [2, 4]] [[1, 3], [2], [4]] [[1, 4], [2, 3]] [[1], [2, 3, 4]] [[1], [2, 3], [4]] [[1, 4], [2], [3]] [[1], [2, 4], [3]] [[1], [2], [3, 4]] [[1], [2], [3], [4]]
●集合の分割 (2)
k 個の要素をもつ集合 ls を要素数が等しい m 個の部分集合に分割することを考えます。部分集合の要素数 n は k / m になります。今回は分割の仕方をすべて求める高階関数 group_partition() を作ります。
リスト : 集合の分割 (2)
function group_part_sub(f, ls, n, m, a)
if length(ls) == 0
f(a)
else
for i = 1 : length(a)
if length(a[i]) < n
push!(a[i], ls[1])
group_part_sub(f, ls[2:end], n, m, a)
pop!(a[i])
end
end
if length(a) < m
push!(a, [ls[1]])
group_part_sub(f, ls[2:end], n, m, a)
pop!(a)
end
end
end
group_partition(f, n, m, ls) = group_part_sub(f, ls[2:end], n, m, [[ls[1]]])
group_partition(println, 2, 2, [1,2,3,4])
group_partition(println, 2, 3, [1,2,3,4,5,6])
group_partition() は partition_of_set() を改造するだけで簡単に作成することができます。生成する部分集合の大きさを n に、部分集合の個数を m に制限するだけです。i 番目の部分集合に要素を追加する場合、length(a[i]) が n 未満であることをチェックします。新しい部分集合を追加する場合、length(a) が m 未満であることをチェックします。これで集合をグループに分けることができます。
[[1, 2], [3, 4]] [[1, 3], [2, 4]] [[1, 4], [2, 3]] [[1, 2], [3, 4], [5, 6]] [[1, 2], [3, 5], [4, 6]] [[1, 2], [3, 6], [4, 5]] [[1, 3], [2, 4], [5, 6]] [[1, 3], [2, 5], [4, 6]] [[1, 3], [2, 6], [4, 5]] [[1, 4], [2, 3], [5, 6]] [[1, 5], [2, 3], [4, 6]] [[1, 6], [2, 3], [4, 5]] [[1, 4], [2, 5], [3, 6]] [[1, 4], [2, 6], [3, 5]] [[1, 5], [2, 4], [3, 6]] [[1, 6], [2, 4], [3, 5]] [[1, 5], [2, 6], [3, 4]] [[1, 6], [2, 5], [3, 4]]
●カークマンの 15 人の女生徒
15 人の女生徒が毎日 3 人ずつ 5 組に分かれて散歩をするとき、1 週間 (7 日) のうちに、どの女生徒も他のすべての女生徒と 1 回ずつ同じ組になるような組み合わせを作ってください。
出典 : 大村平 (著), 『数理パズルの話』, 日科技連出版社, 1998
「カークマンの 15 人の女生徒」の解法プログラムは group_partition() を改造することで簡単に作成することができます。
リスト : カークマンの 15 人の女生徒
const check_table = Vector{Int}[
[], [], [], [], [],
[], [], [], [], [],
[], [], [], [], [],
]
function check_student(xs, y)
for x = xs
if y in check_table[x]; return false; end
end
true
end
function add_student(xs, y)
for x = xs
push!(check_table[x], y)
push!(check_table[y], x)
end
end
function del_student(xs, y)
for x = xs
pop!(check_table[x])
pop!(check_table[y])
end
end
function kirakman_sub(ls, a, b)
if length(ls) == 0
push!(b, a)
if length(b) == 7
println(b)
throw("found!!!")
else
kirakman_sub(collect(2:15), [[1]], b)
end
pop!(b)
else
for i = 1 : length(a)
if length(a[i]) < 3 && check_student(a[i], ls[1])
add_student(a[i], ls[1])
push!(a[i], ls[1])
kirakman_sub(ls[2:end], a, b)
pop!(a[i])
del_student(a[i], ls[1])
end
end
if length(a) < 5
push!(a, [ls[1]])
kirakman_sub(ls[2:end], a, b)
pop!(a)
end
end
end
function kirkman()
try
kirakman_sub(collect(2:15), [[1]], Vector{Vector{Int}}[])
catch e
println(e)
end
end
@time kirkman()
$ julia kirkman.jl [[[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]], [[1, 4, 7], [2, 5, 10], [3, 6, 13], [8, 11, 14], [9, 12, 15]], [[1, 5, 14], [2, 4, 15], [3, 8, 12], [6, 9, 11], [7, 10, 13]], [[1, 9, 13], [2, 7, 12], [3, 4, 11], [5, 8, 15], [6, 10, 14]], [[1, 8, 10], [2, 11, 13], [3, 5, 9], [4, 12, 14], [6, 7, 15]], [[1, 6, 12], [2, 9, 14], [3, 10, 15], [4, 8, 13], [5, 7, 11]], [[1, 11, 15], [2, 6, 8], [3, 7, 14], [4, 9, 10], [5, 12, 13]]] found!!! 7.981730 seconds (34.62 M allocations: 3.647 GiB, 3.64% gc time, 5.05% compilation time) 実行環境 : Julia ver 1.6.4, Windows 10, Intel Core i5-6200U 2.30GHz