2016年11月30日水曜日

予定表

徐々に年末が近付きつつある. 年末独特の心配事は他にもあるが, まずは来年の予定の記帳だ. 通常の手帳は12月(かせいぜい1月かそれも第1週だけとか)で終りであり, 次の会合が翌年の始めだったりすると, 書き込む場所がない.

また例年に送られてきて使っている企業の手帳も年末にならないと届かない.

そういうわけで, この数年, もしかしたら十数年は10月ころになると, 次のような簡単な予定表をA4の紙に印刷し, 手帳に挟んで利用していた.

下に示すのは2017年版である. 左上の01S, 02M, ...は1月1日(日曜), 1月2日(月曜), ... を表す. 右端の列の31Fは3月31日(金曜)に対応する.

1日分が非常に小さいが, この予定表を作るころは, 新春の予定も少いのであまり困ることはなかった. 絵から分るように日曜だけは太字にしてある.

もちろん新しい手帳が届くまでの使命で, 手帳が来るとそちらへ転記し, お役御免になる.



ところが, 世の中では皆さんがスマートフォンで予定を管理するようになり, 学生のころから毎年使っていた手帳は, 企業でも今年(2016)の分から遂に作るのを止めた.

今年はとりあえず市販の手帳に予定を記入してきたが, やはり使い易いとは思わなかった. そこでこの暮は, 来年の予定表を自作することにした. その方が楽しくもある.



上の写真がこうして出来た個人用予定表の表紙である. 右側が表紙, 左側が裏表紙である. 大きさはB7版になる. (こうして開いた状態はB6版である.)

裏表紙はすぐ分るように, 2017年の七曜表である. しかし, 最上段が左から1月2月3月, 次の段が4月5月6月, ... ときて, よく見ると5段あり, 最下段は実は2018年の1月2月3月なのである.

この予定表の特徴は, 2017年1月から2018年3月までになっていることだ. つまり3ヶ月のバッファーを準備している.

従って, 表紙の次は



のように2017年1月になり, 裏表紙の前は



のように2018年3月である.

1日分が狭いのではないかといわれそうだが, 最近は予定表に記入する事項も少くなり, たとえ沢山書き込む日があっても, 欄外に書いて線で結ぶリンク構造にすれば楽々収まる.

そういう次第で来年はこの予定表を活用しようと考えている.

いい遅れたが, 表紙の西暦年号の下は曜日(上段)と月齢(下段)を計算する定数(枠内は左から1月2月3月 | 4月5月6月 | 7月8月9月 | 10月11月12月)である.

どちらも私のブログの話題になったもので, m月d日の曜日は左欄外の4にm月の定数とdを足し, 7で除した法をとる. 月齢は左欄外の1に足し, 30で除した法をとる. 右欄外は2018年の定数である.

曜日の定数が1月と2月で赤になっているのは, 閏年には1を引くことの注意である. 閏年には欄外の定数も赤で表示するようにプログラムしてある.

こういう予定表を作るには紙の中央で綴じる必要がある. 中綴じホチキスを探したが, 「マックス多機能スケール ナカトジール」というものを見付けそれを購入して使った.

ただ, 上の写真で分るように, 丁度折り目のところで留めるのは結構むずかしい.

予定表を吐き出すPostScriptのプログラムを書くのは, これまでカレンダーのプログラムを何度も書いたことがあるので, 簡単である. 来年, 再来年の祝日はインターネットで探すとすぐに見付かる.

2016年10月5日水曜日

ループ検出

前回のこのタイトルのブログで, Gosper流のループ検出法の説明をしたが, これでループが発見できる理由があれだけではいまいちである.

私がどう納得しているかも書いておく必要があろう. xiが次々と置かれていく様子を添字だけで再掲する.

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
これらの表中のxiに, 次のxがあるかどうかを調べるのだが, そこで添字の差だけの表にすると
iT0 T1 T2 T3
1 -1
2 -2  -1
3 -1  -2
4 -2  -3 -1
5 -1  -4 -2
6 -2  -1 -3
7 -1  -2 -4
8 -2  -3 -5  -1
9 -1  -4 -6  -2
10 -2  -1 -7  -3
つまり1つ前2つ前に同じものがあると, それらはT0で見付かり, 3つ前4つ前にあると, それらはT1で見付かり,...というわけだ. Tの添字が大きいほど, 繰り返し開始時のμの値の曖昧さが大きいことがこの表から判明する.

TAOCPでは演習問題3.1-7にこの話があり, その解答には

xnが表Tρnに格納されると, それはその後xn+1, xn+2, ..., xn+2ρn+1と比較される.

と記載してあるが, 上の説明の方が分り易いと私は思っている.

Gosperの元の記事はHAKMEM 132 (http://home.pipeline.com/~hbaker1/hakmem/flows.html#item132)にある.

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でシミュレートしてみたが, これで合っているようであった.