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日に日食がある.

2008年11月29日土曜日

Knuth先生の小切手

8月の米国出張の際, オークランドのFranz社を訪問した. 鉄道マニアの私のことゆえ, 往きは当然Bartで行ったが, 帰りはFranzの小俣さんがサンフランシスコまで車で送ってくれた. その途中サブプライム問題による売り家の看板をいくつも見かけ, 金融危機をいささか実感した.

ところでこの度, Knuth先生がTAOCPのエラー発見者に対する小切手発行を止めるというので, 金融危機が急遽身近かになった. 完全主義者のKnuth先生は, 自著のエラーを最初に報告した人に十進法で2.56ドル(十六進法で1.00ドル)の小切手を贈ることにしてきた. 小切手は通常は1ヶ月以内に振出人に戻って来るが, Knuth先生の小切手は殆んど現金化されないので(「cashedは僅かで残りはcached」だそうだ), 多くの小切手が世界中に残留しており, 口座番号が知られる危険性があり, 10010年に1度の金融危機の現状では放置出来ない仕儀と相なった.

その代りとして, Pincus惑星に支店のある, San Serriffe銀行に, 小切手受取人の口座を開き, そこに支払い証明(certificate)を積み立てることにしたというのだ. もちろんこれは架空の銀行で, われわれLisperならfoo銀行のbar支店というところである.

PincusとかSan Serriffeは, Knuth先生は好きらしく, TAOCP分冊0には, Pincus惑星の論理学者がどうのこうの(演習問題7.1.1-2)とか, San Serriffe島にはm羽の鳩がいる(演習問題7.1.1-41)とか, 出てくる. 私は漫画やSFには暗いので, これらの名詞がそういうものと何か関係があるかは存ぜぬ.

さて, その新システムの支払い証明がコンピュータ屋には面白い. いろいろ新機軸が見らる.

米国で暮した人は知っているように, ある受取人に2.56ドル渡す場合, 小切手には

Pay to the order of の次に受取人を書き, その右の方に $2.56 と書く.

中央の枠には綴りで Two and 〜〜〜〜〜〜〜〜 56/100 Dollars

のように書く. 従ってHundredとかThousandの綴りも知らなければならない. ドル以下は100分のなんとかと書くことになっている.

新システムにはKnuth先生の趣味がふんだんに見られ,

まず受取人の右の金額は十六進法で 0x$ 1.00 となる.

中央の枠は

One and 〜〜〜〜〜〜〜〜〜 no/256 Hexadecimal Dollars

と変った.


これは実は私として, 大変気持ちが悪い. Knuth先生のホームページなどに, $2.56の由来が書いてあるが, 「256 pennies is one hexadecimal dollar」 がその理由だ.

私に言わせれば, 1は十進でも十六進でも1であり, one hexadecimal dollarはすなわち one decimal dollar である. 10016¢なら25610¢ = $2.5610で, これが理由なら納得出来る. 10016¢ !=10010¢なのだ.

エラー発見は$2.56であるが, 改良案を示すと32¢頂ける. つまり0x$0.20 である. 改良案2件なら 0x$0.40, 5件なら 0x$0.a0, 7件なら 0x$0.e0というわけだ. そこで新システムの中央の枠だが, 改良案1件の場合

Zero and 〜〜〜〜〜〜〜〜 32/256 Hexadecimal Dollarsとなるのだろうか?

Zero and 〜〜〜〜〜〜〜〜 20/100 Hexadecimal Dollarsではないかと思うのだが.

さらに仮りにエラーを1210件見つけたら,

Twelve and なのかな. まさか C and とか書かないであろう.

問題は十六進の桁の読み方, A,B,...,Fの綴り方がないことである. 私は私なりの提案を持つが, それはまたいずれ.

追記 その後Knuth先生から32¢の支払い証明が送られてきた. 確かに0x$0.20とあり, 中央の枠には案の定

Zero and 〜〜〜〜〜〜〜〜 32/256 Hexadecimal Dollars

と書いてあった. 32¢を送る航空便には94¢の切手が貼ってあった.

2008年11月17日月曜日

萬葉集物語

1月ほど前, 新聞の出版広告を見ていたら, 思いがけず「森岡美子著 萬葉集物語」が復刊されたとあった.

実はこの本は, 私が小学生の頃, 繰り返し読んだ思い出深い書物である. その頃の本は, 家が空襲で焼けた時, 一緒に焼けてしまったはずだ.

しかし本書を読んだお蔭で, 萬葉集に対する興味が深まり, 多くの和歌や少なからぬ長歌も諳じた.

ところで左様な昔の本が, 2/3世紀もたった今頃, 復刊されたのは, 驚きである. 実は私の読んだ本とは別に, 多少の手直しを経て, 戦後の昭和27年に冨山房から新たに刊行されていたらしい. 今回復刊されたのは, それであった. 早速求めて読み直したのはいうまでもない.

私が愛した本は, ところどころに色刷りのページがあり, そこには「あおによし」「田子の浦に」など有名な和歌が万葉かなで書いてあった. それが冨山房版にはなかった.

この本は, 言葉使いが大変ていねいである. 「みなさまは, どうお考えになりますか」という調子である. 私は地の文の文体までは覚えていなかったが, こうい調子の本であったらしい.

小学生のころ, 著者が森岡美子という方だとは知っていたが, 女子学習院で教鞭をとられていたことは, 今回知った. そして90歳を超されてご健在であった.

考えてみると, 昔の先生は言葉が一般にていねいだったようにも思う. 私が中学1年の時, 「みなさまが将来論文をお書きになるなら」という先生もいらして, われわれも大人になったのだと実感した.

それはとにかく, 昔の優れた本が, 復刊されるのは, 嬉しいことである.

2008年11月7日金曜日

酉の市

11月5日は一の酉だった. したがって今年は三の酉まである. 酉の市で人出が多いのは, 大鷲神社らしい. 以前, 何回か酉の市の晩に大鷲神社に繰り出したが, 熊手を求める人と, 冷やかしと, もちろんお参りの人とでごった返し, 身動きもままならない様子である.

今年は6月に竜泉の一葉記念館を訪れたので, その帰りに大鷲神社に立ち寄ってみたら, 誰もいなかった.

三の酉といえば, 三の酉まである年は火事が多いといわれる. しかし酉の日は12日毎にあり, 11月は30日なので, 三の酉のある確率は0.5になり, 2年に1度は火事が多いことになる.

世の中もそういう解釈らしいが, 三の酉はあるとすれば11月の25日から30日なので, かなり寒さが身にしみ, 火事と結び付いたのであろう. 三の酉がない年の二の酉は11月の19日から24日までになければならず, 酉の市の寒さの印象は異るのであろう.

酉の市の平均が2.5回あるのと同様なのが, 土用の丑の日だ. 夏の土用は立秋の前の18日間なので, 12日毎の丑の日確率は1.5になる. 二の丑のある夏は暑いといわれないのはなぜか.

ところで曜日も干支も, 閏や大小の月に関係なく, 正確に繰り返される. したがって計算は簡単だ. 一の酉であった(今年の)11月5日のユリウス日数は, 例の個人用電卓によると, 2454776で, そのmod 12は8である. 従って剰余と干支は

子丑寅卯辰巳午未申酉戌亥
11 0 1 2 3 4 5 6 7 8 910

と対応することが分かる. (ユリウス日に1を足して12の法をとると, 0,...,11が子,...,亥に対応すると覚えた方がよいかも.)

同様なことは, 年の干支と西暦の間にもあり, 12で割りきれるのは申年, 60で割りきれるのは庚申である. これを知っていると壬申の乱(672年)は確に申年だと了解出来る.

昔のことだが, 三島由紀夫が自決したのは, 1970年11月25日. この日のユリウス日は2440916, mod 12はやはり8で, 実は三の酉であった.

2008年10月31日金曜日

ユリウス日

数日前に書いた個人用電卓にユリウス日のキーがあった. これは特に週日を求めるためのものである. ユリウス日に1を足し, 7で法をとると日曜を0, 月曜を1, ..., 土曜を6とする週日が得られるのである. またこの電卓は, 除算のとき, スタックトップには商が, そのすぐ下には剰余が出るので, 20081020 [JD] 1 [+] 7[/] [Pop] -> 1 となり, 本年10月20日は月曜と分かる. -4712年1月1日も月曜であったわけだ. ユリウス日は, -4712年1月1日 世界時正午を0.0とし, それからの通日である. 2008年10月20日のユリウス日が2454760というのは, 正確にはその日の世界時正午, JSTで午後9時のユリウス日が2454760.0であるということだ. それを省略してその日のユリウス日という. -4712年1月2日のユリウス日は1である. 理科年表や天文年鑑にも計算用の表がある.
A表              B表    *   C表      *
-1000 1355808     0  0,-1    1   0   0
 -900 1392333     1   365    2  31  31
 -800 1428858     2   730    3  59  60
 ...              3  1095    4  90  91
 1300 2195883     4  1460   ...
 1400 2232408     5  1826   11 304 305
 1500 2268933    ...        12 334 335
 1500 2268923    98 35794
 1600 2305448    99 36159
...
 2000 2451545
例えば2004年4月23日のユリウス日ならA表から2000年の2451545, B表から4年の1460, C表から4月の右の閏年用の91をとり, それに23を足して(+ 2451545 1460 91 23) -> 2453119が得られる. このA表は, 例えば2000年なら, -4712年1月1日から1999年12月31日までの日数である. 従ってこれは2001年1月1日のユリウス日のになりそうだが, B表の最初0年に, 閏年なら-1, 平年なら0とあるのが重要で, 2000年は閏年なので, 1引くから結局1月0日のユリウス日になる. 変換の式もサーチエンジンで探せるが, 私が情報処理学会誌のHaskell Programmingに書いた,
mon0, mon1 :: [Int]
mon0 = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]
mon1 = [0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335]

gleap, jleap :: Int -> Bool -- Gregorian, Julian暦のうるう年
gleap y = if y `mod` 100 == 0 then y `mod` 400 == 0
                              else y `mod` 4 == 0
jleap y = y `mod` 4 == 0

julianDate :: Int -> Int -> Int -> Int
julianDate y m d =                   -- aからgは途中の値
 let a = (y + 4712) * 365
     b = (y + 4712 + 3) `div` 4
     c = if y > 1601 then y' `div` 400 - y' `div` 100 else 0
           where y' = y - 1601
     e = if [y,m,d] >= [1582,10,15] then -10 else 0
     f = (if leap y then mon1 else mon0) !! (m - 1)
           where leap = if y > 1600 then gleap else jleap
     g = d - 1
  in a + b + c + e + f + g           -- 最後の結果
のようにも書ける. インターネットなら, https://www.onlineconversion.com/julian_date.htmにJulian Date Converterがあり, Year, Month, Day(, Hour, Minute, Second)を入れ, Computer Julian dateを押すと結果が返ってくる. 個人用電卓では結局上の式を使った. 電卓は実はやっとiPod touchに実装した. 電車の中でもテストが出来て楽しいが, 実用化にはまだ改良の余地は多々ある.

2008年10月29日水曜日

多面体描画道楽

TAOCPにこういう絵がある. (演習問題7.2.1.2-60)

{0,1,2,3}の24通りの全順列を, 隣り同士の交換で実現しようというものである. 天辺にある0123を最初の順列とする. 左の2個を交換すると1023になり, それは天辺から左へ稜線を辿ると出会う1023である. 次は右の2個を交換し1032が得られる. このようにオレンジ色の線を一筆書きでたどれば全順列を通過するというわけである. これはこの立体のHamilton閉路(つまり全頂点を回って元へ戻る経路)になっている.

この立体は英語では"truncatedoctahedron", 日本語では「切頭八面体」というらしいが, いまいちな感じだ. これは正八面体の頂点を切り落としたものである.

赤い線で示すのが正八面体で, 黒線のように切り取る.

ところでこの立体は, 正八面体の頂点を削れば出来るが, 立方体(正六面体)をの頂点を削っても出来る.


従って「切頭六面体」でもある. ただ削り取る量ははるかに多い. もともとこの立体には, 正方形が6枚, 六角形が8枚あるのだが, 正八面体の6つの頂点が6枚の正方形に辿り着き, 立方体の8つの頂点が8枚の六角形に辿り着く.

同様なものには, サッカーボールともいわれる, 「切頭二十面体」がある. (truncatedicosahedron)



これも, 名前が示すように, 正二十面体の12個の頂点から削り, 12個の正五角形に辿り着くのと, 正十二面体の20の頂点から削り, 20個の正六角形に辿り着くのがある.



赤線で示す正二十面体の頂点から削って, 「切頭二十面体」を作る.



この方は赤線で示すのが正十二面体であり, その頂点から削り始めると, 正五角形の面が段々と削られ, 小さい正五角形が現れる.

さらに削り続けるとどうなるか. 正五角形は消滅して頂点になる. 一方正六角形は正三角形になる. 要するに正二十面体になってしまう. これが正十二面体と正二十面体は相対の関係にあることの証明である. 一方から他方へ移るちょうど途中にサッカーボールがあったのだ.

このことは, 最初の立体でも同じで, 正六面体と正八面体は相対の関係にあることが分かる. 正四面体は自分同士で相対の関係にある.

ではその途中はどういう立体か. その「切頭四面体」
(truncatedtetrahedron)を書いてみた.



あまり見かけない立体ではあるが, 4個の頂点から現れた正三角形が4個の正三角形を六角形にして浸食し, ついには頂点になり, 正三角形の方が面になって, 面と頂点の機能を逆転した正四面体になる. この六角形がちょうど正六角形になった瞬間が上の図である. つまり削られて短くなりつつある元の稜と, 削りながら長くなりつつある新しい稜の長さはどこかで等しくなるのである. この点の座標の計算は簡単である. 逆転した正四面体を青線で示す.

2008年10月28日火曜日

個人用電卓

Hewlett Packardの電卓16Cを愛用していた. 自分でもそういう電卓を設計したいと思い, 何年か前にHappy Hacking Calculatorのシミュレータを書いた. これはパソコンの中で動いているので, 計算するならパソコンを使えばよく, お遊びの域を出なかった.

それが最近, iPod touch用のシミュレータがあるので, iPod用に書き直してみた. それについて書いてみたい. (HHCalcの手引き)

見かけは下の通りである.



左半分が入力キー, 右半分が演算キーである. 入力キーの上部に結果を表示する.

演算の機能は次の通り:

[/] y/x
[*] y*x
[-] y-x
[+] y+x
[Hex] 基数を16進に
[Dec] 基数を10進に
[Oct] 基数を8進に
[Ent] 入力終り
[Pop] スタックトップを捨てる
[Dup] スタックトップを複製
[Exc] x<->y
[Chs] -x
[Sqt] √x
[Fct] xの素因数分解
[Pow] xのy乗
[JD] ユリウス日 2008年10月20日のJDは 20081020 [JD] ->2454760
-4712年1月1日のJDは 47120000[Chs]101[+][JD] -> 0

今のところはシミュレータで動いているだけだ. [JD]のようなマイナーな機能は別画面で指定するようにし, [JD]の場所は別画面起動用にするのもいいかと考えている.

2008年10月2日木曜日

単調関数

The Art of Computer Programming の先の方を読んでいたら, monotone-function function というのに出会った. もちろん, monotoneな関数とは, 単調な, 俗語でいえば「右肩上がり」か「右肩下がり」なものである.

しかし単調関数の関数となると, どうもそう簡単なものではないらしい. そこでの説明は以下の通り. 式が多いので, texで書いて張り付けてある.



上のμnの式の右辺のbigwedgeの下にある, 0≤ij≤2nは, 「nビットのiの各ビットがjの各ビットより大きくない対について」という意味である.

TAOCP流にいうと, the bit string i=i0...in-1 is regarded as contained in or equal to the bitstring j=j0...jn-1 if and only if ikjkfor all k.


これでやってみるとmonotone functionは
0変数で2種類(0, 1),
1変数で3種類(0, x, 1),
2変数で6種類(0, xy, x, y, x&ory, 1)あり,
3変数の20種類は以下のとおりである.
0, (∧ x (∧ y z)), (∧ y z), (∧ z x), (∧ (∨ x y) z), z, (∧ x y), (∧ (∨ z x) y), y, (∧ (∨ y z) x), (∨ (∧ x y) (∨ (∧ y z) (∧ z x))), (∨ (∧ x y) z), (∨ (∧ z x) y), (∨ y z), x, (∨ (∧ y z) x), (∨ z x), (∨ x y), (∨ x (∨ y z)), 1

3変数関数はSchemeでチェックしたので, Cambridge Notationになっている.
monotone functionは否定なしで, ∧, ∨ だけで書けるというのはこれで分かる.



この図では, 各キューブが1つのmonotonic functionを示す. 座標軸はxが手前, yが右, zが上である. 角の黒丸が関数値が1になることを示す. 左上が3変数関数の0に, その右隣りがxyzに, のように対応している. 右下が最後の1である. 大学にいた頃, monotonic functionは, このキューブを羊羮だと考え, 黒丸のコーナーと黒丸なしのコーナーが包丁で一度に切り分けられるもの, といったことがある.その様子は見て取れよう.

ついでだが, この図で黒丸の数は0から8まである. 0個の関数は1, 8個のは1である.
1個のは xyz, 7個のは xyzである.
2個のは xy 型, 6個のは xy 型だ.
3個のは x ∧ (yz)型, 5個のは x ∨ (yz)型.
4個のは x 型か (xy) ∨ (yz) ∨ (zx)=(xy) ∧ (yz) ∧ (zx) これは自己相対で<x, y, z> (多数決)である.

μnの真理値表は以下の通り:

μ0=(1,1)
μ1=(1,0,1,1)
μ2=(1,0,0,0,0,0,0,0,1,0,1,0,1,0,1,1)

最後の真理値表から, 2変数関数が単調になる真理値は(一番左を0番として)
0,0,0,0(0番)
1,0,0,0(8番)
1,0,1,0(10番)
1,1,0,0(12番)
1,1,1,0(14番)
1,1,1,1(15番)
だから, 0,0,0,0は0に, 1,0,0,0はxyに, 1,0,1,0はyに, 1,1,0,0はxに, 1,1,1,0はxyに, 1,1,1,1は1に対応する.

この出てきた文脈は, BDD(binary decision diagram)の合成である. これについてはまたにしよう.

2008年9月26日金曜日

パラメトロン計算機完成50年

本年(2008年)3月26日のパラメトロン計算機PC-1の完成50年記念から半年がたつ.

記念の日の講演会会場で, NIIの橋爪君たちのパラメトロン加算回路のデモがあったが, それ以外にもこの半年にはパラメトロンのニュースがいくつかあったので記録しておく.


NTTの研究所でもパラメトロンの試作をしているというニュースを読んだ.

「パラメトロン・コンピュータよ再び」NTT研が低消費の論理素子を開発 EE Times Japan 2008年6月号

長さ250μm, 幅85μm, 厚さ1.4μmの長さの両側から圧電で正弦波の力を加えると, 上下に曲がりながら振動する. それが0相, π相になるという.


9月になると, 逗子の井村さんから, パラメータ発振に成功したというメイルを貰った. 励振と発振の様子NOTらしいものが見える.

ただ加算回路までは出来ていないらしい.

パラメータ励振に興味を持つ人がいるのは嬉しいことだ.

2008年9月19日金曜日

多面体描画道楽

小学生の頃, 算術か図画か工作の教科書に正五角形の書き方が載っていた. 正三角形, 正方形, 正六角形ならコンパスで描けるが, 正五角形は分度器を使うものと思っていたから, この方法には驚かされた. その後, 1辺の長さ1の正五角形の対角線の長さが(√5+1)/2であることを知ると,なんのことはなく, 驚きは霧散した.

この描き方は正五角形の書き方をサーチエンジンで探せばあるだろうが, 図のように底辺ABを引く. 辺の長さを1とする. その中点Cから垂線を立て, CD=1の点Dをとり, ADをDの方向へ延長してDE=1/2に点Eをとる. AEに等しくCDの延長線上にAFをとると, Fは正五角形の頂点になる.



正五角形をPostScriptで描くのは簡単だ.
/l 200 def     % 辺の長さ
0 0 moveto % 左下の位置へ移動
0 72 216 {dup cos l mul exch sin l mul rlineto} for
% 0度から72度おきに216度まで
closepath stroke


描いて楽しいのは正多面体である. 高木貞治先生の「数学小景」には, 正十二面体と, 正二十面体のEuclid幾何学原本による描き方が載っている.

どちらも立方体から始める. 正十二面体では, 面を構成する正五角形の対角線を1辺とする立方体をまず描く. これが内接立方体になる. 対角線を1とすると, 辺の長さは(√5-1)/2である. これををaとする. 立方体の上のの面には, 中央に前後に長さaの線を引く. 左右の面には, 中央に上下に, 前後の面には左右にそれぞれ長さaの線を引く. 図ではこれらの線を多少太めに示した.

さて上の面の線の両端U,Vから上方へ長さa/2の線UP, VQを立てる. 右の面の線の上端から右へ長さa/2の線WRを立てる. こらP,Q,Rと元の立方体の頂点A,Bとが正五角形になる.

中央の線は各面にあるから計6本で, そのそれぞれから2枚の屋根のような正五角形が作られるから, 十二面体になる.




一方正二十面体は同じような立方体の内部に作る. これが今度は外接立方体になる. 先ほどの各面の中央の線の両端U,V,W,X,Y,Zの隣り同士を結べばよい. すると三角形の屋根が12枚の他, 立方体の頂点に対応する三角形が8枚あり, 合わせて二十面体になる.

ここからはご用とお急ぎのないお立ち会いに. 私のかような図の描き方をご披露する.

上の図で, 私は手前をx, 右手をy, 上をzと見ることにしている. そして手前を左に寄せ, 上を手前に引き, x軸方向から見た位置に各点を移して描く. 左に寄せる角度をθ, 手前に引く角度をφとすると, 画面の左右方向Xと上下方向Yについては

X=- x sin θ + y cos θ+ 300
Y=(-x cos θ - y sin θ) sin φ+ z cos φ + 400
(300と400は原点を適当に移動)

だから, x, y, zを変換(map)し, X, Yに移動し(moveto), または線を延す(lineto)命令をmm, mlとすると

/phi 15 def /theta 15 def
/map {3 dict begin /z exch def /y exch def /x exch def
x theta sin mul neg y theta cos mul add 300 add
x theta cos mul neg y theta sin mul sub phi sin mul
z phi cos mul add 400 add end} def
/mm {map moveto} def /ml {map lineto} def


基盤の立方体は直ぐ描けて(稜の長さの半分lは多面体ごとに決める)
/ln l neg def
l l ln mm l l l ml l ln l ml l ln ln ml l l ln ml ln l ln ml
ln l l ml ln ln l ml l ln l ml l l l mm ln l l ml stroke

でよい.

正十二面体では
/l 1 5 sqrt add 2 div 60 mul def
/a 5 sqrt 1 sub 2 div l mul def /an a neg def

とした後
a 0 l a add mm %P
l l l ml %A
0 a l add a ml %R
ln l l ml %B
an 0 a l add ml %Q
a 0 l a add ml stroke %P

正二十面体では
/l 1 5 sqrt add 2 div 80 mul def
/a 5 sqrt 1 sub 2 div l mul def /an a neg def

とした後
an 0 l mm %V
a 0 l ml %U
l an 0 ml %Y
l a 0 ml %Z
a 0 l ml %U
0 l a ml %W
l a 0 ml %Z
0 l a mm 0 l an ml stroke %W X

で描くことが出来る.

数学小景には『正十二面体と正二十面体とは, 自然界が産み能わざる結晶で, それは純然たる「脳産物」である. 』の記述がある.

2008年9月14日日曜日

八進法算盤

下の図は, 私が1960年代の始め頃使っていた, 八進法の算盤である.



当時の計算機はIBM7040で, これはIBM704と同様のアーキテクチャゆえ, 命令語もデータも, 1語36ビットであった. 命令語はその36ビットを左からprefix(3ビット), decrement(15ビット), tag(3ビット), address(15ビット)の部分に分けた. このaddressとdecrementがLispでいうcarとcdrの語源になっているのは周知のとおり. (car=content of address part of register, cdr=content of decrement part of register) 各部分が3の倍数のビット数なので, 命令語の表示は八進法だったのである.

従って, コンソールの前で, 八進法の計算(加減算だけだが)をすること, しばしばであった.

その時考えたのは, 市販の算盤は十進法だが, 下の1の珠を1個外せば, 八進法の算盤として使えるということだ. 早速12桁の算盤を売っていないかと探したが, 市販の算盤はすべて奇数桁であり, 12桁のものはなかった. やむを得ず13桁の算盤を買い求め, 糸鋸などで縮めたのが, 上掲の算盤である.

通常は上が5の珠, 下が1の珠だが, 八進法算盤では, 上が4の珠となる. それ以外は十進法での計算と余り違わない.

この算盤は多いに実用的であった.

引き出しにまだ残っていたので, 実物の写真も載せておく.


八進法算盤

私が小学校で習ったときは, 下の1の珠は4個であったが, 私の母は1の珠が5個の算盤を使っていた. そのずーっと前は, (中国だけかも知れないが) 上の5の珠も2個あったらしい. 1の珠5個で十進法の計算をしていたことを考えると, なにも下の珠を外さなくても, 八進法で計算出来たわけだが, 実際にやってみると改造算盤には及ばない.

余談だが, 昔は病院にいくと, 4号室とか19号室とかはなかった. 4と9は「死」と「苦」を連想するので, 使わなかったのだ. 0,1,2,3,5,6,7,8の8個の数字を使うから, これも一種の八進法である. 私の算盤で, 上の珠を5と解釈すると, この八進法になる. 私はこの八進法を通常のと区別するため, 迷信法といったことがあるが, 東京女子大学の水谷先生がこの命名に興味を持たれたらしい.

最近のアーキテクチャは, 十六進法が多いので, 八進法算盤はお蔵入りになっている. 十六進法算盤はどういう構造にすべきか, 上下それぞれ3個というのもありかも知れないが, 操作出来るだろうか. 目もくらくらしてくる.



ところで, Hewlett PackardのHP16Cが市販されてから, 八進法, 十六進法の計算はもっぱらHP16Cを利用している.

HP16C

2008年9月3日水曜日

和綴じ本

8月にSilicon Valleyに行った折, Stanford大学のDon Knuthさんの家に寄った. 私は数年前からThe Art of Computer Programming(TAOCP)の監訳をしていて, エラーを見つけてはメイルで報告し, 先方からは航空便でコメントが送られてきたりしていたが, 会うのは初めてであった

TAOCPの第4巻は, ペーパーバックの分冊(Fascicle)の形で, これまで分冊2, 3, 4, 0が刊行されている. これらの分冊は刊行前に, Pre-Fascicleとしてインターネットからダウンロード出来る. 私はPre-Fascicleを和綴じにしていたので, 今回その1冊をDonに見せた. Donは奥さんのJillが和綴じが大好きだと喜び, 階下に伝えにいった. やがてJillが和綴じの本を持って2階へ来, しばらく和綴じの話題になった.




上の図は和綴じを示す. 細い線は本の輪郭, 太い線は綴じ糸を示す. 本の方は隠線消去してある. 糸の実線は外から見える部分, 破線は本の中である. 本には綴じる側に穴を4個所開け, 図のように糸を通して綴じる. この糸の線はEuler Circuitになっている. つまり一筆で書け, 出発点に戻れる. だから1本の糸で綴じられる.

しかしトポロジーの問題ならそこまでだが, 和綴じの場合は糸が節点から外れないように通す必要がある.

和綴じは別に一筆書きの問題ではないから, 糸はそれぞれの穴の中を3回通す. つまり穴に向かって3本の糸が集まるが, どの糸も穴に入る. Eulerの件は, 線が2本増えるから, 各節点が偶点であることに変りはなく, やはり1本の糸で綴じることが出来る.

本の縦の長さをl, 綴じ目から背までの幅をm, 本の厚さをhとすると, 糸の長さは2l+8m+18hで, B5版の場合lは257ミリ, 私はmを10ミリとしているから, hを5ミリとして, 糸は(+ (* 257 2) (* 10 8) (* 5 18))=684ミリ要る. ぎりぎりでは作業がし難いから, まず1メートルくらい用意するのが常識だ.

和綴じのもう1つの技法は, 糸の端が本の中に隠れていることである. 穴のどれかの中に, 始点と終点が潜んでいる.

私流の糸の通り方を描いたのが, 下の図である. 上の図と見比べて欲しい. こういう図をPostscriptでプログラムするのは楽しい.


2008年8月26日火曜日

San Franciscoのケーブルカー

今回はSan Franciscoで少し時間があったので, 久しぶりにケーブルカーに乗った. 前回乗ったのは, 50年前, 最初の訪米の帰路, Los AngelesからSan Francisco経由で羽田へ向かおうとしたが, JALの機材に問題があり, 1日San Franciscoに滞在することになった. 早速ケーブルカーに乗りにいった.

その後San Franciscoは数回訪れたが, ケーブルカーに乗るような時間はなかった.

さて, 以前からケーブルカーにつき, いくつかの疑問があった. ケーブルカーは線路の下を動いているケーブルを, 車両の下に出ているグリップで掴み, 車両を動かす仕掛けである. 駅ではグリップを弱め, ブレーキを掛ければ止る. 動くときはグリップを強める.

問題はPowell通りとCalifornia通りのように, 2本の路線の交差はどうなっているか, であった. しかし案の定, 一方の線では, ケーブルを一旦放し, 交差点を過ぎてからケーブルをグリップしなおすのであった. ケーブルを放すのはPowell通りの方である.



この写真で横方向がPowellで, 左の先がMarcket通りである. このPowell線の向うの線, つまり右側通行でMarcket行きの線路が, 多少下がっているのが分かる. ケーブルカーは交差点でケーブルを放し, 少し過ぎた地点で, 車体が下がってグリップがケーブルの位置になり, そこで再びケーブルをグリップしなおすのである.

手前の線には, 交差点の直前に黄色の線が引いてある. これはグリップマン, すなわち運転手にここでケーブルをリリースせよという目印である.

坂の多いところでは, ケーブルカーは水平になった交差点内で停車し客扱いをする. Powell線は, 交差点を過ぎて, 線路が下がった辺りで停車し, 乗降させるとともに, ケーブルをグリップするようだ.

もう1点は登坂の下の話だ. 水平な交差点の下を通るケーブルは, 登坂の上から引かれると, 何もしなければ持ち上がってしまう. スキーリフトの鉄塔を見ると, ケーブルは大体は鉄塔についているプーリーの上を通るが, 谷間の鉄塔だと, プーリーの下を通っている. リフトの腕は水平に出ているから, ケーブルがプーリーの上でも下でもあまり問題はないが, San Franciscoのケーブルカーは, 真上にグリップがあるから, その辺はどうなっているか.

答えはDepression Beamという仕掛けであった.



左はケーブルカーの進行方向から見た摸式図で, 黒はグリップである. 灰色はグリップされているケーブルの断面である. 赤はDepression Beamでプーリーでオレンジ色のケーブルを押さえている. ケーブルカーが遠いときは, ケーブルはオレンジ色の位置で, 上がろうとするのをプリーが押さえている.

ケーブルカーが近づくと, ケーブルは徐々に押し下げられ, プーリーから離れる. 一方ケーブルカーのグリップが近づくにつれDepression Beamを横に押しだす. 赤の点線のような位置に移動し, グリップは無事に通過出来る.

ケーブルカーが通りすぎると, プーリーは元の位置に戻り, またケーブルも上がってきて,プーリーで押さえられる.

Depression Beamが入っていると思われる部分は, 鉄板で覆われ, 保守しやすいようになっていた.

2008年8月25日月曜日

Lermerの篩機械

先日のことだが, Mountain ViewにあるComputer History Museumに行ってきた. 思い込みの深い多くの計算機に対面できた. ただ残念ながらLehmerの初代の自転車のチェーンによる篩機械はドイツに出張中とかで, 見られなかった. これについて私は以前こう書いた. (bit 1997年7月号)

『「先日のことだが, ...」計算整数論者でカリフォルニア大学バークレイ校(以下UCB)教授だったDerrick Henry Lehmer(以下Dick)は昔のエピソードでもこう語り始めるのを好んだ. 私もひとつ真似をしよう.

先日のこと, Dickは集めてきた自転車のチェーンの部品でチェーンの輪を19個作った. それぞれの輪のリンク数は67までの素数だが, 13までの素数に対する輪は小さ過ぎるので64, 27, 25, 49, 22, 26とした. これをそれぞれ同軸に固定した10枚歯のスプロケットに掛け, 回転数をカウンターで数えながら軸をモーターで回転した. 輪は零の位置(または素数の倍数の位置)にピンを立てた時, カウンターが大きな整数Nになるまで軸を回転し, 停止時にどのピンも真上になければNには67以下の素因数はないわけだ.

また探すべき数が, いろいろな素数での法がいくつかという条件さえわかれば, その位置にピンを立て, ピンと接点, リレー, モーターを組合せ, 条件を満たす数が見つかった時点で回転を停止する機構も工夫した.

軸は300rpmで回転したので, 1日に432万個(=10×300×60×24)の整数を調べることができた. だが, この篩計算機を作った目的のMersenne数, 2257-1の素数性を調べるは$1070年かかるらしい. さよう, Dickがこの篩を試作したのは, 彼がUCBの学部学生であった1926年のことである. この機械の複製はボストンのコンピュータミュージアムにあるが, 展示はしていないらしい. 』

Computer History Museumは1990年代にBostonからCaliforniaへ移動している. インターネットで探すとcomputer history museumに写真は見つかった.



どういう問題を解くか. Lehmerのような整数論者は高尚な問題を解くだろうが, 簡単な例はChinese Remainder Theoremのようなものである. つまり3で割れば1余り, 5で割れば2余り, 7で割れば3余る数はなにか.

3での剰余がp, 5での剰余がq, 7での剰余がrなら70p+21q+15r mod 105が解で, 今の70+42+45=157, 157 mod 105 = 52 となる. この程度なら篩機械はいらない.

2008年7月24日木曜日

楕円アルゴリズム

MITのHakmemには意外な話題がたくさんある. プロッタで「円を描くアルゴリズム」もその1つだ.

newX=oldX-epsiron*oldY
newY=oldY+epsiron*newX

下の式で, newXを使うところがなんとなく怪しげなところ. さっそくPostscriptでやってみる.

200 200 translate
/x 100 def /y 0 def /eps 0.01 def
x y moveto
1 1 628 {pop
/x x y eps mul sub def %x=x-eps*y
/y y x eps mul add def %y=y+eps*x
x y lineto
} for stroke

とやるとちゃんと円が描ける.

一方, yの計算に昔のxの値を使うには, Postscriptのスタックを活用し

/x x y eps mul sub /y y x eps mul add def def

とすると, おかしいことになる. やはり新しいxが正解である.

epsの値を徐々に大きくすると, 円ではなく, 45度傾いた楕円を描くプログラムだったことが分かる. そこで楕円の長軸と短軸の長さを計算し重ねてみた.

200 200 translate
/x 100 def /y 0 def /eps 0.5 def
x y moveto
1 1 6.28 eps div {pop
/x x y eps mul sub def
/y y x eps mul add def
x y lineto
} for stroke
gsave
1 0 0 setrgbcolor %赤にする
45 rotate %45度傾ける
1.154700 0.894427 scale %長軸と短軸の長さにscaleする
100 0 moveto 0 0 100 0 360 arc stroke
grestore
0 0 1 setrgbcolor
-100 0 moveto 100 0 lineto
0 -100 moveto 0 100 lineto stroke


(100,0)から出発したときのepsに対する長軸と短軸の長さと離心率を数値計算で求めると

eps 長軸 短軸 離心率
0.5 115.47 89.44 0.6325
0.1 102.60 97.59 0.3086
0.01 100.25 99.75 0.0997

のようになる. 要するに「円を描くアルゴリズム」ではなく, 「楕円を描くアルゴリズム」であった.

2008年7月14日月曜日

3 not problem

私が3notの話を最初に聞いたのは, 1960年代の始め頃, サバティカルで東大に滞在していたイリノイ大学のDavid Muller先生からだ. 研究室のお茶の時間に「アメリカの計算機科学科でよく話題になる問題」の1つとして紹介された.

この話が記述されているのはあまり見た記憶がないが MITのHakmemには出ている. item 19の2-NOTs problemがそれである.

x,y,zの3入力, x',y',z'の3出力を持つブラックボックスがある. 入出力の関係は

x'=not x
y'=not y
z'=not z

である. ブラックボックスには, andとorは好きなだけ使われているが, notは2つしかないことが分かっている. 内部はどうなっているか.

Muller先生からこの話を聞いたわれわれは, その日は他の仕事がまったく手につかなかった. 解けてみるとなるほどと目から鱗の問題であった. 要は3は2ビットで表現できるということ.

Marvin Minsky先生の「Computation: finite and infinite machines」には閾値素子によるヒントが書いてある. 最近こういう話題は流行らないみたいで残念だ.

2008年7月9日水曜日

箱根登山電車のスイッチバック

「箱根の山は天下の険」へ行ってきた. あいにく停滞中の前線のせいで, 「雲は山をめぐり 霧は谷をとざす」という状態であり, 景観は今一つであった. まぁ議論の合宿なので, 景観はどうでもよい.

合宿からの帰路, 箱根登山電車で下山した. あじさい電車としてこの季節に有名で, 1万株といわれるあじさいは, 時として車体に触れるほどであった.

この線は名前通りに急勾配の登山電車で, 全線単線. またスイッチバックでも知られている. スイッチバックは昔は中央本線にも数多くあった. 甲府より東京側と甲府以遠とで, 構内配線の基本構造が違っていたが, 電車列車が多くなり, スイッチバックの駅も激減した.

スイッチバックと単線の駅構内配線をおさらいする.


スイッチバックには通過型と折り返し型がある. 図(A)は通過型で, 細い横線は上が低く, 下が高い等高線を示す. 太いのが線路だ. 通過型では, 坂を登り降りする通過列車はスイッチバックを無視し, LからHへ, またはHからLへ直進する. この駅(または信号所)に停車する列車は, Lから来た場合, 図の右方の行き止まり線へ進入し, 左方の行き止まり線へ逆行する. いずれかの行き止まり線にホームがある. その後H方向へ出発する. この型のスイッチバックは, 停車中の列車を引き出すのに, 水平部分の線区が必要ということで設けられた. 折り返し線は等高線に沿っている. 行き止まり線に入ったり, 逆行したりで面倒そうにも思えるが, ボッボッボーと逆行の汽笛を鳴らし, いとも簡単に戻り始める列車を楽しめた.

中央線のはこの型であった. 残存するのの1つに篠ノ井線の姥捨駅がある.

図(B)は折り返し型で, 全列車が, 折り返し線を使って登り降りする.

箱根登山鉄道とともに山岳鉄道の代表であった草軽電鉄の上信両州の境, 国境平駅のすぐ北に二度上という駅があった. 坂が急なので, 荷物を半分にし, 2回に分けて上げたという地名であるが, 草軽電鉄も二度上で折り返し型スイッチバックした.

この度の箱根登山線では, 金魚ばちといわれる運転席の後に陣取り, 箱根登山の線路の様子を見てきた. この線には, 両終点を除き, 信号所も含めて, 通過型の駅と折り返し型の駅がある. 折り返しがスイッチバックを兼ねる. それらは図に示すと, 通過型は図(A), 折り返し型は図の(B), (C)であった. (B)と(C)は左右が逆なだけで, 機能的には同じと見てよい. (B)は上大平台と出山の信号所, (C)は大平台であり, 彫刻の森の他は, 仙人台信号所を含めてすべての通過型の駅の構造は(A)である. (D)については後述する.

Hは山の上(強羅)側, Lは下(湯本)側である. 図にはポイントがいくつも示してあり, 脇にSとあるのは, スプリング型のポイントで, 定位の方の線がつなげて描いてある. つまり対向に進むとき(進行方向へ分岐している)は, 実線がつながっている方向へ進行する. 背向に進むとき(左右からの線路が合流する)は, 線が切れていても(ポイントが反対を向いていても), 列車は車輪のフランジで無理に押し分けて進む. 列車が去った後, スプリングにより, ポイントは定位に戻る.

線路の横にある短い実戦はホームのつもり.

さて(A). 山から降りてきた列車は, 対向のスプリングポイントで左へ分岐し, 上り列車用のホームへ進入する. 下から登ってきた列車も, 同様に左へ分岐し, 下り列車用のホームへ進入する. 駅の構造は, 進入してくる列車が入れる方を定位とする.

出発すると, 再び単線区間に入るので, 通常は出発方向の先に安全側線が設けられるが, どの駅も下向きにしか安全側線はなかった. あえぎながら登ってくる列車は, 暴走して通過する心配はないと考えたのであろうか.

折り返し型の駅の構造も面白い. 原則として, 進入する列車は渡りを通らず直行する. 上り下りの列車が同時に進入できるためであるのはいうまでもない. この構造だと, 駅から出発する両列車は, 同一区間を同一方向に走ることがあるので, その先のポイントの転轍操作は必要になる.

2つのわたり線はなぜダイアモンドクロスになっていないかと不思議に思うのは当然だ. そうすると図(B)を描きなおした図(D)のように4つのポイントはすべてスプリングにすることが出来るからである. つまりHからの列車は最初の対向ポイントSを直進, 次の背向ポイントでも直進, ホームに入る. 逆行してからは, 最初の対向ポイントで渡り, 次の背向ポイントでL方向へ入る. Lからの列車も同様である.

このことは昔石本祐吉君と話したことがあり, その後石本君が出版した「鉄のほそ道--写真で綴る線路の話--(アグネス技術センター 1996)」でも触れられている. 石本君の解釈は, ダイアモンドクロスは, 揺れがひどく, それを避けたのかも知れないというのである.

2008年7月2日水曜日

Lisp 50歳

織田信長の 人間五十年 下天のうちを比ぶれば 夢まぼろしの如くなり は有名だ. ところで最近計算機のまわりで50年を迎えるものが目立つ.

パラメトロン計算機 PC-1 は今年3月26日に完成50年を祝った. 来年1月のプログラミング・シンポジウムは第50回だ(49年ということ). 同じく1960年に設立された, 情報処理学会も50周年の事業を準備中と聞く.

ところで今度は本年10月20日に, Lisp50歳の誕生日を祝うというニュースが来た. 計算機に比べて, プログラム言語はいつが誕生日か判定が難しいところだ. 上のニュースによると, 1958年10月にMcCarthyが発表した論文にLISPという名前が始めて使われたとある. 20日というのは, 丁度そのころOOPSLAが開催されているかららしい.

McCarthyのLispの論文といえば, Recursive Functions of Symbolic Expressions and Their Computation by Machine, Part Iが知られていて, CACM, 1960年4月号に載ったもので, McCarthyのウェブページにも, This was the original paper on LISPと書いてある. これはそのウェブページからダウンロード出来る. しかしこれが1958年10月に出た論文と同じかは不明である. 似たような論文がいくつかあるからだ.

1958年の11月頃, 私はMITを訪れた. その折, MITに滞在中の藤村先輩から面白いからといってMcCarthyの論文を渡された. それがひょっとして1958年10月の論文かも知れないと思っている. (CACMに載った論文の紹介は, 情報処理学会誌44巻3号にある.)

最近あまり聞かなくなった, Algolも最初に発表されたのは, 1958年であったから, 同様にお祝いされてもいいはずだが, すぐに改良版Algol 60としてより広く知られるようになり, Algol 58の名前を覚えている人は少ない. そのAlgol 60も今や風前の灯火で, AlgolはJISからも抹殺された. しかしAlgol 60の精神はScheme(の特にRevised Report)に受け継がれている.

OOPSLAの会議では「next 50 years of Lisp」についても議論するという. 50年後も夢幻にならずに使われているだろうか.

2008年6月30日月曜日

高橋先生

堀口大学の詩に

「待つ間の長さ 会う間の短さ 時のおなかは蛇腹です」

のようなのがあったと思う.

今日は6月30日 茅の輪くぐりの日だ.

ことしは月曜だが, 1985年と2002年の6月30日はともに日曜であった.

1985年6月30日は高橋秀俊先生が, 2002年6月30日は高橋延匡先生が他界され, 私にとっては忘れられぬ日だ.

どちらもつい先日のようだが, それからそれぞれ23年と6年が経った.

高橋秀俊先生は1985年1月末に体調をくづされ, 入院生活を送られた. 始めの頃はこの病室からは電車がよく見えてということだったが, 段々と重体になられ, 亡くなる前日, 私を病院に呼ばれた. 小さい声でいろいろ話された. 翌日夕方, 逝去されたと電話を受け, 雨の中, 病院に向かった.

その辺のことは, コンピュータソフトウェアの2巻4号に書いた.

高橋延匡先生の方は急であった. 3月の情報処理学会の全国大会でお目にかかったと記憶するが, あっという間のご逝去であった. HITAC5020やアクレディテーションの活動など, 一緒だったので, 残念なことであった.

ちょうど学会誌の編集長だったので, 会誌にこう書いた.

元副会長 拓殖大学教授 高橋延匡先生, 二豎膏肓に入って療養ご専念中とかねてより聞き及びおりしところ, 紫陽花雨に打たるる6月30日, 不帰の客となられる. 無念千万. 翌7月1日は満69歳の誕生日のはずであった. OSの開発に, 情報の教育に, またアクレディテーションの準備にと終始真摯に立ち向かわれた先生のお姿は永く語り継がれよう.

今日は梅雨の晴れ間, 今年の紫陽花はすでに色あせているが, 毎年6月30日になると, 故人を偲ぶことにしている.

2008年6月25日水曜日

夏至の日に

今年の夏至ももう過ぎた.2至2分の中では私の一番好きな日だ. (2至2分 = 夏至, 冬至, 春分, 秋分 これらの英語は難しい. summer solstice, winter solstice, vernal equinox, autumnal equinox) 今の家では, 6月になるとどこからか郭公が来て, 近くの高い木やアンテナの上で鳴く. この2年くらい郭公を聞かなかったが, 今年は何日か鳴き声が聞けた. 「藤の花咲けども鳴かぬほととぎす雲のいずこに迷いけるらむ」を思い出す季節でもある.

6月12日の朝日新聞の夕刊に, 昭和基地で撮影した「幻の日の出」なる写真が出ていた. 南極は白夜と反対で極夜の最中だが, 地平線の直ぐ下にいる筈の太陽が蜃気楼で地上に見えたという図である. 地平線の下ですれすれだとそういうことも十分ありうる.

遥か昔「南極のスコット」という映画を見た. Scott隊が南極からの帰途, 全員遭難した話である. (http://en.wikipedia.org/wiki/Robert_Falcon_Scott) その後極夜になり, 冬至が近づいて再び太陽が戻り, 捜索隊によって遭難の様子が判明する話だ. アナウンスでは「太陽が戻ってきた」といったと記憶するが字幕には「夏が来た」とあった.

北半球の夏は南半球の冬だというのか, 南半球の冬は暑く, 夏は寒いというのか浅学にして知らない. (私は後者の説を支持する.)

夏至は24節気の1つ. 森口繁一先生の「数理つれづれ」に24節気間の時間の図がある. 先生が描かれた方法は分からぬが, 同じようなものを描いてみた. 本年の24節気の月日時分を天文年鑑から取り込んだのが以下のデータである.


( 0 22 15 8 '冬至) ( 1 6 8 25 '小寒) ( 1 21 1 44 '大寒)
( 2 4 20 0 '立春) ( 2 19 15 50 '雨水) ( 3 5 13 59‘啓蟄)
( 3 20 14 48 '春分) ( 4 4 18 46 '清明) ( 4 20 1 51 '穀雨)
( 5 5 12 3 '立夏) ( 5 21 1 1 '小満) ( 6 5 16 12 '芒種)
( 6 21 8 59 '夏至) ( 7 7 2 27 '小暑) ( 7 22 19 55 '大暑)
( 8 7 12 16 '立夏) ( 8 23 3 2 '処暑) ( 9 7 15 14 '白露)
( 9 23 0 45 '秋分) (10 8 6 57 '寒露) (10 23 10 9 '霜降)
(11 7 10 11 '立冬) (11 22 7 44 '小雪) (12 7 3 2 '大雪)
(12 21 21 4 '冬至)

これをユリウス日を計算するプログラムに入れ, 差を取ると(単位は日)

14.720139 14.721527 14.761112 14.826389 14.922916 15.034028
15.165278 15.295138 15.425001 15.540277 15.632639 15.699306
15.727778 15.727777 15.681251 15.615277 15.508334 15.396528
15.258333 15.133333 15.001388 14.897917 14.804167 14.751389

となり, これをpostscriptで図にする.




地球の軌道が真円でないので, 太陽が(地球が)黄道上を移動する時間が変わるからだ. ところで夏至では昼間の時間が一番長いが, 日出, 日入が一番早く, また遅いのは夏至の日ではない. この説明は情報処理学会誌45巻3月号の編集系独白に書いたので, 読まれた方もあろう.

2008年5月30日金曜日

天文時計の鐘

N.J.A.Sloaneの書いたMy Favorite Integer Sequencesという記事を見つけた. 面白いことがいろいろある中で, あ!と思ったのは

1, 2, 3, 4, 32, 123, 43, 2123, 432, 1234, 32123, 43212, 34321, 23432, 123432, 1234321, 2343212, 3432123, 4321234, 32123432, 123432123, 43212343, 2123432123, 432123432, 1, 2, 3, 4, 32, 123,

の数列であった.

説明によると, これはプラハの天文時計(astronomical clock)の鐘の鳴り方だそうだ.

数列の先頭からこのように, 1時には1; 2時には2,...,と鳴り, 正子(0時)には24時として432123432と鳴って, 1時には再び1と鳴る.

この数列の特徴は,

a) n 時に鳴る数列の和はちょうど n になる.

b) カンマを外してみると, 数列は基本パターン"123432"の繰り返しである.

ということだ.

1,2,3,4と書いてあるが, 鐘は1つで, KnuthのTAOCP, v4f2, 7.2.1.2にあるCambridge Forty-Eightのような, 4つの鐘をある順列で鳴らすのではない.

1は1回, 2は2回鳴り, 5時の32はまず3回鳴り, 次に2回鳴る. それを足しながら数えれば時刻が分かる仕掛けだ.しかし24まで間違えずに数えるには, 忍耐もいるであろう.

なぜこういうことが出来るか.

実は, この不思議な数列は24まで続くだけでなく, どこまでも作れる. 計算機が手元にあれば, お茶の子だ.

(1 (1))
(2 (2))
(3 (3))
(4 (4))
(5 (3 2))
(6 (1 2 3))
(7 (4 3))
(8 (2 1 2 3))
(9 (4 3 2))
(10 (1 2 3 4))
(11 (3 2 1 2 3))
(12 (4 3 2 1 2))
(13 (3 4 3 2 1))
(14 (2 3 4 3 2))
(15 (1 2 3 4 3 2))
(16 (1 2 3 4 3 2 1))
(17 (2 3 4 3 2 1 2))
(18 (3 4 3 2 1 2 3))
(19 (4 3 2 1 2 3 4))
(20 (3 2 1 2 3 4 3 2))
(21 (1 2 3 4 3 2 1 2 3))
(22 (4 3 2 1 2 3 4 3))
(23 (2 1 2 3 4 3 2 1 2 3))
(24 (4 3 2 1 2 3 4 3 2))
(25 (1 2 3 4 3 2 1 2 3 4))
(26 (3 2 1 2 3 4 3 2 1 2 3))
(27 (4 3 2 1 2 3 4 3 2 1 2))
(28 (3 4 3 2 1 2 3 4 3 2 1))
(29 (2 3 4 3 2 1 2 3 4 3 2))
(30 (1 2 3 4 3 2 1 2 3 4 3 2))

そこでどこかに繰り返しパターンがないかと, 図を描いてみた. すると1から15までのパターンと16から30までのパターンが殆んど同じになる. 16からは基本パターン1,2,3,4,3,2を1サイクル余計に持っているだけである. 従って31からも同じになると判明した.



10,11,12,13,14はそれぞれ15のcomplementの5,4,3,2,1と相補のパターンだし, 6と9, 7と8も相補パターンなことが分かる.

24がちょうどパターンの最後で終るのも幸いしている. 1から24までの和が300で基本パターンの和15の倍数なので, 当然だが.

1,2の繰り返しでも同様なことが出来ることが判明した.

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

しかし, 1,2,3,4,3,2ほど劇的ではない.

2008年5月28日水曜日

火星

火星探査機がパラシュートで降下するのを, 他の火星周回衛星から撮影したニュースには驚いた.

いかにも画面を眺めて, 今だ! とシャッターを切ったかのようなシーンだが, そんなことは可能だろうか.

火星は昨年12月に衝の位置(火星, 地球, 太陽の順に並ぶ)にあり, 地球, 火星間の距離は0.6天文単位, 今年12月に合の位置(地球, 太陽, 火星の順に並ぶ)になり, 2.6天文単位だ. 今はその中間で, 多分2天文単位くらいであろう. 1天文単位は光で8分ちょっとの時間なので, そう考えると画面を眺めて地上から制御するのはまったく無理と分かる.

NASAのページを見ると, やはりフェニックスの軌道情報を送り, プログラムで制御して撮影したとあった. それにしても相当な離れ業である.

Camera pointing for the image from HiRISE used navigational information about Phoenix updated on landing day. The camera team and Phoenix team would not know until the image was sent to Earth whether it had actually caught Phoenix.

さらに普段は火星に向けて下を向いているカメラを, 斜め方向を撮影するため, 衛星の向きも変えたという. 見事な写真が送られてきたときの, 担当者の喜びが想像できる.

2008年5月27日火曜日

火星

太陽系の中で, 地球と火星は双子の惑星といわれる. 公転周期は地球が1年, 火星が2年. 公転半径は地球が1天文単位, 火星が1.6天文単位. (2年の方を知っていれば, Keplerの法則: 公転周期2/公転半径3=一定 から

(expt (expt 2 2) (/ 1 3)) => 1.5874010519681994)

だが, Bodeの法則を知っていれば距離は一発だ.)

この火星には火星人がいると思われていた. 運河らしいものが見えるからだ. 運河かどうかは分からないが, 火星の表面には多数の筋が見え, それが十字形で交差している. 一方, ガラスの割れ目など自然に出来る筋は, 新しい割れ目が既に割れている割れ目を横切らないから丁字形になる.

というわけで, 火星人襲来というSFはいくつかあるらしい. 私が小学生のころ愛読したのに, 海野十三の「火星兵団」がある. ある晩, 千葉県印旛沼の辺りに火星からのロケットが着陸. 火星人が地上に降り立つ.

H.G.Wellsによると, 火星人は, おむすびのような頭から, 糸のような足が生えている, 細い足を持つ蛸のようだと思われていた. 重力が小さいから細い足で支えられる. 食べ物は苔みたいだから, それをなめる口はおおきい. 空気が薄いので音声は伝わらず, 目が口のようにものを言うので, 目もおおきい. 耳は退化した.Wellsはそういう火星人を想像した.

私が小学生のころ, 上野の科学博物館で, こういう火星人を描いた絵はがきを売っていた.

「火星兵団」によると, 火星にくらべて空気の濃い地球では, 潜水服か宇宙服のような容器を着て, 地上で活躍した.

火星人は「ヒュウーヒュウープクプク」というような言葉を喋る. そのままでは日本人には理解出来ぬが, 彼らはすでに日火両語の音声翻訳機を持っていて, 一方から火星語を入れるともう一方から日本語が出る. 逆ももちろん可能.

まぁそういう話であった.

さて今回, NASAの火星探査機フェニックスが火星の北極近くに着地(着火?)したというニュースだ. 送られてきた火星の地上(火上?)の光景は, 運河は見えねど, 想像とあまり違わない. (いちいち「地」を「火」と言い直すのも面倒だねぇ.)

こう描きながら想像するのは, 火星人はフェニックスをみて, 「これが地球人か」と思ったかということだ. 「地球人襲来」というニュースが流れているかと夏の夢を楽しんでいる.

2008年4月24日木曜日

同一直線上の3点

久野君たちの努力により, Beautiful Codeの翻訳がでた. 以前 三省堂の洋書の棚にあるのを見たことはあったが, その時はパスした.

翻訳をみると, なにしろ多くの人がそれぞれのプログラム言語で書いた自分のプログラムを(それもかなり大きい部分を)自讚しているから, 読むのが大変そうである.

短くて面白かったのは, 33章「『本』のためにプログラムを書く」であった. 要するに平面上の3点A, B, Cの座標が与えられた時, その3点が同一直線上にあるかを判定するプログラムを書くのだ.

私がやっても多分こういうアプローチになるであろうという風に話は展開していく.

まずA,Bの2点を通る直線の式を決め, 点Cがそれに乗っているかを問うもの. これは最初の2点がy軸と平行な線上にあるときの始末が面倒.

次はABを通る直線の勾配と, ACを通る直線の勾配を計算し, それらが一致するかを見るもの. 勾配が無限大になるときはnilを返すようにすると, nil同士もeqで比べればtになるので便利なようだが, 気持は悪い.

3番目はBCの距離a, CAの距離b, ABの距離cを計算し, a,b,cを大きさの順にソートし, 最大の距離が残りの2つの距離の和になるかどうかを見る. この難点はPythagorasの定理で距離を求めるのに関数sqrtを使うことだ. 誤差は避けられない.

このようにして, 最後に辿り着いたのは, 3点で構成する三角形の面積を求め, それが0なら三角形はつぶれて3点は同一直線上にあることが分かるというもの.

たしかに高校生のとき, 三角形 ABC の面積は行列式

|Ax Ay 1|
|Bx By 1|
|Cx Cy 1|

の1/2と習った. 面積が0かとうかを見るのだから, 1/2はいらない, 正負も問題にならないから, A,B,Cを反時計まわりにおくことにこだわることもない.

(Ax*By+Bx*Cy+Cx*Ay)-(Ax*Cy+Bx*Ay+Cx*By)

の計算だから誤差も何もない. これで決りだ.

今年の夏のプロシンのテーマは「プログラムの品格」のようなものだ. 本書もひとしきり話題になりそうな予感.

2008年4月14日月曜日

dancing links

前回のdancing linksの続きである.

TAOCP 4巻分冊0が4月18日に発売になるとアマゾンからメイルが来た. その分冊0は7章の最初の部分で, そこにGraeco-Latin squareの話題が登場する.

Graeco-Latin squareの作り方は, 互いにorthogonalな2つのlatin squareを用意し, それを合わせればよいのだが, あるlatin square Lからそれとorthogonalなlatin square Mの作り方として, Lからtransversalという組を探す. これは直接orthogonalなlatin squareを作るより容易であると本文に書いてある.

その続きは次のようだ.

Once the transversals are known, we're left with an exact cover problem, of 10 stages, which is much simpler than the original 90 stage problem (6). All we need to do is cover the square with ten transversals that don't intersect --- because every such set of ten is equivalent to a latin square M that is orthogonal to L.

10とか90とかの数値は, ここでは10x10のsquareを求めているからである.

しかしこれだけの説明では, よく分からなかったので, 考えた結果がこのblogの趣旨である.

大きいsquareは面倒なので, 4x4を対象とする.
例えば Lとして

L=((0 3 1 2)
(2 1 3 0)
(3 0 2 1)
(1 2 0 3))

は各行 各列に0,1,2,3が1回ずつ現れるから, latin squareになっている. これとorthogonalなlatin squareを求めるべく, Lのtransversalを探す. transversalはLの各行から1個 各列から1個 各文字から1個をとる組で,

A=((0 ) B=(( 3 ) C=(( 1 ) D=(( 2)
( 1 ) (2 ) ( 0) ( 3 )
( 2 ) ( 1) (3 ) ( 0 )
( 3)) ( 0 )) ( 2 )) (1 ))

E=((0 ) F=(( 1 ) G=(( 2) H=(( 3 )
( 3 ) (2 ) ( 1 ) ( 0)
( 1) ( 0 ) (3 ) ( 2 )
( 2 )) ( 3)) ( 0 )) (1 ))

の8個が存在する.

このうちA,B,C,Dを重ねるとちょうどLになり, またE,F,G,Hを重ねてもちょうどLになる.

定義が後回しになったが, LとMがorthogonalというのは, LとMから対応する位置の要素をとって組にすると, 同じ値の組が複数回は出来ないということである.

したがって各transversalの0,1,2,3の位置には0,0,0,0のように同じ値をおけば, orthogonalなlatin squareが得られる. つまり

A,B,C,Dから
M0=((0 1 2 3)
(1 0 3 2)
(2 3 0 1)
(3 2 1 0))

E,F,G,Hから
M1=((0 3 1 2)
(1 2 0 3)
(2 1 3 0)
(3 0 2 1))

ができる (LM)の組を作ると M0とM1のそれぞれから

((00 31 12 23) ((00 33 11 22)
(21 10 33 02) (21 12 30 03)
(32 03 20 11) (32 01 23 10)
(13 22 01 30)) (13 20 02 31))

が得られorthogonalなことが分かる.

Lからdancing linksでtransversalを得る方法, transversalを使わず, 直接orthogonalなlatin square(Mateという)を得る方法はTAOCP ex7-17参照のこと.

2008年4月9日水曜日

dancing links

KnuthのTAOCP, 7巻の分冊が出回り始めた. その最初の方にexact cover problemを解くのにdancing linksが適してると書いてある.

dancing linksについてはKnuthが楽しんで書いた論文があり, 講義のビデオも存在する. exact cover proble (日本語では「敷き詰め問題」か) は条件をならべた行列から解の候補を選び, それにより使えなくなった行や列を外し, またback trackするときは, 外した行や列をもとにもどす.

この行列はスパースなので, 上下左右に双方向リンクで実装すると, うまくいくというのがdancing linksの話の味噌である.

私としては面白そうなアルゴリズムを知ったら, さっそくコーディングして例題を走らせ, 体験によってアルゴリズムを理解, 記憶することが多い.

dancing linksについてもさっそくプログラムを書いた. googleで探すとjavaなんかで書いたプログラムも見つかるが, 自分で書くのが楽しいのでschemeで書くことにした.

Knuthの論文の最初の簡単な例題は, 直ぐにうまくいったが, 当然dancing linksを使えばうまく解けると書いてある他の例は, 結構てこずった.

8×8-2×2の60区画にペントミノのピースを置く問題である. ピースXを左上の(3,3)に置くというのに挑戦た. ぜんぜん解が出てこないのである. 当然プログラムがおかしいかといろいろ見てみたが, 合っているとしか思えない.

実は候補の列を選ぶのに, ブランチ数の少ないものを選ぶと書いてあったが, そうしなくても解が出るようにも読めたので, 残っている列の最左端から選んだのが失敗であった. あるとき思いついて, ブランチ数最小の列を選んでみたら, たちまち最初の解が出た. 超簡単な例では, 最左端でも問題なく解が得られたが, ペントミノ程度になると, ブランチ数を考慮する必要があった.

dancing linksの双方向リンクによる削除挿入は, 野下浩平君とその学生の一松宏君が8クイーン問題を解くのに考えたといわれている. しかし縦と横はすべてにクイーンを置くので, exact cover problemであるが, 斜め方向は半分くらいしか埋めないので, この辺の扱いは別になる. Knuthの例の論文をよく読むと, こういうのはgeneralized exact cover problemというので, 列にはprimaryとsecondaryを用意する. primaryの列名は双方向リンクで接続するが, secondaryの列名は自分自身で循環するリンクにすると書いてあった. そういう風にプログラムをなおすと, なるほどうまく行った.

数独にも適しているといわれるdancing linksである. 数独ではいくつかの解が与えられているので, 解を探し始める前に, 与えられた解について列や行を削除する(coverする)必要がある. 解に相当する行を探す手間が必要で, この辺はプログラムでかなり工夫が必要であった.

最後に, 上下左右に双方向リンクにしないでも, 実装できるのではないかとプログラムをしてみた. ノードはlispの対で, carに下へのリンク, cdrに右へのリンクを置き, 循環するようにリストを実装すると, たしかにこれでもうまく出来ることがわかった.

2008年3月23日日曜日

中山道

以前から少しずつ歩いていた中山道533Kmも遂に京都三条大橋へ辿り着いた. (詳しくいえば草津宿から大津宿を経て三条大橋までは旧東海道) 一緒に歩いてくれた多田君, 寺田君, 岩崎君に感謝したい. 実はまだ日本橋から本郷東大前を通り埼京線板橋駅まで8Kmが残っているが, これは車やバスで幾度も通り, 勝手知ったる道なので, もう済んだようなものでもある.

中山道は人気があり, 解説書やウェブページも多く, 歩こうと思えば問題は少ない. われわれは「ボクたちが歩く中山道」を手引きにした. また歩き出してから気づいたが, 道中処々に赤いシールが貼ってあり, 分岐点などで重要な道案内になっていた. 誰が貼ったのか知らないが有り難かった.

トラックがびゅんびゅん走る国道を歩くこともあったが, だいたいは静かな田舎道で, スローライフそのものである. それと引換えに山道では, 宿泊場所が左程ない; 食堂は国道を横切るときに探すしかない; という制約があることも判明した.

それにしても春秋のよい季節にのんびり歩けるのは最高であった. 昔の人は16日で歩き通したらしいから, 1日あたり33Km程度の速さである. これではのんびり歩く気分ではなかったを思われる. (中山道最大のイベント和宮東下でも25日) われわれは半日行程では10~15Km. 1日行程では20~25Kmであった.

昔と違うのは電車線と併走している場合, 草臥れると最寄駅から電車で宿泊地へ向かい, 翌日また電車で戻って歩き出せたことであった. 追分~下諏訪, 大井~太田のような山道区間を除けば行程はかなりフレキシブルであった.

自治体により中山道の周知は様々であった. 「歴史の道 中山道」の道標が立っているところもあれば, 全く無関心な町もあった. 中山道の保存にもっと尽力して欲しいものである. 中山道を歩く人がさらに増えるように.

2008年3月12日水曜日

3シリンダ機関車

前回書いた3シリンダ機関車の中央シリンダの制御シャフトの続きで, リンク機構はベクトルで考えると明瞭なことが分かった.

120度づつ離れたベクトル a, b, cがある. aに2 to 1 leverを接続すると, 相手の端は逆向きで大きさが半分のベクトル a'になる. これはちょうどbとcの中間なので, a'を中心にしたbの動きはちょうとcになるわけだ.

2008年3月8日土曜日

3シリンダ機関車

時間があったので鉄道博物館へいった. いつも何か発見があり楽しい.

国産の3シリンダ機関車C53の走り装置が置いてあった. 3シリンダ機関車は, 両側のシリンダは普通の蒸気機関車のとそっくりだが, 外から見えない台車の中央に3つめのシリンダーがある.

その中央のシリンダからはクランク軸になっている動輪に力を伝達する. 問題は中央の制御シャフトがどうなっているかであるが, 実に意外な構造であった. 両側の滑り弁を動かす制御シャフトが機関車の先頭の方へ突き抜け, それをリンクでつないで中央の制御シャフトを動かすのであった.

googleで探したら,
http://www.watercressline.co.uk/tw/pics/bitn2to1.jpg
にその図解があった.

こういう計算をしてみた. 一方の制御シャフトの動きをa, 他方のそれをb, 中央のそれをcとする.

a=sin x
b=sin (x+2π/3) = sin x cos 2π/3 + cos x sin 2π/3=-1/2 sinx + √3/2 cos x
c=sin (x+4π/3) = sin x cos 4π/3 + cos x sin 4π/3=-1/2 sinx - √3/2 cos x

これを眺めると a + b = -c とわかる.

たしかにベンツのマークのような, 120度ずつはなれたベクトルをa, b, cとすると, a+bはcと反対向きのベクトルになるから, 納得できる.

さて上の図解によると, 下のvalve spindle linkの動きをfulcrum pinを中心にして2 to 1 leverで上へ伝達する. その点を中心として, 上のvalve spindle linkの動きをequal leverで中央のvalve spindle linkへ伝えている.

従って, 下の動きをaとすると, 2 to 1 leverの上の動きは -1/2 aになる. 上の動きをb, 中央の動きを x とすると, bとxの平均が-1/2 aなのだから (x + b)/2 = -1/2 a. したがって x = - (a + b)となり, ちょうどcの動きとなっている.

すばらしい.

2008年2月22日金曜日

銀河

京都へ出張した帰路, 3月14日で運転をやめる寝台急行銀河を利用した. 以前銀河で下ったときの乗客は半分以下で, 車掌といつまで運転するのかなぁと話したりした. 運転打ち切りの直前ともあって今回は満員であった.

寝ながら聞く鉄橋を渡る響き, すれ違う列車の音, まれに停車する駅の物憂いアナウンス, 「おはようございます. 何月何日何曜何時何分です」という車掌の声で朝を知る. 独特の寝台列車の朝の雰囲気である.

子供のころから寝付かれない夜は今寝台列車で旅していると考えるとすぐ眠れた. 「銀河鉄道の夜」も何度も読んだ.

最近利用する寝台列車は函館からの帰りの北斗星ばかりだが, 利用回数の多いのは銀河であった.

新幹線を併走する路線で寝台列車を利用するのは, 物好きといわれても仕方ない. 能登は残っているし, これからは金沢へ向かう能登と函館から戻る北斗星をときどき利用するだけになりそうだ.

ついでに
英字新聞を読んだ人からアメリカでは洪水で何万という人が寝ていて流されたと聞いたがこのsleeperは寝ている人ではないらしい. sleeperは寝台車であると同時に線路の枕木のことでもある.

2008年2月14日木曜日

knuthのGPS

電通大での講義の時に,GPS(general problem solver)の話をした.Algol 60の名前呼びの機能を最大限活用するものだ.


real procedure gps(i,n,z,v); real i,n,z,v;
  begin for i:=1 step 1 until n do z:=v; gps:=1 end


つまり4つの引数i,n,z,vはvalueと宣言していないからすべて名前呼びである.

これを使うプログラムは

i:=gps(i, if i=0 then -1.0 else i,p,
 if i=1 the 1.0 else
  if gps(a,i,z,if a=1 then 1.0 else
   entier(a)×(entier(i)÷(entier(a))=entier(i)∧a<i then 0.0
    else z)=z then
     (if p<m then=p+1 else i×gps(a,1.0,i,0.0)) else p)

プログラムを実行するとpにm番目の素数が入る

最初の行はiを制御変数にしてfor文をまわす.iは1になっているので終点はiということ,すなわちiはどこまでも増え続ける.forループで実行するのは変数pへの代入で,まずはiが1なので1.0が入る.
次はiが2になって代入文へくる.今回からpへ代入する値はelse以降3行目からになる.
そこでgpsを呼ぶ.今回は制御変数aで1からiまで増える.代入先はzで初回は1.0が入る.その後aは1ずつiまで増える.代入される値は4行目で計算する.ここでentierはAlgol 60の基本関数で小数点以下を切り捨てる. ÷の記号は商の小数点以下を切り捨てる除算. したがってここではiがaで割り切れるかを見ている. 割り切れてしかもaがiより小さければiには約数がある,素数ではない. そうならzにoを入れ,そうでなければzのまま. zには最初1を入れたが途中1つでも約数があると最後は0になっている.
一方このgpsから脱出してきたときの値はgpsの宣言の最後にあるように1である. zがそれに等しいというのはiが素数であったことである. そこで最後の行のelseより前を実行してその値をpに代入する. pがmより小さければp+1だからpはmになるまで素数のたびに1ずつ増える. 素数でなければp, つまり不変.
ということでm番目の素数が見つかったところで中ほどのelseの次に来る. iには今m番目の素数が入っている. gpsの値1を掛けてpに入れる. gpsの仕事はiに0を代入して外側のgpsのループを停止させることである. 問題は最後の乗算でiを先に評価すればよいが, gpsを先に評価すると大切なiの値が変わっていることだ.
まだ他にも微妙な点はあるが,いちおうこのようなプログラムだ.

2008年2月5日火曜日

計算機プログラムの構造と解釈

この冬学期,電気通信大学の大学院で記号処理特論という講義をしました.この学期には「計算機プログラムの構造と解釈」の4章,5章を説明しました.90分の講義15回でカバーするのは結構大変です.

私は実装に興味があるので,使い方の例題はほとんど無視してしまいました.また問題にも触れませんでした.

相当なスピードで進んだので,消化不良の学生さんもいたと思われます.

希望的にはMITのサイトからプログラムをダウンロードし,さまざまな例題でプログラムを走らせ,体験を通して理解して欲しかったと思います.(そうした人がいなかったとは言い切れません)

私個人は4章,5章が1350分で講義できたという経験が貴重です.講義の機会をいただいた電気通信大学に感謝しています.

2008年1月30日水曜日

パラメトロン計算機完成50年

2008年3月26日は,パラメトロン計算機 PC-1が東京大学物理学教室高橋研究室において完成して50年目にあたります.それを記念し,当時の関係者,利用者,PC-1に関心をお持ちの方をお誘いして,当日の午後に講演会,晩にパーティーを開催することに致しました.

この他、PC-1の回路図,ライブラリ,関係書籍の展示も計画しております.

詳細はプログラムページにあります.