2015年3月1日日曜日

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

2月21日のブログで疑似素数発生のプログラムが出来たから, 次は素因数分解のプログラムを書こう. 私がSchemeの環境で使っているfactorizeはこういうものだ.
 (factorize 32) => (2 2 2 2 2)
 (factorize 12) => (3 2 2)
 (factorize 108) => (3 3 3 2 2)
 (factorize 97) => (97)
だから素数ならこのリストの長さが1になるわけで
 (define (prime? n)
  (= (length (factorize n)) 1))

 (map prime? (a2b 2 10)) => (#t #t #f #t #f #t #f #f)
だが, 問題はn=1の時だ.
 (factorize 1) => (1)
と1も素数に判定されてしまう. それはさておき...

個人用電卓Happy Hacking Calculatorでは, n>=2の素因数分解は, pをnの最小の因数としてn=pc*q となるようなp, c, qをスタックに積むようにした. つまり
 n fact -> p, c, q
 2 fact -> 2, 1, 1
 4 fact -> 2, 2, 1
 6 fact -> 2, 1, 3
 8 fact -> 2, 3, 1
12 fact -> 2, 2, 3
のようになる. 計算が終了した時, 最小素因数pが見える. それが最初に入れた数nなら素数だった. そうでないなら, その最小素因数が何回あったかcがスタックのすぐ次で分かる. さらに次に最小素因数で割り切った残りqがあるから, 続けてその素因数分解を始める. だから, 108(=2 * 2 * 3 * 3 * 3)は
108 fact -> 2, 2, 27 popを2回
27 fact -> 3, 3, 1
で108=22*33ということが分るのである.

これが便利だったので, HP-16Cのプログラムでも同様なインターフェースとした.

Happy Hacking Calculatorの素因数分解のアルゴリズムをSchemeで示すと
(define (hhcfactorize n)
 (let ((fl 疑似素数差分リスト) (d 2) (c 0) (i 0))
  (define (floop)
   (define (gloop) (set! n (/ n d)) (set! c (+ c 1))
    (if (= (modulo n d) 0) (gloop)))
   (if (< (quotient n d) d) (list n 1 1)
    (if (= (modulo n d) 0) (begin (gloop) (list d c n))
     (begin (if (= n 485) (set! i 5))
       (set! d (+ d (list-ref fl i))) (set! i (+ i 1)) (floop)))))
(floop)))
のように書いたので, 今回はこれをHP-16C用に書き換える. 上のプログラムでは配列flが疑似素数表だが, それに前回のブログのものを使う.
001 43.22.A LBL A
002 44.0    STO 0 ;n->0
003 6       6     ;差分を設定
004 44.01   STO 1
005 44.03   STO 3
006 4       4
007 44.04   STO 4
008 44.06   STO 6
009 44.08   STO 8
010 2       2
011 44.05   STO 5
012 44.07   STO 7
013 44.09   STO 9
014 44.A    STO A
015 44.d    STO D ;2->d
016    1    1
017 44.b    STO B
018 43.35   CLx
019 44.c    STO C ;0->c
020 1       1
021 1       1
022 44.32   STO I ;11->i
023 43.22.0 LBL 0 ;ループバック地点
024 45.0    RCL 0 ;n
025 45.d    RCL D ;d
026 10      /
027 45.d    RCL D ;d in x n/d in y
028 43.01   x<=y ;d<=n/d?
029 22.01   GTO 1
030 45.b    RCL B ;商が除数より小さくなった
031 45.b    RCL B ;1,1,nを積んでもどる
032 45.00   RCL 0 ;n
033 43.21   RTN
034 43.22.1 LBL 1 ;n/d>=d
035 43.06.4 F? 4  ;剰余!=0フラッグ
036 22.03   GTO 3
037 43.22.2 LBL 2 ;divisible
038 45.00   RCL 0
039 45.d    RCL D
040 10      /
041 44.00   STO 0 ;n/d->n
042 45.c    RCL C
043 1       1
044 40      +
045 44.c    STO C ;c+1->c
046 45.0    RCL 0 ;n
047 45.d    RCL D ;d
048 42.09   RMD ;remainder
049 43.40   x=0
050 22.02   GTO 2
051 45.00   RCL 0 ;n
052 45.c    RCL C ;c
053 45.d    RCL D ;d
054 43.21   RTN
055 43.22.3 LBL 3 ;割り切れなかったので次の除数を作る
056 45.d    RCL D
057 45.31   RCL (I)
058 40      +
059 44.d    STO D
060 43.23   DSZ
061 22.0    GTO 0
062 8       8
063 44.32   STO I
064 22.00   GTO 0
017までは1からBのレジスタに差分を入れる. 2の時には最初の除数もDへ入れる(015). 019でカウンタCをクリアする. 022でIレジスタの初期値を11に設定. 023のLBL 0がループバックの場所.

n/d<dを調べ, 終わりなら030から1,1,nをスタックに積んでRTN.

035は026の除算で剰余が0でない, 割り切れなかった時はフラッグ4が立っているから, F? 4で調べ, 立っていればLBL 3へ飛ぶ.

050までは同じ除数で割れるだけ割り(cも増やしつつ), 割れなくなったら051からn, c, dをスタックに積んでRTN.

そこから下は前回の手法による除数の更新である.

このプログラムが消費したメモリーをしらべると, プログラムは64バイトだが, 7の倍数ずつ領域を確保するから, 70バイト使う. レジスタは0からDまで32ビット= 4バイトを14個つかったら56バイト, 合わせて126バイトだったから, 全体で203バイトのメモリーの62パーセントも消費した.

このプログラムをHP-16Cの実機で実行すると, 不思議なことに分解する数の大きさにあまり関係なく30秒ちょっと程度で計算できる. DM-16はさすがに現代の電卓だけあって, 2,3秒で計算する. ただ誰もがいうようにこの電卓のキーの押し心地は最低だ.

前回紹介したシミュレータは大きい数はかなりの時間がかかり, 実用にならぬ. しかし除数が徐々に増えていくのが窓に見えるので待つのも苦にならない.

やはり私の作った個人用電卓の組込み機能が最高だなぁ.

0 件のコメント: