自然数 n を素因数分解する関数 factorization n を定義してください。返り値はリスト [(p, q); ...] で、(p, q) は \(p^q\) を表します。なお、整数 n, p, q はモジュール Num で表すものとします。
val factorization : Num.num -> (Num.num * Num.num) list = <fun>
# factorization (num_of_string "12345678");; - : (Num.num * Num.num) list = [(Int 2, Int 1); (Int 3, Int 2); (Int 47, Int 1); (Int 14593, Int 1)] # factorization (num_of_string "123456789");; - : (Num.num * Num.num) list = [(Int 3, Int 2); (Int 3607, Int 1); (Int 3803, Int 1)] # factorization (num_of_string "1234567890");; - : (Num.num * Num.num) list = [(Int 2, Int 1); (Int 3, Int 2); (Int 5, Int 1); (Int 3607, Int 1); (Int 3803, Int 1)] # factorization (num_of_string "1111111111");; - : (Num.num * Num.num) list = [(Int 11, Int 1); (Int 41, Int 1); (Int 271, Int 1); (Int 9091, Int 1)]
自然数 n の約数の個数を求める関数 divisor_num を定義してください。
val divisor_num : Num.num -> Num.num = <fun>
# divisor_num (num_of_string "12345678");; - : Num.num = Int 24 # divisor_num (num_of_string "123456789");; - : Num.num = Int 12 # divisor_num (num_of_string "1234567890");; - : Num.num = Int 48 # divisor_num (num_of_string "1111111111");; - : Num.num = Int 16
自然数 n の約数の合計値を求める関数 divisor_sum を定義してください。
val divisor_sum : Num.num -> Num.num = <fun>
# divisor_sum (num_of_string "12345678");; - : Num.num = Int 27319968 # divisor_sum (num_of_string "123456789");; - : Num.num = Int 178422816 # divisor_sum (num_of_string "1234567890");; - : Num.num = Big_int <abstr>
自然数 n の約数をリストに格納して返す関数 divisor を定義してください。
val divisor : Num.num -> Num.num list = <fun>
# divisor (Int 12);; - : Num.num list = [Int 1; Int 2; Int 3; Int 4; Int 6; Int 12] # divisor (Int 12345678);; - : Num.num list = [Int 1; Int 2; Int 3; Int 6; Int 9; Int 18; Int 47; Int 94; Int 141; Int 282; Int 423; Int 846; Int 14593; Int 29186; Int 43779; Int 87558; Int 131337; Int 262674; Int 685871; Int 1371742; Int 2057613; Int 4115226; Int 6172839; Int 12345678] # divisor (Int 123456789);; - : Num.num list = [Int 1; Int 3; Int 9; Int 3607; Int 3803; Int 10821; Int 11409; Int 32463; Int 34227; Int 13717421; Int 41152263; Int 123456789] # divisor (num_of_string "1234567890");; - : Num.num list = [Int 1; Int 2; Int 3; Int 5; Int 6; Int 9; Int 10; Int 15; Int 18; Int 30; Int 45; Int 90; Int 3607; Int 3803; Int 7214; Int 7606; Int 10821; Int 11409; Int 18035; Int 19015; Int 21642; Int 22818; Int 32463; Int 34227; Int 36070; Int 38030; Int 54105; Int 57045; Int 64926; Int 68454; Int 108210; Int 114090; Int 162315; Int 171135; Int 324630; Int 342270; Int 13717421; Int 27434842; Int 41152263; Int 68587105; Int 82304526; Int 123456789; Int 137174210; Int 205761315; Int 246913578; Int 411522630; Int 617283945; Big_int] # divisor (num_of_string "1111111111");; - : Num.num list = [Int 1; Int 11; Int 41; Int 271; Int 451; Int 2981; Int 9091; Int 11111; Int 100001; Int 122221; Int 372731; Int 2463661; Int 4100041; Int 27100271; Int 101010101; Big_int ]
完全数 - Wikipedia によると、『完全数（かんぜんすう，perfect number）とは、その数自身を除く約数の和が、その数自身と等しい自然数のことである。』 とのことです。自然数 n 以下の完全数を求める関数 perfect_number を定義してください。
val perfect_number : int -> unit = <fun>
# perfect_number 10000;; 6 28 496 8128 - : unit = ()
友愛数 - Wikipedia によると、『友愛数（ゆうあいすう）とは、異なる2つの自然数の組で、自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。』 とのことです。自然数 n 以下の友愛数を求める関数 yuuai_number を定義してください。
val yuuai_number : int -> unit = <fun>
# yuuai_number 100000;; 220 284 1184 1210 2620 2924 5020 5564 6232 6368 10744 10856 12285 14595 17296 18416 63020 76084 66928 66992 67095 71145 69615 87633 79750 88730 - : unit = ()
整数 n を 1 以上の自然数の和で表すことを考えます。これを「整数の分割」といいます。整数を分割するとき、同じ自然数を何回使ってもかまいませんが、並べる順序が違うだけのものは同じ分割とします。簡単な例を示します。
n = 6
6 分割 : 1 + 1 + 1 + 1 + 1 + 1
5 分割 : 1 + 1 + 1 + 1 + 2
4 分割 : 1 + 1 + 1 + 3
1 + 1 + 2 + 2
3 分割 : 1 + 1 + 4
1 + 2 + 3
2 + 2 + 2
2 分割 : 1 + 5
2 + 4
3 + 3
1 分割 : 6
6 の場合、分割の仕方は 11 通りあります。この数を「分割数」といいます。自然数 n の分割数を求める関数 partition_number を定義してください。
val partition_number : int -> int = <fun>
# partition_number 1;; - : int = 1 # partition_number 2;; - : int = 2 # partition_number 3;; - : int = 3 # partition_number 4;; - : int = 5 # partition_number 5;; - : int = 7 # partition_number 6;; - : int = 11 # partition_number 7;; - : int = 15 # partition_number 8;; - : int = 22 # partition_number 10;; - : int = 42 # partition_number 50;; - : int = 204226
整数 n の分割の仕方をすべて求める高階関数 partition_of_integer fn n を定義してください。
val partition_of_integer : (int list -> unit) -> int -> unit = <fun>
# partition_of_integer print_intlist 5;; 5 4 1 3 2 3 1 1 2 2 1 2 1 1 1 1 1 1 1 1 - : unit = () # partition_of_integer print_intlist 6;; 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 - : unit = ()
リストで表した集合 ls を分割することを考えます。たとえば、集合 (1 2 3) は次のように分割することができます。
1 分割 : ((1 2 3)) 2 分割 : ((1 2) (3)), ((1 3) (2)), ((1) (2 3)) 3 分割 ; ((1) (2) (3))
このように、分割した集合 xs は元の集合 ls の部分集合になります。分割した部分集合の積は空集合になり、分割した部分集合のすべての和を求めると元の集合になります。
ls の分割の仕方をすべて求める高階関数 parititon_of_set fn ls を定義してください。
val partition_of_set : ('a list list -> unit) -> 'a list -> unit = <fun>
# partition_of_set print_intlist2 [1;2;3];; (1 2 3) (1)(2 3) (1 2)(3) (1 3)(2) (1)(2)(3) - : unit = () # partition_of_set print_intlist2 [1;2;3;4];; (1 2 3 4) (1)(2 3 4) (1 2)(3 4) (1 3 4)(2) (1)(2)(3 4) (1 2 3)(4) (1 4)(2 3) (1)(2 3)(4) (1 2 4)(3) (1 3)(2 4) (1)(2 4)(3) (1 2)(3)(4) (1 3)(2)(4) (1 4)(2)(3) (1)(2)(3)(4) - : unit = ()
リストの先頭から順番に添字と要素を関数 fn に渡して畳み込み (縮約) を行う関数 fold_left_with_index fn a ls を定義してください。
val fold_left_with_index : ('a -> int -> 'b -> 'b) -> 'b -> 'a list -> 'b = <fun>
# fold_left_with_index (fun x i a -> (i, x)::a) [] [1; 2; 3; 4; 5];; - : (int * int) list = [(4, 5); (3, 4); (2, 3); (1, 2); (0, 1)]
リストの末尾から順番に添字と要素を関数 fn に渡して畳み込み (縮約) を行う関数 fold_right_with_index fn a ls を定義してください。
val fold_right_with_index : ('a -> int -> 'b -> 'b) -> 'b -> 'a list -> 'b = <fun>
# fold_right_with_index (fun x i a -> (i, x)::a) [] [1; 2; 3; 4; 5];; - : (int * int) list = [(0, 1); (1, 2); (2, 3); (3, 4); (4, 5)]
添字付きのマップ関数 map_with_index fn ls を定義してください。
val map_with_index : ('a -> int -> 'b) -> 'a list -> 'b list = <fun>
# map_with_index (fun x i -> (i, x)) [1; 2; 3; 4; 5];; - : (int * int) list = [(0, 1); (1, 2); (2, 3); (3, 4); (4, 5)]
集合を分割する方法の総数を「ベル数 (Bell Number) 」といい、次の漸化式で求めることができます。
ベル数を求める関数 bell_number n を定義してください。
val bell_number : int -> Num.num = <fun>
# bell_number 0;; - : Num.num = Int 1 # bell_number 1;; - : Num.num = Int 1 # bell_number 2;; - : Num.num = Int 2 # bell_number 3;; - : Num.num = Int 5 # bell_number 4;; - : Num.num = Int 15 # bell_number 5;; - : Num.num = Int 52 # bell_number 10;; - : Num.num = Int 115975 # string_of_num (bell_number 20);; - : string = "51724158235372" # string_of_num (bell_number 30);; - : string = "846749014511809332450147" # string_of_num (bell_number 40);; - : string = "157450588391204931289324344702531067" # string_of_num (bell_number 50);; - : string = "185724268771078270438257767181908917499221852770"
k 個の要素をもつ集合 ls を要素数が等しい m 個の部分集合に分割することを考えます。部分集合の要素数 n は k / m になります。分割の仕方をすべて求める高階関数 group_partition fn n m ls を定義してください。
val group_partition : ('a list list -> unit) -> int -> int -> 'a list -> unit = <fun>
# group_partition print_intlist2 2 2 [1; 2; 3; 4];; (1 2)(3 4) (1 4)(2 3) (1 3)(2 4) - : unit = () # group_partition print_intlist2 2 3 [1; 2; 3; 4; 5; 6];; (1 2)(3 4)(5 6) (1 4)(2 3)(5 6) (1 3)(2 4)(5 6) (1 2)(3 6)(4 5) (1 6)(2 3)(4 5) (1 3)(2 6)(4 5) (1 2)(3 5)(4 6) (1 5)(2 3)(4 6) (1 3)(2 5)(4 6) (1 6)(2 5)(3 4) (1 5)(2 6)(3 4) (1 6)(2 4)(3 5) (1 4)(2 6)(3 5) (1 5)(2 4)(3 6) (1 4)(2 5)(3 6) - : unit = ()
集合を group_partition で分割するとき、その仕方の総数を求める関数 group_partition_number n m を定義してください。引数 n は部分集合の要素数、m は部分集合の個数です。
val group_partition_number : int -> int -> Num.num = <fun>
# group_partition_number 2 2;; - : Num.num = Int 3 # group_partition_number 2 3;; - : Num.num = Int 15 # group_partition_number 3 3;; - : Num.num = Int 280 # group_partition_number 3 4;; - : Num.num = Int 15400 # group_partition_number 3 5;; - : Num.num = Int 1401400
15 人の女生徒が毎日 3 人ずつ 5 組に分かれて散歩をするとき、1 週間 (7 日) のうちに、どの女生徒も他のすべての女生徒と 1 回ずつ同じ組になるような組み合わせを作ってください。
出典 : 大村平 (著), 『数理パズルの話』, 日科技連出版社, 1998
「カークマンの 15 人の女生徒」を解くプログラムを作ってください。
リスト : 素因数分解
let factorization n =
let rec factor_sub n m c =
if (mod_num n m) <> (Int 0) then (c, n)
else factor_sub (quo_num n m) m (succ_num c)
in
let rec iter i n a =
if n =/ (Int 1) then List.rev a
else if n </ i */ i then List.rev ((n, Int 1)::a)
else
let (c, m) = factor_sub n i (Int 0) in
if c =/ (Int 0) then iter (i +/ (Int 2)) n a
else iter (i +/ (Int 2)) m ((i, c)::a)
in
let (c, m) = factor_sub n (Int 2) (Int 0) in
if c >/ (Int 0) then iter (Int 3) m [(Int 2, c)]
else iter (Int 3) n []
素因数分解は素数 2, 3, 5, ... で順番に割り算していけばいいのですが、いちいち素数を求めるのは大変なので、2 と 3 以上の奇数列で割り算していきます。局所関数 factor_sub は n を m で割り算します。このとき、m で割り切れる回数を求めます。factor_sub は m で割った回数と商をタプルに格納して返します。
次に、factor_sub を呼び出して n を 2 で割り算します。それから、局所関数 iter で奇数列を生成します。変数 i は 3 で初期化します。a は結果を格納するリストです。n が 1 になる、または √n < i になったら繰り返しを終了します。そうでなければ、factor_sub を呼び出して n を i で割り算します。奇数列には素数ではないものがありますが、その前に小さな素数で素因数分解されているので、n がその値で割り切れることはありません。
n の素因数分解ができると、約数の個数を求めるのは簡単です。\(n = p^a \times q^b \times r^c\) とすると、約数の個数は \((a + 1) \times (b + 1) \times (c + 1)\) になります。たとえば、12 は \(2^2 \times 3^1\) になるので、約数の個数は 3 * 2 = 6 になります。実際、12 の約数は 1, 2, 3, 4, 6, 12 の 6 個です。
プログラムは次のようになります。
リスト : 約数の個数
let divisor_num n =
List.fold_left (fun a (_, n) -> (a */ (succ_num n)))
(Int 1)
(factorization n)
divisor_num は fold_left を使って n + 1 を a に掛け算していくだけです。
n の素因数分解ができると、約数の合計値を求めるのは簡単です。n の素因数分解が \(p^a\) だった場合、その約数の合計値は次の式で求めることができます。
たとえば、8 の素因数分解は \(2^3\) になり、素数の合計値は 8 + 4 + 2 + 1 = 15 になります。
\(p^a\) の約数の合計値を \(\sigma(p, a)\) で表すことにします。\(n = p^a \times q^b \times r^c\) の場合、\(n\) の約数の合計値は \(\sigma(p, a) \times \sigma(q, b) \times \sigma(r, c)\) になります。たとえば、12 は \(2^2 \times 3\) に素因数分解できますが、その合計値は (4 + 2 + 1) * (3 + 1) = 28 となります。12 の約数は 1, 2, 3, 4, 6, 12 なので、その合計値は確かに 28 になります。
プログラムは次のようになります。
リスト : 約数の合計値
let divisor_sum n =
let rec iter p n a =
if n =/ (Int 0) then (succ_num a)
else iter p (pred_num n) (a +/ (power_num p n))
in
List.fold_left (fun a (p, n) -> a */ (iter p n (Int 0)))
(Int 1)
(factorization n)
局所関数 iter は σ(p, n) を計算します。あとは fold_left で iter の返り値を累積変数 a に掛け算していくだけです。
p が素数の場合、\(p^a\) の約数は次のように簡単に求めることができます。
n の素因数分解が \(p^a \times q^b\) だったとすると、その約数は次のようになります。
たとえば、12 の約数は 24 = (1, 2, 4) と 3 = (1, 3) から、(1, 2, 4) * 1 と (1, 2, 4) * 3 のすべての要素 (1, 2, 4, 3, 6, 12) になります。
プログラムは次のようになります。
リスト : 約数をすべて求める
let divisor n =
let rec divisor_sub p n a =
if n =/ (Int 0) then (Int 1)::a
else divisor_sub p (pred_num n) ((power_num p n)::a)
in
let rec list_product p q a =
match p with
[] -> a
| x::xs -> list_product xs q ((List.map (fun y -> y */ x) q) @ a)
in
let (p, n)::xs = factorization n in
List.sort (fun x y -> compare_num x y)
(List.fold_left (fun a (p, n) -> list_product (divisor_sub p n []) a [])
(divisor_sub p n [])
xs)
局所関数 divisor_sub は pn の約数をリストに格納して返します。局所関数 list_product は 2 つのリスト p, q の要素を掛け合わせたものをリストに格納して返します。あとは fold_left で素因数分解した結果を順番に取り出し、(p . n) を divisor_sub でリストに変換して list_product で累積変数 a のリストと掛け合わせていくだけです。
リスト : 完全数
let perfect_number n =
let rec iter x =
if x <= n then
let (Int y) = divisor_sum (Int x) in
if y - x = x then (print_int x; print_newline ()) else ();
iter (x + 1)
else ()
in
iter 2
完全数を求める perfect_number は簡単です。x の約数の合計値を divisor_sum で求め、その値から x を引いた値が x と等しければ完全数です。print_int で x を表示します。
リスト : 友愛数
let print_yuuai_number n (Int m) =
print_int n;
print_char ' ';
print_int m;
print_newline()
let yuuai_number n =
let rec iter x =
if x <= n then
let m = (divisor_sum (Int x)) -/ (Int x) in
if (Int x) </ m && (divisor_sum m) -/ m =/ (Int x) then
print_yuuai_number x m
else ();
iter (x + 1)
else ();
in
iter 2
友愛数を求める yuuai_number も簡単です。divisor_sum で x の約数の合計値を求め、その値から x を引いた値を変数 m にセットします。m の約数の合計値から m を引いた値が x と等しければ、x と m は友愛数です。print_yuuai_number で x と m を表示します。同じ組を表示しないようにするため、x < m を条件に入れています。
整数 n を k 以下で分割する総数を求める関数を p(n, k) で表します。参考文献『Ｃ言語による最新アルゴリズム事典』によると、p(n, k) は次の式で表すことができるそうです。
r = 1 の場合は簡単ですね。n 個の 1 を選ぶ方法しかありません。同様に n = 1 の場合も、1 を選ぶ方法しかありません。なお、n = 0 の場合は 1 とします。
p(n, k) の場合、まず 1 を選ぶとすると、残りの n - 1 から 1 で分割する方法は p(n - 1, 1) 通りになります。2 を選ぶとすると、残りの n - 2 から 2 以下で分割する方法は p(n - 2, 2) 通りになります。つまり、1 から k までを選んだあとの分割数を計算し、その総和を求めればいいわけです。
簡単な例を示しましょう。次の図を見てください。
p(6, 6) = p(5, 1)
+ p(4, 2) => p(3, 1) + p(2, 2)
=> p(1, 1) + p(0, 2)
+ p(3, 3) => p(2, 1) + p(1, 2) + p(0, 3)
+ p(2, 4) => p(1, 1) + p(0, 2)
+ p(1, 5)
+ p(0, 6)
= 11 通り
p(6, 6) は p(5, 1) + p(4, 2) + p(3, 3) + p(2, 4) + p(1, 5) + p(0, 6) の総和になります。このうち、p(5, 1), p(1, 5), p(0, 6) は 1 になります。p(3, 3) は p(2, 1) + p(1, 2) + p(0, 3) になるので 3 通り、p(2, 4) は p(1, 1) + p(0, 2) になるので、2 通りになります。p(4, 2) はちょっと複雑です。p(4, 2) = p(3, 1) + p(2, 2) になります。ここで、p(2, 2) を求めると p(2, 2) = p(1, 1) + p(0, 2) になるので 2 通りになります。したがって、合計は 11 通りになります。
これをプログラムすると次のようになります。
リスト : 分割数
let partition_number n =
let rec part_num n k =
match (n, k) with
(0, _) -> 1
| (1, _) -> 1
| (_, 1) -> 1
| (_, _) -> if n < 0 || k < 1 then 0
else (part_num (n - k) k) + (part_num n (k - 1))
in
part_num n n
実際の処理は局所関数 part_num で行います。最初の 3 つの節で、分割数が 1 になる場合を定義しています。次の節で n が負または k が 1 未満の場合は 0 を返します。そうでなければ、p(n - 1, 1) + ... + p(n - k, k) を計算します。プログラムでは k の値をひとつずつ減らし、繰り返しを再帰定義で実現しています。なお、このプログラムはナイーブな実装なため、実行速度はとても遅いです。ご注意くださいませ。
リスト : 整数の分割
(* Q55 *)
let make_list x n =
let rec iter n a =
if n = 0 then a
else iter (n - 1) (x::a)
in
iter n []
let print_intlist xs =
List.iter (fun x -> print_int x; print_string " ") xs;
print_newline ()
let partition_of_integer f n =
let rec part_int n k a =
match (n, k) with
(0, _) -> f (List.rev a)
| (1, _) -> f (List.rev (1::a))
| (_, 1) -> f (List.rev ((make_list 1 n) @ a))
| (_, _) -> if n - k >= 0 then part_int (n - k) k (k::a) else ();
part_int n (k - 1) a
in
part_int n n []
基本的な考え方は partition_number と同じです。局所関数 part_int に累積変数 a を追加して、選んだ数値を a に格納していくだけです。n が 0 の場合は f を評価し、n が 1 の場合は a に 1 を追加してから f を評価します。k が 1 の場合は make_list で要素が 1 で長さが n のリストを作成します。そして、それを演算子 @ で a と連結してから f を評価します。
集合を分割するアルゴリズムは簡単です。たとえば、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)) になります。
このように、簡単な方法で集合を分割することができます。実際にプログラムを作る場合、上図を木と考えて、深さ優先で木をたどると簡単です。次のリストを見てください。
リスト : 集合の分割
(* Q58 *)
let rec remove x = function
[] -> []
| y::ys -> if x = y then remove x ys
else y :: remove x ys
let print_intlist2 ls =
let rec print_intlist = function
[] -> ()
| [x] -> print_int x
| x::xs -> print_int x; print_char ' '; print_intlist xs
in
List.iter (fun xs -> print_char '('; print_intlist xs; print_char ')') ls;
print_newline ()
let partition_of_set f ls =
let rec part_set ls a =
match ls with
[] -> f a
| x::xs -> List.iter (fun y -> part_set xs ((x::y)::(remove y a))) a;
part_set xs ([x]::a)
in
part_set (List.rev ls) []
実際の処理は局所関数 part_set で行います。最初の節が再帰呼び出しの停止条件です。次の節の List.iter で、部分集合に要素 x を追加します。匿名関数でリスト a から要素 y を順番に取り出し、a から y を取り除いたリストに x::y を追加します。最後に、新しい部分集合 [x] を a に追加します。ただし、このままでは要素の並び方が逆順になるので、part_set を呼び出す前に List.rev でリスト ls を反転しています。これで集合の分割をすべて求めることができます。
リスト : 添字付き畳み込み (1)
let fold_left_with_index f a ls =
let rec iter i a = function
[] -> a
| x::xs -> iter (i + 1) (f x i a) xs
in
iter 0 a ls
fold_left_with_index は fold_left に添字の処理を追加しただけです。関数 f を呼び出すとき、第 1 引数がリストの要素、第 2 引数が添字、第 3 引数に累積変数が渡されることに注意してください。
リスト : 添字付き畳み込み (2)
let fold_right_with_index f a ls =
let rec iter i a = function
[] -> a
| x::xs -> f x i (iter (i + 1) a xs)
in
iter 0 a ls
fold_right_with_index は fold_right に添字の処理を追加しただけです。関数 f を呼び出すとき、第 1 引数がリストの要素、第 2 引数が添字、第 3 引数に累積変数が渡されることに注意してください。
リスト : 添字付きマップ関数
let map_with_index f ls =
let rec iter i = function
[] -> []
| x::xs -> (f x i)::(iter (i + 1) xs)
in
iter 0 ls
map_with_index も簡単です。実際の処理は局所関数 iter で行います。引数 i が添字を表します。関数 f の第 2 引数に添字を渡すことに注意してください。
リスト : ベル数
let rec comb_num n r =
if n = r || r = 0 then (Int 1)
else (comb_num n (r - 1)) */ (Int (n - r + 1)) // (Int r)
let bell_number n =
let rec iter i = function
x::_ when i = n -> x
| xs -> iter (i + 1) ((fold_left_with_index
(fun x k a -> (comb_num i k) */ x +/ a)
(Int 0)
xs)::xs)
in
iter 0 [(Int 1)]
bell_number は公式をそのままプログラムするだけです。実際の処理は局所関数 iter で行います。第 2 引数は累積変数で、ベル数を逆順で格納します。nＣk は関数 comb_num で求めます。nＣk * B(k) の総和は関数 fold_left_with_index で計算します。累積変数は逆順になっていますが、二項係数は nＣi と nＣn - i の値が同じになるので、そのまま計算しても大丈夫です。
リスト : 集合のグループ分け
let group_partition f n m ls =
let rec group_part ls a =
match ls with
[] -> f a
| x::xs -> List.iter (fun y -> if List.length y < n then
group_part xs ((x::y)::(remove y a))
else ())
a;
if List.length a < m then group_part xs ([x]::a) else ()
in
group_part (List.rev ls) []
group_partition は partition_of_set を改造するだけで簡単に作成することができます。生成する部分集合の大きさを n に、部分集合の個数を m に制限するだけです。部分集合 y に x を追加する場合、length y < n であることをチェックします。新しい部分集合 [x] を追加する場合、length a < m であることをチェックします。これで集合をグループに分けることができます。
グループ分けの総数は次の式で求めることができます。
たとえば、n = 3, m = 5 の場合は次のようになります。
これをそのままプログラムすると次のようになります。
リスト : グループ分けの総数
let rec fact = function
0 -> Int 1
| n when n > 0 -> (Int n) */ (fact (n - 1))
let group_partition_number n m =
let rec group_part_num k a =
if k = 0 then a // (fact m)
else group_part_num (k - n) (a */ (comb_num k n))
in
group_part_num (n * m) (Int 1)
階乗は関数 fact で、組み合わせの個数は関数 comb_num で計算します。要素の個数を変数 k にセットし、累積変数 a に comb_num k n を乗算します。あとは k から n を減算し、k が 0 でなければ処理を繰り返すだけです。最後に a // (fact m) を計算して返します。
「カークマンの 15 人の女生徒」の解法プログラムは group_partition を改造することで簡単に作成することができます。次のリストを見てください。
リスト : カークマンの 15 人の女生徒
exception Kirkman_exit
let check_table : int list array = Array.make 16 []
let rec check_person ls x =
match ls with
[] -> true
| y::ys -> if List.mem x check_table.(y) then false
else check_person ys x
let add_person ls x =
List.iter (fun y -> check_table.(x) <- y::(check_table.(x));
check_table.(y) <- x::(check_table.(y)))
ls
let del_person ls x =
List.iter (fun y -> check_table.(x) <- List.tl check_table.(x);
check_table.(y) <- List.tl check_table.(y))
ls
let iota n m =
let rec iter i a =
if i < n then a
else iter (i - 1) (i::a)
in
iter m []
let kirkman () =
let rec kirkman_sub ls a b =
match ls with
[] -> if List.length b = 6 then
begin
List.iter (fun x -> print_intlist2 x) (List.rev (a::b));
raise Kirkman_exit
end
else
kirkman_sub (iota 2 15) [[1]] (a::b)
| x::xs -> List.iter (fun y -> if List.length y < 3 && check_person y x then
begin
add_person y x;
kirkman_sub xs ((x::y)::(remove y a)) b;
del_person y x
end
else ())
a;
if List.length a < 5 then kirkman_sub xs ([x]::a) b else ()
in
try kirkman_sub (iota 2 15) [[1]] [] with
_ -> ()
15 人の女生徒を 1 から 15 までの数値で表します。変数 check_table は、いっしょに散歩した人を格納する配列です。0 番目はダミーです。たとえば、[1; 2; 3] というグループを作った場合、check_table の 1 番目には [2; 3] を、2 番目には [1; 3] を、3 番目には [2; 3] をセットします。この check_table を使って、同じ女生徒と 2 回以上散歩しないようにグループ分けを行います。
関数 check_person はグループ ls に x を追加するとき、既に散歩した女生徒がいるかチェックします。check_table の y 番目からリストを取り出し、それに x が含まれていれば、y は既に x と散歩をしています。この場合は false を返します。x が ls の女生徒達とまだ散歩していない場合は true を返します。
関数 add_person は check_table にグループ ls と x の関係を追加します。ls の要素を y とすると、check_table の x 番目のリストに y を、y 番目のリストに x を追加するだけです。関数 del_person は ls と x の関係を削除します。ls の要素を y とすると、check_table の x 番目の先頭要素と、y 番目の先頭要素を削除します。
解法プログラム kirkman の実際の処理は局所関数 kirkman_sub で行います。引数 ls が女生徒を格納したリスト、a が作成中のグループ分けを格納するリスト、b が完成したグループ分けを格納するリストです。b の長さが 7 になれば解を見つけたことになります。
プログラムでは ls が空リストになり (a がひとつ完成する)、b の長さが 6 の場合、完成した a を b に追加し、それを List.rev で反転して要素を print_intlist2 で表示します。そうでなければ、a を b に追加して kirkman_sub を再帰呼び出しして次の日のグループ分けを作成します。グループ分けの処理は group_partition とほぼ同じですが、check_person でチェックを行い、add_person で check_table を更新してから、kirkman_sub を再帰呼び出しします。再帰呼び出しから戻ってきたら、del_person で check_table を元に戻します。
それでは実行結果を示します。
# let a = Sys.time () in kirkman (); print_float (Sys.time () -. a);; (15 14 13)(12 11 10)(9 8 7)(6 5 4)(3 2 1) (15 4 3)(14 10 9)(13 11 8)(12 5 2)(7 6 1) (15 12 7)(14 11 1)(13 10 6)(9 4 2)(8 5 3) (15 11 2)(14 7 5)(13 9 3)(12 8 6)(10 4 1) (15 9 6)(14 12 3)(13 5 1)(11 7 4)(10 8 2) (15 10 5)(14 8 4)(13 7 2)(12 9 1)(11 6 3) (15 8 1)(14 6 2)(13 12 4)(11 9 5)(10 7 3) 188.937- : unit = ()
実行時間は 3 分 9 秒 (Windows XP, celeron 1.40 GHz, OCaml ver 3.10.0) でした。けっこう時間がかかりますね。ocamlopt でネイティブコードにコンパイルすると 39.2 秒になりました。興味のある方は高速化に挑戦してみてください。