自然数 n を素因数分解する関数 factorization n を定義してください。返り値はリスト [(p, q), ...] で、(p, q) は \(p^q\) を表します。
factorization :: Integer -> [(Integer, Integer)]
ghci> factorization 12345678 [(2,1),(3,2),(47,1),(14593,1)] ghci> factorization 123456789 [(3,2),(3607,1),(3803,1)] ghci> factorization 1234567890 [(2,1),(3,2),(5,1),(3607,1),(3803,1)] ghci> factorization 1111111111 [(11,1),(41,1),(271,1),(9091,1)] ghci> factorization 11111111111 [(21649,1),(513239,1)]
自然数 n の約数の個数を求める関数 divisor_num を定義してください。
divisor_num :: Integer -> Integer
ghci> divisor_num 12345678 24 ghci> divisor_num 123456789 12 ghci> divisor_num 1234567890 48 ghci> divisor_num 1111111111 16 ghci> divisor_num 11111111111 4
自然数 n の約数の合計値を求める関数 divisor_sum を定義してください。
divisor_sum :: Integer -> Integer
ghci> divisor_sum 12345678 27319968 ghci> divisor_sum 123456789 178422816 ghci> divisor_sum 1234567890 3211610688 ghci> divisor_sum 1111111111 1246404096 ghci> divisor_sum 11111111111 11111646000
自然数 n の約数をリストに格納して返す関数 divisor を定義してください。
divisor :: Integer -> [Integer]
ghci> divisor 12 [1,2,3,4,6,12] ghci> divisor 12345678 [1,2,3,6,9,18,47,94,141,282,423,846,14593,29186,43779,87558,131337,262674,685871, 1371742,2057613,4115226,6172839,12345678] ghci> divisor 123456789 [1,3,9,3607,3803,10821,11409,32463,34227,13717421,41152263,123456789] ghci> divisor 1234567890 [1,2,3,5,6,9,10,15,18,30,45,90,3607,3803,7214,7606,10821,11409,18035,19015,21642, 22818,32463,34227,36070,38030,54105,57045,64926,68454,108210,114090,162315,171135, 324630,342270,13717421,27434842,41152263,68587105,82304526,123456789,137174210, 205761315,246913578,411522630,617283945,1234567890] ghci> divisor 1111111111 [1,11,41,271,451,2981,9091,11111,100001,122221,372731,2463661,4100041,27100271, 101010101,1111111111] ghci> divisor 11111111111 [1,21649,513239,11111111111]
完全数 - Wikipedia によると、『完全数（かんぜんすう，perfect number）とは、その数自身を除く約数の和が、その数自身と等しい自然数のことである。』 とのことです。自然数 n 以下の完全数を求める関数 perfect_number を定義してください。
perfect_number :: Integer -> [Integer]
ghci> perfect_number 10000 [6,28,496,8128]
友愛数 - Wikipedia によると、『友愛数（ゆうあいすう）とは、異なる2つの自然数の組で、自分自身を除いた約数の和が、互いに他方と等しくなるような数をいう。』 とのことです。自然数 n 以下の友愛数を求める関数 yuuai_number を定義してください。
val yuuai_number = fn : IntInf.int -> unit
ghci> :t yuuai_number yuuai_number :: Integer -> [(Integer, Integer)] ghci> 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)]
整数 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 を定義してください。
partition_number :: Integer -> Integer
ghci> partition_number 1 1 ghci> partition_number 2 2 ghci> partition_number 3 3 ghci> partition_number 4 5 ghci> partition_number 5 7 ghci> partition_number 10 42 ghci> partition_number 20 627 ghci> partition_number 50 204226
整数 n の分割の仕方をすべて求める関数 partition_of_integer n を定義してください。
partition_of_integer :: Int -> [[Int]]
ghci> partition_of_integer 5 [[5],[4,1],[3,2],[3,1,1],[2,2,1],[2,1,1,1],[1,1,1,1,1]] ghci> partition_of_integer 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]]
リストで表した集合 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 ls を定義してください。
partition_of_set :: Eq a => [a] -> [[[a]]]
ghci> mapM_ print $ partition_of_set [1,2,3] [[1],[2],[3]] [[1,2],[3]] [[1,3],[2]] [[1],[2,3]] [[1,2,3]] ghci> mapM_ print $ partition_of_set [1,2,3,4] [[1],[2],[3],[4]] [[1,2],[3],[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,4],[3]] [[1,3],[2,4]] [[1],[2],[3,4]] [[1,2],[3,4]] [[1,3,4],[2]] [[1],[2,3,4]] [[1,2,3,4]]
集合を分割する方法の総数を「ベル数 (Bell Number) 」といい、次の漸化式で求めることができます。
ベル数を求める関数 bell_number n を定義してください。
bell_number :: Integer -> Integer
ghci> bell_number 0 1 ghci> bell_number 1 1 ghci> bell_number 2 2 ghci> bell_number 3 5 ghci> bell_number 4 15 ghci> bell_number 5 52 ghci> bell_number 10 115975 ghci> bell_number 20 51724158235372 ghci> bell_number 50 185724268771078270438257767181908917499221852770
k 個の要素をもつ集合 ls を要素数が等しい m 個の部分集合に分割することを考えます。部分集合の要素数 n は k / m になります。分割の仕方をすべて求める関数 group_partition n m ls を定義してください。
group_partition :: Eq a => Int -> Int -> [a] -> [[[a]]]
ghci> group_partition 2 2 [1,2,3,4] [[[1,4],[2,3]],[[1,3],[2,4]],[[1,2],[3,4]]] ghci> mapM_ print $ group_partition 2 3 [1..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]] [[1,6],[2,3],[4,5]] [[1,3],[2,6],[4,5]] [[1,2],[3,6],[4,5]] [[1,5],[2,3],[4,6]] [[1,3],[2,5],[4,6]] [[1,2],[3,5],[4,6]] [[1,4],[2,3],[5,6]] [[1,3],[2,4],[5,6]] [[1,2],[3,4],[5,6]] val it = () : unit
集合を group_partition で分割するとき、その仕方の総数を求める関数 group_partition_number n m を定義してください。引数 n は部分集合の要素数、m は部分集合の個数です。
group_partition_number :: Integer -> Integer -> Integer
ghci> group_partition_number 2 2 3 ghci> group_partition_number 2 3 15 ghci> group_partition_number 3 3 280 ghci> group_partition_number 3 4 15400 ghci> group_partition_number 3 5 1401400
m 個の整数 1, 2, ..., m の順列を考えます。先頭の要素を 1 から数えることとすると、i 番目の要素が整数 i ではない順列を「完全順列」といいます。1 から m までの整数値で完全順列を生成する関数 derangement m を定義してください。
derangement :: Int -> [[Int]]
ghci> derangement 3 [[2,3,1],[3,1,2]] ghci> mapM_ print $ derangement 4 [2,1,4,3] [2,3,4,1] [2,4,1,3] [3,1,4,2] [3,4,1,2] [3,4,2,1] [4,1,2,3] [4,3,1,2] [4,3,2,1]
完全順列の総数を「モンモール数 (Montmort number) 」といいます。モンモール数は次の漸化式で求めることができます。
モンモール数を求める関数 montmort_number を定義してください。
montmort_number :: Integer -> Integer
ghci> montmort_number 1 0 ghci> montmort_number 2 1 ghci> montmort_number 3 2 ghci> montmort_number 4 9 ghci> montmort_number 5 44 ghci> montmort_number 10 1334961 ghci> montmort_number 20 895014631192902121
「ラテン方陣」は数独の枠の条件を無くした方陣です。ラテン方陣の定義を参考文献『数理パズルのはなし』より引用します。 『ラテン方陣を一般的にいうなら、n 行 n 列の正方形の枡に n 種類の記号を n 個ずつ配列して、各行各列に記号の重複のないものを n 次のラテン方陣というのです。』
1 2 3 1 2 3 1 3 2 1 3 2 2 1 3 2 1 3
2 3 1 3 1 2 2 1 3 3 2 1 1 3 2 3 2 1
3 1 2 2 3 1 3 2 1 2 1 3 3 2 1 1 3 2
標準形
2 3 1 2 3 1 3 1 2 3 1 2 3 2 1 3 2 1
1 2 3 3 1 2 1 2 3 2 3 1 1 3 2 2 1 3
3 1 2 1 2 3 2 3 1 1 2 3 2 1 3 1 3 2
図 : 3 次のラテン方陣
このラテン方陣をパズルに応用したものが数独というわけです。3 次のラテン方陣は上に示す 12 通りになります。この中で、最初の行と列の要素を昇順に並べたものを「標準形」といいます。
3 次のラテン方陣の場合、標準形は 1 種類しかありません。ラテン方陣は任意の行を交換する、または任意の列を交換してもラテン方陣になります。3 次のラテン方陣の場合、標準形から行または列を交換することで、残りの 11 種類のラテン方陣を生成することができます。
4 次の標準形ラテン方陣をすべて求めてください。
リスト : 素因数分解
factorization :: Integer -> [(Integer, Integer)]
factorization n =
if c > 0 then iter 3 m [(2, c)] else iter 3 n []
where
(c, m) = factor n 2 0
factor n m c =
if n `mod` m /= 0 then (c, n)
else factor (n `div` m) m (c + 1)
iter i n a
| n == 1 = reverse a
| n < i * i = reverse ((n, 1) : a)
| otherwise =
let (c, m) = factor n i 0
in if c == 0 then iter (i + 2) n a
else iter (i + 2) m ((i, c):a)
素因数分解は素数 2, 3, 5, ... で順番に割り算していけばいいのですが、いちいち素数を求めるのは大変なので、2 と 3 以上の奇数列で割り算していきます。局所関数 factor は n を m で割り算し、m で割り切れる回数を求めます。factor は m で割った回数と商をタプルに格納して返します。
最初に、factor を呼び出して n を 2 で割り算します。それから、局所関数 iter で奇数列を生成します。累積変数 a は結果を格納するリストです。n が 1 になる、または √n < i になったら繰り返しを終了します。そうでなければ、factor を呼び出して 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 個です。
プログラムは次のようになります。
リスト : 約数の個数 divisor_num :: Integer -> Integer divisor_num n = foldl (\a (_, x) -> a * (x + 1)) 1 (factorization n)
divisor_num は foldl を使って x + 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 になります。
プログラムは次のようになります。
リスト : 約数の合計値 sigma :: Integer -> Integer -> [Integer] sigma p n = map (p^) [0 .. n] divisor_sum :: Integer -> Integer divisor_sum n = foldl (\a (p, x) -> a * sum (sigma p x)) 1 (factorization n)
関数 sigma は [1, p, p2, ... pa-1, pa] を求めます。あとは foldl で sum (sigma p x) の返り値を累積変数 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) になります。
プログラムは次のようになります。
リスト : 約数をすべて求める
divisor :: Integer -> [Integer]
divisor n =
let ((p1, x1):xs) = factorization n
product xs ys = [x * y | x <- xs, y <- ys]
in sort $ foldl (\a (p, x) -> product (sigma p x) a) (sigma p1 x1) xs
局所関数 product は 2 つのリスト xs, ys の要素を掛け合わせたものをリストに格納して返します。あとは foldl で素因数分解した結果を順番に取り出し、(p, x) を sigma でリストに変換して product で累積変数 a のリストと掛け合わせていくだけです。
リスト : 完全数 perfect_number :: Integer -> [Integer] perfect_number n = filter (\x -> let y = divisor_sum x in y - x == x) [2 .. n]
完全数を求める perfect_number は簡単です。x の約数の合計値を divisor_sum で求め、その値から x を引いた値が x と等しければ完全数です。
リスト : 友愛数
yuuai_number :: Integer -> [(Integer, Integer)]
yuuai_number n =
foldr (\x a -> let m = divisor_sum x - x
in if x < m && divisor_sum m - m == x
then (x, m):a else a)
[]
[2 .. n]
友愛数を求める yuuai_number も簡単です。divisor_sum で x の約数の合計値を求め、その値から x を引いた値を変数 m にセットします。m の約数の合計値から m を引いた値が x と等しければ、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 通りになります。
分割数を求める関数 partition_number は、関数 p(n, k) を使うと次のようにプログラムすることができます。
リスト : 分割数
partition_number :: Integer -> Integer
partition_number n = part_num n n
where
part_num 0 _ = 1
part_num 1 _ = 1
part_num _ 1 = 1
part_num n k
| n < 0 || k < 1 = 0
| otherwise = part_num (n - k) k + part_num n (k - 1)
関数 part_num は p(n, k) の定義をそのままプログラムしただけです。ただし、このプログラムは二重再帰で何度も同じ値を求めているため実行速度はとても遅くなります。
拙作のページ Algorithms with Python「動的計画法」では、「動的計画法」を使って分割数を高速に求めています。同ページより引用します。
動的計画法を使うと、大きな値でも高速に計算することができます。次の図を見てください。
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 のベクタを用意します。ベクタの添字が 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) の値を求めることができます。
Haskell の場合、リスト (遅延ストリーム) を使って同様の事を行うことができます。次のリストを見てください。
リスト : 別解 (動的計画法)
partition_number' :: Int -> Integer
partition_number' n = iter 2 (replicate (n + 1) 1)
where
iter k xs
| k > n = last xs
| otherwise = let ys = (take k xs ++ zipWith (+) ys (drop k xs))
in iter (k + 1) ys
実際の処理は局所関数 iter で行います。リストの値を更新するとき、新しいリストを ys とすると、先頭の k 個の要素は xs と同じで、それ以降の要素は xs の i 番目の要素と ys の i - k 番目の要素を加算します。この処理は、ys と (drop k xs) を zipWith で加算することと同じです。ys の先頭要素は take k xs で求められているので、zipWith で ys の要素を先頭から順番に取り出していくことができます。
リスト : 整数の分割
partition_of_integer :: Int -> [[Int]]
partition_of_integer n = part_int n n [] []
where
part_int 0 _ xs ys = (reverse xs):ys
part_int 1 _ xs ys = (reverse (1:xs)):ys
part_int n 1 xs ys = (reverse (replicate n 1 ++ xs)):ys
part_int n k xs ys =
let ys' = part_int n (k - 1) xs ys
in if n - k >= 0 then part_int (n - k) k (k:xs) ys' else ys'
基本的な考え方は partition_number と同じです。局所関数 part_int に累積変数 xs, ys を追加して、選んだ数値を xs に格納していきます。n が 0 の場合は reverse xs を ys に格納し、n が 1 の場合は xs に 1 を追加してから ys に格納します。k が 1 の場合はモジュール Data.List の関数 replicate で要素が 1 で長さが n のリストを作成します。そして、それを演算子 ++ で xs と連結してから ys に格納します。
集合を分割するアルゴリズムは簡単です。たとえば、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)) になります。
このように、簡単な方法で集合を分割することができます。実際にプログラムを作る場合、上図を木と考えて、深さ優先で木をたどると簡単です。次のリストを見てください。
リスト : 集合の分割
partition_of_set :: Eq a => [a] -> [[[a]]]
partition_of_set xs = part_set (reverse xs) [] []
where
part_set [] ys zs = ys:zs
part_set (x:xs) ys zs =
part_set xs ([x]:ys) $ foldr (\y a -> part_set xs ((x:y):delete y ys) a) zs ys
実際の処理は局所関数 part_set で行います。累積変数 ys に分割中のリストをセットし、完成したリストを累積変数 zs に格納します。最初の節が再帰呼び出しの停止条件です。次の節の foldr で、部分集合に要素 x を追加します。ラムダ式でリスト ys から要素 y を順番に取り出し、ys から y を取り除いたリストに x : y を追加します。それから、新しい部分集合 [x] を ys に追加します。ただし、このままでは要素の並び方が逆順になるので、part_set を呼び出す前に reverse でリスト xs を反転しています。これで集合の分割をすべて求めることができます。
ところで、拙作のページ「順列と組み合わせ」で作成した関数 select を使うと型クラス制約 Eq a を省くことができます。
リスト : 別解
select :: [a] -> [(a, [a])]
select [] = []
select [x] = [(x, [])]
select (x:xs) = (x, xs) : map (\(y, ys) -> (y, x:ys)) (select xs)
partition_of_set' :: [a] -> [[[a]]]
partition_of_set' xs = part_set (reverse xs) [] []
where
part_set [] ys zs = ys:zs
part_set (x:xs) ys zs =
part_set xs ([x]:ys) $ foldr (\(y, ys') a -> part_set xs ((x:y):ys') a) zs (select ys)
select ys で部分集合 y を一つ選びます。残りの部分集合はリスト ys' に格納されているので、(x:y) を ys' に追加すればいいわけです。、
リスト : ベル数
comb_num :: Integer -> Integer -> Integer
comb_num n r
| n == r || r == 0 = 1
| otherwise = comb_num n (r - 1) * (n - r + 1) `div` r
bell_number :: Integer -> Integer
bell_number n = iter 0 [1]
where
iter i xs
| n == i = head xs
| otherwise = iter (i + 1) $ foldl (\a (x, k) -> (comb_num i k) * x + a) 0 (zip xs [0..]) : xs
bell_number は公式をそのままプログラムするだけです。実際の処理は局所関数 iter で行います。第 2 引数は累積変数で、ベル数を逆順で格納します。nＣk は関数 comb_num で求めます。nＣk * B(k) の総和は関数 foldl で計算します。累積変数は逆順になっていますが、二項係数は nＣi と nＣn - i の値が同じになるので、そのまま計算しても大丈夫です。
リスト : 集合のグループ分け
group_partition :: Eq a => Int -> Int -> [a] -> [[[a]]]
group_partition n m xs = group_part (reverse xs) [] []
where
group_part [] ys zs = ys:zs
group_part (x:xs) ys zs =
let zs' = foldr (\y a -> if length y < n
then group_part xs ((x:y):delete y ys) a
else a)
zs
ys
in if length ys < m then group_part xs ([x]:ys) zs' else zs'
group_partition は partition_of_set を改造するだけで簡単に作成することができます。生成する部分集合の大きさを n に、部分集合の個数を m に制限するだけです。部分集合 y に x を追加する場合、length y < n であることをチェックします。新しい部分集合 [x] を追加する場合、length a < m であることをチェックします。これで集合をグループに分けることができます。
グループ分けの総数は次の式で求めることができます。
たとえば、n = 3, m = 5 の場合は次のようになります。
これをそのままプログラムすると次のようになります。
リスト : グループ分けの総数
fact :: Integer -> Integer
fact 0 = 1
fact n = n * fact (n - 1)
group_partition_number :: Integer -> Integer -> Integer
group_partition_number n m = group_part_num (n * m) 1
where
group_part_num k a
| k == 0 = a `div` fact m
| otherwise = group_part_num (k - n) (a * comb_num k n)
階乗は関数 fact で、組み合わせの個数は関数 comb_num で計算します。要素の個数を変数 k にセットし、累積変数 a に comb_num k n を乗算します。あとは k から n を減算し、k が 0 でなければ処理を繰り返すだけです。最後に a `div` (fact m) を計算して返します。
リスト : 完全順列
derangement :: Int -> [[Int]]
derangement m = perm_sub 1 [1 .. m] [] []
where
perm_sub _ [] ys zs = reverse ys : zs
perm_sub n xs ys zs =
foldr (\x a -> if x /= n
then perm_sub (n + 1) (delete x xs) (x:ys) a
else a)
zs
xs
derangement は簡単です。実際の処理は局所関数 perm_sub で行います。引数 n が順番を表します。引数 xs と ys は累積変数として使います。foldr のラムダ式の中で、数字 x が n と等しくない場合、その数字を選択することできます。等しい場合は選択しません。xs が空リストになったら、reverse で ys を反転して zs に格納します。これで完全順列を生成することができます。
ところで、拙作のページ「順列と組み合わせ」で作成した関数 select を使うと、もっと簡単にプログラムを作ることができます。
リスト : 別解
select :: [a] -> [(a, [a])]
select [] = []
select [x] = [(x, [])]
select (x:xs) = (x, xs) : map (\(y, ys) -> (y, x:ys)) (select xs)
derangement' :: Int -> [[Int]]
derangement' m = perm_sub 1 [1 .. m]
where
perm_sub _ [] = [[]]
perm_sub n xs = [y:zs | (y, ys) <- select xs, y /= n, zs <- perm_sub (n + 1) ys]
局所関数 perm_sub において、select で選んだ数字 y が n と等しくないことを確認するだけです。
リスト : 完全順列の総数
montmort_number :: Integer -> Integer
montmort_number 1 = 0
montmort_number 2 = 1
montmort_number n = (n - 1) * (montmort_number (n - 1) + montmort_number (n - 2))
-- 別解
montmort_number' :: Integer -> Integer
montmort_number' n = iter 1 0 1
where
iter i a b
| n == i = a
| otherwise = iter (i + 1) b ((i + 1) * (a + b))
関数 montmort_number は公式をそのままプログラムしただけです。二重再帰になっているので、実行速度はとても遅くなります。これを繰り返しに変換すると別解のようになります。考え方はフィボナッチ数列と同じです。累積変数 a に i 番目の値を、b に i + 1 番目の値を保存しておきます。すると、i + 2 番目の値は (i + 1) * (a + b) で計算することができます。あとは、b の値を a に、新しい値を b にセットして処理を繰り返すだけです。
リスト : 標準形ラテン方陣を求める
check_latina :: Int -> Int -> [[Int]] -> Bool
check_latina n x xs =
elem x $ map (\ys -> ys !! (n - 1)) xs
latina :: Int -> [[[Int]]]
latina size = solver 1 [1 .. size] [] [[1 .. size]] []
where
solver n [] a b c =
if size - 1 == length b then (reverse (reverse a : b) : c)
else let m = length b + 2
in solver 2 (delete m [1 .. size]) [m] (reverse a : b) c
solver n xs a b c =
foldr (\x z -> if not(check_latina n x b)
then solver (n + 1) (delete x xs) (x:a) b z else z)
c
xs
実際の処理は局所関数 solver で行います。基本的な考え方は完全順列とほぼ同じで、累積変数 a に順列を格納し、完成した順列を累積変数 b に、完成したラテン方陣を累積変数 c に格納します。引数 xs が空リストの場合、順列がひとつ完成しました。b の要素数をチェックして、size - 1 と等しければラテン方陣ができました。(reverse a) を b に追加してから、それを c に追加します。そうでなければ solver を再帰呼び出しします。このとき、先頭要素は b の要素数 + 2 になることに注意してください。
順列を生成する場合、関数 check_latina を呼び出して数字 x を選択できるかチェックします。map で xs に格納されたリストの n - 1 番目の要素を演算子 !! で取り出します。Haskell のリストは 0 から数えることに注意してください。そして、x が map の返り値に含まれているか関数 elem でチェックします。x が含まれていれば、x を選択することはできません。そうでなければ x を選択します。
実行結果は次のようになります。
ghci> mapM_ (mapM_ print) $ latina 3 [1,2,3] [2,3,1] [3,1,2] ghci> mapM_ (mapM_ print) $ latina 4 [1,2,3,4] [2,1,4,3] [3,4,1,2] [4,3,2,1] [1,2,3,4] [2,1,4,3] [3,4,2,1] [4,3,1,2] [1,2,3,4] [2,3,4,1] [3,4,1,2] [4,1,2,3] [1,2,3,4] [2,4,1,3] [3,1,4,2] [4,3,2,1]
ちなみに、標準形ラテン方陣の総数は次のようになります。
I4 = 4 I5 = 56 I6 = 9408 I7 = 16942080
高次の標準形ラテン方陣の総数は、簡単に求めることができない非常にハードな問題だといわれています。