このループの開始位置とループ長を求めるのがループ検出である.
そのひとつの方法にBill Gosperの考案したものがあり, あちこちに説明があるが, あまり真面目に読んだことはなかった. しかし気にはなっていたので, この数日プログラムを書いてみたりしながら, その仕掛けを調べてみた.
いまx1, x2, x3, ... が1,2,3,4,5,6,7,8,9.10,6,7,8,9,10,6,7,8,...であったとしよう.(x0から始める流儀もあり, また私はそれが好きだが, 今回の添字は1から始めることにした.) The Art of Computer Programming(TAOCP)の演習問題3.1,6の記述を借用すれば, 「x1,x2,...,xμ,..., xμ+λ-1は互いに相異なるが, n ≥μの時はxn=xn+λであるμとλがある.」
このμが開始位置, λがループ長で, 上の例ではμ=6, λ=5だ.
記憶場所も時間もふんだんにある時は, xの値をすべて順に記憶しておき, 新しいxについて, それを過去のすべてのxの値のそれぞれと比較する. そこに同じものがあれば, 前の値のあったところがループの開始場所であり, そこから今の値までの長さがループ長である.
しかし, こんなことはとてもやっていられないから, あの手この手を考えることになる. その一つにGosperの提案したものがある. 以下は添字などを私流に多少修正したものである.
これは簡単にいえば, 繰り返しが「現れるや否や」発見するのは諦める, 最新のxの1つ前, 2つ前, 4つ前, 8つ前, ... などを覚えておき, ループに入り同じ値が現れ始めたらなるべく早目にそれを発見しようという方針である.
「1つ, 2つ, 4つ, 8つ,...」から推察されるように, n個のxを調べた時は, log2 n個程度の配列Tを使う.
ところで, n>0を二進法で表した時, もっとも右の1の位置, もっとも左の1の位置をTAOCPに倣い, それぞれρ(n), &lambda(n)とする. ビットの位置は右端を0とする. (上のループ長の定義にもλがあって申し訳ない.)
ρ(1)=0, λ(1)=0; ρ(2)=1, λ(2)=1; ρ(3)=0, λ(3)=1; ρ(4)=2, λ(4)=2; ρ(5)=0, λ(5)=2; ρ(6)=1, λ(6)=2; ρ(7)=0, λ(7)=2; ρ(8)=3, λ(8)=3;ρの方は別名ルーラー関数(物差関数)という.
十進法の物差しとは様子が違い, ひと目盛りが半分ずつに分割されていて, インチの物差しが大体こういう感じである.
ρとλを計算する関数 rhoとlamは次のようだ.
ρは
(define (rho x) (define (r x rh) (if (= (modulo x 2) 0) (r (quotient x 2) (+ rh 1)) rh)) (if (= x 0) 'error (r x 0))) (map rho (a2b 1 16)) => (0 1 0 2 0 1 0 3 0 1 0 2 0 1 0)またλは関数logがあれば
(define (lam n) (inexact->exact (floor (/ (log n) (log 2)))) (map lam (a2b 1 16))=> (0 1 1 2 2 2 2 3 3 3 3 3 3 3 3)ρのように自力で計算するには
(define (lam x) (cond ((= x 0) 'error) ((= x 1) 0) (else (+ (lam (quotient x 2)) 1))))さて配列をT0, T1, T2, ... とする. (こちらは0から始める) そしてxiをTρ(i)に順に置くのである.
xiが次々と置かれていく様子は次のようだ.
i | T0 | T1 | T2 | T3 |
1 | x1 | |||
2 | x1 | x2 | ||
3 | x3 | x2 | ||
4 | x3 | x2 | x4 | |
5 | x5 | x2 | x4 | |
6 | x5 | x6 | x4 | |
7 | x7 | x6 | x4 | |
8 | x7 | x6 | x4 | x8 |
9 | x9 | x6 | x4 | x8 |
10 | x7 | x10 | x4 | x8 |
さて, 上の例のようなループがあった時の同様な配列のx10=10までの変化の様子は,
i | T0 | T1 | T2 | T3 |
1 | 1 | |||
2 | 1 | 2 | ||
3 | 3 | 2 | ||
4 | 3 | 2 | 4 | |
5 | 5 | 2 | 4 | |
6 | 5 | 6 | 4 | |
7 | 7 | 6 | 4 | |
8 | 7 | 6 | 4 | 8 |
9 | 9 | 6 | 4 | 8 |
10 | 9 | 10 | 4 | 8 |
次はx13=8だが, i=8の時に作ったT3に8があるので, 循環が始まっていたのがわかる.
そこで問題はiの行のTkを置いたのは誰かということだ. 最初の表からxを削除した表は次のようだが
i\k | 0 | 1 | 2 | 3 |
1 | 1 | |||
2 | 1 | 2 | ||
3 | 3 | 2 | ||
4 | 3 | 2 | 4 | |
5 | 5 | 2 | 4 | |
6 | 5 | 6 | 4 | |
7 | 7 | 6 | 4 | |
8 | 7 | 6 | 4 | 8 |
9 | 9 | 6 | 4 | 8 |
10 | 7 | 10 | 4 | 8 |
いろいろやってみると, iからiの下からk+1ビットをとったものの, 先頭の0,1を逆にしたものを引くといいらしいことが分かる.
iの下からk+1ビットをとるには, 2k+2-1で1のビット列を作ってマスクする.
i &(2k+2-1)
これを 2kでxorすると先頭のビットが反転する.
それをiから引くのである. Schemeで書くと
(define (m i k) (- i (fix:xor (fix:and i (- (fix:lsh 1 (+ k 1)) 1)) (fix:lsh 1 k))))であって, 実際に計算すると
(do ((i 1 (+ i 1))) ((= i 11)) (display (cons i (map (lambda (k) (m i k)) (a2b 0 (+ (lam i) 1))))) (newline)) (1 1) (2 1 2) (3 3 2) (4 3 2 4) (5 5 2 4) (6 5 6 4) (7 7 6 4) (8 7 6 4 8) (9 9 6 4 8) (10 9 10 4 8)となって上の表と一致する.
従ってi=13の時のx13>=8がi=12, k=3にあり,
プログラムを書いてやってみよう. 最初の数列xiは, μとλから
(define (f i mu lm) (if (< i mu) i (+ (modulo (- i mu) lm) mu))) (map (lambda (i) (f i 6 5)) (a2b 1 19)) => (1 2 3 4 5 6 7 8 9 10 6 7 8 9 10 6 7 8)と作る. プログラム全体と実行結果を下に示す.
(define (iloop i x) (define (kloop k) (if (> k (lam (- i 1))) (iloop i x) (begin (if (= x (vector-ref tab k)) (list i k (m (- i 1) k)) (kloop (+ k 1)))))) (vector-set! tab (rho i) x) (display i) (display (take (+ (lam i) 1) (vector->list tab))) (newline) (set! i (+ i 1)) (set! x (f i mu lm)) (kloop 0)) (define mu 6) (define lm 5) (iloop 1 (f 1 mu lm)) 1(1) 2(1 2) 3(3 2) 4(3 2 4) 5(5 2 4) 6(5 6 4) 7(7 6 4) 8(7 6 4 8) 9(9 6 4 8) 10(9 10 4 8) 11(6 10 4 8) 12(6 10 7 8) => (13 3 8)つまりi=13のとき, k=3の所に同じものがあり, その時のiは8であった. 従ってλは13-8=5であった.
登録:
コメントの投稿 (Atom)
0 件のコメント:
コメントを投稿