本履歴

購入した古本の履歴と時々プログラミング

数学パズル本に敢てプログラムで挑戦(3)

どんどん挑戦していく。明日GCJ本番なのに目が冴えてしまってあかん。本日もAlgorithmic Puzzlesより。

49. Missionaries and Cannibals(宣教師と人食い部族)

3人の宣教師と3人の人食い族が船で川を渡る時(同時に2人まで)、両岸それぞれで人食い族の人数が宣教師よりも真に多くならないようにして渡る方法を求めよ

何処かで見たことのある問題も、結局はグラフの問題に帰着。船がどちらの岸にあるかと宣教師、人食い部族がそれぞれ左岸に何人いるかの情報を頂点とするグラフを幅優先探索すれば良い。具体的には[boat,Missionary,Cannibal]という配列で[0,3,3] -> [1,0,0]のステップ数を求める。結果は11ステップ。実際グラフを書くと相当単純(ほぼ一本道)のグラフになる。やはり数的優位が崩れると食べられてしまうという条件が厳しい。というのも、宣教師、人食い族それぞれ左右の岸にいる状態を考えると、どちらか一方は数的優位が成り立たなくなってしまうから。ちなみに4名以上でやってみたら渡る方法は無かった。

#49. Missionaries and Cannibals
#Three missionaries and three cannibals must cross a river.
#Their boat can only hold two people, and it cannot cross the river by itself
#with no people on board. Each missionary and each cannibal can row the boat.
#If present, missionaries cannot be outnumbered by cannibals.
#How can all six get across the river with the fewest crossings?

N=3
class State
    attr_reader :row_num
    def initialize(boat=0,people=[N,N],row_num=0)
        @boat=boat
        @people=people
        @row_num=row_num
    end

    def legitimate?
        rm,rc=@people
        lm,lc=N-rm,N-rc
        return false if rm>0&&rc>rm
        return false if lm>0&&lc>lm
        true
    end

    def accomplished?
        @boat==1&&@people==[0,0]
    end

    def next_move
        ret=[]
        m,c=@people
        if @boat==0
            ret<<State.new(1,[m-1,  c],@row_num+1) if m>0
            ret<<State.new(1,[m-2,  c],@row_num+1) if m>1
            ret<<State.new(1,[  m,c-1],@row_num+1) if c>0
            ret<<State.new(1,[  m,c-2],@row_num+1) if c>1
            ret<<State.new(1,[m-1,c-1],@row_num+1) if m>0&&c>0
        else
            ret<<State.new(0,[m+1,  c],@row_num+1) if m<N
            ret<<State.new(0,[m+2,  c],@row_num+1) if m<N-1
            ret<<State.new(0,[  m,c+1],@row_num+1) if c<N
            ret<<State.new(0,[  m,c+2],@row_num+1) if c<N-1
            ret<<State.new(0,[m+1,c+1],@row_num+1) if m<N&&c<N
        end
        ret
    end

    def to_code
        @boat.to_s+@people.join
    end
end


init_state=State.new
queue=Array.new
queue.push init_state
visited=Hash.new

until queue.empty?
    state=queue.shift
    next if visited[state.to_code]
    visited[state.to_code]=state.row_num
    if state.accomplished?
        puts state.row_num
        exit
    end

    state.next_move.each{|next_state|
        next unless next_state.legitimate?
        next if visited[next_state.to_code]
        queue.push next_state
    }
end

puts"impossible"

数学パズル本に敢てプログラムで挑戦(2)

いったん書くとこれがなかなか面白いのがblogの中毒的な所、というわけで今日も昨日の続きでAlgorithmic Puzzlesを解いていこう。

45. A Knight’s Shortest Path(ナイトの最短路)

100*100のチェスボード上の左下コーナーにナイトが居る時、右上のコーナーへは最短何手で行けるか?

最短路とあるけれど、ダイクストラアルゴリズムは必要ない。というのも、全ての一手の価値は同じだから幅優先探索でOK。つまり、最初の位置から一手で進める場所を調べて、そこから更に一手で進める場所を・・・というのを、ゴールに辿り着くまで||進める場所がなくなるまで、行う。もちろん同じ場所を何度も行こうとするから、visitedというハッシュにそのマスまでの最短手数を入れてフラグとする。 直感的にはナイトのジャンプはイビツだからうまくたどり着けない可能性やゴール直前で調整する手数が必要かなと思う。

#!ruby

#45. A Knight’s Shortest Path
#What is the minimum number of moves needed for a chess knight
#to go from one corner of a 100 × 100 board to the diagonally opposite corner?

N = 100

def is_pos_regitimate?(x, y)
    ((0...N) === x) && ((0...N) === y)
end

Moves = [
                    [ 1, 2],
                    [-1, 2],
                    [ 1,-2],
                    [-1,-2],
                    [ 2, 1],
                    [-2, 1],
                    [ 2,-1],
                    [-2,-1]
                ]
Init_Pos = [    0,     0]
Goal_Pos = [N - 1, N - 1]

queue = Array.new
queue.push [Init_Pos, 0]
visited = Array.new(N){Array.new(N)}

until visited.dig(*Goal_Pos) || queue.empty?
    pos, move_num = queue.shift
    next if visited.dig(*pos)
    x, y = pos
    visited[x][y] = move_num
    Moves.each{|dx, dy|
        next_pos = [x + dx, y + dy]
        if is_pos_regitimate?(*next_pos)
            queue.push [next_pos, move_num + 1] unless visited.dig(*next_pos)
        end
    }
end

puts visited.dig(*Goal_Pos) || "Impossible!"

ここでは新しく覚えたdigメソッドを使ってみた。多重構造の時にコード数を減らせるし、今回のように*展開で配列も渡せて美味しい。結果は66というなんとなく納得できる数字に。というのもナイトの動きって1行って2行くわけだから、2回動けば対角線に(3,3)に行ける。従って、(100,100)に行くには66回動けばいいでしょうという寸法。他のNでも同じ結果に。N=2のときだけ解がないようだ。

   N | answer
-----------
  10 | 6
 100 | 66
1000 | 666

数学パズル本に敢てプログラムで挑戦(1)

Google Code Jamの季節になると活発になるこのblog。今日(今年w)は数学系パズル本であるAlgorithmic Puzzlesの問題を解いてみましょう。この本、出たばっかりの時にKindleで買ったんですが最近日本語にもなったみたいですね。本屋で見かけた、というか、手にとって見たらどこかで見たことある問題だらけだったので。

120. Penny Distribution Machine (石分配マシーン)

横に長い配列があり、一番左のセルに石が何個か入っている。2個の石が同じセルの中にあった場合は、1個にし右のセルに移動する操作をできなくなるまで行う時;

1)最終結果は石の移動の順番に依存するか?

2)n個の石で始めた時、何個セルが必要か?

3)この操作が止まるまで何回の操作がされるか?

というのが問題の趣旨(超訳)。詳しい解説は本を購入して頂くとして、これをプログラムを組んで解こうというのがこの記事の趣旨。コンピュータに関連する問題だしね 1)、2)は簡単。この操作の最終結果が2進数表現になることが分かれば、1)は依存しないし(厳密には加法の可換性)、2)は2進数にした時のビット数だ。 ただし、3がちょっとすぐには分からなかった。結局何回繰り上げするかだけど、何か法則があるのかどうか? こういう時にサクッとプログラムを作って実験するのが個人的にいいと思う。紙と鉛筆のオールドスクールを否定しているわけじゃないけど。一度プログラムを組むと、コードのパラメータをいじることで、3つの石になったら分配するとかに応用がすぐ効く。また、結構な数の例を見ることが出来るので、洞察を得やすい。

#!ruby
class PDM
    attr_reader :dist_num
    def initialize(n)
        @boxes=Array.new(10,0)
        @boxes[0]=n
        @dist_num=0
    end

    def dist!
        @boxes.each_with_index{|penny,idx|
            if penny>1
                @boxes[idx]-=2
                @boxes[idx+1]+=1
                @dist_num+=1
                return idx
            end
        }
        nil
    end

    def to_s
        @boxes.join(" ")
    end
end

1.upto(20){|n|
    box=PDM.new(n)
    while box.dist!
    end
    puts "%2d: %d"%[n,box.dist_num]
}

早速実行してみると

 1: 0
 2: 1
 3: 1
 4: 3
 5: 3
 6: 4
 7: 4
 8: 7
 9: 7
10: 8
11: 8
12: 10
13: 10
14: 11
15: 11
16: 15
17: 15
18: 16
19: 16
20: 18

増加列なのはそうだけど、パターンがなかなか見えてこない。列? 困った時のオンライン整数列大辞典で検索すると A011371 - OEIS

a(n) = n minus (number of 1’s in binary expansion of n)

nからnの2進数表現に出てくる1の数を引くってどういうことだ??

アハ!

この石分配マシーンの動作をもう一度考えてみよう。2個の石を右に移動する操作のトータル数を求めるんだけど、1操作につき1つ石が減ってる! すなわち石が減った数を数えれば良い、すなわち最初の石の数-残った石の数 それにしても、OEISはすごい。こんな数列まで載っているとは。

こんな感じで数学パズル本の問題をプログラムで解くのは結構面白いし教育的だ。競技プログラミングに特化している本の問題は、アルゴリズミックな問題だらけで初学者が楽しめないんだよね。それにこう言った問題は紙と鉛筆で解くことが前提なので、Straightforwardなプログラムで結構解けちゃったりする。問題のサイズが大きくなると解けないんだけど、そんなものはあまり求められない。いろいろな解法があったりするのも興味深い。というわけで、結構面白いので皆さんも試してみよう!

www.amazon.com

1次元配列中の連続した部分配列で和が最大となるものを求める

UVa: 11951 - Area

事の発端はこの問題(UVa Online Judge)。 長方形の数字の配置が与えられた時にサブ長方形で和が最大となるもの(ただしK以下)を求める問題。DP[x, y] = Sum(1, 1, x, y)を使ってO(N^2*M^2)の解法はTLEっぽいので、最終的に長方形の数字は全て正数であることから、

Sum(x1,y1,x2,y1)<Sum(x1,y1,x2,y1+1)<Sum(x1,y1,x2,y1+2)< ... <Sum(x1,y1,x2,M)

となり、この数列に対しK以下になるようバイナリサーチをしてO(N^2*M*log(M))でACとなり話はこれで終わりだったのですが・・・

Kadane's Algorithm

もっとうまい方法を模索していると1次元バージョンのMaximum subarray problem - Wikipedia, the free encyclopedia を発見。Kadane's algorithmと呼ばれているらしい。うーむ、沈思黙考、アハ、目からうろこが出た。以下私の翻訳。DPなんだけれど、配列Aの連続部分列で和が最大を求めるのに

DP[j]=A[i]+A[i+1]+...+A[j-1]+A[j]の最大値 for i=1 to j

と定義しこのDP[j]、すなわちおしまいの添字がjとなる連続部分列の中で最大和が分かったとすると

DP[1]=max{A[1], 0} //負の場合は取らないを選択
DP[j]=max{DP[j-1]+A[j],0} //負の場合は取らないを選択

が成り立つから、順次DPが求められて、元々欲しかった最大値はDPの中の最大値を返せばOKというアルゴリズム。これがうまく動くのが不思議ですがよく考えるとそうなる。最大の懸念点は、DP[j]=DP[j-1]+A[j]といえるのか、他にjで終わる連続部分列で最大のものがあるのではないかですが、あるとすると

DP[j-1]=A[i]+A[i+1]+...+A[j-1]
DP[ j ]=A[k]+A[k+1]+...+A[j-1]+A[j]

からA[i..k-1](or A[k..i-1])の部分和が負になる必要があり、DPの定義に矛盾するからそんなことはありえないi=kと言えます。 Rubyの実装は簡単。実は直前のdpの値しか使っていないのでdpという配列は省略できまス。それがwikipediaバージョンですね。実際の連続部分配列の箇所を求めたい時にはdpをアップデートする時にindex情報も保存すればOKです。

#!ruby
#O(n)
def kadane(ary)
    n = ary.size
    dp = Array.new(n + 1, 0) #index jで終わる部分配列の中での最大和
    n.times{|i|dp[i] = [dp[i - 1] + ary[i], 0].max}
    dp.max
end

Rubyで任意の順列の次を求める

next_permutation in Ruby

他言語にはnext_permutationがあるのにRubyに無いのはけしからん、ということでnext_permutationを実装しました。 比較可能なオブジェクトからなる配列を渡すだけです。O(n)なので[*1..n]の全順列をiterateする場合はArrayのpermutationの方が圧倒的に速いです。4倍の差がある。あくまでも途中から次のやつ知りたい時専用です。あと、向こうは配列の中のオブジェクトの大小を見ていないのも違う所。 bsearch_indexを使っているので、ちょっと古いRubyだとmethod_missingになりそう。その時はary.each_with_index.bsearch{...}.lastなんかで取り出せばOK。

#!ruby
#O(n)
def next_perm(ary)
    size=ary.size
    return nil if size<2
    #後ろから順に数値が減少する箇所を探す
    idx=nil
    (size-2).downto(0){|i|
        if ary[i]<ary[i+1]
            idx=i
            break
        end
    }
    #条件を満たす箇所が見つからない場合は、後ろから昇順(前からは降順)になっているので最後の順列
    return nil unless idx

    #idxより右側を順列と考えると、左から降順ということは、このpermutationでは大きくなれない
    #したがって、ary[idx],b=ary[idx+1...size]を"繰り上げ"処理する必要がある
    #実際には、bからary[idx]の次に大きいものを見つけ(idxでary[idx]が下がっているのでかならず見つかる)
    #ary[idx]と場所を交換しbをsortすれば良い
    ary[idx+1...size]=ary[idx+1...size].reverse
    jdx=ary[idx+1...size].bsearch_index{|v|ary[idx]<v}
    ary[idx],ary[idx+1+jdx]=ary[idx+1+jdx],ary[idx]
    ary
end

ちなみにGoogle Code Jan Round 1B 2009 B. The Next Numberもこれを使えばほぼ一発です。

def solve(str)
  a=str.chars.map(&:to_i)
  np=next_perm(a.dup)
  if np
    return np*""
  else
    ls="0"
    ls+=a.pop.to_s while a.last==0
    return a.pop.to_s+ls+(a.sort*"")
  end
end

各種Heap処理速度の比較1

The Heap implementation can make significant difference

Google Code Jamも2016 R1BがそろそろということでWEBをさまよっていたらいろいろな人の優先順位付きキューがgemとして公開されていたので、自分の実装と比較してみた。

GitHub - matiasbattocchia/lazy_priority_queue: A priority queue implemented using a lazy binomial heap. It allows change priority operation. を大いに参考にしている、というか比較方法パクっている。本レース競技者はリンク先の方の実装lazy_priority_queueと自作したbinary heap, skew heap, 最後にC言語拡張のPriorityQueue (supertinou)(フィボナッチヒープ)の4実装だ。後者は実行環境がサーバ側では動かない問題点があるので参考招待選手として走ってもらう。PQueue等他との比較はリンク先のlazyとの結果から類推できる。 結果は、C-extentionが速いのは当然として、やはり前回お知らせしたとおりSkew Heapが結構速かった(知ってた)。

API Design : Whether PQ.push(key, value) or PQ.push(value)

今回やってみて感じた気付きは、当初の目的であった速度とは全く別の所、つまりAPIをどう設計するか、だった。自作以外はpush(add, <<)する時に必ず利用者がユニークなキーを割り当てなければならない。つまり、push(key, value)だ。これは今回紹介した以外の他の優先順位付きキューどれでもそうだった。しかし自分的には後で値を更新しないかもしれないのに、単にヒープソートで使いたいだけかもしれないのに(push & pop)ユニークなキーを利用者に求めるのはなにかおかしい気がして、自作の物は内部的にキー管理をし、pushされた時にpushされたオブジェクトに紐づくid(キー)を付与する形にし、それを使うかどうかは利用者任せにしたんだけれど、世界の趨勢とは真逆だったようだ。下のベンチマークの処理で言うと、自作以外はなんか知らないnを追加で引数にとっている。どっちがいいんだろうか、なかなか難しいのお。 面白いのでもうちょっと続く。

#!ruby
require 'benchmark'
require 'lazy_priority_queue'
require 'priority_queue'

def iterator n, push, pop
  n.times do |i|
    (n - i).times do |j|
      push.call i.to_s + ':' + j.to_s
    end

    i.times do
      pop.call
    end
  end
end

N = 1_000

Benchmark.bm do |bm|
  bm.report('Lazy priority queue') do
    q = MinPriorityQueue.new
    iterator N, ->(n){ q.enqueue n, rand }, ->{ q.dequeue }
  end

  bm.report('My Binary Heap') do
    q = MyHeaps::PriorityQueue.new
    iterator N, ->(n){ q.push rand }, ->{ q.pop }
  end

  bm.report('My Skew Heap') do
    q = MyHeaps::SkewHeap.new
    iterator N, ->(n){ q.push rand }, ->{ q.pop }
  end

  bm.report('PriorityQueue (supertinou)') do
    q = CPriorityQueue.new
    iterator N, ->(n){ q.push n, rand }, ->{ q.delete_min }
  end
end
       user     system      total        real
Lazy priority queue       : 15.820000   0.250000  16.070000 ( 16.286343)
My Binary Heap            : 11.440000   0.200000  11.640000 ( 11.942219)
My Skew heap              :  5.970000   0.070000   6.040000 (  6.129232)
PriorityQueue (supertinou):  3.770000   0.050000   3.820000 (  3.848521)

Rubyだとダイクストラ最短路アルゴリズムはバイナリヒープよりスキューヒープ実装の方が速い?!

Skew Heap might be way faster than Binary Heap in Dijkstra Algorithm in Ruby

気になっていたデータ構造Skew heap - Wikipedia, the free encyclopediaを筆のすさびにRubyで書いてみたら、結構速い。もしやとchange_keyも作ってダイクストラの最短経路アルゴリズムに使うと、以前紹介したバイナリヒープ版より2,3割速いようだ。頂点・辺数が小さければ差はほぼないけれど。Rubyは配列アクセス遅いからね、、narrayが標準ライブラリに入ってくれると嬉しいんだが・・・。いずれにせよ、速いことは良いことだ。ただし注意しなくちゃいけないのは、この構造はバランストで無いので、ハックできるようなコンテストだと、スキューヒープが線形な木構造になるような入力を意図的に作られてしまい、簡単に落とされると思う。 SkewHeapNodeクラスのelementに比較可能なオブジェクトを入れるようになっている。ではよろしくお願いします。

#!ruby
#Last Updated at 04/29/2016

#Skew Heap
#ストアするデータは<=>で比較可能でなければならない

class SkewHeapNode
    attr_accessor :element, :left, :right, :parent
    def initialize(element, left = nil, right = nil)
        @element = element
        @left = left
        @right = right   
        @parent = nil
    end
end

class SkewHeap
    attr_reader :size, :root
    def initialize
        @root = nil
        @size = 0
        @position = {}
        @id = 0
    end

    def add(element)
        @id += 1
        @size += 1
        node = SkewHeapNode.new(element)
        @position[@id] = node
        @root = merge(@root, node)
        @id
    end

    def pop
        return nil if @size==0
        @size -= 1
        element = @root.element
        @root = merge(@root.left, @root.right)
        @root.parent = nil if @root
        element
    end

  def change_key(target_id, new_element)
    node = @position[target_id]
    return nil unless node
    if heap_order_property_satisfied?(node, new_element)
        node.element = new_element
    else
        temp = merge(node.left, node.right)
        temp.parent = node.parent if temp
            if node.parent
            if node.parent.left == node
                node.parent.left = temp
            else
                node.parent.right = temp
            end
        else
            @root = temp
        end
        node.element = new_element
        node.left = nil
        node.right = nil
        node.parent = nil

        @root = merge(@root, node)
        @root.parent = nil if @root
    end
    self
  end

  def empty?
    @size == 0
  end

  alias :push :add
  alias :<< :add
  alias :length :size

private
    #p's parent is never changed
    #O(log n)-expected time
    def merge(p, q)
        return p if q.nil?
        return q if p.nil?
        p, q = q, p if p.element > q.element
        q.parent = p
        p.right = merge(p.right, q)
        p.right, p.left = p.left, p.right
        p
    end

    def heap_order_property_satisfied?(node, new_element)
        if node.parent
            return false if node.parent.element > new_element
        end
        if node.left
            return false if node.left.element < new_element
        end
        if node.right
            return false if node.right.element < new_element
        end
        true
    end
end