2011年5月25日水曜日

曜日の計算

1989年12月19日に59歳で他界された基礎論の数学者, 島内剛一(しまうち たかかず)先生を知る人はもう少ないに違いない. 私はプログラミング言語の設計などで, 島内さんには多くの影響を受けた. 基礎論の専門家だけあって, 計算機の見方も一種独特であり, 島内マニュアルは「計算機は順序のあるビットの有限集合である.」のように書き始める.

島内さんの多くの特技の1つが, 曜日の計算であった. 年月日をいうと10秒くらいで何曜と返ってくる.

その秘法は何か. 私の推測では, 島内さんは,

m 1月 2月 3月 4月 5月 6月 7月 8月 9月10月11月12月
h 2 5 5 1 3 6 1 4 0 2 5 0

のような, 月始めの曜日の表h(m)だけを記憶していたらしい.

すると, 21世紀中の 2000+y年m月d日の曜日は, y+floor(y/4)+h(m)+d+4 mod 7 を計算するとそれが曜日になる.

やってみよう. 2011年5月25日の曜日は, 11+floor(11/4)+h(5)+25+4 = 13+3+25+4 = 45, 45 mod 7 = 3 すなわち水曜である.

だれもが知っているように, 平年ならある日の曜日は, 365 mod 7 = 1 ゆえ, 翌年は次の曜日になる. うるう日をまたぐと, 2日増える. これがyとfloor(y/4)を足す理由である.

うるう年は1月2月も2増やしてしまったが, これでは駄目なので, うるう年の1月と2月はh(m)の値から1を引く必要がある.

最後の+4は21世紀中に通用する定数で, 曜日を合わせるためのものだ. そもそもh(m)の表は全体を一斉に増やしても減らしても問題ない. 補正値をそれに応じて変えるだけである.

要はh(m)の表を覚えやすく作ることだ. 島内さんは似たような表を, おそらく無理に覚えてしまったのだろうが, 私は上の表を工夫をして書いた. mod 7を取るのだから, 7を足しても意味は変らないので, 4月を8に, 9月を7にする.

m 1月 2月 3月 4月 5月 6月 7月 8月 9月10月11月12月
h 2 5 5 8 3 6 1 4 7 2 5 0

この表で, 4月から始まる偶数の月の値は, 8, 6, 4, 2, 0である. 一方, 3月から始まる7月までの奇数の月の値は, 5, 3, 1で, その先は負になるから, 心機一転9月は7に, 11月は5にするのである.

残りは, 経験で知っているように, 2月は3月と同じ, 1月は10月と同じと覚える.

3月は5で, 4月は8の出発の値の記憶が肝心である(9月の7も?). そして2か月で2ずつ減る(2か月の間に, 通常は30日の月と31日の月があり, 30+31=61=-2 mod 7だからだ). それと世紀の定数の+4を忘れないこと.

練習問題
2000年1月1日は何曜か. (うるう年に注意)
9.11の2001年9月11日は何曜か.
来年日本で金環食が見られる2012年5月21日は何曜か.
unix timeの31ビットが立つ, 2038年1月19日は何曜か.

今年ももう年の半ばまで来た. 今年中は, 11+floor(11/4)+4=17 そのmod 7は3だから, 3+h(m)+d mod 7で計算する.

曜日に計算では, 無理に誤差をいれる島内の方法や, John Conwayの方法などが面白いが, それらはまたいずれ.

2011年5月23日月曜日

unix time

前回のunix timeのブログに書いたように, Calendricalの方法は, 確かに簡単ではあるが, 1から勘定を始めるプログラムに, 多少の違和感を持っていた私は, やはり0から始めるプログラムが書いてみたく, そのため, 366, 365, 365, 365のように, おなじものが並ぶ後続の前に, それらより1だけ多いか少ない要素が先頭にあるときの処理について考えていた.

つまり例えば1600年1月1日から1999年12月31日の400年間, 146097日には, 最初の100年は1600年がうるう年なので, 36525日あり, 後の3回の100年は, 最初の00年が平年なので, 36524日であって, 36525, 36524, 36524, 36524であり, またその36524日ある100年を4年ごと25組に分けると, 第2組以降は最初の4年が1日少なく, 1460, 1461, 1461, ..., 1461となる.

このような問題提起について, d日目を入力し, その該当する年yと, その年内の日数d'を計算する方法を探した.

何万日を扱うのも憂鬱なので, とりあえず, 先頭が1多い4日, 3日, 3日, 3日のパターンと, 先頭が1少ない, 3日, 4日, 4日, 4日のパターンを考える.

4,3,3,3のパターンでは, 総日数は13日で, 0≤d<13に対し,

d y d'
0,1,2,3 0 0,1,2,3
4,5,6 1 0,1,2
7,8,9 2 0,1,2
10,11,12 3 0,1,2

また, 3,4,4,4のパターンでは, 総日数は15日で, 0≤d<15に対し,

d y d'
0,1,2 0 0,1,2
3,4,5,6 1 0,1,2,3
7,8,9,10 2 0,1,2,3
11,12,13,14 3 0,1,2,3

となるようにしたい.

前者については

(define (foo d)
(let ((y (floor (/ d 3.25))))
(list y (+ d -13 (floor (* 3.25 (- 4 y)))))))

とする. 13は総日数, 4は年の幅, 怪しげな3.25は13/4である.

後者については

(define (bar d)
(let ((y (- 3 (floor (/ (- 14 d) 3.75)))))
(list y (- d (floor (* y 3.75))))))

とする. 14は総日数-1, 3.75は総日数の15を4で割ったものだ.

とりあえず, これでdに対して(y d')を求めてみると

(map foo (a2b 0 13))
=>((0 0.) (0. 1.) (0. 2.) (0. 3.) (1. 0.) (1. 1.) (1. 2.)
(2. 0.) (2. 1.) (2. 2.) (3. 0.) (3. 1.) (3. 2.))

(map bar (a2b 0 15))
=>((0. 0.) (0. 1.) (0. 2.) (1. 0.) (1. 1.) (1. 2.) (1. 3.)
(2. 0.) (2. 1.) (2. 2.) (2. 3.) (3. 0.) (3. 1.) (3. 2.) (3 3.))

とうまく行きそうである.

fooとbarの図を描いてみるとこうなる. 上がfoo, 下がbar.



図でははっきりしないが, dに対するfloorを取る前の値は
(0 .31 .62 .92 1.23 1.54 1.85 2.15 2.46 2.77 3.08 3.38 3.69)

(3.73 3.47 3.2 2.93 2.67 2.4 2.13 1.87 1.6 1.33 1.07 .8
.53 .27 0)
である.

そこで, この3とか4とかを実際の値に変えて関数を書き, クリティカルな値に対してテストすると次の通りだ.

(define (foo d)
(let ((y (floor (/ d 36524.25))))
(list y (+ d -146097 (floor (* 36524.25 (- 4 y)))))))
(map foo '(0 36524 36525 73048 73049 109572 109573 146096))
=>((0 0.) (0. 36524.) (1. 0.) (1. 36523.) (2. 0.) (2. 36523.)
(3. 0.) (3. 36523.))

(define (bar d)
(let ((y (- 24 (floor (/ (- 36523 d) 1460.96)))))
(list y (- d (floor (* y 1460.96))))))
(map bar '(0 1459 1460 2920 2921 4381 4382 5842 33602 35062
35063 36523))
=>((0. 0.) (0. 1459.) (1. 0.) (1. 1460.) (2. 0.) (2. 1460.)
(3. 0.) (3. 1460.) (23. 0.) (23. 1460.) (24. 0.) (24 1460.))

(define (foo d)
(let ((y (floor (/ d 365.25))))
(list y (+ d -1461 (floor (* 365.25 (- 4 y)))))))
(map foo '(0 365 366 730 731 1095 1096 1460))
=>((0 0.) (0. 365.) (1. 0.) (1. 364.) (2. 0.) (2. 364.)
(3. 0.) (3. 364.))


あとはこれを繋げるだけ. 注意すべきは, Calendricalの方は1年1月1日から始まるのに対し, こちらは0年から399年までの400年から始まるので, 0年分の366日を足す; Calendricalは1月1日のfixed dateが1なのに対し, こちらは, 0年1月1日が0なので1日分の修正が必要だ. この修正の結果がd0である. 400年のところは, 常に146097日なので, 計算は簡単だ. 残りは上の関数を利用している.


(define (gr rd)
(let* ((d0 (+ rd 366 -1))
(y400 (quotient d0 146097))
(d1 (modulo d0 146097))
(y100 (floor (/ d1 36524.25)))
(d2 (+ d1 -146097 (floor (* 36524.25 (- 4 y100)))))
(y4 (if (= y100 0) (quotient d2 1461)
(- 24 (floor (/ (- 36523 d2) 1460.96)))))
(d3 (if (= y100 0) (modulo d2 1461)
(- d2 (floor (* y4 1460.96)))))
(y1 (if (and (> y100 0) (= y4 0)) (quotient d3 365)
(floor (/ d3 365.25))))
(d4 (if (and (> y100 0) (= y4 0)) (modulo d3 365)
(+ d3 -1461 (floor (* 365.25 (- 4 y1))))))
(y (+ (* y400 400) (* y100 100) (* y4 4) y1))
(leap (if (= (modulo y 100) 0) (= (modulo y 400) 0)
(= (modulo y 4) 0)))
(e (if leap
'(0 31 60 91 121 152 182 213 244 274 305 335)
'(0 31 59 90 120 151 181 212 243 273 304 334)))
(mon (apply + (map (lambda (f) (if (<= f d4) 1 0)) e)))
(d (- d4 (list-ref e (- mon 1)) -1)))
(list y mon d)))

テスト

(gr 1) => (1. 1 1.)
(gr 734274) => (2011. 5 16.)

こんな具合いだが, まぁいいか.

曜日の計算

森本雅樹君の追悼文集, 「森本さんの宇宙」の最後の年表「森本おじさんの足跡」に「1932年5月14日(何曜かわかりますか) 東京に生まれました」と書いてある. よし, 計算しよう.

その頃の日付と曜日の組で分かっているのは, 1941年12月8日が月曜ということだ. 太平洋戦争開戦の日を私は記憶している. その日, 学校に行ったという他は曜日は怪しい(月曜だったらしいというかすかな記憶はあるが). 一方, ハワイの軍港は, 12月7日の日曜で, 軍艦にいない水兵が多かったという話から, 12月8日は確かに月曜であったろう. その前, 数年の12月8日の曜日を書いてみると,

1941
   40 日 うるう年
   39
   38
   37
   36 火 うるう年
   35
   34
   33
   32

従って, 1932年12月8日は木曜だ. 12月12日と5月9日の曜日は同じだから, 5月14日曜日は木曜の4日後の5日後だから, 土曜か.

そこで, 手元のHHC電卓で確認する. 1932年5月14日を19320514と入力し, JDを押すと2426842. それに1を足して7で割る. 商は346691. 剰余は6. たしかに土曜であった.

ついでに,... 手元にMacBookもあるから, 道具がますます高級になり,

% cal 5 1932
May 1932
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

2011年5月19日木曜日

個人用電卓のプログラミング

またThe Art of Happy Hacking Calculator Programmingの話だ. 2年ほど前のこのブログに, 私のデザインしたHappy Hacking Calculatorのプログラミングの話を何回か書いた. 例えば1から10までの和のプログラムは,

(0Ua"N+U1-"_aG^N)X=

factorialは

("_6_G1+5_"G"1-!*)!=

という具合いである. Gの命令について復習すると, スタックの上から2段目の内容が負なら, プログラム制御をスタック最上段の数だけバックする(左へ進む), 非負ならそのまま進むというのであった. そこで, プログラムを書くのに面倒なのは, 飛び先までの長さの計測である. factorialのプログラムの左から5文字目のGはスタックの最上段に-6があるので, その下が負なら, 右へ(1+5_"G)の6文字を飛び越えることを示す. つまり右のGの次の"のところだ. 当時のプログラムでは, 長さのところを空けておき, 他が出来てから手作業で長さを数えていた.



またこの図の左のようだと, G命令から10文字左へ飛び, 右のようだと9文字左へ行く. それが同じ場所である. 従って, 工夫によってはプログラムが短く出来る.

これをなんとか自動化したい, つまりコンパイラを書きたいとかねがね考えていた.

飛び先までの間に条件付きのbranch命令や無条件のjump命令があると, そのパラメータの文字数が未定であったりするので, 自縄自縛(という表現が適切かな)になり, コンパイラの書き方が決らなかった.

最近, 思いついた方法があり, それでプログラムを生成してみたが, なんとかなりそうである. そのことを書きたい.

コンパイルされるプログラムはこんな形をしている.

(define sum '(0 up 10 'foo dup down + up 1 - dup neg
(branch foo) pop down 'exit))

(define factorial
'(dup neg (branch bar) 1 ent (branch exit) 'bar
dup 1 sub ! * 'exit))

(define ackermann '(
dup 1 - (branch a) exch dup 1 - (branch b)
exch dup 2 - (branch c)
up dup 1 - exch down 1 - A A (jump exit)
'c pop pop 2 ent (jump exit)
'b pop 2 * (jump exit)
'a pop pop 0
'exit
))

sumが10までの和のプログラムで, クォートされているリスト内がプログラム本体である. 0とか10は定数で, そのまま文字に置き換わる. 今の版ではコンパイラの関数によって十進だったり十六進数だったりする. up, dup, down, +, -, negなどが, プログラムの命令だ. (このAckermann関数は, SICPのex1.10にあるものだ.)

(branch foo)はfooへのbranch命令を挿入することを示す. その飛び先は'fooのところである. jump命令では, スタックの判定用に無理に負の数をおいて常に飛ぶのである.



それぞれの図の箱は, 長さ固定のプログラム. その間のbmやその間jmはbranchやjump命令. 矢印はそこから飛ぶ飛び先で, lnは飛ぶ長さである.
bmやjmは, 対応する飛ぶ距離lnの文字に変換した長さlに,

a. 左へのbranchはlの次にGの1文字,
b. 右へのbranchはlの次に_とGの2文字,
c. 左へのjumpは1_の2文字とlの次にGの1文字,
d. 右へのjumpはlの次に_"Gの3文字
をつける. つまり追加文字数は, a, b, c, dのそれぞれで, 1, 2, 3, 3である.

ackermannのプログラムは, 二進十進(あるいは十六進)変換を

(define (f n)
(if (< n 10) 1 (if (< n 100) 2 3)))

と定義し,

(define b0 0)
(define b1 0)
(define b2 0)
(define j0 0)
(define j1 0)
(define j2 0)

(define (iter)
(let* ((l0 (+ 4 b1 4 b2 10 j0 4 j1 3 j2))
(l1 (+ 4 b2 10 j0 4 j1))
(l2 (+ 10 j0))
(l3 (+ 4 j1 3 j2 3))
(l4 (+ 3 j2 3))
(l5 3)
(c0 (+ (f l0) 2))
(c1 (+ (f l1) 2))
(c2 (+ (f l2) 2))
(k0 (+ (f l3) 3))
(k1 (+ (f l4) 3))
(k2 (+ (f l5) 3)))
(if (equal? (list c0 c1 c2 k0 k1 k2) (list b0 b1 b2 j0 j1 j2))
(list b0 b1 b2 j0 j1 j2 l0 l1 l2 l3 l4 l5)
(begin (set! b0 c0) (set! b1 c1) (set! b2 c2)
(set! j0 k0) (set! j1 k1) (set! j2 k2) (iter)))))

のプログラムで(iter)とすると,
 (4 4 4 5 5 4 47 32 15 19 10 3)

が得られる. したがってb047_G, b132_G, b215_G, j019_"G, j110_"G, j23_"G となる.

facotorialは

(define (iter)
(let* ((l0 (+ 2 b1))
(l1 5)
(c0 (+ (f l0) 2))
(c1 (+ (f l1) 2)))
(if (equal? (list c0 c1) (list b0 b1))
(list b0 b1 l0 l1)
(begin (set! b0 c0) (set! b1 c1) (iter)))))

で,
(3 3 5 5)

が得られるから, b05_G, b15_G でよい.

10までの和は,

(2 9)

が得られ, b0はすなわち9Gとなる.

以上は, それぞれの関数用のプログラムを実行した例だが, プログラムの半ばコンパイルした中間出力のリストを使うことも出来るようになった.

ackermann関数の例では, bmやjmをリストにし, プログラムのブロックをその長さlnで置き換えた中間出力

((3 "\"1-") (b a 0 0) (4 ":\"1-") (b b 0 0) (4 ":\"2-") (b c 0 0)
(10 "U\"1-:N1-AA") (j exit 0 0) c (4 "^^2E") (j exit 0 0)
b (3 "^2*") (j exit 0 0) a (3 "^^0") exit))
から

"1-47_G:"1-32_G:"2-15_GU"1-:N1-AA19_"G^^2E10_"G^2*3_"G^^0

を得るプログラムは走るようになっているが, まだ改善の余地だらけで, 公開ははばかられる.

2011年5月18日水曜日

unix time

前回のblog以後, ユリウス日を年月日に変換するプログラムも書きたくなっている.

私の愛読書の1つ, E.M.Reingold, N. DershowitzのCalendrical Calculationsでその辺の計算法を見てみよう. ユリウス日は, Julian暦-4712年1月1日が起点だが, 本書のカレンダーの起点は, Gregorian暦を昔の方へ外挿した1年1月1日を1とし, それからのFixed Day Numberで表わす. この日数をR.D.(Rata Die, fixed dateのラテン語)という.

本書の関数は, 本文では不思議な構文で書いてあり, 付録にはCommon Lispによる実装があるが, Scheme風にすると,

(define gregorian-epoch 1)
(define (gregorian-leap? y)
(and (= (modulo y 4) 0)
(not (member (modulo y 400) '(100 200 300)))))

(define (fixed-from-gregorian y m d)
(+ gregorian-epoch -1 (* 365 (- y 1)) (quotient (- y 1) 4)
(- (quotient (- y 1) 100)) (quotient (- y 1) 400)
(quotient (- (* 367 m) 362) 12)
(cond ((<= m 2) 0)
((gregorian-leap? y) -1)
(else -2)) d))

(fixed-from-gregorian 1 1 1) => 1
(fixed-from-gregorian 2011 5 16) => 734273


この本では, 年内の日数と月の変換は, 2月も30日あるとして計算し, 補正する

(map (lambda (m) (quotient (- (* 367 m) 362) 12)) (a2b 1 13))
=> (0 31 61 92 122 153 183 214 245 275 306 336)

floorをとらないと, 以下のような値である.

(.417 31. 61.583 92.167 122.75 153.333 183.917 214.5
245.083 275.667 306.25 336.833)

隣り同士の差を取ってみると,

(map (lambda (a b) (- a b))
'(31 61 92 122 153 183 214 245 275 306 336)
'(0 31 61 92 122 153 183 214 245 275 306))
=>
(31 30 31 30 31 30 31 31 30 31 30)

補正は, 1月2月はこのまま, 3月から後はうるう年なら1を引き, 平年なら2を引く.

さて, Gregorian暦のy, m, dからfixed dateを得る最初の関数fixed-from-gregorianの解説である.

前年までの経過日数が必要なので, (- y 1)が頻出する.

(* 365 (- y 1)) ;前年までの平日の日数
(quotient (- y 1) 4) ;Julian暦のうるう日の日数
(- (quotient (- y 1) -100)) ;100年の倍数の年はうるうをやめる
(quotient (- y 1) 400) ;しかし400年の倍数なら, やはりうる
う年にする
(quotient (- (* 367 m) 362) 12) ;2月を30日と仮定して, m月
の前月までの日数
(cond ((<= m 2) 0) ;補正
((gregorian-leap? y) -1)
(else -2)) d)) ;dを足す.


最初のgregorian-epochは, 1年1月1日を1にするためである.
最後のテスト例のようにうまく行く.
予想通り, 逆は難しい. プログラムは以下のようだ.


(define (gregorian-year-from-fixed rd)
(let* ((d0 (- rd gregorian-epoch))
(n400 (quotient d0 146097))
(d1 (modulo d0 146097))
(n100 (quotient d1 36524))
(d2 (modulo d1 36524))
(n4 (quotient d2 1461))
(d3 (modulo d2 1461))
(n1 (quotient d3 365))
(y (+ (* 400 n400) (* 100 n100) (* 4 n4) n1)))
(if (or (= n100 4) (= n1 4)) y (+ y 1))))

(define (gregorian-from-fixed rd)
(let* ((y (gregorian-year-from-fixed rd))
(prior-days (- rd (fixed-from-gregorian y 1 1))
  (correction (cond ((< rd (fixed-from-gregorian y 3 1)) 0)
((gregorian-leap? y) 1
(else 2)))
(m (quotient (+ (* 12 (+ prior-days correction)) 373)
367))
(d (+ (- rd (fixed-from-gregorian y m 1)) 1)))
(list y m d)))

(gregorian-from-fixed 1) => (1 1 1)
(gregorian-from-fixed 734273) => (2011 5 16)

と定義しておき, (gregorian-year-from-fixed rd)
で rd に対する年 y を計算する.

(n400 (quotient d0 146097)) :d0に400年の日数が何回あるか見る. 146097は400年の日数.
(d1 (modulo d0 146097)) :その400年内の日数

(n100 (quotient d1 36524)) ;その日数に100年は何回あるか. しかし400年の100年の日数は, 最初の100年の先頭の00年はうるう年なので, 36525日あるから

(d2 (modulo d1 36524)) ;その100年内の日数
(n4 (quotient d2 1461)) ;その日数内の4年の数
(d3 (modulo d2 1461)) ;その4年内の日数
(n1 (quotient d3 365)) ;4年内の1年の数
(n1 (quotient d3 365)) ;4年内の1年の数

基数変換をやっているみたいに簡単なのに驚く. 400年の中の日数d1のいろいろな値から得られるyをみてみる.


(define (gregorian-year-from-fixed rd)
(let* ((n100 (quotient rd 36524))
(d2 (modulo d 36524))
(n4 (quotient d2 1461))
(d3 (modulo d2 1461))
(n1 (quotient d3 365))
(y (+ (* 100 n100) (* 4 n4) n1)))
(if (or (= n100 4) (= n1 4)) y (+ y 1))))

すると, 以下のようになっていることが分かる.

rd y rd y
0- 364 1 365 35794- 36158 99 365
365- 729 2 365 36159- 36523 100 365
730- 1094 3 365 36524- 36888 101 365
1095- 1460 4 366
1461- 1825 5 365 72318- 72682 199 365
1826- 2190 6 365 72683- 73047 200 365
2191- 2555 7 365 73088- 73412 201 365
2556- 2921 8 366
2922- 3286 9 365 108842-109206 299 365
3287- 3651 10 365 109207-109571 300 365
3652- 4016 11 365 109572-109936 301 365
4017- 4382 12 366
4383 4747 13 365 145001-145365 398 365
4748- 5112 14 365 145366-145730 399 365
5113- 5477 15 365 145731-146096 400 366
5478- 5843 16 366


0から勘定を始めるのが基本と思っている私なら, 下の図の左のような関数を書くところだが, 本書の流儀は違う. d1=0ならy=1が返る. d1=146096なら, y=400であった. なるほど! (図の太い横線は366日の年を示す. 横線の左の黒丸は閉区間, 右の白丸は開区間を示す.)




別の図を描くと下のようになる. つまり, 図で網掛けの例外を最後に置くので計算が簡単になっていたのだ.




それなら, 私の流儀でも, 400年紀末から逆に計算すればおなじわけだった.


ところで, これではまだfixedから, 年が得られただけである. 本番のプログラムはこれを下請けに使う(gregorian-from-fixed d)である. その解法はこうだ.

とりあえずこのrdの落ちるyを求める. 次のその年の1月1日のfixed dateを求め, prior-daysとする. その年の3月1日のfixed dateを求め, prior-daysがそれより小さければ, 補正は0, そうでなくて, うるう年なら補正は1, 平年なら補正は2である. 次に前に日数を計算した式で, 月を見つけ, その月の1日のfixed dateとの差から, 日が分かるのである.

fixed-dateを何回も使うが, それだけ分かりやすいアルゴリズムになっている.

こういうプログラムの解読もなかなか面白い.

2011年5月16日月曜日

unix time

unixには1970年1月1日正子(0時0分)からの延秒数を数えている32ビットの時計がある. 2038年1月19日にサインビットが立つといわれ, 2000年問題みたいになにか起きかもしれないが, 私は多分もうこの世にはいず, 状況を知ることはかなわぬ.

この時計の元は, MITのMulticsではないか. Multicsには, 1900年1月1日正子からのマイクロ秒を数える52ビット時計があった. マイクロは10-6だから, 20ビット程度であり, unixの32ビットに対して52ビットなのは分かる. 私がMITに滞在したのは, 1973年9 月から74年7月までだが, その時計のサインビットが立ったのは, その少し前のたしか5月だったとある院生から聞いた.

まず脱線して, それがいつだったか計算してみよう. 例の個人用電卓が活躍する. 251は2251799813685248. これを1日のマイクロ秒864000000000で割る.

商は26062, 剰余は43013685248. つまり1900年1月1日から26062日後を知りたい. それには1900年1月1日のユリウス日2415021に26062を足し, その2441083がユリウス日になる日を知ればよい.

ここから先は電卓から離れ, 理科年表のユリウス日の表による.

すると, 1971年5月11日がその日であることが判明. 日以下を計算すると, 11時 56分 53秒 685248マイクロ秒であった.

jdをユリウス日を計算する関数として, Schemeで検算すると,

(+ (* (- (jd 1971 5 11) (jd 1900 1 1)) 86400000000)
(* 11 3600000000) (* 56 60000000) (* 53 1000000) 685248)
=> 2251799813685248

(factorize 2251799813685248)
=>
(2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2)

(length (factorize 2251799813685248)) => 51


本題へ戻り, 通常のy年mon月d日h時min分s秒からunix timeへの変換の方は, Multicsの場合と同様に簡単だ.

(JD(y,mon,d)-JD(1970,1,1))*86400+h*3600+min*60+s

となる.

一方, unix time tからy,mon,d,h,min,sへの変換は, ユリウス日からy,mon,dへの変換関数があればなんでもないが, まだそういうプログラムは書いたことがない.

カレンダーの変換で一番面倒なのは, Grogorian暦でのうるう年の計算である. しかし, ことunix timeに限れば, 2000年が通常のうるう年なのが幸いだ. そこで書いたのが次のSchemeのプログラムである.

(define (unixtime t)
(let* ((s (modulo t 60)) (m (quotient (modulo t 3600) 60))
(h (quotient (modulo t 86400) 3600)) (d (quotient t 86400))
(y (- (floor (/ (+ d 731) 365.25)) 2))
(c (- d (floor (+ (* y 365.25) 0.25))))
(e (if (= (modulo y 4) 2)
'(0 31 60 91 121 152 182 213 244 274 305 335)
'(0 31 59 90 120 151 181 212 243 273 304 334)))
(n (apply + (map (lambda (f) (if (<= f c) 1 0)) e))))
(list (inexact->exact (+ 1970 y)) n
(inexact->exact (- c (list-ref e (- n 1)) -1)) h m s)))


引数のtがunixtimeである. 最初にs(秒), m(分), h(時)を取り出す. dは通算の日数になる. その後の変数は, yが1970年以降の年数, cがその年内の日数, nが月だ. yはdを365で割ってfloorを取ればよいが, 閏年があるからそうは問屋が卸ろさない.

まず, y年について, 前の年の終りまでの日数は,



欲しい値は,



365の代りに365.25で割ればよさそうに見えるので, 2.25でテストしてみる.

(map (lambda (n) (floor (/ n 2.25))) (a2b 0 20))
=> (0 0. 0. 1. 1. 2. 2. 3. 3. 4. 4. 4. 5. 5. 6. 6. 7. 7.
8. 8.)

4が3個並ぶのがうるう年に対応し, これを2年にしたいから, 3つの0と2つの1の5個をスキップするために, nの代りに(+ n 5)とし, 最後に2を引く.

(map (lambda (n) (- (floor (/ (+ n 5) 2.25)) 2)) (a2b 0 20))
=> (0. 0. 1. 1. 2. 2. 2. 3. 3. 4. 4. 5. 5. 6. 6. 6. 7. 7.
8. 8.)

なるほどうまくいくので, 356.25に修正し, テストする.

(map (lambda (n) (- (floor (/ (+ n 366 365) 365.25)) 2))
'(0 364 365 729 730 1095 1096 1460 1461 1825 1826 2190
2191 2555))
=> (0. 0. 1. 1. 2. 2. 3. 3. 4. 4. 5. 5. 6. 6.)

うまくいく. これを1970に足せばよい.

次に上の表にあった前年までの日数の和を計算するには,

(lambda (y) (+ (* y 365) (quotient (+ y 1) 4)))

(lambda (y) (floor (+ (* y 365.25) 0.25)))

とする. これをdから引くと, その年内の日数cが得られる. 前の年の終りまでの日数の和のように, 前の月の終りまでの日数の和のリストを, うるう年か否かで変数eに用意する.

月nは, このリストで, cが越えるものの数として得る. nが決れば, 月内の日数は, cから先ほどのリストの要素を引いて作る.

これで完成. テストしてみる.

(unixtime 0) => (1970 1 1 0 0 0)
(unixtime (expt 2 31)) => (2038 1 19 3 14 8)

2011年5月10日火曜日

加算の順序

桁数が多い十進数でも, その筆算による足し算は小学校で習う. 最近のことは知らぬが, 昔は算盤も小学校で習った.

ところで, 筆算は下の桁から始め, 算盤は上の桁から足し始める. もちろん, 結果は同じになる. これは一体どういうことかと, 子供のころ不思議であった.

実は, 足し算(引き算もだが)は, 繰上げを正しく処理し, 全部の桁について加算するなら, どの順に桁を選んでやっても構わないのである.

算盤で 8 8 8 8 8 8 8 8 に 1 1 1 2 1 1 1 2 を足すとき, まず最初の3桁で,

8 8 8 に
1 1 1 を足すと
9 9 9 になる.

次に8に2を足すから, 10になるが, 算盤流では8を払ってその桁を0にしてから, 繰上げ処理に入る. つまり, 次の9に1を足す. それには9を払い, さらに繰上げを処理する. 3つの9から繰上げが続き, 9の上の桁に1がでて止る. 都合4回繰り上げた. 繰上げをピリオドで示し, これを 1.0.0.0.0 と表わす.

その次に3回9が得られ, 1.0.0.0.0 9 9 9 となり, 最後の8+2でまた繰上げが発生. 4回繰り上げて1.0.0.0.1.0.0.0.0 が得られる.

以上が算盤であった.

筆算だと右端の8+2から始め, この桁の結果が0, 繰上げがでる. 従って .0 となり, 繰上げがあったことを記憶して, 次の8+1に進む. この和は9だが, 繰上げがあったので, 0になって再び繰上げを記憶, 次へ進む. このようにして, .0.0.0.0 まで来た. 次は8+2なので, 和は10. それに繰上げで11になり, この桁が1となり, さらに繰上げが続く. こうして,1.0.0.0.1.0.0.0.0となる.

左右のどちらから計算しても, 繰上げの出る場所は決っている.

タイガーのような計算機ではどうか.

8 8 8 8 8 8 8 8 に
1 1 1 2 1 1 1 2 を足すと一旦
9 9 9.0 9 9 9.0 となり,

繰上げが発生した場所がマークされる. その後, 下の桁から繰上げセンサーがマークをスキャンし始める. マークがあると, 繰上げピンにより, 次の桁を1増やす. この結果繰上げが生じることがある. 今の場合もそうだ.

9 9 9.0 9 9.0.0 となり, スキャンは今つけたマークのところに来る. そして9を0にし, 次に桁との間にマークする. 繰上げで1増やそうとするとき, その桁が9でなければ繰上げ処理は終る. 0のところが 1.0.0.0.0 となるが, その次にすでにマークがあるので, 続く9を次々を0にして繰上げを伝える.

結局, 上の算盤や筆算と同じ場所で繰上げが生じ, 同じ結果になる.

足し算をランダムな順でやってみよう.

0 3 6 2 7 5 1 4 ← 足す順
1 1 1 2 1 1 1 2 ← 足す数
8 8 8 8 8 8 8 8 ← 元の数

まず順が0の桁で足す. 従って

0 3 6 2 7 5 1 4 ← 足す順
1 1 1 2 1 1 1 2 ← 足す数
9 8 8 8 8 8 8 8 ← 元の数

となる. 繰上げはない. 次に順が1の桁へ行く.

3 6 2 7 5 1 4 ← 足す順
1 1 2 1 1 1 2 ← 足す数
9 8 8 8 8 8 9 8 ← 元の数

順が2の桁へ行く

3 6 2 7 5 4 ← 足す順
1 1 2 1 1 2 ← 足す数
9 8 8.0 8 8 9 8 ← 元の数

と繰上げが生じたからその処理をして

3 6 2 7 5 4 ← 足す順
1 1 2 1 1 2 ← 足す数
9 8 9.0 8 8 9 8 ← 元の数

順が3の桁

3 6 7 5 4 ← 足す順
1 1 1 1 2 ← 足す数
9 9 9.0 8 8 9 8 ← 元の数

順が4の桁, 繰上げも処理する.

6 7 5 4 ← 足す順
1 1 1 2 ← 足す数
9 9 9.0 8 9.0.0 ← 元の数

順が5の桁

6 7 5 ← 足す順
1 1 1 ← 足す数
9 9 9.0 9.0.0.0 ← 元の数

順が6の桁

6 7 ← 足す順
1 1 ← 足す数
1.0.0.0.0 9.0.0.0 ← 元の数

順が7の桁
7 ← 足す順
1 ← 足す数
1.0.0.0.1.0.0.0.0 ← 元の数

これで, お勧めはしないが, どういう順に足してもよさそうだと分かる.

十進法の2つの数を足すとき, 各桁で加算し, 仮の和 0≤t<19を得る.

それと下の桁からの繰上げの有無により, その桁の真の和と, 上の桁への繰上げが得られる.


t<9なら, 繰上げのありなしに拘らず, 上の桁への繰上げは生じない. 真の和は, 繰上げがないならt, あるならt+1.

t=9で, 繰上げがないなら, 繰上げは生じず, 真の和は9. 繰上げがあるなら, 繰上げが生じ, 真の和は0.

10≤tなら, 繰上げのありなしに拘らず, 上の桁への繰上げが生じる. 真の和は, 繰上げがないならt-9, あるならt-10.

上の加算の例でいうと,

8 8 8 8 8 8 8 8
1 1 1 2 1 1 1 2
9 9 9 0 9 9 9 0 ← 仮の和
. . ← 繰上げ
0 0 ← 真の和
. . ← 繰上げ
0 0 ← 真の和
. . ← 繰上げ
0 0 ← 真の和
. . ← 繰上げ
1 1 ← 真の和