Cable master とドント方式と中央値の中央値

はじめに

競技プログラミングで見つけたとある問題が見かけの違う別の問題と本質的に同じで、その最小計算量オーダが意外と小さかったというお話です。

Cable master 問題

その問題は「プログラミングコンテストチャレンジブック」 (通称「蟻本」。表紙に蟻の絵が描いてあるから) の中級編の問題
「Cable master」
（英語版の元問題はこちら）

長さがそれぞれ $ L_i $ であるような $ N $本の紐があります。これらの紐を切って、同じ長さの紐を $ K $本作るときの最長の長さを求めなさい。
答えは小数以下2桁までを切り捨てて出力しなさい。

蟻本では二分探索を使って解の近似値を求めていますが、こういう問題を見たら近似値ではなく厳密な値を求めたくなるのが人間の性（さが）ですよね。

厳密な解はN本の紐のどれかをある本数で割った長さになるはずです。
なぜならちょうど割り切れた長さでなければまだ伸びしろがあるわけですから解が最長の長さであることに矛盾します。
つまり解は

\{ \frac{L_i}{j} \mid 1 \leqq i \leqq N,\, 1 \leqq j \leqq K \}

の中にあります。
この集合を長さ順にソートして順番にその長さで N本取れるか試していけば解が見つかります。
ナイーブにこの方法を行うと計算量オーダは大きいですが、工夫すると O(N * log(N)) の計算量で求められることがわかりました。（K に影響されない！）

Cable master 問題とドント方式

これが計算量の限界かなと思い答えがネット上に載っていないか探したところ、まさにこの問題を取り上げている
情報学会の会誌「情報処理」連載の「プログラム・プロムナード」の「ケーブルマスタ」
を見つけました。
そこには Cable master問題がドント方式で解けると書いてあります。

ドント方式とは Wikipedia の記事によると

政党名簿比例代表において、議席を配分するための最高平均方式（highest averages method）のひとつである。
この方式はベルギーの数学者ヴィクトル・ドント(Victor Joseph Auguste D'Hondt)から名づけられた。

（日本でも比例代表の選挙で使っています）

ドント方式の手順は
（上記 wikipedia の記事の「配分」の「1議席ずつ配分するのではなく次のように考えても同様である。」のところに書いてある方法）
(1) 各政党の得票数を１，２，３･･･の整数で割る
(2) 一人当たりの得票数が多い順（割り算の答えの大きい順）に、1議席ずつ各政党に議席を配分する
です。

例えば 9議席を A党、B党、C党の3党で争っていて各政党の獲得得票数がそれぞれ 120票、90票、60票だったとします。
まず、各政党の得票数を１，２，３･･･の整数で割ると

党	/ 1	/ 2	/ 3	/ 4	...
A党	120 / 1	120 / 2	120 / 3	120 / 4	...
B党	90 / 1	90 / 2	90 / 3	90 / 4	...
C党	60 / 1	60 / 2	60 / 3	60 / 4	...

これらをまとめて値の順に並べて上位9個を取ると（どの党のものかわかるように後ろに党名をつけておきます）

  \{ 120(A), 90(B), 60(C), 60(A), 45(B), 40(A), 30(C), 30(B), 30(A)\}

結果、A党は4議席、B党は3議席、C党は2議席獲得することになります。

ドント方式が Cable master とどう対応するかというと

Cable master	ドント方式
$ K（切って作る同じ長さの紐の本数） $	総議席数
$ N（元々の紐の本数） $	党の数
$ L_i（元々の各紐の長さ） $	各党の得票数
各紐から取れる紐の本数	各党の獲得議席数
答えである紐の長さ	最後に議席獲得したときの割り算の答え （上記の例では 30）

という対応になります。

先ほどの Cable master の解法では

\{ \frac{L_i}{j} \mid 1 \leqq i \leqq N,\, 1 \leqq j \leqq K \}

を長さ順にソートして順番に N本取れるか試していけばよいと書きましたが、実は単に長さ順でN番目の値を取ってくればよいのでした。

Cable master問題がドント方式となぜ本質的に同じ問題なのかについては上記の「プログラム・プロムナード」で見ていただくとして次に進みましょう。

ドント方式を計算量 O(N) で解く

Cable master問題がドント方式と同じだとわかったので今度は 「ドント方式 計算量」 で検索すると、なんと最悪ケースでも計算量 O(N) で解けると書いてある論文を見つけました。

「議席配分法に対する線形時間アルゴリズム」
（ドント方式の計算量 O(N) のアルゴリズムが書いてあるのは、定理 4.2 です。）

なぜ O(N) で計算できるかは論文を見ていただく（ちょっと難しい）ことにして次に進みましょう。

計算量 O(N) で解くプログラム

上記の論文に書いてあるアルゴリズムにしたがってプログラムを書いてみましょう。Ruby で書いてます。

# s_total: 総議席数。論文では S
# vs: 各党の獲得投票数。論文では {v_i}
def dhondt(s_total, vs)
  # 計算しやすいように実数にしておく
  vs = vs.map(&:to_f)
  n = vs.length # 論文でも n
  vsum = vs.sum # 論文では Σ v_i

  # step1
  # 理想配分数
  ss = vs.map do |v| s_total * v / vsum end
  # 初期配分数
  s0s = ss.map(&:to_i)

  # step2
  # 追加候補集合
  cs = vs.map.with_index do |v_i, i|
    s_i = ss[i]
    s0_i = s0s[i]
    # 論文
    # a < 1 + n/S, a = s / s_i なので s < (S + n) * s_i / S
    (s0_i + 1 .. ((s_total + n) * s_i / s_total - 1).ceil).map {|s| s / s_i}
  end

  # step3
  # 不足配分数
  d = s_total - s0s.sum
  # 論文では α_dth
  a_dth = d == 0 ? 1.0: select_kth(cs.flatten, d-1)

  # step4
  # 追加配分数
  dss = cs.map do |as| as.take_while {|a| a <= a_dth}.count end
  # 各党の獲得議席数
  ss = s0s.zip(dss).map do |is, ds| is + ds end

  # Cable master問題では答えの紐の長さ
  z = vsum / s_total / a_dth

  [z, ss]
end

# ソートされていない配列からk番目に大きな数を求める
# k は 0オリジン
def select_kth(xs, k)
  # どう実装するか
end

ただし最後のメソッド select_kth は複雑なので少し解説してから実装することにします。
というのも上記の論文では、ソートされていない N個の数の集合の中から k番目に大きな数を最悪ケースでも計算量 O(N) で得ることを求めているのです。
すぐに思いつく方法は配列をソートして上から k番目の数を求める方法ですが、この方法ではソートするのに O(N*log(N)) かかってしまいます。

最悪ケースでも O(N) の計算量で k番目の数を求めるアルゴリズは見つかっています。（introselect という名前のようです）

パーティションベースの汎用選択アルゴリズムの中の「中央値の中央値」を用いたクイックセレクトにアルゴリズムが書かれています。

クイックソートと似たクイックセレクトという手法を使っていますがクイックソートと同様にピボットの選び方に失敗すると最悪ケースの計算量オーダが大きくなるのでさらに中央値の中央値という手法を使って計算量が大きくならないように工夫しています。

このアルゴリズムを使って select_kth を実装しましょう。

def select_kth(xs, k)
  return xs[0] if xs.length == 1

  pivot = select_pivot(xs)

  lows = xs.find_all do |x| x < pivot end
  lows_len = lows.length
  if k < lows_len
    select_kth(lows, k)
  elsif k == lows_len
    pivot
  else
    highs = xs.find_all do |x| x > pivot end
    select_kth(highs, k - (lows_len + 1))
  end
end

def select_pivot(xs)
  # 中央値の中央値を選ぶ
  median_of_medians(xs)
end

def median_of_medians(xs, i = xs.length/2)
  return median5(xs, i) if xs.length <= 5

  medians = xs.each_slice(5).map do |s| median5(s) end

  pivot = median_of_medians(medians)

  lows = xs.find_all do |x| x < pivot end
  lows_len = lows.length
  if i < lows_len
    median_of_medians(lows, i)
  elsif i == lows_len
    pivot
  else
    highs = xs.find_all do |x| x > pivot end
    median_of_medians(highs, i - (lows_len + 1))
  end
end

def median5(xs, i = xs.length/2)
  return xs.sort[i]
end

これでめでたしと思ったら

と、ここまで書いておいてなんですが、「「中央値の中央値」を用いたクイックセレクト」に

ただし、この方法では最悪時間は確かに線形になるが、平均時間は、実際にはピボット値を無作為に選ぶなどの素朴な方式の方が優れている。

と書いてあるとおり実際にはピボット選択に「中央値の中央値」を使わずに

def select_pivot(xs)
  # 要素をランダムに選ぶ
  xs.sample
end

とするのがよいでしょう。

さらに言えば現実の問題では O(N) も O(N * log(N)) も似たようなものなので Wikipedia のドント方式の記事の「配分」の「次のようにして1議席ずつ議席を配分する」に書いてある方法（最初に紹介した「情報処理」連載の「プログラム・プロムナード」にプログラムが載っています）を使うのが楽でよい思います。

最後はちょっと身も蓋もない結論になってしまいましたがドント方式の理論上の最小計算量にふれて面白かったのではないでしょうか。

これで（多分）本当にめでたしめでたし。