なんとな~くしあわせ?の日記

「そしてそれゆえ、知識そのものが力である」 (Nam et ipsa scientia potestas est.) 〜 フランシス・ベーコン

動的計画法で組み合わせの総数 nCr を求めてみる

もっと早い組み合わせの総数 (combination)の求め方

動的計画法 - DP(Dynamic Programming)

最近AtCoderの問題をちょろちょろ見ているのですが、その中でいわゆるDP(Dynamic Programming)が勝負の鍵を握っているということがわかりました。組み合わせの求め方についてはいろいろやってきましたが、以下の方法では巨大な組み合わせを算出すると時間がかかってしまいます。

nantonaku-shiawase.hatenablog.com
nantonaku-shiawase.hatenablog.com

Combinationのための動的計画法

英語圏のサイトだとDP自体ちゃんと体系化されているっぽいので、以下のサイトからアルゴリズムを学び、やり方を汎化してみます。

www.csegeek.com

Rubyによるサンプルプログラム
  • 要は2次元配列にnCrの答えを予め用意しておき、後で使う時は添字アクセスするのだ
def choose(n, r)
  1 if r == 0 or r == n

  # store C(n,r) in a matrix
  c = Array.new(n+1).map{Array.new(r+1,0)}
  i,j = 0

  for i in 0..n
    # C(i,0) = 1 for i = 0...n
    c[i][0] = 1
  end
  for j in 0..r
    # if n = 0, C(n,r) = 0 for all 'r'
    c[0][j] = 0
  end

  for i in 1..n
    for j in 1..r
      if i == j
        # C(n,n) = 1
        c[i][j] = 1
      elsif j > i
        # case when r > n in C(n,r)
        c[i][j] = 0
      else
        # apply the standard recursive formula
        c[i][j] = c[i-1][j-1] + c[i-1][j]
      end
    end
  end

  return c[n][r]
end

n = 5
r = 2
comb = choose(n,r)

puts "C (#{n},#{r}) = #{comb}"
処理内容を表で表現してみる
  • nCr , n = 3, r = 2 を求める

1. 行 n、列 r の格子を考える

(0,0) (0,1) (0,2)
(1,0) (1,1) (1,2)
(2,0) (2,1) (2,2)
(3,0) (3,1) (3,2)

2. C(i,0) = 1 for i = 0...n

(0,0) = 1 (0,1) (0,2)
(1,0) = 1 (1,1) (1,2)
(2,0) = 1 (2,1) (2,2)
(3,0) = 1 (3,1) (3,2)

3. if n = 0, C(n,r) = 0 for all r

(0,0) = 0 (0,1) = 0 (0,2) = 0
(1,0) = 1 (1,1) (1,2)
(2,0) = 1 (2,1) (2,2)
(3,0) = 1 (3,1) (3,2)

4. j と i、それぞれ 1から最大値まで

  • i == j であれば、1を設定(初期値っぽい)
  • j > i であれば、0を設定(計算に関係ないから)

それ以外は

  • c[i][j] = c[i-1][j-1] + c[i-1][j]
    • 斜め左上のセルと上のセルを合計したものがそのセルの値

(0,0) = 0 (0,1) = 0 (0,2) = 0
(1,0) = 1 (1,1) = 1 (1,2) = 0
(2,0) = 1 (2,1) = 2 (2,2) = 1
(3,0) = 1 (3,1) = 3 (3,2) = 3

実際使う時は

  • nCr の NとRについては、固定値で先にすべての場合を求めておく方がよさそう

参考

http://program-study.hatenablog.com/entry/solve_project_euler_15

どうやらこれを使えばプロジェクトオイラーの15番は解けるようだ。

AOJ - DSL_1_A (Union-Find) を解いてみた

DSL_1_Aは、Coursera Algorithm I で最近習ったばかりのUnion-Findです

AOJ - DSL_1_A

互いに素な集合 Union Find| データ構造ライブラリ | Aizu Online Judge

quick find で求めてみる

自分の回答

時間切れでアウトです

  • 授業でも言われてましたが、quick-findは要素の更新を行う操作がどうしても遅いのです
    id.collect! { |element|
      (element == updated) ? updating : element
    }

quick-unionで求めてみる

自分の回答

これだと余裕で解けました。

追記 2018/07/11

  • 嘘です、07:93秒もかかってるので使い物になりません。。。

quick-unionのミソは、unionにしてもfindの操作にしてもどちらの場合でも最初に木構造の最上位を求めることだと思われます。(それをしないと、配列上で複数の親をもつ構造が表現できない)

# coding: utf-8
class QuickUnion
  attr :id
 
  def initialize(n)
    @id = (0...n).to_a
  end
 
  def root(i)
    while i != @id[i] do
      i = @id[i]
    end
    i
  end
 
  # find, check if p and q has the same root
  def same?(p, q)
    root(p) == root(q)
  end
 
  # union, to merge components containing p and q
  # set id of p's root to the id of q's root
  def unite(p, q)
    i = root(p)
    j = root(q)
    @id[i] = j
  end
end
 
lines = $stdin.read
array = lines.split("\n")
 
n = array[0].split(" ")[0].to_i
q = array[0].split(" ")[1].to_i
 
un = QuickUnion.new(n)
 
for i in 1..q
  com = array[i].split(" ")[0]
  x   = array[i].split(" ")[1].to_i
  y   = array[i].split(" ")[2].to_i
 
  if com == '0'
    # unite(4,3) id[4] = 3, 4の頭を3にする
    un.unite(x, y)
  else
    # same
    puts un.same?(x, y) ? "1" : "0"
  end
end

しかしやっぱりRubyは使いやすい

計算量の削減 O(α(N))

経路圧縮とrankの実装を行うことで処理時間をO(α(N))まで削減できます。

自分の回答

・00:22秒で終了しました。これでようやくAtCoderで使えるレベル

# coding: utf-8
class UnionFind
  attr :id, :rank
 
  def initialize(n)
    @id   = (0..n).to_a
    @rank = Array.new(n,0)
  end
 
  def root(i)
    if @id[i] == i or i == @id[@id[i]]
      i
    else
      @id[i] = root(@id[i])
      @id[i]
    end
  end
 
  # find, check if p and q has the same root
  def same?(p, q)
    root(p) == root(q)
  end
 
  # union, to merge components containing p and q
  # set id of p's root to the id of q's root
  def unite(p, q)
 
    i = root(p)
    j = root(q)
 
    return if i == j
 
    if @rank[i] < @rank[j]
      @id[i] = j
    else
      @id[j] = i;
      @rank[i]+=1  if @rank[i] == @rank[j]
    end
  end
end
 
 
lines = $stdin.read
array = lines.split("\n")
 
n = array[0].split(" ")[0].to_i
q = array[0].split(" ")[1].to_i
 
un = UnionFind.new(n)
 
for i in 1..q
  com = array[i].split(" ")[0]
  x   = array[i].split(" ")[1].to_i
  y   = array[i].split(" ")[2].to_i
 
  if com == '0'
    # unite(4,3) id[4] = 3, 4の頭を3にする
    un.unite(x, y)
  else
    # same
    puts un.same?(x, y) ? "1" : "0"
  end
end

AOJ - ITP1_7_D (行列の積) を解いてみた

ITP1_7_Dは、今自分の中でブームの行列の積の問題です。生で計算しないにしても、機械学習にも大いに関係があります。

AOJ - ITP1_7_D

Matrix Multiplication | Aizu Online Judge
自分の回答

行列の積をプログラムに起こすと?

もし自分がJavaで行列計算をするのであれば、JAMA: Java Matrix Package を使う気がします。JAMAは見た感じちゃんとオブジェクト指向的な気がします。ちなみにJavaの行列計算は全体的にパフォーマンス悪めっぽいので実際の大規模データには向かないかもしれません。

しかし、今回は競技プログラミング的なWEBテストなのでライブラリは使えません、なので生の二次元配列を使います。

問題の内容

問題自体は2つの行列を単純に掛け算するのみです。インプットの例は以下、

1:数値n, m, lが与えられる

3 2 3

2:行列 A[n][m] が与えられる、縦の数がn, 横の数がm

1 2
0 3
4 5

3:行列 b[m][l] が与えられる、縦の数がm, 横の数がl

1 2 1
0 3 2

問題の解法

手計算だとわかるのですが、これをプログラム化するのに結構苦戦しました。

行列計算の参考リンク


行列同士を掛け算すると、外側の要素が答えの要素になります。

[行列の積]
●[m×n型][n×p型]=[m×p型]
積が定義されるためには,左の列数と右の行数が等しくなければなりません。

結果は,[m×p型]となります。

なので、解答を保持する行列C = A*bとおくと、

4:行列 C[n][l] を0で初期化する

0 0 0
0 0 0
0 0 0

5:行列 C[n][l]変数i,j,kを使って埋めます。

変数は以下の範囲で動きます、が最初は行列Cの出来上がりの構造を考えたほうがいいでしょう。

  • 0 <= i < n
  • 0 <= j < m
  • 0 <= k < l