CodeIQ：番号の対応表で作るグループ

私自身が表題の問題を解いた時のプログラムについて解説します。
問題の詳細は「番号の対応表で作るグループ」（CodeIQ）を参照してください。

問題の概要

問題を引用します。
CodeIQ感謝祭では、各参加者には申し込んだ順番に申込番号が付与されていました。
ただし、座席は申し込んだ順番ではなく、会場に到着した順番に着席しました。
このとき、「着席した座席番号」と「申込番号」によってグループを作ることを考えます。

座席番号と申込番号が同じ場合、その人は単独でグループになります。
異なる場合、申込番号に対応する座席番号に座っている人を順に辿り、辿れる人でグループを作ります。

例えば、6人が以下のような座席番号と申込番号だったとします。

座席番号申込番号
12
24
33
41
56
65

このとき、「1, 2, 4」「3」「5, 6」の3つのグループに分けられます。

3人の場合、座席番号と申込番号の対応は全部で以下の6通りが考えられますので、作られるグループ数の期待値は(3+2+2+1+1+2) / 6 = 1.8333…です。

座席番号申込番号
11
22
33
座席番号申込番号
11
23
32
座席番号申込番号
12
21
33
座席番号申込番号
12
23
31
座席番号申込番号
13
21
32
座席番号申込番号
13
22
31

標準入力からグループの数 n が与えられたとき、作られるグループの数の期待値が n を超えるような最小の参加人数 m を求め、その人数 m を標準出力に出力してください。
なお、nは6以下の整数とします。
また、キャンセルした人は考えないため、座席番号と申込番号はともに1～m が一つずつ付与されます。

【入出力サンプル1】
標準入力
2
標準出力
4

私のプログラム

Rubyで解答しています。

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#!/usr/bin/ruby
 
$Memo = {}
 
# A130534のFomulaのT(n,k)を実装したもの
# 高速化のためメモ化
def func(n,k)
    if n==0 && k==0 then return 1
    elsif n==0 then return 0
    elsif n<k then return 0
    else
        v = $Memo[[n,k]]
        if v != nil then return v
        else
            v = func(n-1, k-1) + n*func(n-1,k)
            $Memo[[n,k]] = v
            return v
        end
    end
end
 
# OEIS A130534の数列を作成する
# nは0始まりなので問題に比べて1小さい
def A130534(n)
    ret = []
    for k in 0..n
        ret << func(n,k)
    end
    return ret
end
 
# 問題を解く
def solve(n)
    i = 0
    ptn = 1
 
    # 1<=iの範囲で問題の条件を満たすまでループする
    begin
        i += 1
        arr = A130534(i-1)  # A130534は0始まりなので1小さい
        sum = 0 #全ての組み合わせのグループ数の合計
        ptn *= i    # パターン数
 
        arr.each_with_index{|v, j|
            sum += (j+1)*v
        }
 
#       printf("arr=%s, sum=%d, ptn=%d\n", arr, sum, ptn)
    end while n*ptn >= sum
 
    return i
end
 
# main
while line = gets
    line.strip!
    if line.empty? then next end
 
    p solve(line.to_i)
end

解説

これも★★の問題ですが前回の問題に続いて非常に難しいです。
自力で考えようとするとかなりの数学能力を要求されると思います。私はOEISに頼ったので正しいプログラムを作成することはできていますが、理屈はよくわかっていません。

考え方

コードのA130534()の理屈がわかっていないので説明できません。代わりに、OEISのA130534という数列にたどり着くまでにやったことを説明します。

問題を読んですぐにはコードを思いつかなかったので、なんらかパターンを発見する必要があるだろうと思い、単純な方法で解答を求めるコードを書いてみました。
この時にやったことは、全ての組み合わせを列挙し、それを実際にグループ分けしてみるという方法です。Rubyの場合、Array#permutation()があるので全てのパターンを列挙するのは簡単です。グループ分けは配列の添字を席番号とし、配列に0〜n（実際の番号より1小さい）を入れ、それをたどることでチェックするという何の効率化も考えない方法です。
ただし、入力値はnではなくmにしました。パターンを見つけるためなのでmを1,2,3,…と手動で指定できた方が便利だからです。

とりあえずこのコードを実行してみたらm=10では全く時間的に間に合わず、この時の期待値は3に達していませんでした。
このことから、n=1からn=2よりもn=2からn=3の方がmがずっと大きくなるということは明らかで、n=6となるととんでもなく大きな数になるだろうということがわかりました。つまり、動的計画法にしろ、メモか再起にしろ、どんな手法を使ってもパターンを列挙してチェックするという方法はダメで、計算だけで済ます方法を考えなければならないことがわかりました。

1時間程度数式を考えてみたのですがわかりません。
次に期待値をm=1〜5くらいまで出力して、OEISで検索してみましたがありませんでした。
今度はm=1〜5くらいの「グループ数1〜mの場合のパターン数の合計（1,2,6,24,120,720）を検索してみたらヒットしました（A000142）。
これはm=1,2,3,…のグループ数ごとの数列(1,1,1,2,3,1,6,11,6,1,…)もあるだろうということで探してA130534にたどり着いたというわけです。

A130534()

一般的に良くない関数名の代表のような名前ですが、OEISの番号に一致させています。
中では0〜引数n（問題文のmに相当する）までkを変化させながらfunc()（これもひどい名前です）を呼び出し、結果をretに保持します。このretの添字+1がグループ数で、値がパターン数になります。

func()

全く説明できません。
コメントに書いてある通り、OEISのA130534のFomulaにあった数式をコードにしただけです。
ただし、OEISのコードは当然メモ化などされていない数式だったのでメモ化して高速化しました。メモ化しないとA130534(n=100)くらいで結構な時間がかかってしまいます。

solve()

1から入力値nまで順番にA130534()を呼び出し、期待値がnを超えるまで繰り返します（38〜49行目）。
コメントにも書いてありますがA130534()の引数は0始まり、入力値は1始まりなのでA130534()には1小さい値を渡します。A130534()からグループ数ごとのパターン数が返ってきたら要素ごとにグループ数×パターン数を計算して合計します（44〜46行目）。
変数ptnは全パターン数でi!になりますが、毎回一から計算するとかなり時間がかかりそうだったのでループごとに新たな値を掛けるだけにして計算量をいくらか節約しています（42行目）。
ループの終了条件ですが、問題文のように期待値を求めるのではなく、n×パターン数が合計値を超えたかどうかをチェックするようにしています。多分、浮動小数点でも誤差は問題になりませんが、値があまりにも大きくRubyの浮動小数点の上限を超えてしまうため整数で計算できるようにしました。Rubyは64bitの範囲を超えたら自動的に任意精度整数を使用するので、これなら正しく計算できます。
最終的にループを終了した時のiが答えになります。

雑感

問題の難易度は前回の問題と同じくらいと思いますが、今回の問題は試しに作った単純なコードの結果からパターンを作成して数えるような方法はダメなことが早々にわかったので大分早くできました。OEISで数列を見つけられなかったらできなかったと思いますが。
この解答のコードは多分まだ効率化できると思います。なんでかというとm=2の結果にはm=1の結果、m=3の結果にはm=2の結果が利用できそうだからです。もっとも、メモ化しているのでそれほど効率は良くならないようにも思えますが。
それから、この問題をC言語（などの任意精度整数が無い言語）でできるのでしょうか（任意精度整数の計算を実装すればできますが大変です）？ うまく計算すると64bit整数の範囲に収まるようにできるのかな？