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

本履歴

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

Faster DP Implementation of Integer Partition Function with the Complexity O(N^2)

そこで色々調べた所、蟻本にO(N^2)の方法が載っている事がわかって一件落着。こちらもDPだけど、分割する数を取るのが違う点。

DP[n,m]=nをm個以下で分割するときの総数とするとDP[n,m]=DP[n,m-1]+DP[n-m,m]

nをm個以下の整数の和にするときに0も含めて考えているのがナルホド。そして、0を含むものと含まないもので分けて、うまく次元が低い処理にそれぞれ帰着できているのがすごい(小並感)

#!ruby
M=1000000007
def integer_partition_faster(n)
	dp=Array.new(n+1){Array.new(n+1,0)}
	dp[0][0]=1
	1.upto(n){|j|
		0.upto(n){|i|
			if i-j>=0
				dp[i][j]=(dp[i][j-1]+dp[i-j][j])%M
			else
				dp[i][j]=dp[i][j-1]
			end
		}
	}
	dp[n][n]
end

さて、実際の両者のコードをよく眺めると、結構似てる気がしませんか? 意図的に添字を同じにしたのもあるけど。というわけで、両アルゴリズムで使われている変数のdpを比べてみるとなんと驚き(ドラムロールオン♪)、中身が全く一緒になっているじゃありませんか!! Curiouser and curiouser. Don't you think so?

整数の分割って数論の一分野になっていますが、実は面白い定理が成り立つのです:

整数 n を分割したときの成分の数が k 個以下になるような分割の総数は、成分が k 以下の整数となるような n の分割の総数に等しい。

はい、そのまんまですね(笑)。なんか競プロの練習を通じてこの定理を証明したっぽい感じがしましたので、熱が冷めないうちにblog記事にしてみました。
その名もズバリ整数の分割って本に証明とかヤング図形との対応とかその他諸々興味深いことが詳しく載ってたはず。

整数の分割

整数の分割

Naive DP Implementation of Integer Partition Function with the Complexity O(N^3)

整数を分割したい病に患ってしまった。見る整数はなんでも、4=1+1+1+1=1+1+2=1+3=2+2のように分割してしまう恐ろしい病気だ。ちなみに4の場合の分割数は5である(自身も含める)。Wikipedia(https://en.wikipedia.org/wiki/Partition_(number_theory) )を読んでいた所、この数を計算するアルゴリズムとして次のようなものが載っていた:

補助関数 p(n,m)をm以下の数の和でnを表現する場合の数とするとp(n,n)が求めるもので、p(0,i)=1、p(i,j)=p(i,i) where i

これをそのままコードに落とすとO(N^3)でとても遅い。遅いアルゴリズムは泣きたくなる。

#!ruby
M=1000000007
def integer_partition_naive(n)
	dp=Array.new(n+1){Array.new(n+1,0)}
	dp[0][0]=1
	1.upto(n){|k|dp[0][k]=1}

	1.upto(n){|i|
		1.upto(n){|k|
			if k>i
				dp[i][k]=dp[i][i]
				next
			end
			1.upto(k){|j|
				dp[i][k]=(dp[i][k]+dp[i-j][j])%M
			}
		}
	}
	dp[n][n]
end

Pure Implementation of Dijkstra's Shortest Path Algorithm with Priority Queue in Ruby

プログラム言語 Ruby での純粋ダイクストラ最短経路アルゴリズム(優先順位付きキュー利用)です。O(m*log(n))の実行時間です(n:ノード数、m:辺の数)。
外部ライブラリなどに依存していないため、このままコピペすれば競技プログラミングで使えると思います。Ruby1.9.3..Ruby2.3.0で動きました。1.9.3は2.3.0の1.8倍時間がかかりました。
一応、AOJのGraph II - Single Source Shortest Path II(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ALDS1_12_C)とAtCoderのABC#035 D - トレジャーハント(http://abc035.contest.atcoder.jp/tasks/abc035_d)は通ったので、Rubyのバージョンが高い&&寛容(?!)なオンラインジャッジ(通常時間制限の数倍余裕くれるところ)では通ると思います。しかし、C++なんかと同じタイム制限のところではまず通らないと思います。やはり十数倍..数十時間がかかるようです。

利用上の注意点は、edges, costs共に0-indexedに加工して下さい。問題文で町に番号が1..nと割り振られーとかある場合に注意が必要。またアルゴリズムの性質上負のコストがあるとうまく動きません。プログラムはsから到達できる全ノードを巡った後最小コスト配列を返しますが、目的地がeと明確になっているのなら、visited[e]=trueになったら即終了したほうが速い。競技プログラミングの性質上、凶悪入力にも打ち勝たないといけないのであまり意味はないかもですが。
precursorには、アンコメントの後、最小コストのルートの1個前のノードが入ります。e=precursor[e] while e で後ろから前に辿れます。

以上です。よろしくお願い致します。

#!ruby

#Last Updated at 04/06/2016
#最短経路問題をダイクストラwith優先順位付きキューアルゴリズムで解く
#入力:
#edges: nノードの隣接リスト(0-indexed) e.g. [[1,2,3],[5],[0,1],...] means
# node 0からの辺はnode 1,2,3へ出ている、node 1からはnode 5のみ等々
#costs: コスト(辺の長さなど)、costs[x][y]で正数を返せばOK。
#s: 開始点。0..n-1
#出力:
#sからの最小コストのリスト dijkstra(edges,costs,s)[e]でs-e間の最小コスト
#到達出来ない点はINF
#経路を知りたい場合はprecursorのアンコメントする。precursor[k]でnode kの前のノードがわかる

INF = 1 << 60 #プログラム簡易化のため利用。全コスト和より大きい値にする

#O(m*log n) where m:辺の数
def dijkstra(edges, costs, s = 0)
	#初期化
	n = edges.size #ノード数
	node_to_pq_id = Array.new(n, nil) #nodeとpq.push時に割与えられたkeyの対応を管理
	#precursor = Array.new(n,nil)
	d = Array.new(n, INF) #sから別ノードまでの(これまで分かった)最短コスト
	d[s] = 0
	pq = PriorityQueue.new
	edges[s].each{|u|
		id = pq.push(Node.new(costs[s][u], u))
		d[u] = costs[s][u]
		node_to_pq_id[u] = id
		#precursor[u] = s
	}
	visited = Array.new(n, false) #既に訪れたか?
	visited[s] = true

	#本体処理
	until pq.empty?
		min = pq.pop
		cost,v = min.cost, min.node
		visited[v] = true
		d[v] = cost
		edges[v].each{|w|
			next if visited[w]
			next if d[w] <= cost + costs[v][w]
			#precursor[w] = v
			d[w] = cost + costs[v][w]
			#ノードwまでの冗長だったコストをよりお得なv経由コストにPQ内で変更
			if node_to_pq_id[w]
				pq.change_key(node_to_pq_id[w], Node.new(d[w], w))
			else
				#PQに初めてノードを入れる場合
				id = pq.push(Node.new(d[w], w))
				node_to_pq_id[w] = id					
			end
		}
	end
	d
	#[d, precursor]
end

#PQ用データ構造、コード簡易化のため
class Node
	attr_reader :cost, :node
	def initialize(cost, node)
		@cost = cost
		@node = node
	end

	def <= other
		@cost <= other.cost
	end

	def < other
		@cost < other.cost
	end

	def >= other
		@cost >= other.cost
	end

	def > other
		@cost > other.cost
	end
end

#Last Updated at 04/06/2016

#優先順位付きキュー
#ストアするデータは<=>で比較可能でなければならない
#小さい順か大きい順かは”符号を逆に”コメント箇所を変更
class PriorityQueue
  attr_reader :elements

  def initialize
    @elements = [nil] #ダミー。1-indexedにするため
    @position = {} #idのデータが@elementsのどこにあるかindexを返す
    @id = 0 #付与するid。1ずつインクリメント
  end

  #PQにデータを追加。付与されるidはchange_keyするときに指定する
  #つまり3を2回追加してもidは別になる
  #change_keyしない時は@position,@idの処理は全て不要
  #O(log n)
  def add(element)    
    @id += 1
    @elements << [element, @id]
    @position[@id] = size
    bubble_up(size)
    @id
  end

  #最小(大)値をPQから取り出し返す
  #O(log n)
  def pop
    exchange(1, size)
    max = @elements.pop
    bubble_down(1)
    max.first
  end

  #id(add時に付与される整数)による値の変更
  #値では動作できない(同じ値が入ることがあるため)
  #成功:self、失敗:nil
  #O(log n)
  def change_key(target_id, new_element)
    index = @position[target_id]
    return nil unless index
    if @elements[index].first <= new_element #符号を逆に
      @elements[index][0] = new_element
      bubble_down(index)
    else
      @elements[index][0] = new_element
      bubble_up(index)
    end
    self
  end

  #ストアされているデータ数
  def size
    @elements.size - 1
  end 

  def empty?
  	size == 0
  end

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

  def bubble_up(index)
    parent_index = (index / 2)

    return if index <= 1
    return if @elements[parent_index].first <= @elements[index].first #符号を逆に

    exchange(index, parent_index)
    bubble_up(parent_index)
  end

  def bubble_down(index)
    child_index = (index * 2)

    return if child_index > size

    not_the_last_element = child_index < size
    left_element = @elements[child_index]
    right_element = @elements[child_index + 1]
    child_index += 1 if not_the_last_element && right_element.first <= left_element.first #符号を逆に

    return if @elements[index].first <= @elements[child_index].first #符号を逆に

    exchange(index, child_index)
    bubble_down(child_index)
  end

  def exchange(source, target)
    @position[@elements[source].last], @position[@elements[target].last] = target, source
    @elements[source], @elements[target] = @elements[target], @elements[source]
  end
end

if __FILE__ == $0
  n = 120000
  m = 500000
  bridged = {}
  costs = Array.new(n){Hash.new}
  edges = n.times.map{|i|
    if i == 0
      costs[i][1] = i
      bridged[[0, 1]] = true
      [1]
    elsif i == n - 1
      costs[i][n - 2] = n - 2
      bridged[[n - 1, n - 2]] = true
      [n - 2]
    else
      costs[i][i - 1] = i - 1
      costs[i][i + 1] = i
      bridged[[i, i - 1]] = true
      bridged[[i, i + 1]] = true
      [i-1, i+1]
    end
  }

  m.times{
    i, j = rand(n), rand(n)
    next if i == j
    next if bridged[[i, j]]
    edges[i] << j
    costs[i][j] = rand(n)
  }

  p dijkstra(edges, costs).last
end

ruby-prof

 %self      total      self      wait     child     calls  name
 10.76      7.970     3.695     0.000     4.275  1930414   PriorityQueue#exchange
 10.60      3.642     3.642     0.000     0.000 18678231   Array#[]
  6.42      3.806     2.206     0.000     1.599  3920171   PriorityQueue#size
  5.62      2.239     1.931     0.000     0.309  4960817   Hash#[]=
  5.08      4.941     1.744     0.000     3.197        3   Integer#times
  4.92      2.387     1.689     0.000     0.698  3750646   Node#<=
  4.48      1.540     1.540     0.000     0.000  7557191   Array#first
  2.98      1.023     1.023     0.000     0.000  4469026   Array#[]=
  2.77      4.905     0.953     0.000     3.952   120000   Array#each
  2.68      0.922     0.922     0.000     0.000  4520167   Fixnum#-
  2.48      0.851     0.851     0.000     0.000  4465190   Fixnum#<=
  2.31      0.793     0.793     0.000     0.000  3920172   Array#size
  2.24      0.770     0.770     0.000     0.000  3860829   Array#last

CodeIQ

CodeIQで出題されていたマヨイドーロ問題を解いてみました。出題者は数学ガールでお馴染みの結城先生です!
使ったアルゴリズム動的計画法ですが、どういう風に私がこの解に到達したのかを時系列を追って説明したいと思います。

問題文を読んですぐに思いついたイメージはピンボールで、Xからボールが落ちてきてA,B,Cにあるフリッパーで方向が変わってYに落ちていく映像を想像しました。

最初は定石的に類似の問題を過去解いたことがあるのかどうかを自問してみました。”任意のグラフが与えられた時にグラフ上の点Pから点Qにn ステップで到達するルートの総数は隣接行列のn乗の(P,Q)成分”という事実を知っていたので、使えそうかどうか突っ込んで考えてみました。本質的には、N回反転してYに出るルートの総数が求めるものですから、

N回反転して1ステップでYに出るルートの総数+N回反転して2ステップでYに出るルートの総数+N回反転して3ステップでYに出るルートの総数+・・・

で求められそうです。ポイントは、ステップ数が多くなるにつれ、反転数も必然的に増える必要があるので、この級数の項は途中から零になるということです。このアイデアは見込みがありそうでしばらく考えていましたが、どうも反転というのをうまく表現できないことがわかってきました。というのも、隣接行列は単に繋がっているかどうかの情報しかないので、反転というのをうまく表せないんですね。
そこで、もう一度from scratchで考えてみることにしました。9回反転してYに出るルートの総数を求めることを当面の目標にします。”9回以下の”にしなかったのは、一度に扱う数が1,2,3,4・・・と増えてしまい複雑になると思ったからです。後で”以下の”で考えたほうが扱いやすいことが分かった場合は、ここに戻ってくることにします。

F(N):=N回反転してYに出るルートの総数

と定義すると、F(9)が求めるものですね。これだけでは、ただ言い換えしたにすぎないように見えますが、抽象化は大事です。さて、ここから「F(8)の答えを知っていたらF(9)は分かるだろうか?」と自問してみます。先程までの文章のまま考えていたのではできなかった進展ですね。さて、ここでparity(偶奇性)に気付きます。というのも、偶数回の反転はZに出てしまうので、Yには絶対にたどり着かないからです。そこで、Zに出るルートも一緒に考えると良さそうです。

G(N):=N回反転してZに出るルートの総数

と定義します。自問する内容は「G(8)の答えを知っていたらF(9)は分かるだろうか?」になります。G(8)の全てのルートは最後マヨイCを通りますので、そのCで反転すればYに出ますね。ただし、これがF(9)の全てを網羅しているかというとそうではないですね。というのも、G(8)には最後A->B->C->Zと辿るルートも含まれているので、Cまで来ないでその前のBで反転してYに行くこともありえますから。以上の分析から、GというZに出る総数というだけでは分類が荒くて、最後にBを通ってZに来るかどうかの情報も必要なことがわかりました。逆に、この情報があればF(9)を網羅できるでしょうか? できますね。途中の右往左往は無視すると、究極的には、Yに出るにはBで反転して来るのか、Cで反転して来るのか、に2種類しかありえません。(これを考えている時に気が付いたのは、マヨイB上で反転を繰り返すことが許されるかどうかです。つまり B->B->Bみたいなその場でクルクル向きを変えるのはありかどうか。これは問題のサンプルから無しなのがわかります。) したがって、まとめると、X を出た後の右往左往は置いておいて、途中8回反転した上最後A->Bとなるルートの総数と8回反転した上最後B->Cとなるルートの総数の合計でF(9)は計算できます!(ここでZに出ることは一旦忘れ、マヨイB,Cにとりま辿り着くルートを考えることにします!) では、”8回反転した上最後A->Bとなるルートの総数”を計算するにはどうしたらいいでしょうか? 言い換えると、X->B->◯->・・・◯->A->Bとなるルートで途中8回反転しているものですね。F,Gはもう役に立たないので捨てて、新しい関数を定義します。

H(P,Q,N):=X->B->◯->・・・◯->P->Qとなるルートで途中N回反転しているルートの総数
(P,QはA,B,Cのいずれか)

とします。問いを言い直します。H(A,B,8)はどう計算すべきでしょうか?(この問いの裏には、当然ながら、Hで再帰的に計算できると嬉しいな、という願いがこもっています(笑)) これは簡単ですね。H(B,A,7)そのものですね。つまり、7回の反転でBからAに行き、Aで8回目の反転をするパターンしかありません。それでは、F(9)を計算しようとした時の片割れのH(B,C,8)はどうでしょうか? 必ずH(A,B,8)を含むはずですね。A->Bときて、そのまま直進するパターンです。その他にも、H(C,B,7)も含みますね。C->BときてBで反転するパターンです。逆にこれら以外はありえません。あとは、マヨイドーロの対称性からH(B,A,7)やH(C,B,7)なども再帰的に計算可能です。
以上をまとめると、Hには次の漸化式があります。

H(A,B,N)=H(B,A,N-1)
H(C,B,N)=H(B,C,N-1)
H(B,A,N)=H(A,B,N-1)+H(C,B,N)
H(B,C,N)=H(C,B,N-1)+H(A,B,N)

この順番に計算することで、もれなく計算ができます。初期値はどうでしょうか? 反転0の時は、ただ右に進むだけですので

H(A,B,0)=1
H(C,B,0)=0
H(B,A,0)=0
H(B,C,0)=1

です。H(A,B,0)=1については、XをAとみなしても問題自体に影響はないので無用な複雑化を避けるためそうしています。
出力すべき数はF(N)を足し合わせたもので、F(N)はHを使って

F(N)=H(B,A,N)

です。あとはこれをコードに落とすだけです。Rubyは多次元配列の表記がうざいので、シンボルで次元を減らしています。

最終的にサブミットしたのは下記のRubyプログラムです(コメントは削除&変更)。結果的にnにでかい数があったけれどパスしました。
これがフィボナッチ数列になるのは式の形からわかるけれど、フィボナッチ数列の和が再びフィボナッチ数列になる事実が使えることは考えもしませんでした。フィボナッチ数列の単項を求めるのなら、行列のn乗をO(ln(n))で計算すればもっと速いですからね。

#!ruby

#dp[[:P2Q,i]]は、既にi回反転してある状態でPからQへ向かうルートの総数
dp={}
dp[[:A2B,0]]=1
dp[[:C2B,0]]=0
dp[[:B2A,0]]=0
dp[[:B2C,0]]=1
n=gets.to_i
1.upto(n){|i|
	dp[[:A2B,i]]=dp[[:B2A,i-1]]
	dp[[:C2B,i]]=dp[[:B2C,i-1]]
	dp[[:B2A,i]]=dp[[:A2B,i-1]]+dp[[:C2B,i]]
	dp[[:B2C,i]]=dp[[:C2B,i-1]]+dp[[:A2B,i]]
}

puts (0..n).map{|i|dp[[:B2A,i]]}.inject(:+)

問題と解説
https://codeiq.jp/magazine/2015/12/35521/

CodeIQ

CodeIQで開催されていた@cielavenirさんのRestricted WordsのRubyで提出した自分のソース公開。わーい、評価5頂きました。
方法はmethod_missingでメソッド名を取得する方法。
それだけだとつまんないので、本文を英語の文にしてみました。
後で、スペース出力が必要であることが分かったので、無理やり引数数(?、redundant?)で作りましたw 2**5=32。

#!ruby

def Object.const_missing name
  const_set name,name
end

def method_missing name,*args
  puts args*(args.size**args.first.size).chr unless args.empty?
end

I.am.bored.of Hello, World
How.about.you?

問題
https://codeiq.jp/ace/cielavenir/q431

競技プログラミング

昨日参加できなかったAtcoder#014の問題をRubyで解いたので。全部解けたのは初めてかも。

A. 君が望むなら世界中全てのたこ焼きを赤と青に染め上げよう

偶奇性を問われてる。ワンライナーで。

puts gets.to_i%2<1 ? "Blue" : "Red"

B. あの日したしりとりの結果を僕達はまだ知らない。

しりとりを判定。単語のおしりと次の単語のあたまが等しいか、と既出でないかをハッシュでチェック。

turn=0
words={}
pre_word=""
gets.to_i.times{
  word=gets.chomp
  if words[word] || ((pre_word!="") && pre_word[-1]!=word[0])
    puts turn<1 ? "LOSE" : "WIN"
    exit
  end
  words[word]=true
  pre_word=word
  turn=1-turn
}
puts "DRAW"

C. 魂の還る場所

また偶奇性を問われた。最初の2個の起き方は重要でないのと、どんな4つのRGBパターンが来ても連鎖で消して2個にできるから、それ以降も連鎖で消せるとダメ元で出したら通った。

gets(p)
puts [?R,?G,?B].map{|c|$_.count(c)%2}.inject(:+)

D. grepマスター

x=y=0の時を考えてみて、次に空いている間隔を考えて、それらがx+yのサイズで埋まるかどうかを考えていたら、間隔の集合をソートして二分探索すればよさげに気づいた。x+yより大きい間隔はx+yすべて使えるから。
次に、x+yで間隔全体をカバーできてしまう場合に、間隔の小さい順から足していくことが必要だったので、DPでメモっておいた。

#find the smallest idx s.t. key<=ary[i] for all idx<=i
def bsearch(ary,key,l,r)
  return nil if ary[r]<key
  return nil if l>r
  ret=r
  m=(l+r)/2
  if key<=ary[m]
    temp=bsearch(ary,key,l,m-1)
    ret=temp ? temp : m
  else
    temp=bsearch(ary,key,m+1,r)
    ret=temp if temp
  end
  ret
end

all,n,m=gets.split.map(&:to_i)
l=n.times.map{gets.to_i}
ary=l.each_cons(2).map{|l1,l2|l2-l1-1}.sort
dp=[ary.first]
1.upto(ary.size-1){|i|dp << dp.last+ary[i]}

m.times{
  x,y=gets.split.map(&:to_i)
  pre=[x,l.first-1].min
  post=[y,all-l.last].min
  middle=n
  unless ary.empty?
    idx=bsearch(ary,x+y,0,ary.size-1)
    if idx
      middle+=(x+y)*(ary.size-idx)
      middle+=dp[idx-1] if idx>0
    else
      middle+=l.last-l.first+1-n
    end
  end
  puts pre+middle+post
}

コードゴルフ

そういえば、前に参加していたCodeIQのこっちのゴルフコンペもソースコード公開。最初のCゴルフの2次元連立方程式の問題。
C言語は通常まったく書かないけど(=ゴルフ専用)、頑張ってみたら、なんとか上級超えられた。わーお。
scanfする時に&a,&b,&c,..とかしたくなかったのでアドレスを渡して処理が工夫の箇所。
printfに変な引数付けないと、うまく処理されなかったので、泣く泣くつけている。全く意味不明。激おこ。
CよくCらないからCようがないCー

問題
https://codeiq.jp/ace/ozy4dm/q246

結果発表
http://codeiq.hatenablog.com/entry/2013/04/08/175858

最終版 118バイト

main(n,f,d,c,e,b,a,g){for(;~scanf("%d",&a+n%6);n--%6/5&&printf("%d %d\n",(d*e-b*f)/g,(a*f-c*e)/g,&f,&d,&c))g=a*d-b*c;}