テトロミノビンゴの未解決問題

前回、テトロミノビンゴの問題を解いたが、今回は
テトロミノ＋ビンゴ！の未解決問題
を解いてみた。
一見膨大な計算のように見えるが、まず、カードに書かれていない数字の出方はどうでもいいので、99!通りではなくて、25!通りの計算で十分なのは明らか。
25!( $\sim 10^{25}$ )通りでも手に負えないが、カードのマークの状態が同じものはまとめて計算できるので、 $2^{32} (\sim 3\times 10^7)$ 通りの計算で済む。

まず、マークの数が k個になった時点で、次の条件が満たされる確率 $P_m$ を考える。

カードの状態が m で (m は ${}_{25}\mathrm{C}_k$ 通りのうちのいずれか)、現時点とそれ以前に賞を獲得していない

k+1個のマークを持ち、テトロミノが含まれていない状態を m' とすると、 $P_{m'}$ は
$\Large P_{m'} \,=\, \frac{1}{25-k}\,\sum_m P_m$
で計算できる。右辺の m は、m' からマークをひとつ取り除いた全ての状態を動く。
分数の計算を避けるためプログラム中では $25\times 24 \times \cdots \times (26-k)P_m$ をハッシュ h に記憶している。

最終的なテトロミノ T の出現頻度 $Q_T$ は
$\Large Q_T \,=\,\sum_{m,c}\,\frac{P_m}{25-k(m)}$
となる。ただし状態 m は $2^{25}$ 通りのうちのいずれかで、c は25個のマス目のうちのいずれか。k(m) は m のマークの個数。m のマス目 c にマークをつけたときにテトロミノ T ができる全ての組み合わせ (m,c) について和を取る。
ループカウンタが k のとき、この加算の途中経過に $25\times 24 \times \cdots \times (25-k)$ を掛けたものを、プログラム中ではハッシュ prize_count としている。

カードの状態を表す変数 mark は前回配列だったが、大きなハッシュのキーにするので、簡単なデータ型のほうが良さそうと思って、整数型に変更して要素が1ビットの配列として使っている。
TetrominoBingoクラスは、前回のコードから余分なメソッドを削って、mark を整数型にするために find_tetromino の中身を少しだけいじった。

プログラムはカードの大きさが N×N マスの一般の場合について計算できるようにした。
N=3 の結果は 9!通り総当りの結果と一致している。
実行時間は手元の環境で N=4 のとき 6秒程度、N=5 のとき 4時間程度。Ideone.com 程度の速度なら N=5 で 2〜3時間くらいと思う。
効率化のために、回転、反転で重なり合うカードをまとめて計算するコードも書いてみたが、かえって激重になった。

プログラムの出力は k マス埋まった状態での各賞の出現頻度。
N=5 のときの最終結果は以下のようになった。

 I :  22912778428485 / 405833817138750 =  5.65%
 L : 117184546420835 / 405833817138750 = 28.88%
 O :  28116512702610 / 405833817138750 =  6.93%
 S :  68722249268210 / 405833817138750 = 16.93%
 T :  74170727776841 / 405833817138750 = 18.28%
 - :  94727002541769 / 405833817138750 = 23.34%

class TetrominoBingo
  def initialize
    tetromino_data = [
      [:I, [0,0], [0,1], [0,2], [0,3]],
      [:L, [0,0], [0,1], [0,2], [1,2]],
      [:O, [0,0], [0,1], [1,0], [1,1]],
      [:S, [0,0], [0,1], [1,1], [1,2]],
      [:T, [0,0], [1,0], [2,0], [1,1]]
    ]
    @tetromino = {}
    tetromino_data.each do |name, *coord|
      4.times do
        @tetromino[coord_to_cell(coord)] = name
        @tetromino[coord_to_cell(coord.map {|x, y| [-x, y]})] = name
        coord.map! {|x, y| [-y, x]}
      end
    end
  end

  def coord_to_cell(coord)
    cells = coord.map {|x, y| x + (N + 1) * y} .sort
    top_left = cells.shift
    cells.map {|cell| cell - top_left}
  end

  def group_to_tetromino(group)
    top_left = group.sort!.shift
    @tetromino[group.map {|cell| cell - top_left}]
  end

  def find_tetromino(mark, cell)
    group = [cell]
    k = 0
    while c = group[k] do
      [c-1, c+1, c-N-1, c+N+1].each do |d|
        next if mark[d] == 0 or group.member?(d)
        group << d
        return if group.size > 4
      end
      k += 1
    end
    return group_to_tetromino(group) if group.size == 4
  end
end


# 出力用
def output(k, prize_count)
  puts "k = #{k}"
  denom = (N**2-k+1 .. N**2).inject(:*)
  prize_count[:-] = 0
  prize_count[:-] = denom - prize_count.values.inject(:+)
  gcd = prize_count.values.inject(0) {|g, count| g.gcd(count)}
  prize_count.each do |prz, count|
    puts "#{prz} : #{count/gcd} / #{denom/gcd} = #{100.0*count/denom} %"
  end
  puts
end


# カードは N×N マス
N = 5
tet = TetrominoBingo.new
all_cells = (0 ... N ** 2).map {|pos| pos + pos / N}
prize_count = {:I => 0, :L => 0, :O => 0, :S => 0, :T => 0}
h = {0 => 1}

# このループ1周でマークが k個の状態から k+1個の状態を計算
(N ** 2).times do |k|
  prize_count.each_key {|m| prize_count[m] *= N ** 2 - k}
  new_h = Hash.new(0)
  h.each do |mark, count|
    all_cells.each do |cell|
      next if mark[cell] == 1
      new_mark = mark | 1 << cell
      if prize = tet.find_tetromino(new_mark, cell) then
        prize_count[prize] += count
      else
        new_h[new_mark] += count
      end
    end
  end
  h = new_h
  output(k + 1, prize_count)
end