開平法

手回し計算機は十進である. 二進でも事情は変らない. むしろ簡単といえる.

1000000000 (2^9=512) を筆算で開平する手順は以下のとおり.



左部分は2x+yを作っているから, 最後の平方根はkの行の左部分の半分,　10110(=22)で, 剰余はkの行の中部分の 11100(=28) である. つまり22*22+28=512である.

これをプログラムに書き直したのは, 以下のとおり. 途中経過を印字する行はいらない.

(define (sqrt n)
  (define (ss n s) (if (> (* s 4) n) s (ss n (* s 4))))
  (define (sq n q s)
    (display (list 'n n 'q q 's s)) (newline)
    (if (> s 0)
      (let ((a (+ q s)))
        (if (>= n a) (sq (- n a) (+ (/ q 2) s) (quotient s 4))
             (sq n (/ q 2) (quotient s 4))))
      (cons q n)))
  (sq n 0 (ss n 1)))

下請け関数 sq の引数 n は部分剰余, q は 部分平方根, s は計算中の位置を記憶するストロブである. 各ステップで, 部分剰余と比較する数 a を作り, 引けるなら部分剰余は n-a になり, 次の部分平方根とストロブを作る. 引けなければ部分剰余はそのまま, 平方根は1桁右にずれるだけである.

もう1つの下請け関数 ss は, 1から4倍を繰り返し, n と比較しながら最初のストロブの位置を決める.

512を開平した時, n, q, sの状態は次のとおり.

 (sqrt 512)
(n 512 q 0 s 256)
(n 256 q 256 s 64)
(n 256 q 128 s 16)
(n 112 q 80 s 4)
(n 28 q 44 s 1)
(n 28 q 22 s 0)
=> (22 . 28)

以前のブログに書いた個人用電卓 Happy Hacking Calculatorの開平は, 基本的にこの方法を採用している.


ところで, これまでの二進の開平では, 引けるかどうかを見て, 引ければ引くという方法であった. これに対し二進の開平では, 引放し開平が出来る. つまり引きすぎて部分剰余が負になったら, 足し戻さず, 次のステップでは足すのである. 理屈は面倒なので書かないが, 512をこの方法で開平する様子を示す.



この方法ではまずストロブの位置から1を引く. 部分剰余が正なら, 部分平方根の最後に101をつけたものを引く. 負なら, 011をつけたものを足す. 平方根は最初の無条件の減算を除き, 引いたときは1, 足した時は0として作る.最後の桁を1として追加する. 従ってこの方法では, 平方根は必ず奇数になり, 剰余も正だったり負だったりする.

512を開平すると, 引く, 足す, 引く, 引くとやったので, 1011となり, 最後に1を追加して, 10111 (=23)が得られる. 剰余は -17 である. 23*23-17=512になっている.

平方根が q, 剰余が r で負の時は, 平方根が奇数なので最後の1を0に変る. つまり平方根を q-1にする. 平方根が q の時, 剰余が r なので, 新しい剰余を r' とするとq*q+r=n=(q-1)*(q-1)+r' だからr'=r+2(q-1)+1 とすればよい. 上の例では, 23から1引いた22が平方根, その剰余は-17+2*22+1=28 と得られる.

開平法

手回し計算機で開平するには, 手回し計算機のことも説明する必要があろう.

1950年頃からリヒテンシュタインで製造されていた, Curtaという小さい機械式の計算機があった. 設計したのはCurt Herzstarkで, 第二次世界大戦中, ドイツ軍に抑留されていたときに設計したといわれている.

サーチエンジンでCurta Calculatorを探すとサイトがいくつも見つかるが, 精巧なシミュレータはここ, マニュアルはここにある.

シミュレータの説明用に私が描いたCurtaの絵は下の通り.






上が正面図, 中が平面図である. 下は説明用の図(表?)だ.

Curtaは筒状で, 正面下側に8桁のSetting Registerがあり, 赤で示すノブを下げると0から9までの各桁の値が設定できる. シミュレータでは同時に下の表のSETTING DIALにも値が入る. Operation HandleをZero Positionから時計方向に1回転(平面図でHandle から下に出ている矢印をクリック)すると, この値がResult Dialに足される. HandleはZero Positionまで回す.

Curtaの上方には, 左右に5桁分回転できるCarriageがあり, 平面図ではこれが見えている. Carriageの上面には, 11桁のResult Dialと6桁のRevolution Counterがある. Carriageがこの図の位置なら, Setting Registerの値は, Result Dialの最下位に足される. 正面図のCarriageの右の赤い矢印をクリックすると, 平面図のCarriageは反時計方向に1桁分回転し, その位置では下から2桁目に足される.

Operation Handleを1回転すると, Setting Dialの最下位に対応するRevolution Counterに1足される. 従ってSetting Dialに34を置き, Carriageの最初の位置でHandleを6回転し, Carriageを1桁回し, 3回転すると, Result Dialには34*36=1224が入り, Revolution Counterには36が入る.

平面図で11時の方向に出ている輪のようなものは, Clearing Leverで, Revolution CounterとResult Dialの境界に止まり, どちらかに回転すると, そのCounterやDialがリセットされる.

Curtaの正面図で, Setting Registerの右に見えるReversion Leverは, ノブが図の位置では, Handleの回転でRevolution Counterに1足され, ノブを下に下げると, Counterから1引かれる.

Curtaの特色の1つは, 補数計算をするので, 減算でもOperation Handleを加算と同様に時計方向に回すことだ. ただしOperation Handleに右にある上向き矢印をクリックし, 回転軸を上に少し引き抜く. 加算に戻すには, 減算時に下向きになっていた矢印をクリックして, 回転軸を押し下げる.

除算の時は, 被除数をResult Dialに, 除数をSetting Dialにセットし, Reversion Leverを下位置にし, 回転軸を引き抜き, 各桁で引けなくなるまで回転し, Revolution Counterに商を, Result Dialに剰余を得る.

さて, 手回し計算機での開平は, 20x+1, 20x+3, ...,20x+17と引く代りに, 被開平数を最初に5倍し, 100x+05, 100x+15, ..., 100x+85と引く. Setting Dialの最下位を5に固定し, 十の桁を0,1,2,...と変えながら引く. 引ける回数は同じなので, 平方根も5倍せずに得たのと同じである. もちろん剰余は5倍になっている.

その様子を下に示す. まず200000000の5倍をResult Dialに置き, Carriageを右方へ寄せておく. 短い横線は減算を, 長い左向き矢印は, Carriageを1桁左へ移動することを示す.




赤字は十の桁が, 0,1,2,...と変る様子を, 青字はRevolution Counterが増えていく
様子を示す. 最後の剰余は, 前回の開平で得られた3836の5倍である.

開平法

1月13日から15日, プログラミング・シンポジウムに参加してきた.

手回し計算機の話題で, 開平が出来るといったら, 怪訝な顔をする人がいたので, 大きなお世話かとも思うが, 私流に説明したい.

まず筆算の方法のおさらいから. 十進法で200000000の平方根を計算してみる.

開平は基本的には除算と似ている. 除算では被除数を最初の部分剰余とし, 部分剰余から除数の倍数を引いて次の部分剰余を作るというステップを繰り返す. 最後の部分剰余が剰余になる. 各ステップで商が1桁ずつ得られる.

除数の何倍を引けば次の部分剰余が除数より小さくなるかを決めるのが, 除数の難しい点である.

開平も被開平数(radicand)を部分剰余とし, 次々と部分剰余と平方根の1桁を求めながら進行するのは除算と同じである. しかし除数の倍数として商が得られるのでなく, 平方根が除数のような相手なしに徐々に決っていくところが除算と違う.

下の図が筆算の経過で, 左部分と中部分と「=」より右の右部分からなり, 平方根が10000, 14000, 14100, 14140, 14142と順に得られる様子を示す. 右部分は中部分の値の説明で, 筆算ではいらない. この図で, aの行の中部分の200000000が被開平数で, 最初の部分剰余である. 十進の開平には平方根の1桁は, 小数点を中心にした2桁ずつの部分剰余から決める. 最初に検討するのは, 小数点から2桁ずつ左に区切っていった最後の2桁(2だけの)部分である. 2を超えない最大の平方根は1なので, それが左のaとbにある1と1になる. その1と1を掛けた1を中部分の2から引く. cの中部分の10000000は右にあるように, 被開平数から最初の平方根を引いたものである.



最初の平方根1をxとする. 次の桁の平方根をyとすると, そこまでの平方根は10x+y. 部分剰余は被開平数から(10x+y)²=100x²+20xy+y²を引いたものになる. しかしすでに100x²は引いてあるので, 20xy+y²=(20x+y)*yを引くことになる.

cとdの行の左部分の24と4がこの20x+yとyに相当する. cの行の2はaとの行の1と1(つまりxとx)を足したものだ. (20x+y)*yで100を超えない最大のyを探すと4なので, cとdの左に24と4を書き, 中にその積96を書き,100から引いて次の部分剰余400が得られる.

これを繰り返すとkの行で最後の部分剰余3836が得られ, 平方根が14142となる.

以前このブログに書いた個人用電卓 Happy Hacking Calculatorの平方根は剰余も得られるが, この部分剰余がスタックの2段目に得られる.

ところで, 除算でもこの剰余と除数に対し商に何を立てるか決めるのは面倒である. 一方, 手回し計算機では, 除数を引けなくなるまで繰り返し引くことである桁の商がきまる.

手回し計算機での開平法は, 下の図がその説明である. 上の図で24掛ける4を決めた代りにyを1から順に増やすのである. それまでに得た平方根をxとし, まずy=1として(20x+1)*1を引く. 次にy=2として(20x+2)*2を引くのだが,(20x+1)*1はすでに引いてあるので, (20x+2)*2-(20x+1)*1=40x+4-(20x+1)=20x+3を引く. 更にy=3とすると20x+5を引くことになる.



かくして手回し計算機でも, 無事に開平が出来るのである.

sideways addition

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

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



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

x₁₀=x₁⊕x₂⊕x₃,  x₁₁=x₄⊕x₅⊕x₆,  x₁₂=x₇⊕x₈⊕x₉,  x₁₃=x₁₀⊕x₁₁⊕x₁₂,
y₁=⟨x₁x₂x₃⟩,  y₂=⟨x₄x₅x₆⟩,  y₃=⟨x₇x₈x₉⟩,  y₄=⟨x₁₀x₁₁x₁₂⟩

とすると

x₁+…+x₉=x₁₃+2(y₁+y₂+y₃+y₄)

で得られる. 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ビットを右からx₁,...,x₉とする. 3ビットずつに区切り, その区間の和を求めると0から3までになるから, 2ビットで表現出来, それが下のx₁₀,y₁などである. x₁,x₂,x₃の3ビットの二進法の和は, その1の桁(x₁₀)はx₁,x₂,x₃の排他的論理和だし, 10の桁(y₁)はそれらの多数決だから, アルゴリズムにはそう書いてある. 次に3つの区間の3ビットの和が得られたとして, この和を取るわけだが, 和の1の桁はx₁₀,x₁₁,x₁₂の排他的論理和である. 10の桁y₁,y₂,y₃を足せばよく, それにx₁₀,x₁₁,x₁₂からの繰上がり, すなわちそれらの多数決を足す. 10の桁だから2倍する.

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

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

EDSACの語長は35ビットだが, 36ビットと考えて6ビットの区間６つとする. ビットシフトやビット毎ANDや加算して, ６つの区間で並列にその区間の1の数を和を作る. 6までだから3ビットの範囲に収まる. つぎに10⁵10⁵...0⁵1（つまり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重平衡変調器が思ったほどにはうまく働かず, 結局パリティビットなしになった.

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

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が求めるものである.

月齢カレンダー

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

ある日の正午の月齢は, 直前の新月からの時間である. 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.))  

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

月齢カレンダー

これまでカレンダーのプログラムは数えきれないほど書いている. 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日に日食がある.