第30章

クモ紡ぎ入門

(やさしい版) Pearls of Functional Algorithm Design(関数プログラミングによるアルゴリズム設計の真珠)
「ああ、なんと絡み合った網を織ることか、はじめて導出を試みるとき」
(サー・ウォルター・スコットに謝意を込めて)

どんな問題?

グレイ経路とビット列

長さ n のビット列(0 と 1 が並んだもの)をたくさん作ります。ただし、 「i 番目のビットが 1 なら j 番目のビットも 1」のような不等式の制約 (aiaj)が付いています。この制約を満たすビット列だけを、 1回に1ビットだけ変えながら順に生成したい——これがこの章の課題です。

「1ビットずつ変えていく」順番のことをグレイ経路順(Gray path order)といいます。 どのビットを変えたかを並べたものを遷移コード(transition code)と呼びます。

具体例 n = 3 で制約が a1a2a3a2 のとき、たとえば 000 → 010 → 011 → 111 → 110 の順に並べられます。変えたビットの位置を並べると [2, 3, 1, 3] になります。

解けない場合もある

ちょっとやっかいなのは、そもそもグレイ経路順に並べられない場合もあることです。

たとえば n = 4 で a1a2a4 かつ a1a3a4 だと、可能なビット列は 0000, 0001, 0011, 0101, 0111, 1111 の 6 個ですが、これらをグレイ経路順に並べることはできません。 偶数個の 1 を持つ列が 4 個、奇数個の 1 を持つ列が 2 個あり、グレイ経路では 1 ビットずつ変わるので偶数個と奇数個が交互になるはずですが、 数が合わないからです。

クモの巣で考える

制約 aiaj は、ij という有向辺で 表せます。Knuth と Ruskey は、この有向グラフの向きを外した無向グラフが非巡回(totally acyclic、木構造)ならば、 グレイ経路が必ず作れることを示しました。連結な非巡回有向グラフを、彼らはクモ(spider)と呼びました。 上向きの辺と下向きの辺が入り混じった姿がクモに似ているからです。

Fig. 30.1 三本脚のクモ

複数のクモの集まりを(nest)と呼びます。

ポイント Knuth はこの問題を "spider squishing"(クモ潰し)と呼びましたが、著者は「クモは大事にしよう」ということで spider spinning(クモ紡ぎ)と呼び直しています。

色付けで考える

各ノードのビットが 1 のとき、0 のときに塗ると考えます。すると 「上位ノードが白なら下位ノードも白」というルールが自然に出てきます(aiaj は「j が白なら i も白」を意味するので)。

問題は、ある1つの色付けから始めて、毎回ちょうど1ノードだけ色を変えて、合法な色付けを全部たどることです。

この章の目標

目標は、この問題を解くループレスな(ループのない、1ステップが定数時間の)アルゴリズムを、 関数プログラミングの計算によって導くことです。Knuth 自身も自作のループレス版 SPIDERS プログラムを載せていますが、 彼はこう告白しています。

Knuth の告白(原文より) 「このアルゴリズムはかなり微妙で、それを門外漢(dummies)に説明する方法を、私はどうにも思いつかなかった」—— だからこの章のタイトルが「Spider spinning for dummies(門外漢のためのクモ紡ぎ)」なのです。

まずは簡単な場合:木クモによるクモ紡ぎ

木クモとは

クモがすべて木構造(脚がすべて下向き)である場合をまず考えます。この特殊ケースは Koda と Ruskey が調べています。データ型は次の通り。

type Nest = [Spider] data Spider = Node Int Nest
Dart // 木クモ: 根ラベル a と子クモの巣(Nest = List<Spider>) final class Spider { final int a; final List<Spider> nest; const Spider(this.a, this.nest); } typedef Nest = List<Spider>;

Fig. 30.2 二つの木クモからなる巣

クモの巣の遷移コード ncode と、単一クモの遷移コード scode は、往復積(boustrophedon product)関数 boxall を使って次のように書けます。

ncode :: Nest → [Int] ncode = boxall · map scode scode :: Spider → [Int] scode (Node a xs) = a : ncode xs
Dart // 往復積 as □ bs : bs で始め、as の各要素 a について [a] と bs / reverse(bs) を交互に挟む List<int> box(List<int> as, List<int> bs) { final out = <int>[...bs]; var flip = false; for (final a in as) { out.add(a); out.addAll(flip ? bs.reversed : bs); flip = !flip; } return out; } // boxall = foldr (□) [] (右畳み込みで巣全体を □ で連結) List<int> boxAll(Iterable<List<int>> xss) => xss.toList().reversed.fold(<int>[], (acc, xs) => box(xs, acc)); // ncode : 巣の遷移コード List<int> ncode(Nest nest) => boxAll(nest.map(scode)); // scode : 単一クモの遷移コード(根 a のあとに部分クモたちの遷移) List<int> scode(Spider s) => [s.a, ...ncode(s.nest)];
読み方のヒント scode (Node a xs) は「まず根ノード a の色を変える遷移、そのあと部分クモたちの遷移をつなげたもの」です。 短くて美しい定義ですが、まだループレスではありません

ループレス化への道

ループレス化は、次の順序で進めます。

  1. foldrmap fusionfold fusion(畳み込みを合体させる法則)で外側の boxall を消す。
  2. reverse(リストの反転)を、タプリング(結果を組で持ち回す)で消す。
  3. 薔薇木のキューを導入して、全体を実際に定数時間で動くようにする。

boxall = foldr (□) [ ] という事実と map fusion から、まず

ncode = foldr ((□) · scode) [ ]

が得られます。あとは (□) · scode を展開して整理していきます。細かい計算は原文の通りですが、結論として

ncode = foldr op [ ] op (Node a xs) bs = bs ++ [a] ++ foldr op (reverse bs) xs
Dart // 木クモ: fold fusion 後の ncode(まだループレスではないが、外側の □ は消えた) List<int> ncodeOp(Nest nest) => nest.reversed.fold(<int>[], (bs, spider) => op(spider, bs)); List<int> op(Spider node, List<int> bs) { final xs = node.nest; final inner = xs.reversed.fold( bs.reversed.toList(), (acc, s) => op(s, acc), ); return [...bs, node.a, ...inner]; }

となります。ここで reverse を消して薔薇木のキューを使う仕上げをすると、次のループレスな形になります。

ncode = unfoldr step · wrapQueue · fst · foldr op (empty, empty) op (Node a xs) (bs, sb) = (insert bs (Fork a cs), insert sc (Fork a sb)) where (cs, sc) = foldr op (sb, bs) xs
Dart 対訳について この定義は薔薇木のキュー(wrapQueue, Fork, empty, insert, unfoldr step) を前提としており、Dart への完全対訳は本記事の範囲を超えます(キュー自体をまるごと移植する必要があります)。 素朴版 ncodeOp と同じ遷移列を出すこと・毎ステップ定数時間で 1 個ずつ吐き出せることの 2 点だけ押さえておけば十分です。
重要な観察 最初の foldr op (empty, empty) は巣のサイズに対して線形時間で終わります。そして unfoldr step の各ステップは定数時間で 1 個ずつ遷移を吐き出します。だからこの全体がループレスだと言えるのです。

本題:一般のクモによるクモ紡ぎ

脚は上向きにも下向きにも

いよいよ一般のクモ紡ぎに進みます。クモの脚は上向きにも下向きにも他のクモを指せるようになります。 どのノードを「頭」にして持ち上げるかで木の形は変わりますが、表す制約は同じです。

Fig. 30.3 クモとそれを持ち上げてできる木

type Nest = [Spider] data Spider = Node Int [Leg] data Leg = Dn Spider | Up Spider
Dart // 一般クモ: 脚は Up / Dn を区別(sealed class + パターンマッチで表現) final class GSpider { final int a; final List<Leg> legs; const GSpider(this.a, this.legs); } sealed class Leg { final GSpider spider; const Leg(this.spider); } final class Up extends Leg { const Up(super.spider); } final class Dn extends Leg { const Dn(super.spider); } typedef GNest = List<GSpider>;

もう「辺」ではなく「脚」と呼び、Up(上向き)か Dn(下向き)かを区別します。

大事な前提 一般のクモでは、スタート地点が全 0 とは限らないのがやっかいです。たとえば a1a2a3 だと合法なビット列は 000, 001, 100, 101, 111 の 5 つで、 グレイ経路順に並べようとすると奇数重みの列(001, 100, 111)から始めるしかありません。 最初のビット列を返す関数 seed は後回しにします。

白コード・黒コード・そして "cox"

単一クモの遷移コード scode は、頭ノードがのときの遷移列(白コード wcode)と、 のときの遷移列(黒コード bcode)をつなげたものになります。

scode :: Spider → [Int] scode (Node a legs) = wcode legs ++ [a] ++ bcode legs

白コードが終わったところの色付けと、黒コードが始まるところの色付けをうまく揃えるために、□ の変種 ◇(読み方は「コックス」) を導入します。◇ は □ の共役です。

as ◇ bs = reverse ((reverse as) □ (reverse bs))
□ と ◇ の違い たとえば [2, 3, 4] □ [0, 1] = [0, 1, 2, 1, 0, 3, 0, 1, 4, 1, 0]([0,1] で始まる)ですが、 [2, 3, 4] ◇ [0, 1] = [1, 0, 2, 0, 1, 3, 1, 0, 4, 0, 1]([0,1] で終わる)です。

また、□ を ◇ で書き直す便利な等式があります(後で何度も使います)。

as □ bs = if even (length as) then as ◇ bs else as ◇ (reverse bs) (30.1)

脚ごとの白コード・黒コード

coxall = foldr (◇) [ ] と置き、次のように定義します。

wcode, bcode :: [Leg] → [Int] wcode = coxall · map wc bcode = boxall · map bc

脚ごとの白コード wc と黒コード bc は、色付けが正しく繋がるように次の定義に「強制」されます。

wc (Up (Node a legs)) = wcode legs ++ [a] ++ bcode legs wc (Dn (Node a legs)) = reverse (wcode legs) bc (Up (Node a legs)) = reverse (bcode legs) bc (Dn (Node a legs)) = wcode legs ++ [a] ++ bcode legs

簡単な計算で、巣全体の遷移コード ncode は下向き脚として扱った bcode で書けることがわかります。

ncode = bcode · map Dn

Fig. 30.4 ncode の出発点となるプログラム

ループレスなアルゴリズムを組み立てる

ステップ 1:fold fusion で外側を消す(ここが難所)

木クモのときと同じ道筋を辿ります。まず map fusion で次を得ます。

bcode = foldr ((□) · bc) [ ] wcode = foldr ((◇) · wc) [ ]

bcUp の場合、素直に fold fusion を狙うと、次の条件を満たす h が欲しくなります。

reverse (as □ bs) □ cs = h as ((reverse bs) □ cs) (30.2)
つまづきポイント 実は (30.2) を満たす h存在しません。理由は □ が単射(injective)でないからです。 たとえば "abab" □ "aaaba" = "ab" □ "aaabaaaba" のように、違う入力から同じ出力が出ることがあります。

そこで、最後の □ を ◇ に置き換えた次の形なら成立します。

reverse (as □ bs) ◇ cs = h as ((reverse bs) ◇ cs) (30.3)

◇ の結合律を使えば h as bs = reverse as ◇ bs と取れるとわかります。あとは (30.1) を使って □ ↔ ◇ を行き来しつつ、 もう一度 fold fusion を適用して、最終的に

bop (Up (Node a legs)) cs = reverse (foldr bop cs′ legs) where cs′ = if even (length (bcode legs)) then reverse cs else cs

を得ます。wop についても双対な議論で同様の式が得られます。

ステップ 1 のつづき:もう一つの節

bcDn の場合は、++ を □ に分配する規則(第一の真珠で導出した式)を使って計算します。

(xs ++ [y] ++ ys) □ zs = (xs □ zs) ++ [y] ++ (ys □ zs′) where zs′ = if even (length xs) then reverse zs else zs

結果、次の定義が得られます。

bop (Dn (Node a legs)) cs = foldr wop (reverse cs′) legs ++ [a] ++ foldr bop cs′ legs where cs′ = if even (length (wcode legs)) then reverse cs else cs

4 つの節すべてを整理すると Figure 30.5 のプログラムになります。

まとめ(現時点) □ と ◇ は消せました。でも見た目は美しくないし、効率も悪いです。原因は「パリティ(長さの偶奇)を何度も計算し直している」 こと。次のステップでこれを直します。

ステップ 2:パリティクモでパリティを埋め込む

パリティ(wcode legsbcode legs の長さの偶奇)をノードごとにあらかじめ計算して保持しておく ことにします。これがパリティクモ (Spider') です。

data Spider′ = Node′ (Bool, Bool) Int [Leg′] data Leg′ = Dn′ Spider′ | Up′ Spider′

不変条件は次の通り。

w = even (length (wcode legs)) b = even (length (bcode legs))

通常のクモにパリティを取り付ける関数を decorate と呼びます。

decorate :: Spider → Spider′ decorate (Node a legs) = node′ a (map (mapLeg decorate) legs) mapLeg f (Up x) = Up′ (f x) mapLeg f (Dn x) = Dn′ (f x)

スマートコンストラクタ node' は、子どもの脚のパリティ情報からこのノードのパリティを計算します。

node′ a legs = Node′ (foldr op (True, True) legs) a legs op (Up′ (Node′ (w, b) _ _)) (w′, b′) = ((w ≠ b) ∧ w′, b ∧ b′) op (Dn′ (Node′ (w, b) _ _)) (w′, b′) = (w ∧ w′, (w ≠ b) ∧ b′)
なぜこの op で正しいか wcode (leg : legs) の長さの偶奇は、wc leg の偶奇と wcode legs の偶奇の 両方が偶のときにだけ偶になります(as ◇ bs の長さは as の長さ と bs の長さの和だから)。
さらに wc (Up (Node a legs')) = wcode legs' ++ [a] ++ bcode legs' の長さが偶なのは、wcode legs'bcode legs' のパリティが逆のとき。だから w ≠ b が現れます。

これで、都度計算していた length がすべて事前計算済みのフラグに置き換わり、Figure 30.6 のようにスッキリした bopwop が得られます。

bop, wop :: Leg′ → [Int] → [Int] bop (Up′ (Node′ (w, b) a legs)) cs = reverse (foldr bop (revif b cs) legs) bop (Dn′ (Node′ (w, b) a legs)) cs = foldr wop (revif (not w) cs) legs ++ [a] ++ foldr bop (revif w cs) legs wop (Up′ (Node′ (w, b) a legs)) cs = foldr wop (revif b cs) legs ++ [a] ++ foldr bop (revif (not b) cs) legs wop (Dn′ (Node′ (w, b) a legs)) cs = reverse (foldr wop (revif w cs) legs) revif b cs = if b then reverse cs else cs

ステップ 3:reverse とキュー化

残るは reverse を消して、遅延評価が定数時間ずつ動くようにすることです。列 as を組 (as, reverse as) で持ち回すタプリングを使い、連結は

cat a (ws, sw) (bs, sb) = (ws ++ [a] ++ bs, sb ++ [a] ++ sw)

で実装します。reverse は「組を入れ替えるだけ」で終わります。さらに各成分を薔薇木のキューで表せば、Figure 30.7 の 最終ループレスプログラムが完成です。

ncode = unfoldr step · prolog prolog = wrapQueue · fst · foldr bop (empty, empty) · map (Dn′ · decorate) bop (Up′ (Node′ (w, b) a legs)) ps = swap (foldr bop (swapif b ps) legs) bop (Dn′ (Node′ (w, b) a legs)) ps = cat a (foldr wop (swapif (not w) ps) legs) (foldr bop (swapif w ps) legs) wop (Up′ (Node′ (w, b) a legs)) ps = cat a (foldr wop (swapif b ps) legs) (foldr bop (swapif (not b) ps) legs) wop (Dn′ (Node′ (w, b) a legs)) ps = swap (foldr wop (swapif w ps) legs) cat a (ws, sw) (bs, sb) = (insert ws (Fork a bs), insert sb (Fork a sw)) swap (xs, ys) = (ys, xs) swapif b (xs, ys) = if b then (ys, xs) else (xs, ys)

Fig. 30.7 最終的なループレスなプログラム

Dart 対訳の範囲について ここまでの一般クモ・パリティクモ・fold fusion の中間段階・ループレス最終版(薔薇木キュー + unfoldr step)・ seedbsp, wsp による初期/最終状態計算)は、いずれも複数の相互再帰と薔薇木キューを前提とした 純関数的コードで、実装量が大きくなるため本記事の Dart 対訳の範囲外とします。 木クモ版(Spider, box, boxAll, ncode, scode, ncodeOp) で、章の中核となる往復積と fold fusion の考え方は動く形で追えます。
まとめ プロローグ(前処理)はクモ・リスト・キュー・木、そして swap や fold といった役者が入り乱れる四幕劇になりましたが、 それでも巣のサイズに対して線形時間で完了します。以後は毎ステップ定数時間で 1 個ずつ遷移を吐き出せる、 正真正銘のループレスなクモ紡ぎです。

最後に:開始ビット列 seed の計算

状態としての色付け

ラベルが 1〜n のクモの巣に対して、グレイ経路の最初のビット列 a1a2an を返す関数 seed を作ります。 Haskell の Data.Map(有限写像)を使い、「クモの状態」を「ノードラベル → ビット」の写像として表します。

type State = Map.Map Int Bit type Bit = Int

名前衝突を避けるため、次のように別名を付けます。

install = Map.insert union = Map.union start = Map.empty

bseed / wseed で初期・最終状態を持つ

各コードに対して初期状態と最終状態のペアを計算していきます。bseedbcode 用、 wseedwcode 用。パリティ情報が必要なので、パリティクモ上で計算します。

seed = elems · fst · bseed · map (Dn′ · decorate) bseed = foldr bsp (start, start) · map bs wseed = foldr wsp (start, start) · map ws

bs は Figure 30.4 の bc の定義をなぞって、遷移ではなく状態を作ります。

bs (Up′ (Node′ (w, b) a legs)) = (b, install a 1 fs, install a 1 is) where (is, fs) = bseed legs bs (Dn′ (Node′ (w, b) a legs)) = (b, install a 0 is, install a 1 fs) where is = fst (wseed legs) fs = snd (bseed legs)
読み方のヒント bc (Up ...) = reverse (bcode legs) なので、bsUp では初期状態と最終状態が ひっくり返ります(install a 1 fsinstall a 1 is の順に注目)。
黒コードを見ているので、ラベル a のビットは 1。Dn のときは wcode legs ++ [a] ++ bcode legs なので、開始時は a が 0、終了時は 1 になります。

bsp と wsp で状態をまとめる

bspwsp は □ と ◇ の状態版です。□ が「始点は as、終点は as の偶奇次第で bsreverse bs」だったことを思い出すと、次の定義に納得できます。

bsp (b, ia, fa) (ib, fb) = (union ia ib, union fa (if b then fb else ib)) wsp (w, ia, fa) (ib, fb) = (union ia (if w then ib else fb), union fa fb)

Fig. 30.8 関数 seed

Dart // 章末の動作確認: 木クモ版 ncode / scode を素朴定義と fold fusion 版の両方で走らせる void main() { // 本文冒頭の例: n = 3, a1 ≤ a2, a3 ≤ a2 // ラベル 2 が根、1 と 3 が下の脚 final spider = const Spider(2, [Spider(1, []), Spider(3, [])]); print(scode(spider)); // => [2, 3, 1, 3] … 000→010→011→111→110 print(ncode([spider])); // => [2, 3, 1, 3] print(ncodeOp([spider])); // => [2, 3, 1, 3] (fold fusion 版も同じ結果) // Fig. 30.2 のような「二つの木クモからなる巣」 final nest = <Spider>[ const Spider(1, [Spider(2, [])]), const Spider(3, []), ]; print(ncode(nest)); // 巣全体の遷移コード print(ncodeOp(nest)); // 素朴版と一致 }

まとめ

最終的な感想 Knuth が「dummies に説明できない」と嘆いたループレスなクモ紡ぎ。関数プログラミングの計算則(fusion、タプリング、 スマートコンストラクタなど)を積み重ねると、こんなにも長い道のりですが、一歩ずつ機械的に導出できるのがこの章の見どころです。

1半順序集合 Sイデアルとは、S の部分集合 I で、 xI かつ xy ならば yI であるようなもののこと。

2ちなみに "to box and cox" は「交互にする」という意味で、これがまさに □ と ◇ の役割です。 この語はジョン・マディソン・モートンの喜劇 Box and Cox(1 幕の喜劇)に由来し、部屋を共有する二人の下宿人ボックスとコックスが 昼夜で交代して使っていた話にちなみます。

参考文献

Filliâtre, J.-C., and Pottier, F. (2003). Producing all ideals of a forest, functionally. Journal of Functional Programming 13 (5), 945–56.

Knuth, D. E. (2001). SPIDERS: a program downloadable from www-cs-faculty.stanford.edu/~knuth/programs.html.

Knuth, D. E. and Ruskey, F. (2003). Efficient coroutine generation of constrained Gray sequences (aka deconstructing coroutines). Object-Orientation to Formal Methods: Dedicated to The Memory of Ole-Johan Dahl. LNCS 2635. Springer-Verlag.

Koda, Y. and Ruskey, R. (1993). A Gray code for the ideals of a forest poset. Journal of Algorithms 15, 324–40.