第27章

整列を保った挿入

(やさしい版) Pearls of Functional Algorithm Design(関数プログラミングによるアルゴリズム設計の真珠)

どんな問題?

本棚を思い浮かべてください。棚には N 個の仕切り(スロット)があり、最初は全部空です。

床には本の山があります。本を1冊ずつ取り上げて棚に置いていくのですが、棚の上の本は常にアルファベット順に並んでいなければならないというルールがあります。新しい本を置く場所を空けるには、すでに並んでいる本を左右にずらす必要があります。

目標 本を全部並べ終わるまでの移動の合計回数を最小にしたい。移動の距離は問いません(1マス動かすのも100マス動かすのも「1回」と数えます)。

例:PEARLS を並べる

たとえば「PEARLS」の6文字を大きさ6の棚に並べる例です。

123456移動
P0
EP0
AEP2
AEPR1
AELPR1
AELPRS2

この例では合計6回の移動で済んでいます。

この章のゴール

本章では、合計 Θ(N log3 N) 回の移動で仕事を終えるアルゴリズムを作ります。この計算量が最良だと予想されています。アルゴリズムそのものよりも、「なぜこの計算量になるのか」を丁寧に解析するのがこの章の見どころです。

素朴なやり方(そしてなぜ遅いか)

まず最小移動回数を見てみる

N = 1 から 12 までについて、必要な最小移動回数 m(N) を表にすると次のようになります。

N = 1 2 3 4 5 6 7 8 9 10 11 12 m(N) = 0 1 2 4 6 8 11 14 17 21 24 29
豆知識 「PEARLS」が6回で並べられるからといって、あらゆる6文字の単語が6回で並べられるわけではありません。実は、6文字なら最悪でも8回で並べられる戦略を見つけるのは面白い脇道の問題です。m(N) の閉じた式はまだ知られていませんが、m(N) = Θ(N log3 N) だという強い証拠があります。

素朴な案:詰めて挿入

いちばん単純な案は、「新しい要素は必要なだけ既存の要素をずらして挿入する」というものです。もし k 番目の要素がそれまで挿入した全ての値より小さいと、それまでの k−1 個を全部1つずつずらすはめになります。

悪い場合 最悪の場合は N(N−1)/2 回の移動が必要になります。これは N2 のオーダーで遅すぎます。

やや改善:真ん中に挿入

「空いている領域の真ん中に置く」やり方もあります。最初の ⌊log N⌋ 個は移動なしで置けます。でも半分ぐらい埋まってくると結局移動が増え、また N2 オーダーに落ちてしまいます。

もう少しマシな案:特殊要素で均等再配置

ある整数 M を選び、iM 番目(i = 1, 2, …, k)を特殊要素と呼ぶことにします。特殊要素を挿入するときは、これまで並べた要素を全部動かして、空きスロットが均等になるように配り直します。

i=1k (iM − 1) = Θ(N2/M)
Θ(∑i=1N/M Ci M2) = Θ(MN log N)

ここで Ci = iM/(NiM) は特殊要素挿入後の「空きスロット1個あたりの周りの要素数」の目安です。

両方の合計をバランスさせるように M ≈ √(N/log N) を選ぶと、総コストは Θ(√(N3 log N)) となります。これは Θ(N1.5) より少し悪い程度で、まだ目標の Θ(N log3 N) には届きません。

改良版アルゴリズム

2つの改良を組み合わせる

目標に届くためには、次の2つを両方改良する必要があります。

改良1:特殊要素は「後ろほど密に」

特殊要素の番号を次のように選びます。

n1 = ⌊N/2⌋, n2 = ⌊3N/4⌋, n3 = ⌊7N/8⌋, …, nk = ⌊(2k−1)N/2k

ここで k = ⌈log N⌉ です。N = 1000 なら特殊要素は次の10個。

500, 750, 875, 937, 968, 984, 992, 996, 998, 999

特殊要素の挿入にかかる総コストはたかだか Θ(N log N) です(それぞれ最悪 N 回、それが k = ⌈log N⌉ 回だから)。

改良2:非特殊要素は「木」と「密度条件」で挿入

N 個のスロットを、大きさが均衡した二分木の葉として並べたと想像します。N = 11 のときの例:

[0 … 10] / \ [0 … 4] [5 … 10] / \ / \ [0,1] [2,3,4] [5,6,7] [8,9,10] /\ /\\ /|\ / \ [0][1] [2][3,4] [5][6,7] [8][9,10] │ │ │ /\ │ /\ │ /\ [0][1] [2][3][4] [5][6][7] [8][9][10]

挿入したい位置 p から根までの道は、p を含むだんだん大きな区間の列 v0, v1, …, vk になります。ここで、次の記号を使います。

フェーズ ini−1ni の間の要素を挿入している時期)における挿入戦略は、次の密度条件を満たす最小の j の区間 vj に挿入することです。

S(vj) < δ(i, j) L(vj)

ここで密度関数 δ は次のとおり:

δ(i, j) = ((2i − 1) / 2i)j/k

挿入時には、要素と S(vj) 個の既存要素を vj 全体に均等に配り直します。この操作のコストは S(vj) 回です。

密度条件はいつでも満たせる δ(i, 0) = 1 なので、葉 v0 でも空いていれば条件を満たします。また根 vk では S(vk) < ni ≤ δ(i, k)L(vk) が成り立ちます。つまり常に条件を満たす区間が経路上に存在するので、必ず挿入場所が見つかります。

DANGEROUS の挿入例

「DANGEROUS」の9文字を大きさ9の棚に挿入した結果は次のようになります(図 27.1)。

123456789移動
D0
AD0
ADN0
ADGN2
ADEGN0
ADEGNR3
ADEGNOR2
ADEGNORU7
ADEGNORSU0
図 27.1 DANGEROUS の挿入

総コストの上界

フェーズ i に挿入される非特殊要素は CN/2i 個です。最悪の場合、これらは全部同じ挿入点にぶつかります。要素 p が区間 vjp でシャッフル(再配置)を引き起こしたとします。

ここで、vjp の直前の再配置直後の占有数を S0(vjp) とし、Δ(vjp) = S(vjp) − S0(vjp) と置きます(前回配り直してから、新たに何個増えたか)。

次節で証明する 2 つの不等式が鍵になります。

S(vjp) < 4k 2i (Δ(vjp) + 1)(27.1)
p=1C Δ(vjp) ≤ kC(27.2)

これらを組み合わせるとフェーズ i の非特殊要素の総コストは 4k(k+1)N 以下となり、フェーズ全体を足すと

i=1k 4k(k+1) N = 4k2(k+1) N = Θ(N log3 N)
まとめ 特殊要素分 Θ(N log N) と非特殊要素分 Θ(N log3 N) を合わせて、総コストは Θ(N log3 N) に収まります。さらに詳しく見ると、最初の N/2 個は 4N log2 N 以下、最初の 3N/4 個は 8N log2 N 以下、…と分かります。もし棚のサイズを (1+ε)N と少し余裕を持たせれば、Θ((N log2 N)/ε) 回にまで下がります。

証明

(27.1) の証明

jpj と略します。示したいのは S(vj) < 4k 2i (Δ(vj) + 1)。使う性質は次の4つ:

これらから、次のように計算していきます。

2Δ(v_j) ≥ {Δ(v_j) ≥ Δ(v_{j−1}) と Δ の定義より} 2S(v_{j−1}) − 2S_0(v_{j−1}) ≥ {2S_0(v_{j−1}) ≤ S_0(v_j) + 1 ≤ S(v_j) + 1 より} 2S(v_{j−1}) − S(v_j) − 1 ≥ {S(v_{j−1}) ≥ δ(i, j−1) L(v_{j−1}) より} 2 δ(i, j−1) L(v_{j−1}) − S(v_j) − 1 ≥ {2L(v_{j−1}) ≥ L(v_j) − 1 より} δ(i, j−1) (L(v_j) − 1) − S(v_j) − 1 > {δ(i, j) L(v_j) > S(v_j) と δ(i, j) ≤ 1 より} (δ(i, j−1) / δ(i, j) − 1) S(v_j) − 2

したがって:

Δ(vj) + 1 > (1/2)(δ(i, j−1)/δ(i, j) − 1) S(vj)

δ の定義から

δ(i, j−1)/δ(i, j) − 1 = (2i/(2i−1))1/k − 1

x1/k = 2(1/k) log x、0 < y < 1 で 2y − 1 ≥ y/2 を使うと

δ(i, j−1)/δ(i, j) − 1 ≥ (1/(2k)) log(2i/(2i−1))

最後に 2i log[2i/(2i−1)] ≥ 1 を用いて、4k 2i (Δ(vj) + 1) > S(vj) が示せます。これが (27.1) です。

(27.2) の証明

より一般に、1 ≤ abC について ∑p=ab Δ(vjp) ≤ k(ba + 1) を示します。この左辺を P(a, b, k) と書き、k(木の高さ)への依存を明示します。

Δ(vjp) の値は「vjp を含む区間で最後に再配置した時(それを jq とする)以降に vjp に入った要素数」で、これは qp の間にある点の数に等しくなります(図 27.2 の最初の図)。

Δ(vjb) = d(0 ≤ dba)とすると、P は 2 つの和に分解できます。後ろ側の和は最大 j の値が jb − 1 以下だと分かるので:

P(a, b, k) ≤ max0≤d≤b−a (P(a, b−d−1, k) + P(b−d, b−1, jb−1) + d)

jbkP は第3引数について単調なので:

P(a, b, k) ≤ max0≤d≤b−a (P(a, b−d−1, k) + P(b−d, b−1, k−1) + d)

ba ≤ 1 か k = 0 なら P(a, b, k) = 0 なので、あとは単純な帰納法で P(a, b, k) ≤ k(ba + 1) が示せます。

実装

データ表現

Haskell では配列を「(インデックス, 値) のペアのリスト」で表します。

type Array a = [(Int, a)]
Dart // 配列は (インデックス, 値) のペアのリスト。 // インデックスは昇順、値も昇順(挿入対象)である前提。 typedef Arr<A> = List<(int, A)>;

ルール:各 i は 0 ≤ i < n、値は増加順。

ラベル付けと挿入の全体像

主要な挿入関数はこう定義します。

insertAll n = scanl (insert n) [] · label n
Dart // scanl: 空配列から始めて label 付き入力を順に insert していき、 // 各ステップの中間結果を全部集めたリストを返す。 List<Arr<A>> insertAll<A extends Comparable<A>>(int n, List<A> xs) { final labeled = label(n, xs); final steps = <Arr<A>>[<(int, A)>[]]; var acc = <(int, A)>[]; for (final ix in labeled) { acc = insert(n, acc, ix); steps.add(acc); } return steps; }

特殊要素を求める

specials :: Int → [Int] specials n = scanl1 (+) (halves n) halves n = if n == 1 then [] else m : halves (n − m) where m = n div 2
Dart // 半分ずつ削っていく: halves 11 = [5, 3, 1, 1] List<int> halves(int n) { final result = <int>[]; while (n > 1) { final m = n ~/ 2; result.add(m); n = n - m; } return result; } // 累積和: specials 11 = [5, 8, 9, 10] List<int> specials(int n) { final hs = halves(n); final result = <int>[]; var acc = 0; for (final h in hs) { acc += h; result.add(acc); } return result; }

例:specials 11 = [5, 8, 9, 10]。

ラベル付け

label :: Int → [a] → [(Int, a)] label n xs = replace 1 (zip [1..] xs) (specials n) replace i [] ns = [] replace i ((k, x) : kxs) ns | null ns = [(0, x)] | k < n = (i, x) : replace i kxs ns | k == n = (0, x) : replace (i+1) kxs (tail ns) where n = head ns
Dart // 各要素にフェーズ番号を付ける。 // 特殊要素(specials に含まれる位置)は 0、その間の要素は 1, 2, ... と番号付け。 List<(int, A)> label<A>(int n, List<A> xs) { final ns = specials(n); final result = <(int, A)>[]; var i = 1; var nsIdx = 0; for (var k = 1; k <= xs.length; k++) { final x = xs[k - 1]; if (nsIdx >= ns.length) { result.add((0, x)); } else if (k < ns[nsIdx]) { result.add((i, x)); } else { // k == ns[nsIdx] result.add((0, x)); i += 1; nsIdx += 1; } } return result; }

たとえば label 11 [1..11] は次を返します:

[(1,1),(1,2),(1,3),(1,4),(0,5),(2,6),(2,7),(0,8),(0,9),(0,10),(0,11)]

insert 本体

insert :: Ord a ⇒ Int → Array a → (Int, a) → Array a insert n as (i, x) = if i == 0 then relocate (0, n) x as else relocate (ℓ, r) x as where (ℓ, r) = ipick n as (i, x)
Dart // i == 0(特殊要素)なら (0, n) 全体を配り直す。 // そうでなければ ipick が返す区間 (l, r) で再配置。 Arr<A> insert<A extends Comparable<A>>(int n, Arr<A> as, (int, A) ix) { final (i, x) = ix; if (i == 0) { return relocate((0, n), x, as); } final (l, r) = ipick(n, as, ix); return relocate((l, r), x, as); }

relocate と distribute と spread

relocate :: Ord a ⇒ (Int, Int) → a → Array a → Array a relocate (ℓ, r) x as = distribute (add x (entries (ℓ, r) as)) (ℓ, r) as entries (ℓ, r) as = [x | (i, x) ← as, ℓ ≤ i ∧ i < r] add x xs = takeWhile (< x) xs ++ [x] ++ dropWhile (< x) xs
Dart // 区間 (l, r) にある既存要素を取り出し、x を挿入したうえで // 区間全体に配り直す。 Arr<A> relocate<A extends Comparable<A>>((int, int) lr, A x, Arr<A> as) { final (l, r) = lr; return distribute(add(x, entries(lr, as)), lr, as); } // 区間 [l, r) 内にある値だけを昇順のまま取り出す。 List<A> entries<A>((int, int) lr, Arr<A> as) { final (l, r) = lr; return [for (final (i, x) in as) if (l <= i && i < r) x]; } // 昇順リストに x を昇順を保って挿入。 List<A> add<A extends Comparable<A>>(A x, List<A> xs) { final result = <A>[]; var inserted = false; for (final y in xs) { if (!inserted && y.compareTo(x) >= 0) { result.add(x); inserted = true; } result.add(y); } if (!inserted) result.add(x); return result; }

distribute は「区間の外の要素はそのまま、区間内は spread で並べ直す」関数:

distribute :: [a] → (Int, Int) → Array a → Array a distribute xs (ℓ, r) as = takeWhile (λ(i, x) → i < ℓ) as ++ spread xs (ℓ, r) ++ dropWhile (λ(i, x) → i < r) as
Dart // 区間外(左と右)はそのまま、区間内は spread で並べ直す。 Arr<A> distribute<A>(List<A> xs, (int, int) lr, Arr<A> as) { final (l, r) = lr; final left = [for (final p in as) if (p.$1 < l) p]; final right = [for (final p in as) if (p.$1 >= r) p]; return [...left, ...spread(xs, lr), ...right]; }

spread はリストと区間を半々に分けて再帰的に配ります:

spread :: [a] → (Int, Int) → Array a spread xs (ℓ, r) | null xs = [] | n == 0 = [(m, head xs)] | n > 0 = spread ys (ℓ, m) ++ spread zs (m, r) where (n, m) = (length xs div 2, (ℓ + r) div 2) (ys, zs) = splitAt n xs
Dart // 要素列を区間 (l, r) に均等に配る。 // 中央 m を境に、リストを半々に分けて左右へ再帰。 Arr<A> spread<A>(List<A> xs, (int, int) lr) { final (l, r) = lr; if (xs.isEmpty) return <(int, A)>[]; final n = xs.length ~/ 2; final m = (l + r) ~/ 2; if (n == 0) return [(m, xs.first)]; final ys = xs.sublist(0, n); final zs = xs.sublist(n); return [...spread(ys, (l, m)), ...spread(zs, (m, r))]; }

ipick:挿入区間の選択

ipick :: Ord a ⇒ Int → Array a → (Int, a) → (Int, Int) ipick n as (i, x) = if p < q then (p, q) else head [(ℓ, r) | (j, (ℓ, r)) ← zip [0..] (ipath n p), let s = length (entries (ℓ, r) as), densityTest i j s (r − ℓ)] where (p, q) = ipoint n x as
Dart // 挿入区間を決める。まず ipoint で候補点 (p, q) を求め、 // p < q なら空きがあるのでそのまま使う。 // p == q(空きなし)なら葉から根までの区間列 ipath を歩き、 // 密度条件を満たす最小の区間を選ぶ。 (int, int) ipick<A extends Comparable<A>>(int n, Arr<A> as, (int, A) ix) { final (i, x) = ix; final (p, q) = ipoint(n, x, as); if (p < q) return (p, q); final path = ipath(n, p); for (var j = 0; j < path.length; j++) { final (l, r) = path[j]; final s = entries((l, r), as).length; if (densityTest(i, j, s, r - l, n)) return (l, r); } throw StateError('no interval satisfies density condition'); }

ipoint / ipath / intervals

ipoint :: Ord a ⇒ Int → a → Array a → (Int, Int) ipoint n x as = search (0, n) as where search (p, q) [] = (p, q) search (p, q) ((i, y) : as) = if x < y then (p, i) else search (i+1, q) as ipath n p = reverse (intervals (0, n) p) intervals (ℓ, r) p | ℓ + 1 == r = [(ℓ, r)] | p < m = (ℓ, r) : intervals (ℓ, m) p | m ≤ p = (ℓ, r) : intervals (m, r) p where m = (ℓ + r) div 2
Dart // ipoint: x を挟む隣接スロットの位置 (p, q) を返す。 // p == q なら「空きなし」(隣接する2つの値の間に x が入る)。 (int, int) ipoint<A extends Comparable<A>>(int n, A x, Arr<A> as) { var p = 0; var q = n; for (final (i, y) in as) { if (x.compareTo(y) < 0) return (p, i); p = i + 1; } return (p, q); } // 葉から根に向かって p を含む区間の列。 List<(int, int)> ipath(int n, int p) => intervals((0, n), p).reversed.toList(); // (0, n) から始めて、p を含む半分に降りていく区間の列。 List<(int, int)> intervals((int, int) lr, int p) { final result = <(int, int)>[]; var l = lr.$1; var r = lr.$2; while (true) { result.add((l, r)); if (l + 1 == r) return result; final m = (l + r) ~/ 2; if (p < m) { r = m; } else { l = m; } } }

密度テストは任意精度で

浮動小数点だと丸め誤差が怖いので、任意精度整数で計算します:

densityTest i' j' n s' w' = 2 ↑ (i * k) * s ↑ k < (2 ↑ i − 1) ↑ j * w ↑ k where (i, j, s, w) = convert toInteger (i', j', s', w') k = toInteger (ceiling (logBase 2 (fromIntegral n)))
Dart // 密度条件 s < ((2^i - 1)/2^i)^(j/k) * w を、任意精度整数(BigInt)で // s^k * 2^(i*k) < (2^i - 1)^j * w^k と等価変形して判定する。 bool densityTest(int iArg, int jArg, int s, int w, int n) { final k = (log(n) / log(2)).ceil(); final i = BigInt.from(iArg); final j = BigInt.from(jArg); final bs = BigInt.from(s); final bw = BigInt.from(w); final bk = BigInt.from(k); final two = BigInt.two; final lhs = two.pow((i * bk).toInt()) * bs.pow(k); final rhs = (two.pow(iArg) - BigInt.one).pow(jArg) * bw.pow(k); return lhs < rhs; }

convert f (a, b, c, d) = (f a, f b, f c, f d) です。

動作確認

Dart import 'dart:math'; void main() { // 特殊要素の位置: specials(11) = [5, 8, 9, 10] print(specials(11)); // "DANGEROUS" の 9 文字を大きさ 9 の棚に挿入する。 final n = 9; final xs = 'DANGEROUS'.split(''); final steps = insertAll<String>(n, xs); for (final step in steps) { final row = List<String>.filled(n, '-'); for (final (i, ch) in step) { row[i] = ch; } print(row.join(' ')); } // 期待される最終行: A D E G N O R S U }

おわりに

関連する問題

この問題は、より一般的なオンラインリストラベリング(Bender et al., 2002)の制限版です。オンラインリストラベリングでは、動的な集合(要素は挿入も削除もされる)から 0 〜 N−1 の整数ラベルへの写像を保ちながら、ラベルの順序を元の順序と一致させます。

用途としては順序保守(order maintenance)が代表的です。順序保守は、リストに対して insert(x, y)(y の直後に x を挿入)、delete(x)、query(x, y)(xy の前?)を扱う問題です。順序保守自体はラベルを必ずしも必要としませんが、多くの解ではオンラインラベリングを部品として使います。

この章の意義

見どころ このアルゴリズムの面白さは Haskell 実装そのものではなく、密度関数の巧みな選び方と、なぜ Θ(N log3 N) が達成できるのかの解析にあります。

Ω(N log3 N) が下界なのかは今も未解決ですが、Zhang (1993) は「空きスロットが常に均等になるように移動を制限するアルゴリズム」については、これが下界であることを証明しています。本章のアルゴリズムはこの意味で「なめらか」であり、なめらかでない方式でこれより少ない移動を達成できるとは考えにくいものの、可能性を排除するのは難しいのです。

この真珠は Bird and Sadnicki (2007) を元にしています。詳しい参考文献はそちらへ。

参考文献

Bender, M. A., Cole, R., Demaine, E. D., Frach-Colton, M. and Zito, J. (2002). Two simplified algorithms for maintaining order in a list. Lecture Notes in Computer Science, Volume 2461. Springer-Verlag, pp. 139–51.

Bird, R. S. and Sadnicki, S. (2007). Minimal on-line labelling. Information Processing Letters 101 (1), 41–5.

Zhang, J. (1993). Density control and on-line labeling problems. Technical Report 481 and PhD thesis, Computer Science Department, University of Rochester, New York, USA.