2016年9月22日木曜日

ループ検出

数列xiが例えば32ビットの語の範囲の値しかとれなくて, xi+1=f(xi)を繰り返していくと, いつかは昔のある値が現れ, 以後同じ数列を辿ることになって, ループを回り始める.

このループの開始位置とループ長を求めるのがループ検出である.

そのひとつの方法にBill Gosperの考案したものがあり, あちこちに説明があるが, あまり真面目に読んだことはなかった. しかし気にはなっていたので, この数日プログラムを書いてみたりしながら, その仕掛けを調べてみた.

いまx1, x2, x3, ... が1,2,3,4,5,6,7,8,9.10,6,7,8,9,10,6,7,8,...であったとしよう.(x0から始める流儀もあり, また私はそれが好きだが, 今回の添字は1から始めることにした.) The Art of Computer Programming(TAOCP)の演習問題3.1,6の記述を借用すれば, 「x1,x2,...,xμ,..., xμ+λ-1は互いに相異なるが, n ≥μの時はxn=xnであるμとλがある.」

このμが開始位置, λがループ長で, 上の例ではμ=6, λ=5だ.

記憶場所も時間もふんだんにある時は, xの値をすべて順に記憶しておき, 新しいxについて, それを過去のすべてのxの値のそれぞれと比較する. そこに同じものがあれば, 前の値のあったところがループの開始場所であり, そこから今の値までの長さがループ長である.

しかし, こんなことはとてもやっていられないから, あの手この手を考えることになる. その一つにGosperの提案したものがある. 以下はそれを私流に多少修正したものである.

これは簡単にいえば, 繰り返しが「現れるや否や」発見するのは諦める, 最新のxの1つ前, 2つ前, 4つ前, 8つ前, ... などを覚えておき, ループに入り同じ値が現れ始めたらなるべく早目にそれを発見しようという方針である.

「1つ, 2つ, 4つ, 8つ,...」から推察されるように, n個のxを調べた時は, log2 n個程度の配列Tを使う.

ところで, n>0を二進法で表した時, もっとも右の1の位置, もっとも左の1の位置をTAOCPに倣い, それぞれρ(n), &lambda(n)とする. ビットの位置は右端を0とする. (上のループ長の定義にもλがあって申し訳ない.)

ρ(1)=0, λ(1)=0; 
ρ(2)=1, λ(2)=1; 
ρ(3)=0, λ(3)=1; 
ρ(4)=2, λ(4)=2; 
ρ(5)=0, λ(5)=2; 
ρ(6)=1, λ(6)=2; 
ρ(7)=0, λ(7)=2; 
ρ(8)=3, λ(8)=3; 
ρの方は別名ルーラー関数(物差関数)という.



十進法の物差しとは様子が違い, ひと目盛りが半分ずつに分割されていて, インチの物差しが大体こういう感じである.

ρとλを計算する関数 rhoとlamは次のようだ.

ρは
(define (rho x)
 (define (r x rh)
   (if (= (modulo x 2) 0) (r (quotient x 2) (+ rh 1)) rh))
 (if (= x 0) 'error (r x 0)))

(map rho (a2b 1 16)) => (0 1 0 2 0 1 0 3 0 1 0 2 0 1 0)
またλは関数logがあれば
(define (lam n)
 (inexact->exact (floor (/ (log n) (log 2))))

(map lam (a2b 1 16))=>
(0 1 1 2 2 2 2 3 3 3 3 3 3 3 3)
ρのように自力で計算するには
(define (lam x)
 (cond ((= x 0) 'error)
       ((= x 1) 0)
       (else (+ (lam (quotient x 2)) 1))))
さて配列をT0, T1, T2, ... とする. (こちらは0から始める) そしてxiTρ(i)に順に置くのである.

xiが次々と置かれていく様子は次のようだ.

iT0 T1 T2 T3
1x1
2x1 x2
3x3 x2
4x3 x2x4
5x5 x2x4
6x5 x6x4
7x7 x6x4
8x7 x6x4 x8
9x9 x6x4 x8
10x7 x10x4 x8
注意すべきは, この表のうち見えているのは最後の行だけということである. つまりi=8の時には, x1, x2, x3, x5はすでに消えている.

さて, 上の例のようなループがあった時の同様な配列のx10=10までの変化の様子は,
iT0 T1 T2 T3
1 1
2 1  2
3 3  2
4 3  2 4
5 5  2 4
6 5  6 4
7 7  6 4
8 7  6 4  8
9 9  6 4  8
10 9 10 4  8
のようだ. ここでx11=6になった時, この配列は9, 10, 4, 8であり, この中に6はないから, 6を(ρ(11)=0だから)T0のところに置く. したがって配列は6, 10, 4, 8になる. 次, x12=7も配列にないから, 7をT2に置き, 配列は6, 10, 7, 8になる.

次はx13=8だが, i=8の時に作ったT3に8があるので, 循環が始まっていたのがわかる.

そこで問題はiの行のTkを置いたのは誰かということだ. 最初の表からxを削除した表は次のようだが
i\k 0  1  2  3
11
21 2
33 2
43 24
55 24
65 64
77 64
87 64 8
99 64 8
107 104 8
これで例えばi=7, k=2ならこの4を置いたのはi=4であったということが知りたい.

いろいろやってみると, iからiの下からk+1ビットをとったものの, 先頭の0,1を逆にしたものを引くといいらしいことが分かる.

iの下からk+1ビットをとるには, 2k+2-1で1のビット列を作ってマスクする.

i &(2k+2-1)

これを 2kでxorすると先頭のビットが反転する.

それをiから引くのである. Schemeで書くと
(define (m i k)
   (- i (fix:xor (fix:and i (- (fix:lsh 1 (+ k 1)) 1)) 
    (fix:lsh 1 k))))
であって, 実際に計算すると
(do ((i 1 (+ i 1))) ((= i 11))
(display (cons i 
 (map (lambda (k) (m i k)) (a2b 0 (+ (lam i) 1)))))
(newline))

(1 1)
(2 1 2)
(3 3 2)
(4 3 2 4)
(5 5 2 4)
(6 5 6 4)
(7 7 6 4)
(8 7 6 4 8)
(9 9 6 4 8)
(10 9 10 4 8)
となって上の表と一致する.

従ってi=13の時のx13>=8がi=12, k=3にあり, (m 12 3)=8が分ったから, λ=13-8=5であった.

プログラムを書いてやってみよう. 最初の数列xiは, μとλから
(define (f i mu lm)
 (if (< i mu) i (+ (modulo (- i mu) lm) mu)))

(map (lambda (i) (f i 6 5)) (a2b 1 19)) => 
(1 2 3 4 5 6 7 8 9 10 6 7 8 9 10 6 7 8)
と作る. プログラム全体と実行結果を下に示す.

(define (iloop i x)
 (define (kloop k)
  (if (> k (lam (- i 1))) (iloop i x)
  (begin 
   (if (= x (vector-ref tab k)) (list i k (m (- i 1) k))
     (kloop (+ k 1))))))
(vector-set! tab (rho i) x)
(display i) (display (take (+ (lam i) 1) (vector->list tab))) 
(newline)
(set! i (+ i 1))
(set! x (f i mu lm))
(kloop 0))

(define mu 6) (define lm 5)
(iloop 1 (f 1 mu lm))

1(1)
2(1 2)
3(3 2)
4(3 2 4)
5(5 2 4)
6(5 6 4)
7(7 6 4)
8(7 6 4 8)
9(9 6 4 8)
10(9 10 4 8)
11(6 10 4 8)
12(6 10 7 8)

=> (13 3 8)
つまりi=13のとき, k=3の所に同じものがあり, その時のiは8であった. 従ってλは13-8=5であった.

2016年8月13日土曜日

複素数用計算尺

複素数用の計算尺があるという話をブログに書いたのは, 2009年8月のことだ.

最近すこし図を描き直し, Webでも楽めるようにしたので, その話をしたい.



上の図は対数の複素数を示す. 点Aから右へ延びる実線は1より大きい正の実数軸. Aは(虚数部が0の)複素数1. Bは複素数2, ...で, 2,3,4,...と対数的に間隔が狭くなっている.

Aで実数軸と交わる実線の曲線は, 実数部が1で虚数部が何かの複素数である.

一方, 点Cから右へ延びる破線はiより大きい正の(?)虚数軸. Cは(実数部が0の)複素数i, D は複素数2i, ... で, こちらも2i,3i,4i,...と対数的に間隔が狭くなっている.

Cで虚数軸と交わる破線は, 実数部が何かで虚数部がiの複素数である.

従って点Eは, 複素数1+iを示す.

実はこの上の図は下の真数のGauss平面を対数化したもので, 真数の図では実数軸と虚数軸の交点に0があるが, 対数では0は無限に遠くにあるので, 実数軸と虚数軸が平行になってしまう. 対数の図で, 実数軸と虚数軸の間にある部分は, 真数のGauss平面の第一象限に対応する.



この図には描いてないが, 正の実数軸の反対側には負の実数軸があるわけで, それは最初の対数の図では, Cから右に伸びる虚数軸のさらに上に現る. Cから上へACの距離だけ行ったところに, 負の実数の出発点 -1があって, 正の実数軸と同様に=1, -2, -3,...が右に出てくるのである.

さらにその上というか, Aから出発する正の実数軸の下というかに, 負の虚数軸があるのだが, 図がごちゃごちゃするから省略してある.

省略はそれだけではない. 10, 10i -10などの右には, これと同じパターンがあって, 20,30,... 20i,30i,... -20,-30,...と続くのである.

また左には0.9, 0.8,...があるのは理解できるであろう. つまり上の図はほんの一部を示しているに過ぎない.

もっと広範に描いて見ると, 下のような図になる.



前置きが長くなった. この図を使って複素数の乗算をするにはどうするかを述べよう. 通常の計算尺では, 2つの対数尺A, Bが次の図の上のように上下に並んでいて,

2掛ける何かを計算したければ, 下の図のようにB尺の1をA尺の2に合わせ, 例えばB尺の3に相対するA尺の目盛り6を見て, 2×3は6なることを 知るのである.

対数でも同様で, 1+iに何かを掛けようとするなら, B尺に相当する赤い図の1(点A)をA尺に相当する黒い図のEに合わせ, B尺の何か(例えば1+i(E))に 相対するA尺の値(Dつまり2i)を知る.



(1+i)×(1+i)=12+i2+2iで, 1の2乗とiの2乗の和は消えるから2iで正解だ.

このwebプログラムは http://www.iijlab.net/~ew/complexsr.htmlに置いてある. B尺のAを置きたいところをクリックすると, A尺が薄い色になり, 濃い色のB尺の1がその場所へ移動するから, B尺の乗数に相対するA尺を読めばよい.

2016年7月7日木曜日

プラニメーター

何日か前に書いたプラニメーターのブログに機能をもうひとつ追加した.

シミュレーターの画面にsquare, circleなどクリックすべき箱があったが, その隣に無名の箱があった. 実はここには曲線で囲まれた図形の例を置くつもりであった.

それが出来たので, シミュレーターを更新した. そのcurveの箱をクリックすると次のような図になる. この水銀整流管のようなものは「公立はこだて未来大学」のロゴである. その面積を求めることができる.



startをクリックするとメーターに



134.179が現れる. 従って面積は(* 134.179 400)=>53671.6である.

ところでこの自動追尾の方法だが, もとのロゴのビットマップを次のように作り,

0000000001FFFFF800000000
000000003FFFFFFFC0000000
00000001FFFFFFFFF8000000
0000000FFFFFFFFFFF000000
0000003FFFFFFFFFFFC00000
000000FFFFFFFFFFFFF00000
000003FFFFFFFFFFFFFC0000
00000FFFFFFFFFFFFFFF0000
00003FFFFFFFFFFFFFFFC000
0000FFFFFFFFFFFFFFFFE000
その輪郭線をとりだすプログラムにかけると
0000000001FFFFF800000000
000000003E000007C0000000
00000001C000000038000000
0000000E0000000007000000
000000300000000000C00000
000000C00000000000300000
0000030000000000000C0000
00000C000000000000030000
00003000000000000000C000
0000C0000000000000002000
のようなものが得られる. この線を辿るデータを作るプログラムをSchemeで書いた. シミュレーターではそのデータを追いながらプラニメーターの先端を動かしているのである.

ところで始めのビットマップに1がいくつあるか数えてみた. 54003個あった. プラニメーターの結果は53671.6だったから, まぁ大体うまくいっていると思ってよかろう.

2016年7月4日月曜日

プラニメーター

小学生の時, 算数の時間にいろいろな図形の面積の求め方を習った. その折, 先生から三角形とか台形とかではない複雑な図形の面積はどう求めるか聞かれた.

私の当然の答は方眼紙に写しとり, 図形内部の1ミリ方眼を数えることであった. もうひとつの答は, 厚紙に写しとり, 図形を切り抜いて重さを測り, それと同じになる矩形の厚紙を作り, その面積を求めるのであった.

後者の答は象の重さを測るのに, 象を船に乗せ, その時の吃水線の位置を船端に印し, 次に同じ吃水線に船が沈むまで小石を船の乗せてその小石の総重量を測るという中国の子供の話にヒントを得ていたかもしれない.

その後大学生くらいになって, プラニメーター(面積計)なるものがあることを知った. しかし, 自分で使う機会はなかった.

東大工学部計数工学科に勤めるようになった時, その学科では以前の学生実験でプラニメーターを使うものがあったと聞いたが, すでにプラニメー ターやタイガー計算機を使う実験はなくなっていた.

プラニメーターは機械式アナログ計算機である. プラニメーターの図は当今ではインターネットで探すといくらでも見付かる.

その一つが次である.

とあるページから借用させていただいた.



これは詳しくいうと Amsler polar planimeter という形式のものだ.

使い方は右上のOは下にピン, 上に重りのあるもので, 机のある場所に重りで固定する.

そして右下のAのピンで測りたい図形の周囲を一巡する. それに従い, 2本の腕は回転したり腕の間の角度を変たりしなが動く. 腕の蝶番BのあたりにCという車輪がみえるが, この車輪の軸はBAの腕に平行で, 特定の方向の回転だけを蓄積する.

すると不思議にも車輪の移動した長さにBAの腕の長さを乗じたものが面積になるのである. 長さ×長さだからディメンジョン的には合っている.

今回はプラニメーターの原理はかつ愛し, 最近プログラムしたプラニメーターのシミュレーターのことを書きたい.

そのプログラムは
http://www.iijlab.net/~ew/planimetersim.html
に置いてある. そこへ行くと



のような図が現れる. 左下のcircleをクリックすると, 今度は



のような図が現れる.

左端の赤の点が固定点(最初の図のO)で, 四角や円の上の青の点が図形をトレースする点(同A)である. その間の2本の線が, プラニメーターの2本の腕で, 折れ曲っているところが腕の蝶番(B)に対応する. 折れ曲りの近くの直交している線分が, 計測用の車輪の位置のつもりだが, その回転が左上のメーターみたいなもにに示される.

たとえば四角の図を表示している場合, 左下のstartをクリックすると, 青の点が移動し, 計測が始まる. 一周すると停止する. その時のメーターの読みは399.99になっている.



この図の四角は一辺が400の正方形で, プラニメーターの二の腕(BA)の長さも400になっている. メーターの読みが約400なので, それに腕の長さ400を乗ずると, 正方形の面積160000が得られる.

計測点Aの移動はプログラムで制御しているとはいえ, 凄い精度である.

次に円をやってみよう. 円は直径が400, 半径200なので, 面積はπ×200×200である.

それを腕の400で割ると, メーターの読みはπの100倍になるはずである. 実際にやってみると314.156が得られて安心する.



しかし, 手動でプラニメーターを操作すると, もちろんこんな精度は得られない. だからといって, 正方形の場合, 計測点を移動する際に物差しを当ててずらしてはならないと学生実験の手引書に注意がある. つまり線から一方に外れているとその誤差がそのまま集積されるからで, 手動でふらふら動かす方が誤差が左右に揺れてキャンセルされるので望ましいと書いてあった.

2016年5月1日日曜日

Prefix adder

60年くらい前, 東大物理学科高橋研究室でパラメトロン計算機PC-1を作っていたころ, 36ビットの並列加算器で, 高速に繰上げ(carry)したいという問題があった.

後藤さんはパラメトロンの多数決素子を巧妙に使い, 桁数に対して対数時間で繰上げを検出する回路を考案して, それを取り入れた.

先日出版されたDigital Design and Computer Architecture Arm Editionを眺めていたら, 同じような考えによる, Prefix Adderという回路の話があった. 今回は それが話題である.





この上の図が16ビットの二進数Ai, Bi (i=0,1,...,15)と下から来る繰上げCinを足すときのCarry-Lookaheadの回路だ.

上の図の上段の白い四角, 中頃の黒い四角, 下段の角のとれた白い四角の説明が下の図である.

ある桁 i が, 下の桁 i-1 からの繰上げの有無に関係なく, 上の桁 i+1 へ繰上げを生じるとき, 繰上げをgenerateするということで, Gi=1, そうでないとき, Gi=0とする. AiとBiが共に1なら繰上げがあるから, Gi=AiBi.

ある桁 i が, 下の桁 i-1 からの繰上げを上の桁 i+1 へ伝えるとき, 繰上げをpropagateするということで, Pi=1, そうでないとき, Pi=0とする. AiとBiのいずれかが1なら繰上げがあるから, Pi=Ai+Bi.

そうすると, i からi+1へ出る繰上げCi

Ci=Gi+PiCi-1

になる. このGやPをPrefixという.

これは単一の桁についての話であったが, iからjの桁のブロックについても同じことがいえる. ただGi, Piの代わりにGi:j, Pi:jのように書く.

そうだと思って上の図を見てみると, 左下のあたりに黒四角で左下隅に14:-1と書いたものがある. つまり15の桁への下全体のブロックからの繰上げはG14:-1 (と, P14:-1 と, Cin) で決るということを示す.

GiとPiがAiとBiとで書けたように, Gi:jとPi:jがそのブロックを途中で分けたGi:kとGk+1:j, Pi:kとPk+1:jで同じように書くと, Gi:jはi,k間で繰上げが出るか, k+1,j間で繰上げが出て, i,k間で繰上げが伝わるかだから, Gi:j=Gi:k+Pi,kGk-1,j.

Pi:jはi,k間でもk-1,j間でも繰上げが伝わるかだから, Pi:j=Pi:kPk-1:j.

ここまではDigital Designの本に書いてあることだ. 上の図のように半分, 半分と区間を分るにはi,jからkをどう決めればよいかは, この本には書いてない. そこがこの本とThe Art of Computer Programmingの違うところである. TAOCPならk(i,j)が詳しく記述してある筈だ.

しかしそう難しいことではなくて, λx(x>0)をxを二進法で表したとき, 最も左にある1のビット番号とすると, (λ1=0, λ2=1,λ3=1, λ4=2,...)

k=2λ(i-j)+j

と簡単に得れれる. Schemeでシミュレートしてみたが, これで合っているようであった.

2016年3月1日火曜日

菱形六十面体

菱形六十面体についてはこれまで何度もこのブログに書いた. 夫々の菱形が黄金比であることは, MathWorldにそう書いてあったのを鵜呑みにしたきらいがある.

今回はそれを確かめてみたい.

下の図は正十二面体である. 中心を通る赤, 緑, 青の線はx, y, z軸を示す. それぞれの軸は正十二面体の30本ある稜線のうち6本の中央を通過する. 次の図との関係で偶数だけがついているが, その頂点を持つ正五角形で座標を考える.


頂点間の距離, 稜線0, 2の長さを2とすると, 正五角形の対角線8, 4の長さは2φ. 従って4の頂点の座標は(φ, φ, φ). それを基準に他の座標も分かって, 1は(φ2, -1, 0);2は(φ2, 1, 0); 6は(1, 0, φ2); 8は(φ, -φ, φ)である. (1のx座標は中心から0までの距離と, 中心から4までの距離が等しいことから計算できる.) この正五角形に5個の菱形を作ったのが下の図の菱形六十面体だ.



これをx軸の方向から撮影した写真がこれである.



写真では視差があるから, 注目している正五角形のついて正確な図にしたのが次の図だ.



ちょうど1のところをx軸がこちらへ向かっている. その両側へ1離れたところに02が来る. また1から真上に行ったところに中心10があるわけだが, その高さは現時点では未定.

さて菱形は真上から見れば菱形であるが, 少し横から見れば菱形に見える方向も皆無ではないものの, 一般には菱形に見えない. しかし平行四辺形であるのは確かである.

従って0,1,10,9の菱形, 1,2,3,10の菱形を考えると, 0,91,3の辺は1,10の辺と平行である. また9,100,1と平行なども分かる.

すると8,75,40,1などに平行である. 一方 10,5, 6,7の菱形は対角線の方向から見ているから菱形のままであり, 10,55,6の長さは等しい. 5, 6の高さの差(φ2-φ=1)は既知だから, 10の高さも確定し, 計算するとφ-1=1/φである. 5,4の幅が1だから, 5のy座標も1/φと決まる.

ところで残りは奥行きである. 2, 3 の辺は, 横から見ると, 正面図の6, 5の辺に相当していることが分かる. 従って3の奥行きが分かり, 3, 4 がz軸に平行なことも分かる. 従って10, 5, さらに5, 6もz軸に平行だったことも分かった.

つまり正面から見た菱形10,5, 6,7は真の形であった.

という次第で菱形の縦横比は1:φであったわけである.

2016年2月24日水曜日

MathematicaのCal

いろいろなプログラム言語でCalのプログラムを書いてみた. 今回はMathematicaである.

Mathematicaというと, 一発処理のような印象である. 例えば Prime[25] と入力すると, 25番目の素数 97が出力される. しかしなんとかプログラムを書いてみようと思い, 例によってCalを書いてみた. 結局いろいろなノウハウを覚えつつ試行錯誤の結果完成したのが次のプログラムである. その見どころを説明しよう.



0行目. 仮引数はyr_のように後に下線をつける. 本体で使うときはその下線はいらない. システムの関数手続きは大文字で始めるが, ユーザの関数は小文字で始める. 本体の定義は:=の後に書く.

1行目. Module[{ の次にローカル変数を書く. ローカルな関数名も書く.

2,3行目. printd[d_]は, その月の最後の日付をzとして, 1≤d≤zの時は出力し, それ以外は3文字の空白にする. さらにdが1桁なら前に2文字の空白, 2桁なら1文字の空白をつける.

4,5行目. prline[d_]は, dからd+6までをprintdで文字列にしたものを, StringJoin(<>を使う)し, Printで出力する.

6,7行目. 年と月を出力. また曜日名を出力.

8,9行目. うるう年なら, ローカル変数lを1に, そうでなければ0にする.

10行目. 各月の朔日の曜日を決めるため, 次の年の初めまでの日数を計算.

11行目. 1,2月の日数と, 月初めから年末までの日数(mod 7)を設定.

12,13行目. 3月以降の月に上のような値を計算する.

14行目. 朔日の曜日wを確定する.

15行目. dを第1週の日曜に相当する日(1-w)から月末日まで, 7ずつ増やしてprlineする.

こうして出来たのが次だ.



なれない言語でプログラムを書くのはしんどい.