2019-06-29

グレッグ・イーガン経由で知った順列生成アルゴリズム

SF作家のグレッグ・イーガンさんは最小超置換文字列の分散探索プロジェクトを主催しています。それ自体も興味深いんですが、解説ページで紹介されていた副産物の順列を列挙するアルゴリズムが面白かったので紹介&実装してみます。

Local Ruleのセクションで解説されている内容の紹介です。)

(追記:順列生成アルゴリズムをライブラリ化しました deltam/perm: Permutation generator based on group theory written in Go.

順列生成アルゴリズム

まず順列の順番を入れ替える写像sigma、deltaを次のように定義します。

sigma: (2,3,4,...,n-1,n,1)
delta: (3,4,5,...,n-1,n,2,1)

始点の順列を(n-2, n-3, ..., 3, 2, n, 1)として、それに対して次のルールAを繰り返し適用します。

ルールA
以下の条件をすべて満たすなら現在の順列にdeltaを適用し、そうでないならsigmaを適用する。

  1. 左端の数字がnではない。
  2. 順列中のnの右隣の数字をrとする(nが右端なら左から2番めをrとする)。左端の数字をfとし、r が 1+(f-2) mod (n-1) と等しい。
  3. 現在の順列が(1, n, n-1, ..., 3, 2)ではない。

(2番めの条件が分かりにくいので具体例を挙げると、順列(3, 1, 2, 4) ならf=3かつr=1で、1+(3-2) mod (4-1) = 1+1 = 2なので r = 1 ≠ 2となります)

始点からルールAを繰り返し適用した系列はすべての順列が一回だけ登場するものになっています。

実験

n=3の順列で実験してみます。まず始点は(2,3,1)です。

(2,3,1) : ルールA 1.OK 2.OK 3.OK -> delta
(1,3,2) : ルールA 1.OK 2.OK 3.NG -> sigma
(3,2,1) : ルールA 1.NG -> sigma
(2,1,3) : ルールA 1.OK 2.OK 3.OK -> delta
(3,1,2) : ルールA 1.NG -> sigma
(1,2,3) : 終了

全順列が出てきました。
4以上だと手計算が面倒なのでGoで書いてみました。

https://play.golang.org/p/hjt_iTAJWPo

$ go run perm_gen.go |head
[3 2 4 1]
[2 4 1 3]
[1 3 4 2]
[3 4 2 1]
[2 1 4 3]
[1 4 3 2]
[4 3 2 1]
[3 2 1 4]
[1 4 2 3]
[4 2 3 1]
$ go run perm_gen.go | wc -l
      24
$ go run perm_gen.go |sort|uniq|wc -l
      24

うまくいってますね。(4,3,2,5,1)から始めると5の全順列120個が出てくるので適当に書き換えて試してみてください。

イーガン先生が書いたJavaScriptによる実装も解説ページにリンクがあります(Supermutate.jsというやつ)。ページトップの文字がぐるぐる動いているのはそれで文字の全順列を表示しているみたいです。

なぜうまくいくのか?

グレッグ・イーガン先生は最小超置換の長さの上界を証明するためにAaron Williams氏の論文の結果を使っています。

[1307.2549] Hamiltonicity of the Cayley Digraph on the Symmetric Group Generated by σ = (1 2 … n) and τ = (1 2)

この論文では順列と2つの写像による群をケイリーグラフというグラフ構造に変換したものの上に、ハミルトン路を構成しそれが必ず存在するということを証明しています。
ハミルトン路とはグラフのすべての頂点をそれぞれ一回だけ通る経路のことです(終点→始点にならない場合もある)。頂点=順列なのでそのような経路を生成するアルゴリズムは順列生成アルゴリズムになります。上で紹介したアルゴリズムはイーガン先生が論文から少しアレンジしたものです。

上記の論文ではハミルトン路の存在以外にも次のことが証明されています。

  • 論文が対象としているケイリーグラフにはちょうど2つのサイクルカバーが存在し、すべての頂点はそのどちらか一方に含まれる
  • サイクルカバーを生成するアルゴリズムが存在する

サイクルカバーとはグラフ上の各頂点を一回だけ通るループのことです(すべての頂点を巡るものはハミルトンサイクルと呼ばれる)。この順列生成アルゴリズムはまず片方のサイクルカバーを生成し、始点に戻る直前の頂点でもう一方のサイクルカバーに移動して全頂点をめぐります。

サイクルカバーの次の頂点を生成するためのルールがルールAの1.2.です。3.は最初のサイクルカバーの終点を判別し、もう一方のサイクルカバーに移動するためのルールです。

普通の順列生成なら(1,2,3)から始めたいところですが、サイクルの終点を明確にするために(2,3,1)から始めてるんですね。ルールAを調整すれば(1,2,3)を始点に出来るかもしれないですが、出力するときに変換をかますほうが楽かも。

ケイリーグラフ上のハミルトン路の存在証明は上記の論文が初めてではないようですが、すでに知られているものよりかなりシンプルになったそうです。Goで書いたものも70行ほどなので素晴らしいですね。

最小超置換との関わり

元論文ではdeltaは先頭二文字を入れ替える(2,1,3,4,...,n-1,n)写像になっています。イーガン先生がさらに二文字を後ろに持っていくようにしたのは超置換文字列を作るためです。
超置換はnの順列がすべて部分列として含まれる文字列のことです。sigma、deltaで順列を変換すると前の順列と先頭n-1、n-2文字が共通になります。変換するごとに1,2文字を追加していけばそれは超置換文字列になります。

上記のGoプログラムのmainをこんなふうに書き換えれば超置換を出力します。

func main() {
 p := []int{2, 3, 1}
 for i := 0; i < len(p); i++ {
  fmt.Print(p[i])
 }
 add := 0
 for i := 0; i < 6-1; i++ {
  if ruleCheck(p) {
   delta(p)
   add = 1
  } else {
   sigma(p)
   add = 2
  }
  for i := len(p) - add - 1; i < len(p); i++ {
   fmt.Print(p[i])
  }
 }
}
// Output: 2313232121312123

この方法で作った超置換は部分列に重複が多く、残念ながら最小にはなりません。n=3の最小超置換は123121321です。
ただしこれによって数学的に証明された超置換を生成するアルゴリズムを作ることができて、超置換の長さの上界を示すことができます。

最小超置換問題の全体的な説明はこちらの記事が参考になります。イーガン先生の貢献も解説されています。
ハルヒ問題(最小超置換問題)の紹介 - Qiita

(おまけ)イーガン先生による最小超置換の分散探索プロジェクト紹介

Distributed Chaffin Method Results

いまアマチュア数学者業界で最小超置換が熱い! 単純な予想を覆す複雑さを持ち、なんか新発見できるんではという希望があります。
イーガン先生は現在n=6の超置換を全探索してすでに知られたものが本当に最小か追試しています。

最小超置換の探索アルゴリズムのChaffin Method(これも考え方が興味深いので解説書きたい)を複数台のコンピュータで実行できるようにしたのがDistributed Chaffin Methodです。サーバもクライアントもイーガン先生が書いてます。

superperm/DistributedChaffinMethod at master · superpermutators/superperm

クライアントはCとJavaのものがあります。チーム名を決めて実行するとプロジェクトページに実績が載るので誘い合わせて参加するといいかも(私はぼっちチームですが…)。

現在のところ最小超置換を見つけることの実用的意味はまったくありませんが、すこしでも面白そうと思ったならクライアントを動かしてみませんか?

イーガン先生の熱情がどこから来るのか不明ですが、このような著作があるのと関係あるのかもしれません。(もしくは短編『ルミナス』のように数論の<不備>を発見しようとしているのかも?)




コード全体

// permutation generator
// author @deltam

package main

import "fmt"

func main() {
 p := []int{3, 2, 4, 1}
 for i := 0; i < 24; i++ {
  fmt.Println(p)
  if ruleCheck(p) {
   delta(p)
  } else {
   sigma(p)
  }
 }
}

// 2 3 4 ... n-1 n 1
func sigma(p []int) {
 n := len(p)
 f := p[0]
 copy(p[0:n-1], p[1:n])
 p[n-1] = f
}

// 3 4 5 ... n-1 n 2 1
func delta(p []int) {
 n := len(p)
 f0 := p[0]
 f1 := p[1]
 copy(p[0:n-2], p[2:n])
 p[n-2] = f1
 p[n-1] = f0
}

func ruleCheck(p []int) bool {
 n := len(p)
 // 1.
 if p[0] == n {
  return false
 }
 pos := 0
 for i := 0; i < n; i++ {
  if p[i] == n {
   pos = i
   break
  }
 }
 pos++
 if pos == n {
  pos = 1
 }
 // 2.
 if p[pos] != 1+(p[0]-2+n-1)%(n-1) {
  return false
 }
 // 3.
 for i := 0; i < n; i++ {
  if p[(i+1)%n] != n-i {
   return true
  }
 }
 return false
}

おわり。


(訂正 2019-06-30:ハミルトンサイクルとサイクルカバーの用語の意味を取り違えてたので修正)

2018-10-10

カスタマイズしやすい重み付きレーベンシュタイン距離ライブラリをGoで書きました

deltam/go-lsd-parametrized: Weighted Leveshtein Distance and its extended interfaces written in Go.

重み付きレーベンシュタイン距離(編集距離)のライブラリを書きました。文字ごとに追加・削除・置換のコストが指定できます。文字列の長さによる正規化や複数の文字列から近いものを選ぶ関数も含まれてます。

エラーメッセージをラフに分類したいという動機から作り始めました。類型的で機械学習が必須なわけでない分類タスクにお役立てください。
詳しい使い方はGitHubのREADMEを見てください。

go get -u github.com/deltam/go-lsd-parametrized

カスタマイズ


重みの調整でだいたいの用途は足りるかなーと思います。編集方法ごとに重みを付ける場合はlsdp.Weights構造体、さらに文字ごとに重みを付けたい場合はlsdp.ByRune()を使います。

    // weighted
    wd := lsdp.Weights{Insert: 0.1, Delete: 1, Replace: 0.01}
    fmt.Println(wd.Distance(a, b))
    // Output:
    // 0.22

    // weighted by rune
    wr := lsdp.ByRune(&lsdp.Weights{1, 1, 1}).
        Insert("g", 0.1).
        Insert("h", 0.01).
        Replace("k", "s", 0.001).
        Replace("e", "i", 0.0001)
    fmt.Println(wr.Distance(a, b))
    // Output:
    // 0.1111

調整するときに分類精度を評価するのが面倒なのでtools.Evaluate()というのも作りました。tools.EvaluateCSV()ならCSVでテストデータを用意できるので便利。

もっと細かく調整したい場合はlsdp.DistanceMeasurer インターフェイスを実装する必要があります。
ただロジック的な変更のみでパラメータが無い場合はいちいち構造体を定義するのも面倒なので、こちらの記事を参考に直接関数で指定できるようにしてみました。

Goのメソッドは構造体以外にでも定義できるしそれが便利なこともよくある - Qiita

// 長さの差を距離とする
func lenDiff(a, b string) float64 {
    d := utf8.RuneCountInString(a) - utf8.RuneCountInString(b)
    return math.Abs(float64(d))
}

func main() {
    var d lsdp.DistanceFunc = lenDiff
    fmt.Println(d.Distance("kitten", "shitting"))
    // Output:
    // 2

    group := []string{"", "a", "ab", "abc"}
    s, dist := lsdp.Nearest(d, "xx", group)
    fmt.Println(s, dist)
    // Output:
    // ab 0
}


Goで書いた感想


Goでライブラリを書く経験はそれほどなかったのでいろいろ勉強になりました。

interfaceの使い方


Effective Goでこんな一節があったのでNormalized()の実装で実践してみました。

Generality
If a type exists only to implement an interface and will never have exported methods beyond that interface, there is no need to export the type itself.

Effective Go - The Go Programming Language

構造体を作ってInterfaceとして返すようにしておけば実装の詳細を隠せるのであとで変更しやすいです。Go moduleだと後方互換性が大事なのでこうしておくとのちのち楽です。

ベンチマークテストたのしい


Goの標準ライブラリのコミットログを見てるとベンチマークテストの結果が貼り付けられたものがあります(例:string)。効果が分かりやすいし、なんかかっこいいので真似することにしました。

Big Sky :: golang でパフォーマンスチューニングする際に気を付けるべきこと

単純な比較ならbenchcmpというのもあります。入力をランダム生成して複数回試行するようなベンチマークならbenchstatのほうが良さそう。


アルゴリズム的なやつ


コスト計算を工夫してメモリ削減


レーベンシュタイン距離のアルゴリズム(解説)では文字列1、2の先頭からの部分列に対して段階的にコストを計算する方法を使います。空行から文字列全体なので文字列1,2がそれぞれN文字、M文字ならば N+1行 M+1列の表を作ります。

詳しい説明は省きますが、現在の行を計算するのに必要なのは一つ前の行の値だけなので実際に必要な記憶領域は(N+1)*2だけです。レーベンシュタイン距離の実装を見てみるとN+1の配列2つを互い違いに使うことで表の代用をする方法が多いです。

ベンチマークテストを追加した影響でもっと効率化できないか考えるようになりメモリ削減のアイデアを思いつきました。表のあるコスト(i, j)の計算に必要なコストは左(i-1, j)、左上(i-1, j-1)、上(i, j-1)の3つだけです。コスト(i, j)の計算が終わったら次は(i+1, j)に移るので、実は左上(i-1, j-1)の値はもう参照しません。だったら左(i-1, j)の値で上書きしちゃえばいいのでは?というが思いついたアイデア。これによって記憶領域はN+2で済みます。

たぶん既出アイデアだと思うけどそこそこ高速化できたので気分がいい。

Performance improvement for accumulateCost() · deltam/go-lsd-parametrized@61be4fe

ベンチマークテスト万歳!

重み付き〜の定義


最初かんぜんに勘違いしてた、はずかしい…

アルゴリズム - 重み付きレーベンシュタイン距離の定義について - スタック・オーバーフロー

勘違いの副産物で削除・挿入・置換ごとの最短回数を集計する関数ができました(通常のレーベンシュタイン距離はコストは全部合計するため編集種別ごとの回数は分からない)。
lsdp.CountEdit()





なんかあったらissue、PullReq投げてくださいー。

おわり。

2018-03-13

FitDesk X-2.0を修理した

FitDesk X-2.0というエアロバイクを愛用してます。



デスクが付いているので運動しながら読書をしたりPC作業をしたり便利です。
もう2年以上ほぼ毎日使っているのですが、すこし前からペダリングに違和感を覚えるようになりました。音がうるさくなりペダル回転も引っかかるような感触がします。

エアロバイクの故障について調べてみたら故障の原因はだいたいがベルトかベアリングの破損のようです。私のFitDeskの保証期間はすでに切れていたので何とか自分で修理できないか試すことにしました(安くない買い物だったし)。

検索してみるとFitDesk製造元のYoutubeチャンネルにギアボックスを分解するための手順の解説動画があります。








どうやらクランクを外さないとギアボックスは開けない、クランクを外すには自転車の専用工具がいるらしい。
しょうがないので必要な工具を用意しました。


修理に必要な工具リスト

  • クランク抜き工具
    • この商品はFitDeskのクランクに合いました(共通規格なのかな?)

  • ゴム製ハンマー

  • ソケットレンチ
    • クランク抜きに必要なのはサイズ13のみ
    • 固く締まっていたので力が込めても痛くなさそうな柄のやつがいいです

  • モンキーレンチ
    • ペダルも外すのに力がいるので柄が長いほうがいいです

  • プラスドライバー
    • カバーを外すのに必要




ギアボックスを開く


クランクを外す→前の接地部を外す→カバーを外すの手順でいきます(バランスが悪いのでカバーを外した後に接地部を戻して軽く固定してます)。




ギアボックスを開くと中に黒い金属の玉が落ちてました。ベアリングが壊れて中のボールが落ちてきたみたいです。



実際ギアボックスも分解してホイールやベルトを取り出してみると、ベルトを裏から抑えているベアリングがこんなふうに壊れてました。



他の部分に異常はなさそうなのでこのベアリングを交換すればOKだと判断しました。



ベアリングを交換して修理完了


保証期間過ぎてるけどダメ元で販売代理店の株式会社PODに連絡してみると送料のみで譲っていただけることに。感謝です!


さっそく交換して元通りになりました。


自前でベアリングを用意して交換する場合


ベアリングは消耗品らしいので今度また壊れたらこんなふうに直そうと思ってます。
  • サイズを調べて同じベアリングを購入して差し替える
    • マニュアルの部品リストには「Bearing 6200ZZ」と書いてあるのでこれで大丈夫そう

    • ボルトにローラーベアリングがくっついた構造なので取り外すのにまた専用工具がいる……



故障を予防するために


マニュアルには「30分毎に20分の休憩を」と書いてありますが、健康上の理由だけでなくベアリングやベルトが熱を持ちすぎないようにとの理由があるそうです。ダイエットのためとしても連続で漕ぎ続けるのは膝にもエアロバイクにも良くないみたいですね。




【参考リンク】

プログラミングしながら運動できるエアロバイク FitDesk X-2.0を買った - orangain flavor
 私が購入したきっかけはこの記事。しかし振動があるので入力作業はあまりしてません。難しいドキュメントやソースコードを読むときなどは眠くならなくて超捗る。

2015-08-27

SICP問題1.45の回答と説明(SICP1.3.4, ex1.45)

この記事の対象読者:『計算機プログラムの構造と解釈』を読んだことがあって、問題1.45に納得がいってないひと



1.3.4 値として返される手続き (計算機プログラムの構造と解釈 第二版)

SICPの読んでたときにけっこう苦労した練習問題。
n乗根を求めるための平均緩和法の適用回数を「実験して確かめよ」って書いてあるけど、たぶん実験だけでは分からない。回答自体はカンニング(検索)して見つけたけど、なぜそうなのかをちゃんと説明しているのを見つけられなかったので、この記事で説明してみる。


平均緩和法、不動点探索の復習


回答のまえに平均緩和法や不動点探索についておさらいしておきます。

不動点探索(fixed-point)ではx=f(x)となる点をみつけます。y=f(x)を計算してさらにf(y)を計算して、というのを繰り返す。図解すると次の赤線のような軌道を描きます。



ところがこの方法では上手くいかない場合がある。次のようにy=xを軸として対称な曲線だと同じ軌道をグルグルしてしまい不動点に行き着かない。



このようなときに平均緩和法を適用するといい具合に非対称に歪めてくれる。この変形のおかげで不動点探索がうまくいくようになります。





四乗根を平均緩和法で求める


平均緩和の不動点として立方根を見つけることが出来る. 困ったことに, この方法は四乗根には使えない.

問題文にはこう書いてあるけど、実際やってみると四乗根を見つけることはできます。見つけられないのは(上で説明したとおり)y=xの軸に対称な場合ですが、平均緩和法を適用すれば非対称になるはずなので見つけられるということです。
でもこれをコーディングして実行すると三乗根に比べてだいぶ遅いです。

;;; 安全に計算の動きをおうため、計算回数の上限を指定できるように改造した
 (define (fixed-point-limited f first-guess limit)
   (define (close-enough? v1 v2)
     (< (abs (- v1 v2)) tolerance))
   (define (try guess count)
     (let ((next (f guess)))
       ;; 計算過程を表示する
       (display count) (display ": ") (display next) (newline)
       (if (or (= count limit) (close-enough? guess next))
           (list count next) ; 計算回数も一緒に返す
           (try next (+ count 1)))))
   (try first-guess 1))

 ;;; 平均緩和法1回での4乗根計算手続き
 (define (root4-ad1 x limit)
   (fixed-point-limited
     (average-damp (lambda (y) (/ x
                                  (* y y y))))
    1.0
    limit))

 (root4-ad1 16 10)
 ; (10 1.8597959434977196)

 (root4-ad1 16 100)
 ; (100 1.9364452032204889)
 (root4-ad1 16 1000)
 ; (1000 1.9781656271972445)

 ;(time (root4-ad1 81 0))
 ; real 176.408
 ; user 175.530
 ; sys    0.380
 ; (449999982 3.000050001249936)

81の四乗根を求めるのに4億回も計算して176秒も掛かっている。



なぜ四乗根の計算は遅い?


四乗根の不動点探索の様子を見てみましょう。平均緩和法を一回適用したグラフはこのようになります。その上での不動点を求める軌道はこのように螺旋を描きます。



比べて、平均緩和法を2回適用した場合は、このようにまっすぐに不動点に向かう軌道になります。



直感的に上より下のほうが早く不動点に行き着きそうですよね。


不動点探索に向いた曲線とその求め方


不動点探索に向いた曲線、早く探索が収束するための曲線ってどんなのでしょう。上のふたつのグラフを比較してみると、グラフの最下点の位置が重要そうに思えます。
最下点が不動点の位置よりだと、探索軌道が螺旋を描いてしまう。最下点が不動点よりだとまっすぐな探索軌道を描く。
最下点と不動点のx軸上の位置関係が、最下点<=不動点となっていると良さそうです。


ということで平均緩和法を複数回適用した場合に、その最下点がどうなるか調べてみれば良さそうですね。
調べやすくするために平均緩和法を適用した数式を変形してみます。

・準備
定数mのn乗根を求めるための不動点探索用関数を以下のように定義する。
\[g_{n}(x)=\frac{m}{x^{n-1}}\]

$g_{n}(x)$に平均緩和法をk回適用した関数をつぎのように定義する
\[
\begin{eqnarray}
f_{n,k}(x) & = & \frac{1}{2} (x+ \frac{1}{2} (x+\ldots \frac{1}{2} (x+ g_n(x) )\ldots)) \nonumber \\
& = & \frac{1}{2} (x+\frac{1}{2} (x+ \ldots \frac{1}{2} (x+ \frac{m}{x^{n-1}} ) \ldots )) \nonumber \\
& = & (\frac{1}{2} + \frac{1}{2^2}+ \frac{1}{2^3}+ \ldots + \frac{1}{2^k})x + \frac{1}{2^k} \frac{m}{x^{n-1}} \nonumber \\
& = & \frac{2^k-1}{2^k}x + \frac{1}{2^k} \frac{m}{x^{n-1}} \nonumber \\
& = & \frac{1}{2^k}((2^k-1)x+ \frac{m}{x^{n-1}})
\end{eqnarray}
\]

$f_{n,k}$の最下点を調べたい。そのため式(1)を微分して0になる点を調べる。
\[
Df_{n,k}(x)= \frac{1}{2^k}(2^k-1+(1-n) \frac{m}{x^n})
\]
\[
\begin{eqnarray}
Df_{n,k}(x) & = & 0 \nonumber \\
\frac{1}{2^k}(2^k-1+(1-n) \frac{m}{x^n}) & = & 0 \nonumber \\
2^k-1+(1-n) \frac{m}{x^n} & = & 0 \nonumber \\
\frac{m}{x^n} & = & \frac{1-2^k}{1-n}
\end{eqnarray}
\]

さて$f_{n,k}$の不動点は$\sqrt[n]{m}$である。不動点にあるとき式(2)の左辺は1になる。
\[
\begin{eqnarray}
1 & = & \frac{1-2^k}{1-n} \nonumber \\
1-n & = & 1-2^k \nonumber \\
n & = & 2^k \nonumber \\
\log_2 n & = & \log_2 2^k \nonumber \\
& = & k \log_2 2 = k \nonumber \\
\log_2 n & = & k
\end{eqnarray}
\]


回答


ようやく求める平均緩和回数がわかりそうだ。式(3)の k は整数にしないと回数にならないので、左辺を切り上げるか切り捨てるかしないといけない。

  • 切り上げ:式(2)の右辺は1より大きくなる。左辺のxは不動点より小さい値でなくてはならない
  • 切り捨て:式(2)の右辺は1より小さくなる。左辺のxは不動点より大きい値でなくてはならない

不動点と最下点xの位置関係は「最下点x<=不動点」が良いってことはグラフを描いて分かってました。なので式(3)は切り上げのほうが良さそうです。
まとめると平均緩和法の適用回数kはつぎの式になります。

\[ k= \lceil \log_{2}{n} \rceil \]



;;; n乗根を探す手続き
 (define (root-n x n)
   (fixed-point
    ((repeated average-damp
                 ;; 回数:2を底とするlog nの小数点以下切上げ
                 (ceiling (/ (log n) (log 2))))
     (lambda (y) (/ x
                    (expt y (- n 1))))) ; exptはべき乗
    1.0))

 ;; test
 (root-n 2 2)
 ; 1.4142135623746899
 (root-n 27 3)
 ; 3.0000234260125787
 (root-n 256 8)
 ; 2.0000000000039666
 (root-n 65536 16)
 ; 2.000000000076957

参考リンク


切り捨てが多いけど理由説明は特にないんですよね。上のグラフの議論から切り上げのほうが良いと私は考えてます。

2015-05-26

『コンピュータシステムの理論と実装』第5章 CPUを作ったよ




個人的に進めてる『コンピュータシステムの理論と実装』の実装メモを書きますよ。
今回はCPUを作るのが山場でした。

「第5章 コンピュータアーキテクチャ」ではHackコンピュータというこの本オリジナルのコンピュータを作ります。前章まででALU、レジスタやプログラムカウンタなんかは作ってあるので、それらを組合せてなんとかします。

命令メモリ、CPU、データメモリを合体させるとHackコンピュータになるらしい(P.105、5.3.3)。メモリは前章のRAMを利用すればすぐ作れるので、問題はCPUを作ることです。


P.102に概略図が載ってるのでそのとおりに作るんですが、詳細が省かれてるので実際にHDLで書こうとすると大変でした。この本は実践に重きを置いているので、詳細を省いているのはたぶん意図的です。実装を通じて学んで欲しいということでしょう。

下に私の考えた実装コードを載せますけど、本の思想に忠実に自分で考えたいというひとはネタバレ注意です。(GW中にやったけどハマった部分があって3日ぐらいかかってしまったw)










HDLコード
nand2tetris/CPU.hdl at master · deltam/nand2tetris




苦労したところはジャンプ命令のパースでした。ぜんぶ論理回路でやらないといけないので慣れないとどこから手を付ければいいのか分からないです。
ジャンプ命令パース部分の論理回路を部分的に書いてみました。



まず命令コードのジャンプ命令部分をデマルチプレクサのselに入力して何の命令なのかはっきりさせます(該当するジャンプ命令のビットだけ1が立つ)。
ALUなどの論理回路を作ったときの経験から学んだコツとして"必要な値は最初に全部作っておいて、必要なときに組み合わせる"というのがあります。
ジャンプ命令はALUの出力値の大小によって分岐するので、条件判断に必要な次の値とその否定を全部作っておきます。

  • zr − ALUの出力が0と等しい
  • ng − ALUの出力が0未満
  • gt − ALUの出力が0より大きい(zrとngから作る)

つぎにそれぞれのジャンプ命令に必要な論理回路を作って条件判断させておきます。
最後にデマルチプレクサの出力値とAndして(該当命令の条件判断だけが残る)、全体をORでまとめればジャンプ命令がパースできたことになります。

このロジックというか、考え方を思いつくのに苦労しました。CPU作るのって大変ですねー。



苦労したけど、エミュレータでサンプルコードがちゃんと動くのを確認できたときはかなりの達成感でした。「やったぜ! 俺はコンピュータを作ったぞ!」www


次の第6章はこのHackコンピュータ用のアセンブラを実装します。ここから先はソフトウェアの世界になるので少し楽になるんじゃないかなーと思ってます(コンパイラ、OSの実装とかあるけど)。

2015-04-19

AND、ORだけではNOT回路は作れないのか

「コンピュータシステムの理論と実装」読書会のときの疑問を追っかけて考えてみました。

NotとAnd/Orからならすべての論理回路が作れることは分かりました。ではAnd、OrのみからNot回路は作れないのでしょうか?
直感的に「作れない」って感じるんですが、きちんと論理的に説明するのはけっこう難しい気がします。

ではどうやったら説明できるか。

Notは1入力1出力の回路です。まずは入力をaとおきます。最初の入力の可能性としてはa、定数0、定数1のどれかしかありません。
最初の入力から次の出力の可能性は次の表の通りです。やっぱりa,0,1のどれかになってしまいます。

入力1 入力2 And出力 Or出力
a 0 0 a
a 1 a 1
a a a a
0 0 0 0
0 1 0 1
1 1 1 1


つぎに1出力ってことを考えます。And、Orのみで回路を作るならケーブルをまとめる最後の素子はAndまたはOrのどっちかになります。そして回路の最初の入力も次の出力もa,0,1になってるので、最後の出力もやっぱりa,0,1になってしまう。
雑に図解するとこんな感じ。


途中の回路がどんな複雑でもAnd/Orのみなら中間出力も変わらないので、最後の論理素子の入力も変わらない。だから最後までNot aを作ることはできない。


……という説明を考えてみたけど、なんだかすっきりしない。背理法か何かでさっぱりと証明する方法はないのかなー。




追記:ツイートしたら@chiralさんが速攻で返信くれました。だいぶすっきりした!

2015-04-17

#nandris コンピュータシステムの理論と実装 読書会(2) に参加してきました

「コンピュータシステムの理論と実装」は論理回路から順を追ってコンピュータの仕組みを知ろうという教科書です。
Nand回路から初めて論理回路、加算器から機械語、OSや高級言語のコンパイラを経てテトリスを作っちゃおうという、デアゴスティーニでやりそうな構成。

この本の内容はNand2Tetrisというオンライン学習コースで、私はこのTED Talksで知りました。

独習コンピュータ教程



Nandからテトリスまでつなぎ目無く作っていくというのは、なんとなく冒険心を刺激されるものがありますね。
学習コースでは資料全部、ハードウェアシミュレータのようなツールも全部公開されてるので本当に独習できるんですが、全部英語なのが難点でした。
いつかやろうと思って放置してましたが日本語訳が出たというなら、今でしょッ! ということでちょうどよく開催されていた読書会に参加!

コンピュータシステムの理論と実装 読書会(2) - connpass



第2回は主に論理回路とブール代数でした。ここらへんは情報系工業高校に通った身なので死ぬほどやったところ。すごく懐かしかったです。
論理回路のキモはNand回路が一種類があればAnd、Or、Notがすべて構成できるってこと。さらにそれらを使ってつくるマルチプレクサなども作れます。

書いたHDLをハードウェアシミュレータで動かしてみることもできます。ハードウェアシミュレータ+今後必要なツールはここからダウンロードできます。

The Elements of Computing Systems / Nisan & Schocken



NandかNorから他の論理素子を作れることは分かったけど、別のものを元に論理素子を作ることは出来ないのかな? ということでマルチプレクサからAndとNotを作ってみました。
Andは $out = b \cdot sel$。

Notは $out = \overline{sel} $。

HDLのコードで書くとこんな感じ(HDLでは1,0はtrue,falseで表すらしい)。

CHIP Not {
    IN in;
    OUT out;

    PARTS:
    Mux (a=true, b=false, sel=in, out=out);
}

CHIP And {
    IN a, b;
    OUT out;

    PARTS:
    // Put your code here:
    Mux (a=false, b=a, sel=b, out=out);
}

元になる論理素子からNotが作れればなんとかなりそうですね。
And、OrだけじゃNotは作れないだろうけど、そのことってどうやって証明するんだろうなー。




ブール算術の章は2の補数のところまでで今回は終了。2の補数は2進数でマイナスの数を表す約束事ですが、これを採用すれば足し算回路だけで引き算も表現できるので超重要です。使わないと引き算回路も別に作らないといけなくなって不経済。
1の補数というのもあるのですが、これだと−0と+0が出来てしまって効率的じゃないのです。途中で@yshigeruさんが紹介してくれた式変形が、1の補数から2の補数への変化をわかりやすく表現していて良かったです。

\[ a + \overline{a} = -1 \]
\[ \overline{a} + 1 = -a \]



次回はちょっと長めですがそのぶんゆっくりできるでしょう。興味ある人はどうぞ。
コンピュータシステムの理論と実装 読書会(3) 休日編 - connpass







【おまけ】Nandのみで作ったAnd,Or,NotのHDLコード
NandOnly.hdl