本日のお題
フィボナッチ数列を一般項から求めるプログラムを Scheme で書いてみます。
フィボナッチ数列の一般項は、
面倒くさい数式を頑張って書きました。褒めてくださいｗ
ポイント
このお題を解くにあたってのポイントは、一般項に出てくる無理数（√５）をいかに誤差が出ないように計算するかです。
プログラム
汎用の手続きから。一般項に出てくるｎ乗の処理には逐次平方を使用して高速化します。
(define nil '()) ;; general (define (square x) (* x x)) (define (divisible? x y) (= (remainder x y) 0)) (define (negated x) (* -1 x)) (define (enumerate-interval low high) (if (> low high) nil (cons low (enumerate-interval (+ low 1) high)))) (define (filter predicate sequence) (cond ((null? sequence) nil) ((predicate (car sequence)) (cons (car sequence) (filter predicate (cdr sequence)))) (else (filter predicate (cdr sequence))))) (define (square-sqrt x) (mul x x)) (define (fast-expt b n) (fast-expt-iter b n (make-sqrt nil 1))) (define (fast-expt-iter b n a) (cond ((= n 0) a) ((even? n) (fast-expt-iter (square-sqrt b) (/ n 2) a)) (else (fast-expt-iter b (- n 1) (mul b a)))))
無限ストリームの手続き。後述の素数判定で使用します。
;; streams (define the-empty-stream nil) (define (stream-null? stream) (null? stream)) (define (stream-car stream) (car stream)) (define (stream-cdr stream) (force (cdr stream))) (define (cons-stream a b) (cons a b)) (define (stream-filter pred stream) (cond ((stream-null? stream) the-empty-stream) ((pred (stream-car stream)) (cons-stream (stream-car stream) (delay (stream-filter pred (stream-cdr stream))))) (else (stream-filter pred (stream-cdr stream))))) (define (integers-starting-from n) (cons-stream n (delay (integers-starting-from (+ n 1)))))
素数判定の処理です。無限ストリームにエラトステネスの篩を使用して素数の無限ストリームを生成します。
;; primes stream (define primes (cons-stream 2 (delay (stream-filter prime? (integers-starting-from 3))))) (define (prime? n) (define (iter ps) (cond ((> (square (stream-car ps)) n) #t) ((divisible? n (stream-car ps)) #f) (else (iter (stream-cdr ps))))) (iter primes)) (define (fact-prime-nums num) (define (iter n lis) (if (= n 1) lis (let ((pn (stream-car (stream-filter (lambda (x) (divisible? n x)) primes)))) (iter (/ n pn) (cons pn lis))))) (iter num nil))
後述で定義するパッケージで手続きを登録するテーブルを用意します。
;; tagging (define (attach-tag type-tag contents) (if (eq? type-tag 'scheme-number) contents (cons type-tag contents))) (define (type-tag datum) (cond ((pair? datum) (car datum)) ((number? datum) 'scheme-number) (else (error "Bad tagged datum -- TYPE-TAG" datum)))) (define (contents datum) (cond ((pair? datum) (cdr datum)) ((number? datum) datum) (else (error "Bad tagged datum -- CONTENTS" datum)))) ;; table (define *table* nil) (define (get op type) (let ((item (filter (lambda (x) (and (eq? (car x) op) (equal? (cadr x) type))) *table*))) (if (null? item) #f (caddr (car item))))) (define (put op type item) (set! *table* (cons (list op type item) (filter (lambda (x) (not (and (eq? (car x) op) (equal? (cadr x) type)))) *table*)))) (define (apply-generic op . args) (let ((type-tags (map type-tag args))) (let ((proc (get op type-tags))) (if proc (apply proc (map contents args)) (error "No method for these types -- APPLY-GENERIC" (list op type-tags)))))) ;; generic arithmetic operations (define (add x y) (apply-generic 'add x y)) (define (sub x y) (apply-generic 'sub x y)) (define (mul x y) (apply-generic 'mul x y)) (define (div x y) (apply-generic 'div x y)) (define (equ? x y) (apply-generic 'equ? x y)) (define (irrational? x) (apply-generic 'irrational? x))
平方根を扱うパッケージを用意します。平方根用の四則演算を実装しています。このパッケージを使用することで、C 言語風に言うと「平方根型」が使えるようになります。
(define (install-sqrt-package) ;; private (define (irrationals x) (car x)) (define (constant x) (cdr x)) (define (roots x) (map car (irrationals x))) (define (coefs x) (map cdr (irrationals x))) (define (root ir) (car ir)) (define (coef ir) (cadr ir)) (define (make-irrat root coef) (list root coef)) (define (make-sqrt irrats const) (define (iter irs new-irs c) (if (null? irs) (cons new-irs c) (let ((a (coef (car irs))) (x (root (car irs)))) (cond ((or (= a 0) (= x 0)) (iter (cdr irs) new-irs c)) ((= x 1) (iter (cdr irs) new-irs (+ c a))) ((integer? (sqrt x)) (iter (cdr irs) new-irs (+ c (* a (sqrt x))))) (else (let ((sf (square-factors x))) (iter (cdr irs) (append new-irs (list (make-irrat (root sf) (* a (coef sf))))) c))))))) (iter irrats nil const)) (define (square-factors x) (define (iter facts p sf) (if (null? facts) (make-irrat (/ x (apply * (map square sf))) (apply * sf)) (if (= p (car facts)) (iter (cdr facts) 0 (append sf (list (car facts)))) (iter (cdr facts) (car facts) sf)))) (iter (fact-prime-nums x) 0 nil)) (define (add-sqrt x y) (define (add-iter x-irs y-irs) (cond ((null? x-irs) y-irs) ((null? y-irs) x-irs) ((= (root (car x-irs)) (root (car y-irs))) (append (list (make-irrat (root (car x-irs)) (+ (coef (car x-irs)) (coef (car y-irs))))) (add-iter (cdr x-irs) (cdr y-irs)))) ((< (root (car x-irs)) (root (car y-irs))) (append (list (car x-irs)) (add-iter (cdr x-irs) y-irs))) (else (append (list (car y-irs)) (add-iter x-irs (cdr y-irs)))))) (make-sqrt (add-iter (irrationals x) (irrationals y)) (+ (constant x) (constant y)))) (define (sub-sqrt x y) (add-sqrt x (negated-sqrt y))) (define (mul-sqrt x y) (define (mul-cst c sqrt) (define (iter irs) (if (null? irs) nil (append (list (make-irrat (root (car irs)) (* c (coef (car irs))))) (iter (cdr irs))))) (make-sqrt (iter (irrationals sqrt)) 0)) (define (mul-y-irs x-irrational y-irrationals) (define (iter x-ir y-irs new-irs c) (if (null? y-irs) (list new-irs c) (iter x-ir (cdr y-irs) (append new-irs (list (make-irrat (* (root x-ir) (root (car y-irs))) (* (coef x-ir) (coef (car y-irs)))))) c))) (apply make-sqrt (iter x-irrational y-irrationals nil 0))) (define (mul-x-irs x-irs) (if (null? x-irs) (add-sqrt (add-sqrt (mul-cst (constant x) y) (mul-cst (constant y) x)) (make-sqrt nil (* (constant x) (constant y)))) (add-sqrt (mul-y-irs (car x-irs) (irrationals y)) (mul-x-irs (cdr x-irs))))) (mul-x-irs (irrationals x))) (define (div-sqrt x y) (if (constant? y) (mul-sqrt x (make-sqrt nil (/ 1 (constant y)))) (if (rationalize? y) (mul-sqrt x (reciprocal-sqrt y)) (/ (sqrt->float x) (sqrt->float y))))) (define (negated-sqrt x) (define (iter irs) (if (null? irs) nil (append (list (make-irrat (root (car irs)) (* -1 (coef (car irs))))) (iter (cdr irs))))) (make-sqrt (iter (irrationals x)) (* -1 (constant x)))) (define (constant? x) (= (length (irrationals x)) 0)) (define (rationalize? x) (if (= 0 (constant x)) (<= (length (irrationals x)) 2) (<= (length (irrationals x)) 1))) (define (reciprocal-sqrt x) (define (simple-reciprocal s) (let ((irrat (car (irrationals s)))) (let ((a (coef irrat)) (x (root irrat))) (make-sqrt (list (make-irrat x (/ 1 (* a x)))) 0)))) (define (rationalized-reciprocal s) (let ((irrat (irrationals s)) (c (constant s))) (if (= c 0) (let ((a (coef (car irrat))) (x (root (car irrat))) (b (coef (cadr irrat))) (y (root (cadr irrat)))) (let ((d (- (* (square a) x) (* (square b) y)))) (make-sqrt (list (make-irrat x (/ a d)) (make-irrat y (/ (negated b) d))) 0))) (let ((a (coef (car irrat))) (x (root (car irrat)))) (let ((d (- (* (square a) x) (square c)))) (make-sqrt (list (make-irrat x (/ a d))) (/ (negated c) d))))))) (if (and (= (constant x) 0) (= (length (irrationals x)) 1)) (simple-reciprocal x) (rationalized-reciprocal x))) (define (sqrt->float x) (define (iter irs sum) (if (null? irs) (+ sum (constant x)) (iter (cdr irs) (+ sum (* (coef (car irs)) (sqrt (root (car irs)))))))) (iter (irrationals x) 0)) (define (sqrt->integer x) (inexact->exact (floor (sqrt->float x)))) ;; interfaces (define (tag x) (attach-tag 'sqrt x)) (put 'add '(sqrt sqrt) (lambda (x y) (tag (add-sqrt x y)))) (put 'sub '(sqrt sqrt) (lambda (x y) (tag (sub-sqrt x y)))) (put 'mul '(sqrt sqrt) (lambda (x y) (tag (mul-sqrt x y)))) (put 'div '(sqrt sqrt) (lambda (x y) (tag (div-sqrt x y)))) (put 'equ? '(sqrt sqrt) (lambda (x y) (and (equal? (irrationals x) (irrationals y)) (= (constant x) (constant y))))) (put 'irrational? '(sqrt) (lambda (x) (not (eq? (irrationals x) nil)))) (put 'make 'sqrt (lambda (irrats const) (tag (make-sqrt irrats const)))) (put 'float 'sqrt (lambda (x) (sqrt->float (contents x)))) (put 'integer 'sqrt (lambda (x) (sqrt->integer (contents x)))) 'done) (define (make-sqrt irrats const) ((get 'make 'sqrt) irrats const)) (define (sqrt->float x) ((get 'float 'sqrt) x)) (define (sqrt->integer x) ((get 'integer 'sqrt) x)) ;; install packages (install-sqrt-package)
以上です。これで計算の準備が整いました。
軽くテスト
√1.5 を計算してみます。
gosh> (div (make-sqrt '((3 1)) 0) (make-sqrt '((2 1)) 0)) (sqrt ((6 1/2)) . 0)
これでは人間にはわかりにくいですね。sqrt->float 手続きを使用して浮動小数点数に変換してみます。
gosh> (sqrt->float (div (make-sqrt '((3 1)) 0) (make-sqrt '((2 1)) 0))) 1.224744871391589
来春発売予定のニーアレプリカントのバージョンアップ版は 1.22474487139... です。合っていますねｗ
(make-sqrt '((3 1)) 0)
は 1 * √3 + 0 を表しています。つまり、√3 です。いろいろ手抜きなので係数の 1 と定数項の 0 は省略できませんｗ また、いきなり分数は書けない仕様なので div 手続きを使って表現します。先の√1.5 は、√3 / √2 です。有理化すると √6 / 2 です。このように計算中では平方根型で表現するので誤差が出ません。
gosh> (add (make-sqrt '((2 1)) 0) (make-sqrt '((3 1)) 0)) (sqrt ((2 1) (3 1)) . 0)
この例では √2 と √3 を足していますが、平方根型では計算を無理に行わずに保留することで誤差がでない仕組みになっています。人間が筆記で計算するときと同じ方法です。
gosh> (sqrt->float (add (make-sqrt '((2 1)) 0) (make-sqrt '((3 1)) 0))) 3.1462643699419726
平方根型を数値（浮動小数点数）に変換すると、無理数の値を永遠に表示することは不可能ですので誤差が生じます。
一般項の実装
さて、いよいよ本題の一般項を実装してフィボナッチ数列を生成してみます。
(define (fib n) (let ((a (fast-expt (div (make-sqrt '((5 1)) 1) (make-sqrt nil 2)) n)) (b (fast-expt (div (sub (make-sqrt nil 1) (make-sqrt '((5 1)) 0)) (make-sqrt nil 2)) n))) (sqrt->integer (div (sub a b) (make-sqrt '((5 1)) 0)))))
実行結果
フィボナッチ数列の 10 番目まで求めてみます。
gosh> (map fib (enumerate-interval 1 10)) (1 1 2 3 5 8 13 21 34 55)
問題なさそうです。フィボナッチ数列10000 番目を求めてみます。ついでに処理時間も計ります。
gosh> (time (fib 10000)) ;(time (fib 10000)) ; real 0.001 ; user 0.000 ; sys 0.000 33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088830269235673613135117579297437854413752130520504347701602264758318906527890855154366159582987279682987510631200575428783453215515103870818298969791613127856265033195487140214287532698187962046936097879900350962302291026368131493195275630227837628441540360584402572114334961180023091208287046088923962328835461505776583271252546093591128203925285393434620904245248929403901706233888991085841065183173360437470737908552631764325733993712871937587746897479926305837065742830161637408969178426378624212835258112820516370298089332099905707920064367426202389783111470054074998459250360633560933883831923386783056136435351892133279732908133732642652633989763922723407882928177953580570993691049175470808931841056146322338217465637321248226383092103297701648054726243842374862411453093812206564914032751086643394517512161526545361333111314042436854805106765843493523836959653428071768775328348234345557366719731392746273629108210679280784718035329131176778924659089938635459327894523777674406192240337638674004021330343297496902028328145933418826817683893072003634795623117103101291953169794607632737589253530772552375943788434504067715555779056450443016640119462580972216729758615026968443146952034614932291105970676243268515992834709891284706740862008587135016260312071903172086094081298321581077282076353186624611278245537208532365305775956430072517744315051539600905168603220349163222640885248852433158051534849622434848299380905070483482449327453732624567755879089187190803662058009594743150052402532709746995318770724376825907419939632265984147498193609285223945039707165443156421328157688908058783183404917434556270520223564846495196112460268313970975069382648706613264507665074611512677522748621598642530711298441182622661057163515069260029861704945425047491378115154139941550671256271197133252763631939606902895650288268608362241082050562430701794976171121233066073310059947366875
ほぼ瞬時に計算結果が出力されました。
検算
普通にフィボナッチ数列を求める処理で検算してみます。「メモ化」という技法を使用して高速化しています。
(define nil '()) (define (make-table) (let ((local-table (list '*table*))) (define (assoc key records) (cond ((null? records) #f) ((equal? key (caar records)) (car records)) (else (assoc key (cdr records))))) (define (search keys) (define (iter k t) (let ((r (assoc (car k) (cdr t)))) (cond ((equal? r #f) t) ((not (list? (cdr r))) r) (else (iter (cdr k) r))))) (iter keys local-table)) (define (lookup keys) (define (iter k t) (cond ((equal? t #f) #f) ((null? k) (cdr t)) (else (iter (cdr k) (assoc (car k) (cdr t)))))) (iter keys local-table)) (define (make-elem keys value) (if (null? (cdr keys)) (cons (car keys) value) (list (car keys) (make-elem (cdr keys) value)))) (define (last keys value) (define (iter r c) (if (not (list? r)) c (iter (cadr r) (+ c 1)))) (make-elem (tail keys (iter (s keys) 0)) value)) (define (head l c) (if (or (null? l) (<= c 0)) nil (cons (car l) (head (cdr l) (- c 1))))) (define (tail l c) (reverse (head (reverse l) c))) (define (insert! keys value) (let ((r (search keys))) (if (list? r) (if (equal? (car r) '*table*) (set-cdr! r (cons (make-elem keys value) (cdr r))) (let ((rr (last keys value))) (set-cdr! r (cons rr (cdr r))) )) (set-cdr! r value))) local-table) (define (print-table) local-table) (define (dispatch m) (cond ((eq? m 'lookup) lookup) ((eq? m 'insert!) insert!) ((eq? m 'print-table) print-table) (else (error "Unknown operation -- TABLE" m)))) dispatch)) (define (lookup k t) ((t 'lookup) (list k))) (define (insert! k v t) ((t 'insert!) (list k) v)) (define (print-table t) ((t 'print-table))) (define (memoize f) (let ((table (make-table))) (lambda (x) (let ((previously-computed-result (lookup x table))) (or previously-computed-result (let ((result (f x))) (insert! x result table) ; (print (print-table table)) result)))))) (define (fib n) (cond ((= n 0) 0) ((= n 1) 1) (else (+ (fib (- n 1)) (fib (- n 2)))))) (define memo-fib (memoize (lambda (n) (cond ((= n 0) 0) ((= n 1) 1) (else (+ (memo-fib (- n 1)) (memo-fib (- n 2))))))))
fib 手続きは愚直な実装です。memo-fib がメモ化版です。
fib 手続きのほうは、n=100 の時点でいつ応答が返るのかわからないぐらい遅いです。
gosh> (time (memo-fib 10000)) ;(time (memo-fib 10000)) ; real 2.206 ; user 2.188 ; sys 0.016 33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088830269235673613135117579297437854413752130520504347701602264758318906527890855154366159582987279682987510631200575428783453215515103870818298969791613127856265033195487140214287532698187962046936097879900350962302291026368131493195275630227837628441540360584402572114334961180023091208287046088923962328835461505776583271252546093591128203925285393434620904245248929403901706233888991085841065183173360437470737908552631764325733993712871937587746897479926305837065742830161637408969178426378624212835258112820516370298089332099905707920064367426202389783111470054074998459250360633560933883831923386783056136435351892133279732908133732642652633989763922723407882928177953580570993691049175470808931841056146322338217465637321248226383092103297701648054726243842374862411453093812206564914032751086643394517512161526545361333111314042436854805106765843493523836959653428071768775328348234345557366719731392746273629108210679280784718035329131176778924659089938635459327894523777674406192240337638674004021330343297496902028328145933418826817683893072003634795623117103101291953169794607632737589253530772552375943788434504067715555779056450443016640119462580972216729758615026968443146952034614932291105970676243268515992834709891284706740862008587135016260312071903172086094081298321581077282076353186624611278245537208532365305775956430072517744315051539600905168603220349163222640885248852433158051534849622434848299380905070483482449327453732624567755879089187190803662058009594743150052402532709746995318770724376825907419939632265984147498193609285223945039707165443156421328157688908058783183404917434556270520223564846495196112460268313970975069382648706613264507665074611512677522748621598642530711298441182622661057163515069260029861704945425047491378115154139941550671256271197133252763631939606902895650288268608362241082050562430701794976171121233066073310059947366875
同じ値が出力されます。大丈夫そうです。
処理時間ですが、普通に求める方法ではメモ化しても一般項から求める方法に勝てませんでした。一般項で求める方法はとても速い事がわかります。
まとめ
・フィボナッチ数列を一般項から求める処理を Scheme で実装して、フィボナッチ数列の 10000 番目の値を求めてみました。
・新しく「平方根型」を実装しました。平方根型では有理数を扱うため誤差は生じません。
・平方根型を数値（浮動小数点数）に変換したい場合は計算の最後で行いますが、無理数の小数部を永遠に表示することは不可能ですので誤差が生じます。
