読者です 読者をやめる 読者になる 読者になる

本履歴

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

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

「の」の字形の無限に循環する配列のk番目

K-th Value of Japanese Hiragana Character "の"-like Infinitely Circulating Array

コンピュータは有限しか扱えないので前の状態から次の状態が決まる系は最終的に必ず循環します。そんな循環する状態列のk番目をここでは求めたいと思います。例えば、みんな大好きフィボナッチ数列1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, ...で純粋数学的にはこのまま増大していくのですが、有限しか扱えないコンピュータでは10^9+7で割った余りとかにして32bitに収まるよう妥協します。ここでは試しに11で割った余りにして最初の項から見てみると1, 1, 2, 3, 5, 8, 2, 10, 1, 0, 1, 1, 2, 3, 5, 8, 2, ...となり、すぐに循環してしまいました。フィボナッチ数列は定義から前の2項の和で次の項が決まるので、実際は数列内に1, 1パターンが2回目に出てきた所で計算(2項を足して11で割る)を打ち切って、それまで計算した結果を使ってk番目の値は簡単に分かります(これを後で見ます)。一般的には鳩の巣原理 - Wikipedia から、11*11+1番目まで計算すれば必ず既に循環していると言えます。ただし、10^9+7で割る例の場合は、最悪(10^9+7)^2+1番目まで計算しないと循環しない可能性があるので、有名な行列のpowmodを使うべきです。前述のフィボナッチ数列の例では一番最初に戻ってしまいましたが(円環)、一般的には途中の状態に戻ることが多く、その場合、自分は配列が「の」の形に似ているなあと常々感じています(一番最後に戻る=同じ状態が永遠と続く場合は「の」の円が潰れて「つ」になってしまいますが(笑))。英語でなんと表現するんでしょう? 似たような形のアルファベットはbやd, e,, p, qって結構ありますね。実際のプログラムは簡単で、どこに戻るのかで場合分けをすればOKです。ret_ptが「の」の字の合流点で、円周長がsize - ret_ptです。

#!ruby
#”の”の字型無限循環配列のk番目を求める

#ary: 循環配列
#ret_pt: 配列最後まで行った時にどこに戻るか
#k: 求める番目、当然配列のサイズより大きいことが予想される
#0-indexed
#O(1)

def no_at(ary, ret_pt, k)
    size = ary.size
    return ary[k] if k < size
    return ary[k % size] if ret_pt == 0
    return ary.last if ret_pt == size - 1
    k -= ret_pt
    r = k % (size - ret_pt)
    ary[ret_pt + r]
end

フィボナッチ数列の例で使うときはno_at([1, 1, 2, 3, 5, 8, 2, 10, 1, 0], 0, 1000)でもいいですしno_at([1, 1, 2, 3, 5, 8, 2, 10, 1, 0, 1], 1, 1000)でもOK。ではどうやって循環していることを調べるかについてですが、これはRubyだと現在の状態をハッシュに入れておけば二度目に出てきた時判別できます。オススメなのがハッシュの値を配列のindexにする方法。Ruby likeなpseudo codeだと

S = 系の初期状態
H = Hash.new
A = Array.new
Until H[S]
  A << S
  H[S] = A.size-1 #配列A内の状態Sの位置
  S = NextState(S)
End

Print no_at(A, H[S], k)

で計算量はO(|A|)すなわち状態の総数になります。したがってこれを使うときは状態の場合の数がどのくらいになるか事前に見積もることが必要です。 実際、mod 11フィボナッチ数列の最初の15項目を上記に当てはめて求めてみると、

#!ruby
M = 11
s = [1, 1]
h = Hash.new
a = Array.new
until h[s]
    a << s
    h[s] = a.size - 1
    s = [s.last, (s.first + s.last) % M]
end

15.times{|k|
    p no_at(a, h[s], k)
}

# [1, 1]
# [1, 2]
# [2, 3]
# [3, 5]
# [5, 8]
# [8, 2]
# [2, 10]
# [10, 1]
# [1, 0]
# [0, 1]
# [1, 1]
# [1, 2]
# [2, 3]
# [3, 5]
# [5, 8]

となりうまく行ってます。ちなみに、今回は純粋に状態を返していますのでフィボナッチ数列の実際の項を求めるには、リターン値の加工が必要です。このようなのの字ループですが、応用例は他にも 線形合同法 - Wikipediaなど枚挙にいとまがないくらいです。 以下、この考え方が使えるオンラインジャッジ問題の二例です。

Infinity Maze

例えば、Infinity Maze | Aizu Online Judgeでは、ロボットが与えられた迷路上で直進し壁にぶつかったら右に向きを変えてまた直進ということを繰り返していく時のL歩進んだロボットの最終位置と向きが問われますが、ロボットの迷路上の位置と向きが与えられればそれ以降の動きが完璧に決まるので、そのうちロボットは循環することが分かります。つまり「の」の字メソッドが使えます。状態数は迷路のサイズ(100*100)とロボットの向き(4通り)の組み合わせで高々40,000となり、Rubyでも十分間に合います。

Google Code Jam Qualification Round 2010 C. Theme Park

Dashboard - Qualification Round 2010 - Google Code Jam

こちらの問題も同じ考えが適用できます。問題の趣旨は、ジェットコースターに複数グループ(Max 1,000)がキューになって並んで居て、乗り終わったグループは再び整列することを繰り返す時に一日に何人がジェットコースターに乗るかを答えます。ジェットコースターは一日108回走りキャパシティは109で、各グループはグループメンバー全員が乗れない場合は次のを待ちます。これもキューの最初のグループが決まればその後の状態は決まるということと状態数が1000であることから、高々1000+1回シミュレーションすれば「の」の字無限ループに入るはずです。実際、問題のサンプルにもあるグループ[1, 4, 2, 1]キャパ6、4回運転の例を取ってみると、

[1, 4, 2, 1] =>5人
[2, 1, 1, 4] =>4人
[4, 2, 1, 1] =>6人
[1, 1, 4, 2] =>6人
[2, 1, 1, 4]

となり無事循環しました(ちなみに総搭乗者数は21)。この問題が少し難しいのは108回目の状態(どのグループが一番最初に並んでいるか)ではなく、各状態に付随する数値の総和を求めることです。これもno_atメソッドを少し変形するだけで求められます。スタートからのの字の合流点までとクルクル回転するところと最後のところの3パターンです。実際はスタートから合流点と最後のところは一緒にできるので2パターンの場合分けですみます。メソッドno_sumは各配列に対し1回しか呼ばれないのでinjectしていますが、何回も呼ばれる場合は総和計算をDPにすべきです。

#!ruby
#”の”の字型無限循環配列のk番目までの総和を求める
def no_sum(ary, ret_pt, k)
    size = ary.size
    return ary[0..k].inject(:+) if k < size
    return (k / size) * ary.inject(:+) + ary[0..k % size].inject(:+) if ret_pt == 0
    return ary.inject(:+) + (k - size + 1) * ary.last if ret_pt == size - 1
    k -= ret_pt
    q, r = k.divmod(size - ret_pt)
    ary[0..ret_pt + r].inject(:+) + q * ary[ret_pt..size - 1].inject(:+)
end

これを用いてTheme Parkを解いてみると、下記のようになりました。ラージも十分通ります。注意点は合計用の配列bを用意しているのと、状態がどのグループが先頭にいるかで決まるためindex付きの配列groupを用意しているのと、全部のグループが乗車したらそれ以上乗車を打ち切るようにrotate_numがnを超えないように監視しているのと、0-indexedにするためr-1にしている所になります。本質的にはpseudo codeそのままなのがお分かりいただけるはずです。

#!ruby
def solve(r,k,n,g)
    group = g.each_with_index.to_a
    s = 0
    h = Hash.new
    a = Array.new #for state
    b = Array.new #for summation
    until h[s]
        a << s
        h[s] = a.size - 1
        passenger = 0
        rotate_num = 0
        while passenger + group.first.first <= k && rotate_num < n
            passenger += group.first.first
            group.rotate!(1)
            rotate_num += 1
        end
        b << passenger
        s = group.first.last
    end

    no_sum(b, h[s], r - 1)
end

以上をまとめると、前の状態から次の状態が決まる系では「の」の字型循環になり、そのk番目の状態(やk番目までの和)はループする箇所をdivやmodすることにより状態総数の計算量で求められます。 個人的には、循環は仏教の輪廻的な考え方なので好きです。この宇宙も有限個の原子で出来ているならばそのうち「の」の字循環することが分かります。えっ?もう何回も循環してる!?

補足

知り合いから「の」でなく「6」では?と言われた。確かに書き順からしてそうだ(笑)

Identifies mirrored and rotated one

したがって面白くするため、回転・鏡像を含めて同一視した場合にどのくらいになるか手で数えてみた。なんとなく偶数になると思っていたが、N=4の場合は7という数字が出てきて不思議な感じ。これは配置によって回転・鏡像変換8通り全てが異なる順列になる場合や1234=4321のように2つしか同一視できないなど結構複雑だ。そこで早速プログラムを組んでカウントしてみた。1,1,2,7,23,115まで求めて後はOEISに突っ込めばやはり誰かが既に求めていた(https://oeis.org/A000903
この数列に関する面白い結果は特になさそうでちょっと残念。

#飛車を相互にアタックしないように配置するパズル
#回転鏡像含め同一視した中でいくつの場合があるか数える
#0-indexed
#O(N!)

def rotate(a)
	a.each_with_index.sort_by(&:first).reverse.map(&:last)
end

def mirror(a)
	table=[*0...a.size].reverse.each_with_index.to_a
	a.map{|v|table[v]}
end

N=6
count=0
h={}
[*0...N].permutation(N).each{|perm|
	next if h[perm]
	count+=1
	until h[perm]
		h[perm]=true
		perm=rotate(perm)
	end
	perm=mirror(perm)
	until h[perm]
		h[perm]=true
		perm=rotate(perm)
	end	
}

p count

N Rooks = N!

チェスにハマったので、8 Queensのカウンターパート的8 Rooksを考えてみた。チェスでRook(ルック、lookと発音は違う)は将棋で言うところの飛車の動きをする駒で、それを8*8のチェスボード上に互いに攻撃しあわないよう置くことを考えることになる。これは少し考えれば明らかなように順列そのままである。すなわち、最初の列に置けるパターンは8通り。次の列は最初に置いた駒の行を外して7通り(以下同)となって順列だ。クイーンの斜めの効きがないため非常にシンプルになる。