自然数 n を素因数分解する関数 factorization n を定義してください。返り値はリスト [(p, q); ...] で、(p, q) は \(p^q\) を表します。なお、整数 n, p はモジュール IntInf で表すものとします。
val factorization = fn : IntInf.int -> (IntInf.int * int) list
- factorization(12345678); val it = [(2,1),(3,2),(47,1),(14593,1)] : (IntInf.int * int) list - factorization(123456789); val it = [(3,2),(3607,1),(3803,1)] : (IntInf.int * int) list - factorization(1234567890); val it = [(2,1),(3,2),(5,1),(3607,1),(3803,1)] : (IntInf.int * int) list - factorization(1111111111); val it = [(11,1),(41,1),(271,1),(9091,1)] : (IntInf.int * int) list
自然数 n の約数の個数を求める関数 divisor_num を定義してください。
val divisor_num = fn : IntInf.int -> int
- divisor_num(12345678); val it = 24 : int - divisor_num(123456789); val it = 12 : int - divisor_num(1234567890); val it = 48 : int - divisor_num(1111111111); val it = 16 : int
自然数 n の約数の合計値を求める関数 divisor_sum を定義してください。
val divisor_sum = fn : IntInf.int -> IntInf.int
- divisor_sum(12345678); val it = 27319968 : IntInf.int - divisor_sum(123456789); val it = 178422816 : IntInf.int - divisor_sum(1234567890); val it = 3211610688 : IntInf.int - divisor_sum(1111111111); val it = 1246404096 : IntInf.int
自然数 n の約数をリストに格納して返す関数 divisor を定義してください。
val divisor = fn : IntInf.int -> IntInf.int list
- divisor(12); val it = [1,2,3,4,6,12] : IntInf.int list - divisor(12345678); val it = [1,2,3,6,9,18,47,94,141,282,423,846,...] : IntInf.int list - divisor(123456789); val it = [1,3,9,3607,3803,10821,11409,32463,34227,13717421,41152263,123456789] : IntInf.int list - divisor(1234567890); val it = [1,2,3,5,6,9,10,15,18,30,45,90,...] : IntInf.int list - divisor(1111111111); val it = [1,11,41,271,451,2981,9091,11111,100001,122221,372731,2463661,...] : IntInf.int list
完全数 - Wikipedia によると、『完全数（かんぜんすう，perfect number）とは、その数自身を除く約数の和が、その数自身と等しい自然数のことである。』 とのことです。自然数 n 以下の完全数を求める関数 perfect_number を定義してください。
val perfect_number = fn : IntInf.int -> unit
- perfect_number(10000); 6 28 496 8128 val it = () : unit
友愛数 - Wikipedia によると、『友愛数（ゆうあいすう）とは、異なる2つの自然数の組で、自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。』 とのことです。自然数 n 以下の友愛数を求める関数 yuuai_number を定義してください。
val yuuai_number = fn : IntInf.int -> unit
- 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 val it = () : unit
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
図 : 整数の分割
整数 n を 1 以上の自然数の和で表すことを考えます。これを「整数の分割」といいます。整数を分割するとき、同じ自然数を何回使ってもかまいませんが、並べる順序が違うだけのものは同じ分割とします。
たとえば 6 の場合、分割の仕方は 11 通りあります。この数を「分割数」といいます。自然数 n の分割数を求める関数 partition_number を定義してください。
val partition_number = fn : int -> int
- partition_number(1); val it = 1 : int - partition_number(2); val it = 2 : int - partition_number(3); val it = 3 : int - partition_number(4); val it = 5 : int - partition_number(5); val it = 7 : int - partition_number(6); val it = 11 : int - partition_number(7); val it = 15 : int - partition_number(8); val it = 22 : int - partition_number(9); val it = 30 : int - partition_number(10); val it = 42 : int - partition_number(20); val it = 627 : int - partition_number(50); val it = 204226 : int
整数 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 = fn : (''a list list -> unit) -> ''a list -> unit
- partition_of_set print_intlist2 [1,2,3]; ( 1 2 3 ) ( 1 )( 2 3 ) ( 1 2 )( 3 ) ( 1 3 )( 2 ) ( 1 )( 2 )( 3 ) val it = () : 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 ) val it = () : unit
リストの先頭から順番に添字と要素を関数 fn に渡して畳み込み (縮約) を行う関数 fold_left_with_index fn a ls を定義してください。
val fold_left_with_index = fn : ('a * int * 'b -> 'b) -> 'b -> 'a list -> 'b
- fold_left_with_index (fn(x,i,a) => (x, i)::a) [] [1,2,3,4,5]; val it = [(5,4),(4,3),(3,2),(2,1),(1,0)] : (int * int) list
リストの末尾から順番に添字と要素を関数 fn に渡して畳み込み (縮約) を行う関数 fold_right_with_index fn a ls を定義してください。
val fold_right_with_index = fn : ('a * int * 'b -> 'b) -> 'b -> 'a list -> 'b
- fold_right_with_index (fn(x,i,a) => (x, i)::a) [] [1,2,3,4,5]; val it = [(1,0),(2,1),(3,2),(4,3),(5,4)] : (int * int) list
添字付きのマップ関数 map_with_index fn ls を定義してください。
val map_with_index = fn : ('a * int -> 'b) -> 'a list -> 'b list
- map_with_index (fn(x, i) => (x, i)) [1,2,3,4,5]; val it = [(1,0),(2,1),(3,2),(4,3),(5,4)] : (int * int) list
集合を分割する方法の総数を「ベル数 (Bell Number)」といい、次の漸化式で求めることができます。
ベル数を求める関数 bell_number n を定義してください。
val bell_number = fn : IntInf.int -> IntInf.int
- bell_number(0); val it = 1 : IntInf.int - bell_number(1); val it = 1 : IntInf.int - bell_number(2); val it = 2 : IntInf.int - bell_number(3); val it = 5 : IntInf.int - bell_number(4); val it = 15 : IntInf.int - bell_number(5); val it = 52 : IntInf.int - bell_number(10); val it = 115975 : IntInf.int - bell_number(20); val it = 51724158235372 : IntInf.int - bell_number(30); val it = 846749014511809332450147 : IntInf.int - bell_number(40); val it = 157450588391204931289324344702531067 : IntInf.int - bell_number(50); val it = 185724268771078270438257767181908917499221852770 : IntInf.int
k 個の要素をもつ集合 ls を要素数が等しい m 個の部分集合に分割することを考えます。部分集合の要素数 n は k / m になります。分割の仕方をすべて求める高階関数 group_partition fn n m ls を定義してください。
val group_partition = fn : (''a list list -> unit) -> int -> int -> ''a list -> unit
- group_partition print_intlist2 2 2 [1,2,3,4]; ( 1 2 )( 3 4 ) ( 1 4 )( 2 3 ) ( 1 3 )( 2 4 ) val it = () : 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 ) val it = () : unit
集合を group_partition で分割するとき、その仕方の総数を求める関数 group_partition_number n m を定義してください。引数 n は部分集合の要素数、m は部分集合の個数です。
val group_partition_number = fn : IntInf.int * IntInf.int -> IntInf.int
- group_partition_number(2,2); val it = 3 : IntInf.int - group_partition_number(2,3); val it = 15 : IntInf.int - group_partition_number(3,3); val it = 280 : IntInf.int - group_partition_number(3,4); val it = 15400 : IntInf.int - group_partition_number(3,5); val it = 1401400 : IntInf.int
リスト : 素因数分解
fun factorization(n : IntInf.int) =
let
fun factor_sub(n, m, c) =
if n mod m <> 0 then (c, n)
else factor_sub(n div m, m, c + 1)
fun iter(i, n, a) =
if n = 1 then rev a
else if n < i * i then rev ((n, 1)::a)
else
let
val (c, m) = factor_sub(n, i, 0)
in
if c = 0 then iter(i + 2, n, a)
else iter(i + 2, m, (i, c)::a)
end
val (c, m) = factor_sub(n, 2, 0)
in
if c > 0 then iter(3, m, [(2, c)])
else iter(3, n, [])
end
素因数分解は素数 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 個です。
リスト : 約数の個数
fun divisor_num(n : IntInf.int) =
foldl (fn(x, a) => a * (#2(x) + 1)) 1 (factorization n)
divisor_num は foldl を使って 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 になります。
プログラムは次のようになります。
リスト : 約数の合計値
fun divisor_sum(n : IntInf.int) =
let
fun iter(p, n, a) =
if n = 0 then a + 1
else iter(p, n - 1, a + IntInf.pow(p, n))
in
foldl (fn(x, a) => a * iter(#1(x), #2(x), 0)) 1 (factorization n)
end
局所関数 iter は σ(p, n) を計算します。あとは foldl で 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) になります。
プログラムは次のようになります。
リスト : 約数をすべて求める
fun divisor(n : IntInf.int) =
let
fun divisor_sub(p, n, a) =
if n = 0 then 1::a
else divisor_sub(p, n - 1, IntInf.pow(p, n)::a)
fun list_product([], _, a) = a
| list_product(x::xs, q, a) =
list_product(xs, q, (List.map (fn y => y * x) q) @ a)
val ((p, n)::xs) = factorization(n)
in
ListMergeSort.sort
(op >)
(foldl (fn((p, n), a) => list_product(divisor_sub(p, n, []), a, []))
(divisor_sub(p, n, []))
xs)
end
局所関数 divisor_sub は pn の約数をリストに格納して返します。局所関数 list_product は 2 つのリスト p, q の要素を掛け合わせたものをリストに格納して返します。あとは foldl で素因数分解した結果を順番に取り出し、(p . n) を divisor_sub でリストに変換して list_product で累積変数 a のリストと掛け合わせていくだけです。
リスト : 完全数
fun perfect_number(n) =
let
fun iter(x) =
if x <= n then
let
val y = divisor_sum(x)
in
if y - x = x then print(IntInf.toString(x) ^ "\n") else ();
iter(x + 1)
end
else ()
in
iter(2)
end
完全数を求める perfect_number は簡単です。x の約数の合計値を divisor_sum で求め、その値から x を引いた値が x と等しければ完全数です。print で x を表示します。
リスト : 友愛数
fun yuuai_number(n) =
let
fun iter(x) =
if x <= n then
let
val m = divisor_sum(x) - x
in
if x < m andalso divisor_sum(m) - m = x
then print(IntInf.toString(x) ^ " " ^ IntInf.toString(m) ^ "\n")
else ();
iter(x + 1)
end
else ()
in
iter(2)
end
友愛数を求める yuuai_number も簡単です。divisor_sum で x の約数の合計値を求め、その値から x を引いた値を変数 m にセットします。m の約数の合計値から m を引いた値が x と等しければ、x と m は友愛数です。print で 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 通りになります。
これをプログラムすると次のようになります。
リスト : 分割数
fun partition_number(n) =
let
fun part_num(0, _) = 1
| part_num(1, _) = 1
| part_num(_, 1) = 1
| part_num(n, k) =
if n < 0 orelse k < 1 then 0
else part_num(n - k, k) + part_num(n, k - 1)
in
part_num(n, n)
end
実際の処理は局所関数 part_num で行います。最初の 3 つの節で、分割数が 1 になる場合を定義しています。次の節で n が負または k が 1 未満の場合は 0 を返します。そうでなければ、p(n - 1, 1) + ... + p(n - k, k) を計算します。プログラムでは k の値をひとつずつ減らし、繰り返しを再帰定義で実現しています。なお、このプログラムはナイーブな実装なため、実行速度はとても遅いです。ご注意くださいませ。
リスト : 整数の分割
(* 問題26 *)
fun make_list(x, n) =
let
fun iter(0, a) = a
| iter(n, a) = iter(n - 1, x::a)
in
iter(n, [])
end
(* 問題68 *)
fun print_intlist([]) = print("\n")
| print_intlist(x::xs) = (print(Int.toString(x) ^ " "); print_intlist(xs))
fun partition_of_integer f n =
let
fun part_int(0, _, a) = f (rev a)
| part_int(1, _, a) = f (rev (1::a))
| part_int(n, 1, a) = f (rev (make_list(1, n) @ a))
| part_int(n, k, 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, [])
end
基本的な考え方は 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)) になります。
このように、簡単な方法で集合を分割することができます。実際にプログラムを作る場合、上図を木と考えて、深さ優先で木をたどると簡単です。次のリストを見てください。
リスト : 集合の分割
(* 問題27 *)
fun remove(x, []) = []
| remove(x, y::ys) =
if x = y then remove(x, ys) else y :: remove(x, ys)
fun print_intlist2([]) = print("\n")
| print_intlist2(x::xs) = (
print("( ");
List.app (fn y => print(Int.toString(y) ^ " ")) x;
print(")");
print_intlist2(xs)
)
(* 問題69 *)
fun partition_of_set f ls =
let
fun part_set([], a) = f a
| part_set(x::xs, a) = (
List.app (fn y => part_set(xs, (x::y) :: remove(y, a))) a;
part_set(xs, [x]::a)
)
in
part_set(rev ls, [])
end
実際の処理は局所関数 part_set で行います。最初の節が再帰呼び出しの停止条件です。次の節の List.app で、部分集合に要素 x を追加します。匿名関数でリスト a から要素 y を順番に取り出し、a から y を取り除いたリストに x::y を追加します。最後に、新しい部分集合 [x] を a に追加します。ただし、このままでは要素の並び方が逆順になるので、part_set を呼び出す前に rev でリスト ls を反転しています。これで集合の分割をすべて求めることができます。
リスト : 添字付き畳み込み (1)
fun fold_left_with_index f a ls =
let
fun iter(_, [], a) = a
| iter(i, x::xs, a) = iter(i + 1, xs, f(x, i, a))
in
iter(0, ls, a)
end
fold_left_with_index は foldl に添字の処理を追加しただけです。関数 f を呼び出すとき、第 1 引数がリストの要素、第 2 引数が添字、第 3 引数に累積変数が渡されることに注意してください。
リスト : 添字付き畳み込み (2)
fun fold_right_with_index f a ls =
let
fun iter(_, [], a) = a
| iter(i, x::xs, a) = f(x, i, iter(i + 1, xs, a))
in
iter(0, ls, a)
end
fold_right_with_index は foldr に添字の処理を追加しただけです。関数 f を呼び出すとき、第 1 引数がリストの要素、第 2 引数が添字、第 3 引数に累積変数が渡されることに注意してください。
リスト : 添字付きマップ関数
fun map_with_index f ls =
let
fun iter(i, []) = []
| iter(i, x::xs) = f(x, i) :: iter(i + 1, xs)
in
iter(0, ls)
end
map_with_index も簡単です。実際の処理は局所関数 iter で行います。引数 i が添字を表します。関数 f の第 2 引数に添字を渡すことに注意してください。
リスト : ベル数
(* 問題32 *)
fun comb_num(n, r) =
if n = r orelse r = 0 then 1 : IntInf.int
else comb_num(n, r - 1) * (n - r + 1) div r
(* 問題73 *)
fun bell_number(n) =
let
fun iter(i, ls) =
if i = n then hd ls
else iter(i + 1,
(fold_left_with_index
(fn(x, k, a) => comb_num(i, Int.toLarge(k)) * x + a)
0
ls)::ls)
in
iter(0, [1])
end
bell_number は公式をそのままプログラムするだけです。実際の処理は局所関数 iter で行います。第 2 引数は累積変数で、ベル数を逆順で格納します。nＣk は関数 comb_num で求めます。nＣk * B(k) の総和は関数 fold_left_with_index で計算します。累積変数は逆順になっていますが、二項係数は nＣi と nＣn - i の値が同じになるので、そのまま計算しても大丈夫です。
リスト : 集合のグループ分け
fun group_partition f n m ls =
let
fun group_part([], a) = f a
| group_part(x::xs, a) = (
List.app (fn y => if length(y) < n
then group_part(xs, (x::y)::remove(y, a))
else ())
a;
if length a < m then group_part(xs, [x]::a) else ()
)
in
group_part(rev ls, [])
end
group_partition は partition_of_set を改造するだけで簡単に作成することができます。生成する部分集合の大きさを n に、部分集合の個数を m に制限するだけです。部分集合 y に x を追加する場合、length y < n であることをチェックします。新しい部分集合 [x] を追加する場合、length a < m であることをチェックします。これで集合をグループに分けることができます。
グループ分けの総数は次の式で求めることができます。
たとえば、n = 3, m = 5 の場合は次のようになります。
これをそのままプログラムすると次のようになります。
リスト : グループ分けの総数
fun fact(0) = 1
| fact(n : IntInf.int) = n * fact(n - 1)
fun group_partition_number(n, m) =
let
fun group_part_num(k, a) =
if k = 0 then a div fact(m)
else group_part_num(k - n, a * comb_num(k, n))
in
group_part_num(n * m, 1)
end
階乗は関数 fact で、組み合わせの個数は関数 comb_num で計算します。要素の個数を変数 k にセットし、累積変数 a に comb_num k n を乗算します。あとは k から n を減算し、k が 0 でなければ処理を繰り返すだけです。最後に a div (fact m) を計算して返します。