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