HP-16Cのプログラム技法

HP-16Cにあって私の個人用電卓にない機能のひとつが実数演算である. なにしろHP電卓設計チームのひとりがIEEE 754の代表のWillian Kahanなのだからあるのは当たり前か.

しかしHP-16Cの内部表現は56ビットということを別として, 私には分かっていることもまだ分からないこともある.

整数から実数に演算を切り換えると, y*2^xなる値がxレジスタに入る. 反対に実数から整数に切り換えると, その時のxの値がy*2^xなるようにxとyが設定される.

たとえば円周率3.141592...は

(* 1686629713. (expt 2 -29)) => 3.1415926534682512

だから, 1686629713 ENTER 29 CHS FLOAT とすると窓に 3.141592 00 が現れる.

ここでまたDECにもどすと, xが-30, (normalizeされた)yが 3373259425 となり

(* 3373259425. (expt 2 -30)) => 3.1415926525369287

実数演算でなにかやってみよう. 円周率の連分数展開などが面白そうだ.

3.141592 これから 整数部の3を引く.
1.415926 -01 1/xをとる
7.062513 00 7を引く
6.251333 -02 1/xをとる
1.599658 01 15を引く
9.965869 -01 1/xをとる
1.003424 00 1を引く
3.424719 -03 1/xをとる
2.919947 02 ...

という次第で連分数は3,7,15,1,291,...

したがってπの近似値はまず3, 次は3+1/7=22/7, 次は3+1/(7+1/15) = 333/106, 次は有名な113/355. 捨てたのが1/292という項だから精度は抜群である.

さてこれをプログラムに書き直そう. 1/xはそういう命令があるから簡単だ. 問題は逆数をとったあと, 整数部を引くことである. この電卓にはfloorとかroundの命令は見当たらない.

いろいろ考えた末, 最初の例にあるように, yにある整数に, 2のxにある整数乗を掛けたのが実数だから, yにある整数をxにある整数の絶対値だけ右シフトすれば整数部が得られることが分った.

ということで書いた連分数展開のプログラムが以下のである.

001 43.22.C  LBL C
002 44.0     STO 0 ;0に格納しておく
003 24       DEC   ;整数計算に
004 44.32    STO I ;羃数をIに於く
005 33       R↓  ;yを取り出す
006 43.22.0  LBL 0 ;ループバック地点
007 42.b     SR    ;右1ビットシフト
008 43.24    ISZ   ;インデックス-1 結果0?
009 22.0     GTO 0
010 31       R/S   ;一旦停止
011 0        0     ;羃を0にして
012 42.45.48 float.;実数演算に
013 45.0     RCL 0 ;元の数を取り出し
014 34       x<>y  ;x,yを交換して
015 30       -     ;引く
016 43.26    1/x   ;逆数をとる
017 22.C     GTO C ;次の桁の計算に戻る

FLOATにしてからレジスタに3.141592....を入力し, GSB Cとするとしばらくrunningした後3が現れて停止する. 値を読み取り, R/Sを押すとまた走りだす. 後は7, 1, 291と続く.

TAOCPの第2巻に連分数の例がある.

π=3+//7,15,1,292,1,1,1,2,....//
e=2+//1,2,1,1,4,1,1,6,...//
φ=1+//1,1,1,1,.....//

8/29=//3,1,1,1,2//
√8/29=//1,1,9,2,2,3,2,2,9,1,2,1,9...//
など

これらを計算してみると楽しいが, 始めの円周率の例で見たように, 292のところがそろそろ限界で, HP-16Cの桁数(56ビット)では破綻し始める.