2011年11月7日月曜日

素因数探し

昨日のブログ(篩を使った素因数探し)は, TAOCPのアルゴリズムを理解するのが目的であったので, 同書のアルゴリズムの特徴であるgoto文をそのまま反映していた.

しかし, もっとSchemeらしくするにはどうするか.

前のアルゴリズムでは, 絶えずmoduloを取っているのが問題であった. あのアルゴリズムでは, 篩が2重の配列になっているので, 添字を配列の範囲に収めるためにmoduloを取るのである.

しかし, Scheme風にすると, 配列はリストになり, リストとなれば, 自転車のチェーンのように無限リストが作れる.

すなわち, (define foo '(0 1 2 3)) と設定し(図の上), (set-cdr! (list-tail foo 3) foo)
とすると, (list-tail foo 3)でfooのcdrを3回とり, 3のセルに達する. そのcdrのnilをfooに書き換えると(図の下), 3の次が0になり, 無限ループが出来る.



第2の改良点は, 篩の要素を0と1にせず, #fと#tにする. そうするとandが一発でとれる. TAOCPでは[述語]というIverson blacketを多用していて, これは述語が真のとき1, 偽のとき0になるものなので, 前回のプログラムもそうなっていた. これを(1 1 1... 1)とequal?で真偽を判定した.



lispのandは先頭からみて, 偽をみつけると途端に終了するから, この方が早いのである.

第3は, 無限リストを次々とcdr downするのに, いちいち代入するのではなく, 引数として末尾再帰で渡すことである.

このようにして書直したのが, 次のプログラムである.


(define (algorithm454d n)
(define (makesieve m n)
(let* ((x (a2b 0 m))
(x2 (map (lambda (a) (modulo (* a a) m)) x))
(s (map (lambda (b)
(if (member (modulo (- (* b b) n) m) x2) #t #f)) x)))
(set-cdr! (list-tail s (- (length s) 1)) s)
s))
(define s3 (makesieve 3 n))
(define s5 (makesieve 5 n))
(define s7 (makesieve 7 n))
(define s11 (makesieve 11 n))
(define s13 (makesieve 13 n))
(define s17 (makesieve 17 n))
(define s19 (makesieve 19 n))
(call-with-current-continuation
(lambda (throw)
(define (next x s3 s5 s7 s11 s13 s17 s19)
(if (and (car s3) (car s5) (car s7) (car s11)
(car s13) (car s17) (car s19))
(let ((y (sqrt (- (* x x) n))))
(if (integer? y) (throw (cons (+ x y) (- x y))))))
(next (+ x 1) (cdr s3) (cdr s5) (cdr s7) (cdr s11)
(cdr s13) (cdr s17) (cdr s19)))
(let ((x (inexact->exact (ceiling (sqrt n)))))
(next x
(list-tail s3 (modulo x 3))
(list-tail s5 (modulo x 5))
(list-tail s7 (modulo x 7))
(list-tail s11 (modulo x 11))
(list-tail s13 (modulo x 13))
(list-tail s17 (modulo x 17))
(list-tail s19 (modulo x 19)))))))


解が見つかったとき, 脱出するのにcall-with-current-continuationを使っているが, これが結局一番簡単なようである.

引数をぞろぞろ引き摺っていくのは, 素数の個数を変更するのに困るわけだが, とりあえずはこれでさくさく動く.

篩全体をリストにするプログラムも書いてはみたが, mapをとったりするので, 上のプログラムより遅かった. プログラムを動的に生成するという考えもあるが, 分かり難くもなり, 今はためらっている.

0 件のコメント: