曜日の計算

前回のブログで島内先生の計算尺方式の万年七曜表を話題とした.

余談ながら現在東京理科大学の近代科学資料館で開催中の企画展「計算する器械たち---アナログコンピュータ展」には 円筒形の計算尺も展示されている. 今回はそういう円筒方式の万年七曜表の話だ.

下の図を見てほしい. ちょっと分り難いがこれが円筒式のもの(展開図ともいうべきもの)だ.



前回のブログを見た人なら思い出すだろうが, 最上段の右から0,1,...,15 少しあけて15,16,...,20とあるのは西暦の上2桁である. 右側の15まではユリウス暦, 左側の15からはグレゴリオ暦用のものだ.

その下, 5段位に左下がりに並ぶのは西暦の下2桁である. 段々と下るのは円筒の卷いたとき, 21に22を自然に繋げるためである. なかなか凝っている(次の図を見ると分る). 44は丁度切れ目になってしまった.

次は2, 8, 5,...と現れて, 4と7が重さなっているから, これは月名である. 展開図の範囲に4回出てくる.

最後は, 今度は右下がりだが, これは一箇月のカレンダーである.

もうひとつ青色の枠の図がある. 週日名が書いてあるほか. 枠の所が3カ所くり抜いてある.

島内さんの記事では, 西暦上2桁と月名が1枚の紙でCといい, 下2桁とカレンダーがもう1枚の紙でBといい, 青色の枠の紙はAである. 同じ幅に描いてあるが, A,B,Cの順の幅がすこしずつ短くなり, Aが一番外側, Bが真中. Cが内側でそれぞれ円筒状になって外中内の順で差し込まれている.

勿論底の場所は同じで, 展開図では左右に移動するように, 円筒は回転できる.

前回の例と同様に2014年7月の曜日を見るには,



のようにCの20の真下にBの14を合せる. そしてCの7月がAの月名を入れる青色の小さい枠に入るように合せると, 図のように青色の大きい枠がカレンダー部分に重さなり, 7月1日が火曜と判明するのである.

この図では, Bの紙が円筒になっているのが見えるであろう.

難しいのは例の丸めのところである. この図でも青色枠内の4と7が枠の左右の中央にないのが見てとれる. 4と7の中央は1.125という位置にあり, 枠の中央が1なのである. どこか中央かは, 下のカレンダー部分の枠の左右がカレンダーの数字の丁度真中に來るように離散的に合せることで分るのだ.

もうひとつの例



この例は2000年1月で, 閏年であり, 1月, 4月, 7月が同じ曜日の配置になる(Aの同じ枠に入る)ことが見てとれる.

さてこの計算尺の仕掛けの説明だが, mod 7の範囲だけの図にしたのが次だ.



最上段の線が西暦の上, その下の線が西暦の下で, 2014年になるようにその20と14を対応させた図になっている. つまりBの紙が右に移動している. (輪だから左かもしれない.) 移動していることを示すため赤で描いてある.

その下の線が月名で, これは最上段と同じ紙, 同じ位置にある.

さて赤の線の基準に対して7月の位置を計算すると, g(14)とf(20)が合っているから, 黒の基準位置はg(14)-f(20). 従って7月の位置はg(14)-f(20)+b(7). ところがfのスケールは反対に取ってあるから, g(14)+f(20)+b(7)になって, 実際の値を入れると
17+5.875+6.25=29.125=1.125 mod 7
従ってこれを丸め, 7月0日の値が1. すなわち7月1日は2(火曜)となる.

7月のすぐ右に9月がある. だから9月0日の値は0, 9月1日は1(月曜)だ.

従って青色の枠は下のように作ればよい.



上の小さい枠に7月があると, その1日は火曜. 上の枠がひとつ右に移ると, 曜日の名前が右にひとつずれて, 9月1日が月曜になるという仕掛である.

曜日の計算

私が2011年5月25日や6月5日のブログに書いた故島内剛一(しまうちたかかず)先生の曜日の計算は興味深いものであった. 先生の数学セミナーの記事(1978年7月号)には驚くようなアルゴリズムが沢山登場する.

その中で, しばらく私が読み飛ばしてはいたものの, 気になっていたのは, 「計算尺方式の万年七曜表」であった.

先日来, ちょっと読み直してみたが, もとの図を島内さんがどのようにして書いたか不明である. 雑誌に掲載された図の正確さもよく分らないので, 使い方が分かった時に, 自分で図を書き直してみた. (元の記事には「日のひと目盛の1/8=0.125までの精度はどうしても必要」と書いてある.)

ざっといえば, 下のようなものである.



上下2枚の図があり, どちらも半分が上下逆になっているのは, 中央で折って表と裏にするのである. 島内さんは上をB, 下をCといった.

Bの上半分は横に2倍にしたカレンダーである. 左側のどの列から始めても一週間分があるようになっている. 下半分は西暦の下2桁の表だ. こちらもかなり冗長である.

Cの下半分は月名である. 左の方, 3と11が縦に書いてあるのは, 3月と11月の曜日が同じことをしめす.

上半分は西暦の上2桁になっている. 0から15までの部分と15から22までの部分があるが, 前者はユリウス暦, 後者はグレゴリオ暦に対応する.

私は島内さんの記事のこれらの図をコピーし, 切り抜き, BとCを計算尺のように滑らせながら操作してみて, 使い方が分ってきた.

逆さまの図は見にくいから, 正立させた図にしよう. 上下の余白も省いた.



この上の図では上と下が表面と裏面に相当する. それぞれの面の途中の横線が計算尺の滑る位置だ. 上のカレンダーの中央付近に緑色の枠が見えるのは, 通常の計算尺にあるカーソルのようなもので, その内部がその年その月のカレンダーになる.

前の説明の繰り返しだが, 上から月, 日, 年の下2桁, 年の上2桁である.

さて, 2014年7月のカレンダーを見るには下の図のようにする.



下の西暦上2桁の20と, 上の西暦下2桁の14が合うように下の目盛を右にずらす. すると上の月名の目盛は(上下反転したので)同じ長さだけ左にずれる.

日の表の上端に曜日の区間の境界を示す青い線があり, 7の下の対応する区間の中央にM(monthのつもり)がくるように緑の枠を滑らすと, 7月1日が火曜になり, 7月のカレンダーが得られる仕掛けである.

目盛合わせが正確にできるよう, この図には多くの縦線が引いてある. ユリウス暦とグレゴリオ暦の間には微妙な差があるので, グレゴリオ暦の部分は青線で表示する.

この計算尺でカレンダーが出来る魔法のような理由は次のとおり.

式が導出された経過は省略するが, 島内方式で西暦y年m月d日の曜日は

round(f(y idiv 100) + g (y % 100) + b (m) + d) % 7

で計算する.

y idiv 100は西暦の上2桁. y % 100は下2桁である.

f(y)= - floor(3(y+1)/4) - y + 6 + (y % 4 == 0 ? -0.125 : +0.125)
(ユリウス暦に対するfの式は省略)
g(y)= y==0? -0.25 : (y + floor (a/4) - y % 4 ==0 ? -0.5 : 0)
bは1月から12月について
6.75, 2.75, 3.25, 6.25, 1.25, 4.25, 6.25, 2.25, 5.25, 0.25, 3.25, 5.25
である.

y=15から22までのfの値は
0.125, 5.875, 4.125, 2.125, 0.125, 5.875, 4.125, 2.125

y=0から15までのgの値は
-0.25, 1, 2, 3, 4.5, 6, 7, 8, 9.5, 11, 12, 13, 14.5, 16, 17, 18

アルゴリズムを解明すべく, これを次のような目盛に記入する. gは7の法をとってある. 目盛を下からf, g, m, dということにする. dの線の上の1〜7はカレンダーの一番上の数である.



いま2014年の計算をしようとして, 下の図のようにfの目盛の20とgの目盛の14を合わせると(赤い線), gの0はfの1.875に合う.

これはf(20)=5.875, g(14)=3の和 5.875+3=8.875=1.875(mod 7) である.



通常の計算尺でA尺とB尺を使い, log x + log y を計算するには, A尺のlog xのところにB尺の1(原点)を合わせ, B尺のlog yに対応するA尺を見るとlog x + log yになっている.

しかし今の計算では, gの目盛が逆向きにとってあるので, 通常の計算尺の減算の要領で加算が出来るようになっている. 結局fが右に1.875ずれているわけで, 同時に一番上のmも左1.875(=f+g)ずれている. mでの7月のb=6.25に対応するdの位置(緑の線)は, f+g+bの値f(20)=5.875, g(14)=3, b(7)=6.25の7を法とした和1.125)になっており, 丸めると1になる.

7月はこの値が1. ということは1日はd=1を足すから2になって火曜になる. 従ってカレンダーの中央の水曜を2日の上に置くことになる. dの値の1の上に2とある理由だ. 丸めた値が2だと1日は3だから水曜で, カレンダーの中央を1日に置く. dの目盛のカレンダーの値はそのように決めてある.

この説明も回りくどいが, 実は私も理解に到達するのに結構苦労した. 島内さんの発想に驚くばかりである.

微分解析機

東京理科大学近代科学資料館で企画展「計算する器械たち」を開催中だ.

そこで毎日微分解析機を実演している. デモは例のサークルテストだが, やはり一周すると2パーセントくらいのエラーが出る.

このデモの様子を, 昔微分解析機を使っていらした, 渡辺勝先生にお見せしたら, 「サークルテストでの反転時のバックラッシュが小さいようにお見受けしております」というメイルを頂いた.

ところで手回し計算機や電動計算機が出回ると, Runge Kuttaという計算法で連立m元常微分方程式を解くようになった. 我々がプログラミングを勉強したケンブリッジ大学のEDSACの本の例題にRunge Kutta Gillがあったので, 当時のコンピュータ屋にはお馴染の方法である. この方法については「伊理正夫,松谷泰行,Runge-Kutta-Gill法について,情報処理,Vol.8,No.2,pp.103-107」に解説がある. 今回はそれを使ってサークルテストを計算してみた. Schemeのプログラムは次のようだ.

(define (runge-kutta-gill m h n)
 ;m 変数の個数, h 分点間隔, n 分点数
(let ((y (make-vector m)) (f (make-vector m))
  (q (make-vector m)) (k 0) (qi 0) (s 0) (r 0)
  (c0 (- 1 (/ 1 (sqrt 2)))) (c1 (+ 1 (/ 1 (sqrt 2)))))
 (define (inity)    ;yの初期値
   (vector-set! y 0 0)
   (vector-set! y 1 0)
   (vector-set! y 2 1))
 (define (calcf)    ;fを計算
   (vector-set! f 0 1)
   (vector-set! f 1 (vector-ref y 2))
   (vector-set! f 2 (- (vector-ref y 1))))
  (define (printy)  ;yの出力
   (display (list 
     (vector-ref y 0)
     (vector-ref y 1)
     (vector-ref y 2))) (newline))
  (newline) (inity)
  (do ((i 0 (+ i 1))) ((= i m)) (vector-set! q i 0))
   ;qiの初期化
  (calcf) (printy)
  (do ((j 0 (+ j 1))) ((= j n))
  (for-each (lambda (e c) 
   (do ((i 0 (+ i 1))) ((= i m))
    (set! k (* h (vector-ref f i)))
    (set! qi (vector-ref q i))
    (set! r (e))
    (set! s (vector-ref y i))
    (vector-set! y i (+ s r))
    (set! r (- (vector-ref y i) s))
    (vector-set! q i (+ qi (* 3.0 r) (* c k))))
   (calcf))
   (list (lambda () (- (* 0.5 k) qi))
     (lambda () (* c0 (- k qi)))
     (lambda () (* c1 (- k qi)))
     (lambda () (/ (- k (* 2.0 qi)) 6.0)))
   (list -0.5 (- c0) (- c1) -0.5))
  (printy))))

(runge-kutta-gill 3 (/ (atan 1) 2.5) 20)

普通とちょっと違うのはy[0]がxであることだ. 従ってf(というのはdy/dxのことだが)のf[0]=1にしてある.

yの初期値では, y[0](xのこと)=0, y[1](yのこと)=0, y[2](dy/dxのこと)=1. fの計算では f[0]=1, f[1]=y[2], f[2]=-y[1]とする.

Runge Kuttaでは分点を1個進めるのに, x, x+h/2, x+h/2, x+hについて4回fの値を計算するが, そのループはfor-eachで回している. ループごとに変る値については, (lambda (e c) ..)の引数で渡す. 特にeの値は, 呼ばれた時に計算する必要があるので, 引数の方は(lambda () ..)の形にし, 使う時に呼出す.

initiy, calcf, printfなどのサブルーチンが, 関数runge-kutta-gillの中で定義してあるのは, 問題ありだが, yやf読み書きするので, これで我慢している.

最後の行で呼び出しているが, 0から2πまでを20分割で計算する. つまりh=2π/20; πを4*(atan 1)で計算したいからh=(atan 1)/2.5としてある.

たったの20分割なのだが, 一周したときのy[2]が.999868...だから1万分の2の程度の精度だ.

その出力を貰ってサークルを描くPostScriptのプログラムは次の通り.

/ps[[0 0 1]
[.3141592653589793 .30899155257892935 .9510578492071948]
[.6283185307179586 .5877376828378169 .8090352529734782]
[.9424777960769379 .8089575954451165 .5878333484965556]
[1.2566370614359172 .9510010098334781 .30910245672629366]
[1.5707963267948966 .9999670230159169 1.2303914619280843e-4]
[1.8849555921538759 .9510645042444495 -.30886434562367165]
[2.199114857512855 .8090808881734999 -.5876187580148391]
[2.5132741228718345 .5879134969774285 -.8088585919500648]
[2.827433388230814 .3092092738117766 -.9509316169859521]
[3.141592653589793 2.4607017746574233e-4 -.9999340319806838]
[3.4557519189487724 -.30873714204448566 -.951071143410806]
[3.7699111843077517 -.5874998314987273 -.8091265072362017]
[4.084070449666731 -.8087595818584492 -.5879936306339407]
[4.39822971502571 -.9508622132841098 -.30931607883671747]
[4.71238898038469 -.9999010268957623 -3.6909309235748844e-4]
[5.026548245743669 -.951077766707203 .30860994184321267]
[5.340707511102648 -.8091721101619074 .5873809032915227]
[5.654866776461628 -.5880737494657694 .8086605651723112]
[5.969026041820607 -.3094228718001786 .9507927987297933]
[6.283185307179586 -4.921078894070952e-4 .9998680077626148]
] def
280 360 translate
/mv {200 mul exch 200 exch mul moveto} def
/ln {200 mul exch 200 exch mul lineto} def
ps 0 get dup 1 get exch 2 get mv
1 1 20 {ps exch get dup 1 get exch 2 get ln} for stroke

図はこのようだ.

鶴亀鴉算

小学生のころ鶴亀算を習った. 例えば鶴亀合せて10匹. 足は合せて30本. 鶴と亀はそれぞれ何匹か. 全部鶴とすると足は20本. 残りの10本は鶴より2本足が多い亀によるから, 亀は5匹. 従って鶴も5匹. 足の本数を検算すると2*5+4*5=30.

矢野健太郎先生がアメリカの子供にこのような話をしたら, 鶴と亀とは見ただけで分り, 足だってまったく違うという反応だったらしい. そこで矢野先生は2ドルのケーキと4ドルのケーキを合せて10個買い, 合計を30ドルにするには2ドルと4ドルのケーキを何個ずつ買うかと話したそうだ.

小学6年の私は近くに住む中学1年の先輩と話をしていた. 先輩は「鶴をx, 亀をyとすると, x+y=10, 2x+4y=30. これを解くとx=5, y=5が得られる」と得意そうに話した. へぇーすごいなぁと思ったが, 1年後には私もその先輩と同じ中学(詳しくいえば7年制高等学校尋常科)に入学し, 同じように解けるようになった.

さて東京理科大学近代科学資料館で, 6月19日から8月8日まで「計算する器械たち —アナログコンピュータ展—」が開催される.

それに門脇廉君(現 九州大学)が作った機械式アナログ計算機の「三元連立方程式求解機」も展示される. 私は展示準備中のその機械を見たが, 非常に美しく作られていて感動した.

企画展開催中に説明を担当する理科大の学生さんたちと三元連立方程式の例題を作って解いてみたりしていたが, 見学に来た小学生に連立方程式の意味を教えるにはどうするかということになった.

そこで私が考えたのが, このブログの始めにあった鶴亀算であった. でもそれは二元であるが, もう1種類の変わった動物を登場させようと思った.

それがサッカーの3本足の鴉である. ボールを掴んでいるのが3本目の足だ.



これが三元連立方程式求解機用問題.

鶴t, 亀k, 鴉cとする. ただし鴉はサッカーJFAのシンボル 3本足のものである. 頭の数が6, 羽の数が10, 足の数が16のとき, 鶴, 亀, 鴉それぞれ何匹か.

つまり連立方程式
t+k+c=6 (0)
2t+0k+2c=10 (1)
2t+4k+3c=16 (2)
を解くわけだ.

早速三元連立方程式求解機のシミュレータでやってみる.



4枚のプレートがあり, 左からt, k, c, 定数に対応している. また制約するテープは内側から順に上の式(0), (1), (2)に対応している.

プレートの傾斜は右端の定数のを10°にしたので, 左からの傾斜は31.395°, 10.0°, 20.322°,10°である. つまり(t, k, c)=(3, 1, 2).

この図で分かるのは定数項の値が係数に較べて大きいので, 右端の定数のプーリーが分散しているのに対し, 係数の方のプーリーがごちゃごちゃと固まっていることである.

これを解決するには上の連立方程式の変数のそれぞれから1, または2を引くのである.

1を引くと
t+k+c=3 (0)
2t+0k+2c=6 (1)
2t+4k+3c=7 (2)

2を引くと
t+k+c=0 (0)
2t+0k+2c=2 (1)
2t+4k+3c=-2 (2)

2を引たので解いたのが下の図である. 各プーリーが分散している.


傾斜が左から10°, -10°, 0°, 10°になり, (t, k, c)=(1, -1, 0)であるから, 鶴が3, 亀が1, 鴉が2と判明した.

却って難しくなったかしら.

横方向の和の和

昨日のブログではさぼったnusum (2^r-1) から nusum 2^rを計算してみた.

前回のような表を書くと

n	nusum	νn	es
2r-1	nusum 2r-1	r	(r-1 r-2 r-3 ... 0)
2r	nusum 2r	1	(r)

羃の列が(r-1 r-2 r-3 ... 0)のnusum (2^r-1) の式は

nusum (2^r-1) = (r-1)2^r-2 + (r-2+2)2^r-3 + (r-3+4)2^r-4 + ... + (0+2r-2)2^-1

で, これにrを足したものが

nusum 2^r = r.2^r-1

になればよい.

nusum (2^r-1)
={各項の係数の中を計算する}
(r-1)2^r-2+r2^r-3+(r+1)2^r-4+...+(r+(r-2))2^-1
={降順の羃の列で同じ係数を持つものをならべる}
(r-1)2^r-2+(r-1)2^r-3+(r-1)2^r-4+...+(r-1)2^-1  {1行目}
             +2^r-3+2^r-4+...+2^-1                    {2行目}
                    +2^r-4+...+2^-1                    {3行目}
                          +2^r-5+...+2^-1              {4行目}
                          +...
                          +2^-1                          {r行目}
={各行を2^p+2^p-1+... +2^q=2^p+1-2^q で置き換える}
(r-1)(2^r-1-2^-1)
+2^r-2-2^-1
+2^r-3-2^-1
+...
+2⁰-2^-1
={1行目を展開し2行目以下の第1項を足し合せる 2行目以下のr-1個の2^-1を足す}
(r-1)2^r-1-(r-1)2^-1+(2^r-1-2⁰) -(r-1)2^-1
={いろいろ消えて}
r2^r-1-2^r-1-(r-1)+2^r-1-1
=r2^r-1-r.

rを足せばnusum 2^rになる.

存外簡単であった.

横方向の和の和

TAOCPの演習問題713-42にこういうのがある.

「問題42. [M21] e₁ > ... > e_r ≥ 0について, n = 2^e₁ + ... + 2^e_r なら, 和 ∑_k=0^n-1νkを, 冪e₁, ..., e_rを使って表せ.」

ここでνkはkを二進法で表わした時の1の数であり, 横方向の和といったりする.

ν0=0; ν1=1; ν2=1; ν3=2; ν4=1; ν5=2; ν6=2; ν7=3; ν8=1; ν9=2; ...

従って問題はnが与えれた時, n-1までの各整数にあった1の総和nusumを計算するものである.

n=10とすれば, 上のνの値を使いnusum 10 = 0+1+1+2+1+2+2+3+1+2=15

しかし途中のνを全部計算して求めるのはしんどいから, 10を二進法で1010と表わし, 10 = 8 + 2 = 2³ + 2¹. そのe₁=3, e₂=1からnusumを得たいという問題である.

私はこの解答にある J.-P. Allouche, J. Shallitの本 Automatic Sequences (2003)にあるというフラクタルの図が見たかったが, その本が簡単に見られそうもなかったので, PostScriptで描いてみることにした. そこで解答は次のようだ.

「解答. rについての帰納法により,

e₁2^e₁-1+ (e₂+2)2^e₂-1 + ... + (e_r+2r-2)2^e_r-1.

D. E. Knuth, Proc. IFIP Congress (1971), 1, 19--27. この和のフラクタル模様はAlloucheとShallitの本の図3.1と3.2に示してある.] またS_n'(1)も考察せよ. ただし S_n(z) = ∑_k=0^n-1z^νk = (1+z)^e₁+ z(1+z)^e₂+... +z^r-1(1+z)^e_r.」

つまりe₁=3, e₂=1だったから 3*2^3-1 + (1+2)*2^1-1 = 3*2² + 3*2⁰ = 12 + 3 = 15なわけだ.

もう少し計算してみる. 下の表は左端がn. 次がnusum. 各行のnusumは 1行上のnusumに1行上のνnを足したものになっている. その次のかっこ 内はe₁...e_r. さらにその右は上の 式による計算とその和を示す.

n	nusum	νn	es
1	0	1	(0)	0.2-1	0
2	1	1	(1)	1.20	1
3	2	2	(1 0)	1.20 + 2.2-1	2
4	4	1	(2)	2.21	4
5	5	2	(2 0)	2.21 + 2.2-1	5
6	7	2	(2 1)	2.21 + 3.20	7
7	9	3	(2 1 0)	2.21 + 3.20 + 4.2-1	9
8	12	1	(3)	3.22	12
9	13	2	(3 0)	3.22 + 2.2-1	13
10	15		(3 1)	3.22 + 3.20	15

さてこれでrによる帰納法が出来るであろうか.

偶数のnからそれより1大きいn'に進むとき, eの列は 最後にr+1項目として0が追加される. だから

nusum n+1 = e₁2^e₁-1+ (e₂+2)2^e₂-1 + ... + (e_r+2r-2)2^e_r-1 +(e_r+1+2(r+1)-2)2^{e_{r +1}-1} .

e_r+1=0とすると

nusum n+1 = e₁2^e₁-1+ (e₂+2)2^e₂-1 + ... + (e_r+2r-2)2^e_r-1 +2r.2^-1
=nusum n+r = nusum n+νn = nusum n'.

一方n=7から8へのように繰上がりがあるときはどうするか.

31から32の計算をしてみると, eは(4 3 2 1 0) から(5)になる.

このときのnusumは

4.2³+5.2²+6.2¹+7.2⁰+8.2^-1

これを次のように計算してみる.

3(23+22+21+20+2-1)
+(23+22+21+20+2-1)
+(22+21+20+2-1)
+(21+20+2-1)
+(20+2-1)
+(2-1)
=
3(24-2-1)
+(24-2-1)
+(23-2-1)
+(22-2-1)
+(21-2-1)
+(20-2-1)
=
3(24-2-1)
+(25-20)
-5(2-1)
=
3*15.5+31-2.5
=75

nusum 32は(5)だから, 5.2⁴ = 80. さっきの75にν31=5を足して 80になる.

一般的には上の(4 3 2 1 0)が(r-1 r-2 ... 0)なのでr について同様の計算をしなければならない. しかしもう出来そうになったから フラクタル図の方へ急ごう.

nusumの式をPostScriptで書くことになる.

/es {4 dict begin /s exch def /m exch def /n exch def
/l s length def
n 0 eq{s}{n 2 mod 0 eq{n 2 idiv m 1 add s es}
{n 2 idiv m 1 add m s aload pop l 1 add array astore es}
ifelse}ifelse end} def

/nusum {4 dict begin /n exch def
/ee n 0 [] es def /sum 0 def
0 1 ee length 1 sub{/i exch def /e ee i get def
/sum e i 2 mul add 2 e 1 sub exp mul sum add def} for
sum end} def

0.33 setlinewidth
0 1 560{/x exch def
x 0 moveto 0 x nusum 4 div rlineto stroke} for

最初の手続きesはnの冪の列を返す. n 0 [] esのように呼ぶ. 10 0 [] es => [3 1]

手続きnusumが総和の計算である. 32 nusum => 80.0, 560 nusum => 2480.0

n=0から560までを描いたのが下の図だ. フラクタルと言われればそう見えなくもないが.



横軸の目盛は100ごと, 縦軸の目盛は500ごとである.

ビットごとの秘法と技法 から

TAOCPの演習問題7.1.3-20にGosperのハックというのがある.

正の整数xについて

u ← x&-x;
v ← x + u;
y ← v + (((x⊕v)) »2);
としてyを返す.

とにかくやってみよう.

(define (gosper x)
 (let* ((u (band x (- (expt 2 16) x)))
        (v (+ x u))
        (y (+ v (>> (/ (bxor v x) u) 2))))
 y))

(gosper 1) => 2
(gosper 2) => 4
(gosper 4) => 8
(gosper 8) => 16
(gosper 3) => 5
(gosper 5) => 6
(gosper 6) => 9
(gosper 7) => 11
(gosper 11) => 13

という次第で, 1のビットの数が同じの次に大きい整数が得られる.

その理由を考えてみた.

x=α01^a0^b

とする. つまりxは左に適当な0と1の列αがあり, その右に0が1個, その右に1がa >0, その右に0がb ≥0個あるとする.

uの式は有名な最も右の1を取り出すものだから

u=10^b

である.

v=x + uは1の並びの右端に1を足すから

v=α10^a0^b.

x ⊕ v は左のαがキャンセルされるから

x ⊕ v=11^a0^b

になり, これをuで割るから右の0の列が消えて

(x ⊕ v)/u=11^a

つまり1が右端にa +1個ならぶ. これを2ビット右シフトするから右端に1がa -1個ならぶ. これをvに足すから

y = α10^b01^a-1     (a >0に注意)

で次のyが得られたのであった. たしかにハックだね. これはMITのHAKMEMの175番にある.