Yet Another Erlang Problems

●解答66

リスト : 素因数分解

factor_sub(N, M, C) when N rem M =/= 0 -> {N, C};
factor_sub(N, M, C) -> factor_sub(N div M, M, C + 1).

factor_loop(I, N, A) when N < I * I ->
  if
    N > 1 -> lists:reverse([{N, 1} | A]);
    true -> lists:reverse(A)
  end;
factor_loop(I, N, A) ->
  {M, C} = factor_sub(N, I, 0),
  I2 = I + 2,
  if
    C > 0 -> factor_loop(I2, M, [{I, C} | A]);
    true -> factor_loop(I2, N, A)
  end.

factorization(N) ->
  {M, C} = factor_sub(N, 2, 0),
  if
    C > 0 -> factor_loop(3, M, [{2, C}]);
    true  -> factor_loop(3, N, [])
  end.

素因数分解は素数 2, 3, 5, ... で順番に割り算していけばいいのですが、いちいち素数を求めるのは大変なので、2 と 3 以上の奇数列で割り算していきます。関数 factor_sub は N を M で割り算します。このとき、M で割り切れる回数を求めます。factor_sub は M で割った回数と商をタプルに格納して返します。

関数 factor_loop は奇数列を生成します。変数 I は 3 で初期化します。A は結果を格納するリストです。N が 1 になる、または √N < I になったら繰り返しを終了します。そうでなければ、factor_sub を呼び出して N を I で割り算します。奇数列には素数ではないものがありますが、その前に小さな素数で素因数分解されているので、N がその値で割り切れることはありません。あとは、2 で割り算してから factor_loop を呼び出すだけです。

●解答67

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(N) ->
  lists:foldl(fun ({_, X}, A) -> A * (X + 1) end, 1, factorization(N)).

divisor_num は lists:foldl を使って X + 1 を a に掛け算していくだけです。

●解答68

整数 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 になります。

プログラムは次のようになります。

リスト : 約数の合計値

% 整数の累乗
pow(_, 0) -> 1;
pow(X, Y) when Y > 0, Y rem 2 =:= 0 ->
    Z = pow(X, Y div 2),
    Z * Z;
pow(X, Y) when Y > 0 ->
    Z = pow(X, Y div 2),
    X * Z * Z.

prime_sum(_, 0, A) -> A + 1;
prime_sum(P, N, A) -> prime_sum(P, N - 1, A + pow(P, N)).

divisor_sum(N) ->
  lists:foldl(fun({P, M}, A) -> prime_sum(P, M, 0) * A end, 1, factorization(N)).

関数 prime_sub は σ(p, n) を計算します。あとは lists:foldl で prime_sub の返り値を累積変数 a に掛け算していくだけです。

●解答69

p が素数の場合、\(p^a\) の約数は次のように簡単に求めることができます。

n の素因数分解が \(p^a \times q^b\) だったとすると、その約数は次のようになります。

たとえば、12 の約数は 2⁴ = (1, 2, 4) と 3 = (1, 3) から、(1, 2, 4) * 1 と (1, 2, 4) * 3 のすべての要素 (1, 2, 4, 3, 6, 12) になります。

プログラムは次のようになります。

リスト : 約数をすべて求める

% P ^ N の約数を求める
divisor_sub(_, 0, A) -> [1 | A];
divisor_sub(P, N, A) when N > 0 -> divisor_sub(P, N - 1, [pow(P, N) | A]).

divisor_product(Xs, Ys) -> [X * Y || X <- Xs, Y <- Ys].

divisor_sub1([], A) -> A;
divisor_sub1([{P, N} | Xs], A) ->
  divisor_sub1(Xs, divisor_product(A, divisor_sub(P, N, []))).

divisor(N) -> 
  insert_sort(fun(X, Y) -> X < Y end, divisor_sub1(factorization(N), [1])).

関数 divisor_sub は P^N の約数をリストに格納して返します。関数 divisor_product は 2 つのリスト Xs, Ys の要素を掛け合わせた結果をリストに格納して返します。あとは divisor_sub1 で素因数分解した結果を順番に取り出し、{P, N} を divisor_sub でリストに変換して、それを divisor_product で累積変数 a のリストと掛け合わせていくだけです。

●解答70

リスト : 完全数

perfect_number(N, X) when N < X -> ok;
perfect_number(N, X) ->
  M = divisor_sum(X),
  if
    M - X =:= X -> io:write(X), io:nl();
    true -> false
  end,
  perfect_number(N, X + 1).
perfect_number(N) -> perfect_number(N, 2).

完全数を求める perfect_number は簡単です。X の約数の合計値 M を divisor_sub で求め、その値から X を引いた値が X と等しければ完全数です。io:write で X を表示します。

●解答71

リスト : 友愛数

yuuai_number(N, X) when N < X -> ok;
yuuai_number(N, X) ->
  M = divisor_sum(X) - X,
  if
    X < M -> M1 = divisor_sum(M),
             if
               X =:= M1 - M -> io:write({X, M}), io:nl();
               true -> false
             end;
    true -> false
  end,
  yuuai_number(N, X + 1).
yuuai_number(N) -> yuuai_number(N, 2).

友愛数を求める yuuai_number も簡単です。divisor_sum で X の約数の合計値を求め、その値から X を引いた値を変数 M にセットします。M の約数の合計値 M1 から M を引いた値が X と等しければ、X と M は友愛数です。io:write で {X, M} を表示します。同じ組を表示しないようにするため、X < M を条件に入れています。

●解答72

─┬─ 6                           : 6
  │
  ├─ 5 ─ 1                      : 5 + 1
  │
  ├─ 4 ┬ 2                      : 4 + 2
  │     │
  │     └ 1 ─ 1                 : 4 + 1 + 1
  │
  ├─ 3 ┬ 3                      : 3 + 3
  │     │
  │     ├ 2 ─ 1                 : 3 + 2 + 1
  │     │
  │     └ 1 ─ 1 ─ 1            : 3 + 1 + 1 + 1
  │
  ├─ 2 ┬ 2 ┬ 2                 : 2 + 2 + 2
  │     │   │
  │     │   └ 1 ─ 1            : 2 + 2 + 1 + 1
  │     │
  │     └ 1 ─ 1 ─ 1 ─ 1       : 2 + 1 + 1 + 1 + 1
  │
  └─ 1 ─ 1 ─ 1 ─ 1 ─ 1 ─ 1  : 1 + 1 + 1 + 1 + 1 + 1

                    図 : 整数 6 の分割

6 の場合、分割の仕方は上図のように 11 通りあります。分割の仕方を列挙する場合、整数 n から k 以下の整数を選んでいくと考えてください。まず、6 から 6 を選びます。すると、残りは 0 になるので、これ以上整数を分割することはできません。次に、6 から 5 を選びます。残りは 1 になるので、1 を選ぶしか方法はありません。

次に、4 を選びます。残りは 2 になるので、2 から 2 以下の整数を分割する方法になります。2 から 2 を選ぶと残りは 0 になるので 2 が得られます。1 を選ぶと残りは 1 になるので、1 + 1 が得られます。したがって、4 + 2, 4 + 1 + 1 となります。同様に、6 から 3 を選ぶと、残りは 3 から 3 以下の整数を選ぶ方法になります。

6 から 2 以下の整数を選ぶ方法は、残り 4 から 2 以下の整数を選ぶ方法になり、そこで 2 を選ぶと 2 から 2 以下の整数を選ぶ方法になります。1 を選ぶと 4 から 1 以下の整数を選ぶ方法になりますが、これは 1 通りしかありません。最後に 6 から 1 を選びますが、これも 1 通りしかありません。これらをすべて足し合わせると 11 通りになります。

整数 n を k 以下の整数で分割する総数を求める関数を p(n, k) とすると、p(n, k) は次のように定義することができます。

簡単な例を示しましょう。次の図を見てください。

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

分割数を求める関数 partitionNumber は、関数 p(n, k) を使うと次のようにプログラムすることができます。

リスト : 分割数

part_num(0, _) -> 1;
part_num(1, _) -> 1;
part_num(_, 1) -> 1;
part_num(N, K) when N < 0; K < 1 -> 0;
part_num(N, K) ->
  part_num(N - K, K) + part_num(N, K - 1).

part_num(N) -> part_num(N, N).

実際の処理は述語 part_num で行います。最初の 3 つの規則で、分割数が 1 になる場合を定義しています。次の規則は N が負または K が 1 未満の場合を表します。この場合、分割数は 0 になります。最後の規則で、p(n - 1, 1) + ... + p(n - k, k) を計算します。プログラムでは K の値をひとつずつ減らしていることに注意してください。なお、このプログラムはナイーブな実装なため、実行速度はとても遅いです。ご注意くださいませ。

●解答73

リスト : 整数の分割

part_int(F, 0, _, A) -> F(lists:reverse(A));
part_int(F, 1, _, A) -> F(lists:reverse([1 | A]));
part_int(F, N, 1, A) -> F(lists:reverse(yaep04:make_list(1, N) ++ A));
part_int(F, N, K, A) ->
  if 
    N - K >= 0 -> part_int(F, N - K, K, [K | A]);
    true -> false
  end,
  part_int(F, N, K - 1, A).

partition_of_integer(F, N) -> part_int(F, 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 - K が 0 以上であれば part_int を再帰呼び出しします。そのあとで、K の値を -1 して part_int を再帰呼び出しします。

●解答74

リスト : 完全順列

derangement(F, _, [], A) -> F(lists:reverse(A));
derangement(F, N, Xs, A) ->
  lists:foreach(
    fun(X) -> 
      if
        N =/= X -> derangement(F, N + 1, lists:delete(X, Xs), [X | A]);
        true -> false
      end
    end,
    Xs).
derangement(F, N) -> derangement(F, 1, iota(1 , N), []).

derangemet は簡単です。実際の処理は関数 derangement で行います。iota で 1 から N までの数値を格納したリストを生成し、それを引数 Xs に渡します。引数 N が順番を表します。lists:foreach の無名関数で、数字 X が N と等しくない場合、その数字を選択することできます。等しい場合は選択しません。Xs が空リストになったら、reverse で A を反転して F を評価します。これで完全順列を生成することができます。

●解答75

リスト : 完全順列の総数

montmort_number(1) -> 0;
montmort_number(2) -> 1;
montmort_number(N) when N > 2 ->
  (N - 1) * (montmort_number(N - 1) + montmort_number(N - 2)).

% 別解
montmort_number1(N, N, A, _) -> A;
montmort_number1(I, N, A, B) -> montmort_number1(I + 1, N, B, (I + 1) * (A + B)).
mont_mort_number1(N) when N >= 1 -> montmort_number1(1, N, 0, 1).

関数 montmort_number は公式をそのままプログラムしただけです。二重再帰になっているので、実行速度は大変遅くなります。これを繰り返しに変換すると別解のようになります。考え方はフィボナッチ数列と同じです。累積変数 A に i 番目の値を、B に i + 1 番目の値を保存しておきます。すると、i + 2 番目の値は (I + 1) * (A + B) で計算することができます。あとは、B の値を A に、新しい値を B にセットして処理を繰り返すだけです。

●解答76

集合を分割するアルゴリズムは簡単です。たとえば、n -1 個の要素 x₁, ..., x_n-1 を持つ集合を分割したところ、i 個の部分集合 S₁, ..., S_i が生成されたとしましょう。ここに、n 番目の要素 x_n を追加すると、要素が n 個の集合を分割することができます。

新しい要素を追加する場合は次に示す手順で行います。

1. 部分集合 S_k (k = 1 から i まで) に要素 x_n を追加する
2. 新しい部分集合 S_i+1 (要素が x_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]] になります。

このように、簡単な方法で集合を分割することができます。実際にプログラムを作る場合、上図を木と考えて、深さ優先で木をたどると簡単です。次のリストを見てください。

リスト : 集合の分割

part_set(F, [], A) -> F(A);
part_set(F, [X | Xs], A) ->
  lists:foreach(fun(Y) -> part_set(F, Xs, [[X | Y] | lists:delete(Y, A)]) end, A),
  part_set(F, Xs, [[X] | A]).
partition_of_set(F, Xs) -> part_set(F, lists:reverse(Xs), []).

実際の処理は関数 part_set で行います。最初の規則が再帰呼び出しの停止条件です。次の規則の lists:foreach で、部分集合に要素 X を追加します。無名関数 でリスト A から要素 Y を順番に取り出し、A から Y を取り除いたリストに [X | Y] を追加します。最後に、新しい部分集合 [X] を A に追加します。ただし、このままでは要素の並び方が逆順になるので、part_set を呼び出す前に reverse でリスト Xs を反転しています。これで集合の分割をすべて求めることができます。

●解答77

リスト : ベル数

fold_with_index(_, _, A, []) -> A;
fold_with_index(F, I, A, [X | Xs]) ->
  fold_with_index(F, I + 1, F(X, I, A), Xs).

bell_number(N, N, [X | _]) -> X;
bell_number(I, N, Bs) ->
  bell_number(I + 1,
              N,
              [fold_with_index(
                 fun(X, K, A) -> A + comb_num(I, K) * X end,
                 0,
                 0,
                 Bs) | Bs]).
bell_number(N) -> bell_number(0, N, [1]).

bell_number は公式をそのままプログラムするだけです。累積変数 Bs にベル数を逆順で格納します。_nＣ_k は関数 comb_num で求めます。_nＣ_k * B(k) の総和は関数 fold_with_index で計算します。fold_with_index は添字を関数に渡して畳み込みを行います。無名関数の引数 X がリストの要素、K が添字、A が累積変数です。Bs は逆順になっていますが、二項係数は _nＣ_i と _nＣ_{n - i} の値が同じになるので、そのまま計算しても大丈夫です。もちろん、reverse で Bs を逆順にしてから計算してもかまいません。

●解答78

リスト : 集合のグループ分け

group_part(F, _, _, [], A) -> F(A);
group_part(F, N, M, [X | Xs], A) ->
  lists:foreach(fun(Y) ->
      if
        length(Y) < N -> group_part(F, N, M, Xs, [[X | Y] | lists:delete(Y, A)]);
        true -> false
      end
    end,
    A),
  if
    length(A) < M -> group_part(F, N, M, Xs, [[X] | A]);
    true -> ok
  end.
group_partition(F, N, M, Xs) -> group_part(F, N, M, lists:reverse(Xs), []).

group_partition は partition_of_set を改造するだけで簡単に作成することができます。生成する部分集合の大きさを N に、部分集合の個数を M に制限するだけです。部分集合 Y に X を追加する場合、length(Y) < N であることをチェックします。新しい部分集合 [X] を追加する場合、length(A) < M であることをチェックします。これで集合をグループに分けることができます。

●解答79

グループ分けの総数は次の式で求めることができます。

たとえば、n = 3, m = 5 の場合は次のようになります。

これをそのままプログラムすると次のようになります。

リスト : グループ分けの総数

fact(0) -> 1;
fact(N) when N > 0 -> N * fact(N - 1).

group_part_num(0, _, M, A) -> A div fact(M);
group_part_num(K, N, M, A) -> group_part_num(K - N, N, M, A * comb_num(K, N)).

group_partition_number(N, M) -> group_part_num(N * M, N, M, 1).

階乗は関数 fact で、組み合わせの個数は関数 comb_num で計算します。要素の個数を変数 K にセットし、累積変数 a に comb_num(K, N) を乗算します。あとは K から N を減算し、K が 0 でなければ処理を繰り返すだけです。最後に A div fact(M) を計算して返します。

●解答80

参考文献『Ｃ言語による最新アルゴリズム事典』の「式の評価」によると、四則演算の数式は次の構文規則で表すことができます。

式 := 項 (+ | -) 項 (+ | -) 項 ...
項 :- 因子 (* | /) 因子 (* | /) 因子 ...
因子 := 数 | (式)

これをそのままプログラムすると、次のようになります。

リスト : 数式の計算 (中置記法)

% 因子の評価
factor([X | Xs]) when is_number(X) -> {X, Xs};
factor([X | Xs]) when is_list(X) -> {expression(X), Xs}.

% 項の評価
term_sub(V, [X | Xs]) when X =:= '*' ->
  {Y, Ys} = factor(Xs), term_sub(V * Y, Ys);
term_sub(V, [X | Xs]) when X =:= '/' ->
  {Y, Ys} = factor(Xs), term_sub(V / Y, Ys);
term_sub(V, Xs) -> {V, Xs}.

term(Xs) -> {Y, Ys} = factor(Xs), term_sub(Y, Ys).

% 式の評価
expr_sub(V, []) -> V;
expr_sub(V, [X | Xs]) when X =:= '+' ->
  {Y, Ys} = term(Xs), expr_sub(V + Y, Ys);
expr_sub(V, [X | Xs]) when X =:= '-' ->
  {Y, Ys} = term(Xs), expr_sub(V - Y, Ys).

expression(Xs) -> {Y, Ys} = term(Xs), expr_sub(Y, Ys).

関数 expression は「式」を評価します。最初に関数 term を呼び出して「項」を評価します。次に関数 expr_sub を呼び出します。最初の節で、リストが空であれば評価結果 V を返します。演算子が + または - の場合は term を呼び出して評価を行い、返り値を Y と Ys にセットします。そして、V + Y または V - Y を計算します。

関数 term も同様の処理を行います。この場合は最初に関数 factor を呼び出して「因子」を評価します。そして、演算子が * または / の場合は factor を呼び出して評価を続行します。そうでなければ、評価結果 V と残りのリスト Xs をタプルに格納して返します。関数 factor は簡単で、リストの先頭要素 X が数値の場合はそれをそのまま返し、リストであればそれを expression に渡して評価します。

数値の判定は組み込み関数 is_number で、リストの判定は組み込み関数 is_list で行うことができます。それ以外の場合はマッチングしないのでエラーになります。

●プログラムリスト

%
% yaep04.erl : Yet Another Erlang Problems (4)
%
%              Copyright (C) 2011-2024 Makoto Hiroi
%
-module(yaep04).
-export([factorization/1, divisor_num/1, divisor_sum/1, divisor/1, perfect_number/1]).
-export([yuuai_number/1, part_num/1, partition_of_integer/2, derangement/2, montmort_number/1]).
-export([montmort_number1/1, partition_of_set/2, bell_number/1, group_partition/4]).
-export([group_partition_number/2, expression/1]).
-import(yaep01, [iota/2]).
-import(yaep02, [make_list/2, comb_num/2]).
-import(yaep03, [insert_sort/2]).

factor_sub(N, M, C) when N rem M =/= 0 -> {N, C};
factor_sub(N, M, C) -> factor_sub(N div M, M, C + 1).

factor_loop(I, N, A) when N < I * I ->
  if
    N > 1 -> lists:reverse([{N, 1} | A]);
    true -> lists:reverse(A)
  end;
factor_loop(I, N, A) ->
  {M, C} = factor_sub(N, I, 0),
  I2 = I + 2,
  if
    C > 0 -> factor_loop(I2, M, [{I, C} | A]);
    true -> factor_loop(I2, N, A)
  end.

factorization(N) ->
  {M, C} = factor_sub(N, 2, 0),
  if
    C > 0 -> factor_loop(3, M, [{2, C}]);
    true  -> factor_loop(3, N, [])
  end.

divisor_num(N) ->
  lists:foldl(fun ({_, X}, A) -> A * (X + 1) end, 1, factorization(N)).

% 整数の累乗
pow(_, 0) -> 1;
pow(X, Y) when Y > 0, Y rem 2 =:= 0 ->
    Z = pow(X, Y div 2),
    Z * Z;
pow(X, Y) when Y > 0 ->
    Z = pow(X, Y div 2),
    X * Z * Z.

prime_sum(_, 0, A) -> A + 1;
prime_sum(P, N, A) -> prime_sum(P, N - 1, A + pow(P, N)).

divisor_sum(N) ->
  lists:foldl(fun({P, M}, A) -> prime_sum(P, M, 0) * A end, 1, factorization(N)).

% P ^ N の約数を求める
divisor_sub(_, 0, A) -> [1 | A];
divisor_sub(P, N, A) when N > 0 -> divisor_sub(P, N - 1, [pow(P, N) | A]).

divisor_product(Xs, Ys) -> [X * Y || X <- Xs, Y <- Ys].

divisor_sub1([], A) -> A;
divisor_sub1([{P, N} | Xs], A) ->
  divisor_sub1(Xs, divisor_product(A, divisor_sub(P, N, []))).

divisor(N) ->
  insert_sort(fun(X, Y) -> X < Y end, divisor_sub1(factorization(N), [1])).

perfect_number(N, X) when N < X -> ok;
perfect_number(N, X) ->
  M = divisor_sum(X),
  if
    M - X =:= X -> io:write(X), io:nl();
    true -> false
  end,
  perfect_number(N, X + 1).
perfect_number(N) -> perfect_number(N, 2).

yuuai_number(N, X) when N < X -> ok;
yuuai_number(N, X) ->
  M = divisor_sum(X) - X,
  if
    X < M -> M1 = divisor_sum(M),
             if
               X =:= M1 - M -> io:write({X, M}), io:nl();
               true -> false
             end;
    true -> false
  end,
  yuuai_number(N, X + 1).
yuuai_number(N) -> yuuai_number(N, 2).

part_num(0, _) -> 1;
part_num(1, _) -> 1;
part_num(_, 1) -> 1;
part_num(N, K) when N < 0; K < 1 -> 0;
part_num(N, K) ->
  part_num(N - K, K) + part_num(N, K - 1).

part_num(N) -> part_num(N, N).

part_int(F, 0, _, A) -> F(lists:reverse(A));
part_int(F, 1, _, A) -> F(lists:reverse([1 | A]));
part_int(F, N, 1, A) -> F(lists:reverse(make_list(1, N) ++ A));
part_int(F, N, K, A) ->
  if
    N - K >= 0 -> part_int(F, N - K, K, [K | A]);
    true -> false
  end,
  part_int(F, N, K - 1, A).

partition_of_integer(F, N) -> part_int(F, N, N, []).

derangement(F, _, [], A) -> F(lists:reverse(A));
derangement(F, N, Xs, A) ->
  lists:foreach(
    fun(X) ->
      if
        N =/= X -> derangement(F, N + 1, lists:delete(X, Xs), [X | A]);
        true -> false
      end
    end,
    Xs).
derangement(F, N) -> derangement(F, 1, iota(1 , N), []).

montmort_number(1) -> 0;
montmort_number(2) -> 1;
montmort_number(N) when N > 2 ->
  (N - 1) * (montmort_number(N - 1) + montmort_number(N - 2)).

% 別解
montmort_number1(N, N, A, _) -> A;
montmort_number1(I, N, A, B) -> montmort_number1(I + 1, N, B, (I + 1) * (A + B)).
montmort_number1(N) when N >= 1 -> montmort_number1(1, N, 0, 1).

part_set(F, [], A) -> F(A);
part_set(F, [X | Xs], A) ->
  lists:foreach(fun(Y) -> part_set(F, Xs, [[X | Y] | lists:delete(Y, A)]) end, A),
  part_set(F, Xs, [[X] | A]).
partition_of_set(F, Xs) -> part_set(F, lists:reverse(Xs), []).

fold_with_index(_, _, A, []) -> A;
fold_with_index(F, I, A, [X | Xs]) ->
  fold_with_index(F, I + 1, F(X, I, A), Xs).

bell_number(N, N, [X | _]) -> X;
bell_number(I, N, Bs) ->
  bell_number(I + 1,
              N,
              [fold_with_index(
                 fun(X, K, A) -> A + comb_num(I, K) * X end,
                 0,
                 0,
                 Bs) | Bs]).
bell_number(N) -> bell_number(0, N, [1]).

group_part(F, _, _, [], A) -> F(A);
group_part(F, N, M, [X | Xs], A) ->
  lists:foreach(fun(Y) ->
      if
        length(Y) < N -> group_part(F, N, M, Xs, [[X | Y] | lists:delete(Y, A)]);
        true -> false
      end
    end,
    A),
  if
    length(A) < M -> group_part(F, N, M, Xs, [[X] | A]);
    true -> ok
  end.
group_partition(F, N, M, Xs) -> group_part(F, N, M, lists:reverse(Xs), []).

fact(0) -> 1;
fact(N) when N > 0 -> N * fact(N - 1).

group_part_num(0, _, M, A) -> A div fact(M);
group_part_num(K, N, M, A) -> group_part_num(K - N, N, M, A * comb_num(K, N)).

group_partition_number(N, M) -> group_part_num(N * M, N, M, 1).

% 因子の評価
factor([X | Xs]) when is_number(X) -> {X, Xs};
factor([X | Xs]) when is_list(X) -> {expression(X), Xs}.

% 項の評価
term_sub(V, [X | Xs]) when X =:= '*' ->
  {Y, Ys} = factor(Xs), term_sub(V * Y, Ys);
term_sub(V, [X | Xs]) when X =:= '/' ->
  {Y, Ys} = factor(Xs), term_sub(V / Y, Ys);
term_sub(V, Xs) -> {V, Xs}.

term(Xs) -> {Y, Ys} = factor(Xs), term_sub(Y, Ys).

% 式の評価
expr_sub(V, []) -> V;
expr_sub(V, [X | Xs]) when X =:= '+' ->
  {Y, Ys} = term(Xs), expr_sub(V + Y, Ys);
expr_sub(V, [X | Xs]) when X =:= '-' ->
  {Y, Ys} = term(Xs), expr_sub(V - Y, Ys).

expression(Xs) -> {Y, Ys} = term(Xs), expr_sub(Y, Ys).