曜日の計算

私が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番にある.

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

最も左のビットの位置を知る法

TACPにはもっと凄いλ関数の計算法がある. 今回はそれを説明しよう.

xのサイズは知らなければならないが, それがどんなに大きくてもbignumのような2adic integer(2個進数)での処理だ.

まず準備としてあるアルゴリズムを考える. 具体的に4ビットとしよう. 4ビットだから対象は0≤x<16 である. c≤8と最も左のビットだけが1のh=8が登場する.

h|x-cは引かれる方が8≤ <16だから, 0≤x<16について, 差は8-c≤ <16-cの2回繰り返しになる.

たとえばc=3とすれば, 5,6,...,12,5,6,...,12になる. つまりそれぞれの最初のc個だけが<8になる. これにxをビットごとorすれば, 右半分はxの最も左のビットが1なので≥8になる. 従って全体では左端のc個だけが<8, それ以外は≥8になって, hでビットごとandをとると左端のc個が0, それ以外は8になる. c=3の例では0,0,0,8,8,...,8だ. このアルゴリズムで重要なのは, 4ビットの範囲の左から繰下げがないことである.

以下のアルゴリズムにはこれを多用する.

xが如何に大きくても, gを2の羃, nをg²とし, x<2ⁿの⌊lg (x)⌋を計算する.



以下の説明ではg=4, n=16とする. そうすると最初の行にあるlは2¹⁶-1, つまり1111111111111111を2⁴-1, つまり1111で割るから0001000100010001になる. 従ってhはそれを3ビット左へシフトした1000100010001000になる. 最初の行は16ビットを4ビットの区間4つのそれぞれについて, 4ビットがオール0かそうでないかを上のアルゴリズムで調べる. c=1の場合だ. 各区間はその区間がオール0か否かにより, 0か8になる. その結果がt₁だ. だからt₁は各々の区間の結果をp,q,r,sのビットで表すと, p000q000r000s000になる. p,q,r,sの少なくても1つは1である. 次の行はまずaを計算する. 2¹²-1つまり111111111111を2³-1つまり111で割るからa=001001001001である. これを先ほどのt₁に掛けるから, 積は

         p000q000r000s000
      p000q000r000s000
   p000q000r000s000
p000q000r000s000

の和,
p00pq0pqrpqrsqrs0rs00s000
222221111111111
4321098765432109876543210←ビット位置

になり, この下から12,13,14,15ビットをとるからyはpqrsにlを掛けたもの, pqrspqrspqrspqrsになる. 言い換えれば区間の情報を4ビットに詰め込んで4(=g)回コピーしたものが出来た.

今度はt₂の計算だ. bは2²⁰-1を2⁵-1で割るからb=1000010000100001となる. だから式中のbは左端の区間から順に1000, 0100, 0010, 0001で, これらをcとしてまたあのアルゴリズムを使う.

pqrs≥1000>100>10>1ならすべての区間で8になり,
1000>pqrs≥100>10>1なら0888になり,
1000>100>pqrs≥10>1なら0088になり,
1000>100>10>pqrs≥1なら0008になる.
従ってt₂は二進法では
1000100010001000
0000100010001000
0000000010001000
0000000000001000
のいずれかになり, このそれぞれからmは
1111111111111111
0000111111111111
0000000011111111
0000000000001111
になり, 4ビット右シフトしたものとxorするとmとして
1111000000000000
0000111100000000
0000000011110000
0000000000001111
が出来る.

次の行は最も左の1を持つ区間をmでマスクしてとりだし, lを掛けて左へ4個ならべ16ビットに区切ってから12ビット右シフトして右端の4ビットに置き, さらにlを掛けて4区間全体に展開する. 先ほどのpqrsから最も左の区間を選んだように, 今度は最も左のビットの位置がt₃に出来る.

整理すると最も左の1のある区間が右端にあればt₂は0008, その1つ左にあれば0088, 左端の区間なら8888である.

その区間で最も左の1が右端にあればt₃は0008, 左端にあれば8888である.

これからビットの位置λを計算するにはこうする.

t₂を5ビット右シフトするから, 0,4,44,444のいずれかになる.
t₃を7ビット右シフトするから, 0,1,11,111のいずれかになる.

これらを上下から取って足し, lを掛けて12,13,14,15ビット目をとるとそれがλである. たとえば444と1だとすると和は445になり,

   445
  445
 445
445

この4+4+5で13が得られる. 下からの繰り上げはないか心配だが, この計算は十六進であり, 最大の444と111でも5+5+5になって繰り上げはない.

このアルゴリズムは繰り返しがないのでO(1)と書いてあるが, 乗算が5回もあるから実用にはならないとTAOCPはいう.

次はDoループはあるがステップ間のジャンプはないアルゴリズムである.



d=4とするとn=d・2^d=64になる. この値を使ってアルゴリズムを調べてみよう.

B1. λ=0とした後, ⌈lg n⌉=6だから, k=5,4についてxを32ビット, 16ビットの幅に1があるかどうかでxを右シフトし, λを増やす. このステップが終わるとxは最も左のビットを含む16ビットの区間になり,　その区間がどこにあるかに従ってλは0, 16, 32, 48のいずれかになっている.

B2. 16ビットの区間に縮められたxを4つの区間にコピーする.

次にμ_d,_kというのがでてくるが, 右から2^kごとに1と0を繰り返す長さ2^dのビット列である.

μ_d,k=(2^2d-1)/(2^2k+1)

今はd=4だから
μ_4,3=0000000011111111
μ_4,2=0000111100001111
μ_4,1=0011001100110011
μ_4,0=0101010101010101

μ_4,3の0の部分を左半分, 1の部分を右半分といおう. 同様に
μ_4,2の0の部分を左四半分, 1の部分を右四半分といい,
μ_4,1の0の部分を左八半分, 1の部分を右八半分といい,
μ_4,0の0の部分を左十六半分, 1の部分を右十六半分といおう.

すると下の図で灰色の部分が左なんとか半分, 白い部分が右なんとか半分である.


B3. これらのμを繋げてビットごとにnotをとったものとxのandをyとする. つまりxの左なんとか半分の1をとりだす.

B4. x⊕yは右なんとか半分の1になるわけだ. そしてh=8000800080008000としておなじみのアルゴリズムを使う.

そうすると左部分をアルゴリズムのx, 右部分をcとすることになるので, x<cの時に0, x≥cの時に8だったように, 左≥右なら0, そうでなければ8000になる.

上の図で最上段の0〜fはビット番号で, そのビット6に最も左の1があるとすると上から順に
左<右
左>右
左>右
左<右
だから, B4のtはそれぞれの区画で0, 8000, 8000, 0になる.

B5. ⌈lgd⌉=2なので, 0,1のkについて
2⁴-2⁰=15ビット,
2⁵-2¹=30ビット
ずつ左シフトして加える. 8の立つビットを左からp,q,r,sとするとそれらのビット位置はそれぞれp=63, q=47, r=31, s=15である.

まず15ビット左へ移るから, qは62, sは30になりpq, rsは並ぶ. それをさらに30ビット左シフトすると31,30にあったrsは61,60になり, pqrsは並ぶ.

2⁶⁴の剰余をとるから64ビットレジスタの左端に出来る.

先ほどの例ではその4ビットに0110, つまり最も左のビット位置の6ができる. それをn-d=60ビット右シフトし, すでに得られていたλと足すと最も左の1のビット位置が得られるのである.

ところで理解して頂けたであろうか. 私はまずSchemeでプログラムを書いて実行したりして, 様子を見ながらアルゴリズムがどうなっているかを考えた. 分かってしまえば簡単なのだが, 最初に述べたアルゴリズムが何をしているか分からないうちは, なんとも不思議に思うアルゴリズムであった.