2008年12月22日月曜日

sideways addition

アルゴリズムを読んでいると, 簡単なものでも, あれ? と思うのがたまにある.

ビット列のsideways additionというのは, 1のビットの数である.
(1 0 1 0 1 0 1 0 1) なら5である. (Henry WarrenのHacker's Delightではpopulation countという.) あるいは∑xi (xiは0または1) といってもよい.



またまたTAOCPの話題で恐縮だが, x1, x2, ..., x9の9ビットのsideways additionが次のように書いてあった.

x10=x1⊕x2⊕x3,  x11=x4⊕x5⊕x6,  x12=x7⊕x8⊕x9,  x13=x10⊕x11⊕x12,
y1=⟨x1x2x3⟩,  y2=⟨x4x5x6⟩,  y3=⟨x7x8x9⟩,  y4=⟨x10x11x12

とすると

x1+…+x9=x13+2(y1+y2+y3+y4)

で得られる. x⊕y はxとyの排他的論理和; ⟨xyz⟩はx,y,zの多数決(TAOCPではmedianという)を表わす. x⊕y⊕zと⟨xyz⟩で全加算器になる.

非常に美しい式だが, 本当かとも思い, 早速やってみる.

(define base 10)
(define (sidewayadd xs)
(define (xor a b c) (modulo (+ a b c) 2))
(define (med a b c) (quotient (+ a b c) 2))
(let* ((x1 (modulo (quotient xs (expt base 0)) base))
(x2 (modulo (quotient xs (expt base 1)) base))
(x3 (modulo (quotient xs (expt base 2)) base))
(x4 (modulo (quotient xs (expt base 3)) base))
(x5 (modulo (quotient xs (expt base 4)) base))
(x6 (modulo (quotient xs (expt base 5)) base))
(x7 (modulo (quotient xs (expt base 6)) base))
(x8 (modulo (quotient xs (expt base 7)) base))
(x9 (modulo (quotient xs (expt base 8)) base))
(x10 (xor x1 x2 x3)) (x11 (xor x4 x5 x6))
(x12 (xor x7 x8 x9)) (x13 (xor x10 x11 x12))
(y1 (med x1 x2 x3)) (y2 (med x4 x5 x6))
(y3 (med x7 x8 x9)) (y4 (med x10 x11 x12)))
(+ x13 (* 2 (+ y1 y2 y3 y4)))))
(sidewayadd 000000111) => 3
(sidewayadd 111000111) => 6

なるほど. letの中の, 割ったり剰余をとったりしてx1, x2,...を作っているのは, (sidewayadd '(1 1 1 0 0 0 1 1 1))と書くのは面倒なので, (sidewayadd 111000111)としたいからである.

しかし落ち着いて考えてみると, アルゴリズムが正しいことは簡単に分かる.



この図のように9ビットを右からx1,...,x9とする. 3ビットずつに区切り, その区間の和を求めると0から3までになるから, 2ビットで表現出来, それが下のx10,y1などである. x1,x2,x3の3ビットの二進法の和は, その1の桁(x10)はx1,x2,x3の排他的論理和だし, 10の桁(y1)はそれらの多数決だから, アルゴリズムにはそう書いてある. 次に3つの区間の3ビットの和が得られたとして, この和を取るわけだが, 和の1の桁はx10,x11,x12の排他的論理和である. 10の桁y1,y2,y3を足せばよく, それにx10,x11,x12からの繰上がり, すなわちそれらの多数決を足す. 10の桁だから2倍する.

と書いてしまえば, 当然なアルゴリズムであった.

この話題では, 思い出すことが2つある. 1つは英国ケンブリッジ大学のEDSACのライブラリのsideways additionで, 第2版に登場した.

EDSACの語長は35ビットだが, 36ビットと考えて6ビットの区間6つとする. ビットシフトやビット毎ANDや加算して, 6つの区間で並列にその区間の1の数を和を作る. 6までだから3ビットの範囲に収まる. つぎに105105...051(つまり5個の0と1が交互の列)を掛けると中央の6ビットの区間に総和が現れる仕掛けである.

この話はTAOCPのbitwise tricks and techniquesにも書いてある.

もう1つはパラメトロン計算機PC-1の語長が36ビットだったことに関係する. この36なる数は, 9ビット掛ける4なのだ. 9ビット毎にパリティをとり, エラー検出ビットを設ける計画であった. 3個のパラメトロンの信号のパリティは, 3重平衡変調器という回路で簡単に取れるはずで, それを2段用いれば9個のパリティが取れるというもくろみであった. そしてメモリーは1語当り情報36ビット, 検査用4ビット, 計40ビットで構成するはずであった. しかし3重平衡変調器が思ったほどにはうまく働かず, 結局パリティビットなしになった.

ところで, 最初の美しい式はy4まではBoole式だが, 最後に算術演算で足している. これは反則のようにも見えるが, y1+y2+y3+y4が実は(yiは0または1なので)再びsideways additionになり, 再帰的に回路を構成すればよい. Knuthの計算機MMIXには, 64ビットまでのオペランドのsideways addition SADDがあるが, かなりの回路になるであろう.

2008年12月12日金曜日

Gaussの正十七角形

Gaussが正十七角形が描けることを証明したという話は伝説として聞いてはいるが, 高木先生の近世数学史談に詳しく書いてある. 計算を追うのは大変だが, なるほど描けそうだという気分になるには十分である.

角2π/17をφとし, cos φを計算すると

が得られる. 凄い式だが計算機の威力を借りて検算すると,


(+ (/ -1 16)
(/ (sqrt 17) 16)
(/ (sqrt (- 34 (* 2 (sqrt 17)))) 16)
(/ (sqrt (+ 17 (* 3 (sqrt 17)) (- (sqrt (- 34 (* 2 (sqrt 17)))))
(* -2 (sqrt (+ 34 (* 2 (sqrt 17))))))) 8))
-> .9324722294043558
(define pi (* 4 (atan 1)))
(define phi (/ (* 2 pi) 17))
(cos phi) -> .9324722294043558

と驚く.

cos φが分かればsin φも得られ, 下の左の図のように正十七角形が描ける.

まず単位の長さで線分ABを引く. これが正十七角形の1辺になる. 次にBからcos φだけ右へ延長しCをとり, 上へsin φだけ延長してDをとり, BDを結ぶ. これが次の辺になる. 以下同様に15回繰り返せば, 図のように正十七角形が完成する. もちろんGaussが自分でこういう絵を描いたとは思えない. 私に描けたのは, PostScriptのお蔭である.




しかし, ルートなんたらが沢山出てくるから, cos φの長さを取るのはおおごとである. 基本的には右上の図のように, AB=1をとり, 直角にBC=1を延すとAC=√2となる. CからACに直角にCD=1を取ると, AD=√3になり, 同様にして整数の平方根はいくらでも取ることが出来る. Gaussの式には√17が沢山現れるが, それには1辺4と1辺1との直角三角形を作れば斜辺は√17になってくれる.

右下の図は√(a-b)の取り方を示す. ABを√aの長さにとり, ABを直径とする半円を描く. BC=√bを円周上にとれば, ACが求めるものである.

2008年12月9日火曜日

月齢カレンダー

直前のブログ「月齢カレンダー」の月齢計算で, 「ある定数」の決め方がなんとなく杜撰に思われるので, もう少し説得力ある計算法を考えた.

ある日の正午の月齢は, 直前の新月からの時間である. NAOJのウェブページを見ると, 2008年12月27日21時22分が新月だそうだ. すると12月28日正午までの経過時間は14時間38分. 24時間に対しては0.6097...となり, つまりそれが28日の月齢である.

従って望む日の月齢は, 1朔望月を29日499/940とすると, 望む日のRDから2008年12月28日のRDを引き, その29日499/940の剰余プラス0.6097となる.

(define (rd y m d) ;rdの計算
(let* ((p (lambda (n) (= (modulo y n) 0)))
(q (lambda (n) (quotient y n)))
(s `(0 31 ,(if (if (p 100) (p 400) (p 4)) 29 28)
31 30 31 30 31 31 30 31 30 31)))
(- (+ d (* y 365) (q 4) (q -100) (q 400))
(apply + (list-tail s m)))))

(define (ma y m d) ;月齢の計算
(+ (/ (modulo (* (- (rd y m d) (rd 2008 12 28)) 940)
27759) 940.) 0.6097))

(do ((i 1 (+ i 1))) ((> i 12)) ;1日の月齢表の計算
(display (/ (round (* (ma 2009 i 1) 10)) 10)) (newline))

4.6
6.1
4.5
6.
6.5
8.
8.4
9.9
11.4
11.8
13.3
13.8

一方19年7閏法と有理化した朔望月の求め方は次の通り.

1太陽年=365.242198日
1朔望月=29.530589日
1太陽年から12朔望月を引くと10.875166日

これが1朔望月の何分の何に相当するか知りたいので, この比を連分数展開する.

小数xの連分数をn項とるのに

(define (cf x n)
(if (> n 0)
(let* ((y (/ 1 x)) (z (inexact->exact (truncate y))))
(cons z (cf (- y z) (- n 1)))) '()))
として
(cf (/ 10.875144 29.530589) 10)
を計算すると
(2 1 2 1 1 17 2 2 2 2)
が得られる. 17分の1辺りで棄てることにし,
(define (cf2r l)
(if (null? l) 0
(/ 1 (+ (car l) (cf2r (cdr l))))))
で有理数を計算すると
(cf2r '(2 1 2 1 1))→ 7/19

となり, 余りは1朔望月の7/19であることが分かる. 従って19年に7回閏月を置くのである.

次に19太陽年は 365 1/4日の19倍となる. この時間に朔望月は19×12+7回(つまり235回)あるので,

(/ (* (+ 365 (/ 1 4)) 19) (+ (* 12 19) 7))→ 27759/940
(quotient 27759 940)→29
(modulo 27759 940)→499

すなわち29日499/940となる.

追記 上のような月齢の式を書いたが, 最後に0.6097を足すのはこれまたウンシャンである. その回避法は2つあり, まずは2008年12月28日のRDではなく, 2008年12月27日21時22分のRDを基準にすること. つぎは正午ころに新月のある日を見つけることである. たまたま手元の天文年鑑を見ていたら, 2005年8月5日の新月は12時5分であった. これはしたり! 1日は86400秒であり, 5分は300秒なので, 誤差は無視出来よう. 従って

(define (ma y m d) (/ (modulo
(* (- (rd y m d) (rd 2005 8 5)) 940) 27759) 940.))

はそれほどウンシャンではない式だと思っている.

2008年12月8日月曜日

月齢カレンダー

これまでカレンダーのプログラムは数えきれないほど書いている. 2007年の夏のプロシンのお題もカレンダーであった. Unixのコマンドcalのように, 年と月を入力し, 例えば

December 2008
S M Tu W Th F S
1 2 3 4 5 6
7 8 9 10 11 12 13
14 15 16 17 18 19 20
21 22 23 24 25 26 27
28 29 30 31

を出力するものだ. プロシンのときは, mss, scheme, haskell, teco, tex, postscriptで書いた.

最近は不況のせいか, いや私にカレンダーを送ってくれる同僚が殆んど退職したのであろうが, 富士通の「世界の車窓」以外は, カレンダーもあまり貰わなくなった. それでこの数日カレンダーを自作しようかと考えていた.

以前はpcalというプログラムがあったが, まぁそういうものは, 自分で書く方が早い. ところで, そのpcalのカレンダーに, 月齢の絵がついているのを見つけた. lcalともいうらしい. それはどうすれば書けるか, 思いがけなく熱中した.

月齢カレンダーを書くには, 2つの問題を解決しなければならない. まずy年m月d日の月齢の計算法. 次に月相の表示法である.

とりあえずサーチすると月齢の計算法は見つかるが, どの式を見てもウンシャンである. いろいろ考えた末, 私が採用したのは次のものである.

カレンダーに図示するのだから, それほど正確である必要はない. それよりエレガントに計算したい. まずy年m月d日のある起点からの日数がいる. それにはReingold, DershowitzのCalendrical CalculationsにあるRD(Rata Die (fixed date))を使うことにする. これはGregorian暦が昔まで使われたとして, 西暦1年1月1日(月曜)からの通日である.

一方, 朔望月は29.53...日であるが, これはTAOCPの復活祭公式にあるように, 有理数29 449/940にする. つまりRDにある定数を足し, 有理朔望月で除した剰余を月齢とするのである.

ある定数が問題であるが, たとえば自宅にあった2004年と2008年の理科年表の太陽, 月の表の正午の月齢からキャリブレートすれば得られる. こうして計算したのが, 以下の表で, 各年は1月1日から毎月1日の月齢で, 左が私の計算, 右が2004年と2008年は理科年表, 2009年と2010年はどこかのウェブページによるものである. まぁよかろう.


2004 2008 2009 2010
8.5 8.7 | 22.5 22.4 | 4.6 4.6 | 15.2 15.6
10.0 10.2 | 24.0 23.6 | 6.1 5.8 | 16.7 16.8
9.5 9.7 | 23.5 23.0 | 4.6 4.1 | 15.2 15.0
10.9 11.2 | 24.9 24.4 | 6.0 5.5 | 16.7 16.3
11.4 11.6 | 25.4 25.0 | 6.5 6.0 | 17.1 16.6
12.9 12.9 | 26.9 26.6 | 8.0 7.6 | 18.6 18.1
13.3 13.3 | 27.3 27.3 | 8.4 8.3 | 19.1 18.7
14.8 14.7 | 28.8 29.0 | 9.9 10.0 | 20.5 20.3
16.3 16.1 | 0.7 1.3 | 11.4 11.7 | 22.0 22.0
16.8 16.5 | 1.2 1.8 | 11.8 12.3 | 22.5 22.7
18.2 18.0 | 2.7 3.2 | 13.3 13.9 | 23.9 24.3
18.7 18.5 | 3.1 3.4 | 13.8 14.3 | 24.4 24.9


次に月相の図であるが, これは要するに図学の話題で, 図学は教養学部の学生であった頃の得意科目であり, 簡単の単であった.


2つの円があり, 上は月を北の極から見たもの, 下は月を地球から見たものである. 上の図の真上に太陽があるのが新月で, 図はφだけ太陽が右によった, あるいは月が東に進んだことを陰で示す. この陰の境界を地球から見た図に書けばよいわけで, h-h'の水平面で月を輪切りにし, 各断面での陰の位置を下の図へ移せば良い.

これで難問は片付いたので, これも似たようなのがインターネットにあったが, 2009年の各月の奇数日の月相を書いてみた.



7月21日と23日に新月が見える. この中間の22日に日食がある.