整数立方根

平方根に較べ, 立方根はおおいに疎外されている感があるが, 立方根を計算したいこともある.

その場合, 実用的には
(expt 9270 (/ 1 3)) => 21.00680051861007
なことで済ませてしまう.

さて, Warrenの「ハッカーの楽しみ」を眺めていたら, 整数立方根という話題があった(訳書の226ページ). こういうプログラムが書いてある(32ビット用だ).



同書の巻頭の「推薦の辞」に私が書いたように, 本書にはそのアルゴリズムでよい理由が殆んど書いてない.

このアルゴリズムは下の説明のように出来ているのだ.

例によってSchemeに書き直すと




最初のsの設定が30でなく18なのは, そんなに大きい数でテストすることもないからだ. また最後にyだけでなく, xも返すのは, 余りも欲しいからである.

途中の経過を見るdisplayが3行ある. 9270の立方根を計算してみる.



すなわち, 立方根は21で, 剰余が9. 整数の範囲で開立をやめ, 剰余が得られるから, 整数立方根というわけである. 同じものを十進法で書くと,



つまりbの始めの方は, xより小さくなるまで2¹⁸, 2¹⁵, 2¹²,... と減っていき, 4096がxより初めて小さくなった時点でyを1とする. yは最後は二進で10101だが, その最初の1(=2⁴)が決ったのである.

xの方はそれを引いて, 9270-4096=5174になった.

次はyのその下のビットを1にするか, 0にするかを決めることである. y=16の時, その下のビットは8だから(16+8)³=16³+3*16²*8+3*16*8²+8³=13824

しかし, 最初の4096はすでに引いてあるので, 比較するbは, 上の4段目の13824-4096=9728である. こうするよりは, y=2, (y+1)³-y³=3y(y+1)+1を作り, 左シフトする方が楽というのが, プログラムの趣旨であろう.

8のビットは立たなかった. 次は4のビットで, (16+4)³-16³=3904である. これはxより小さいからyは1増えて101になった.

このようにしてyのビットを上から次々を立てていき, 1ビットごとにきめていく. 最後に残ったxが開立の剰余である.

このプログラムの場合は, 方針は最初から見え見えであったが.