導入
フィボナッチ数列の 10000 番目の値は？
考え方
一般的な漸化式による方法では時間がかかりすぎます。「メモ化」という技法を使えばそれなりになんとかなりそうですが、それでも計算が遅いので一般項を使って求めます。
一般項
フィボナッチ数列の一般項は以下の式で表されます。
Fn=\frac{1}{\sqrt{5}}\left\{\Big(\frac{1+\sqrt{5}}{2}\Big)^n-\Big(\frac{1-\sqrt{5}}{2}\Big)^n\right\}
プログラム
まずは基礎的な手続き群を定義します。
(define nil '())
(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 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)))))
(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))
これから実装する「平方根型」手続きはハッシュテーブルに登録するようにしますが、このテーブルに登録するキーを「タグ」、値を「コンテンツ」と呼び、このタグとコンテンツを扱うための手続きを用意しておきます。
(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))))
ハッシュテーブルの操作を行う手続きを定義します。取り出しは get 手続き、登録は put 手続きで行います。テーブル内の手続きの呼び出しは apply-generic 手続きを使います。
;; 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))))))
これから実装する「平方根型」の四則演算の手続きを定義します。それぞれ apply-generic 手続きを使用して、ハッシュテーブル内に登録済みの手続きを呼び出すようにします。
(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))
メインとなる平方根型の定義です。平方根の四則演算に必要な手続きを定義します。公開したい手続きはハッシュテーブルに登録します。
(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))
平方根型に対応した n 乗する手続きを定義します。「逐次平方」というテクニックを使用して計算時間を短縮します。
(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)))))
最後に「平方根型」を使用できるように、install-sqrt-package 手続きを評価します。
(install-sqrt-package)
以上で準備が整いました。
平方根の表現と計算のしかた
平方根型で
\sqrt{3}
は
1*\sqrt{3}+0
と表し以下のように記述します。
(make-sqrt '((3 1)) 0)
平方根の計算してみます。
\sqrt{2}+\sqrt{3}
平方根型を使用すると次のように記述でき、計算結果としてはルート部分をまとめて持つようにします。
gosh> (add (make-sqrt '((2 1)) 0) (make-sqrt '((3 1)) 0))
(sqrt ((2 1) (3 1)) . 0)
平方根型では無理に計算することはせず、人が計算するときのような動作をすることで、誤差を抑える仕組みになっています。掛け算もやってみます。
\sqrt{2}*\sqrt{2}
gosh> (mul (make-sqrt '((2 1)) 0) (make-sqrt '((2 1)) 0))
(sqrt () . 2)
今回の実装では簡略化は行っていないため平方根型のままになっていますが、計算結果はルート部分が無くなるので 2 としたほうが良いでしょう。
一般項を求める手続きを定義する
では、これまで定義してきた手続きを使用して一般項を求める手続きを書いてみます。
(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
結果が見づらいので張っておきます。フィボナッチ数列の 10000 番目の値は、
33644764876431783266621612005107543310302148460680063906564769974680081442166662368155595513633734025582065332680836159373734790483865268263040892463056431887354544369559827491606602099884183933864652731300088830269235673613135117579297437854413752130520504347701602264758318906527890855154366159582987279682987510631200575428783453215515103870818298969791613127856265033195487140214287532698187962046936097879900350962302291026368131493195275630227837628441540360584402572114334961180023091208287046088923962328835461505776583271252546093591128203925285393434620904245248929403901706233888991085841065183173360437470737908552631764325733993712871937587746897479926305837065742830161637408969178426378624212835258112820516370298089332099905707920064367426202389783111470054074998459250360633560933883831923386783056136435351892133279732908133732642652633989763922723407882928177953580570993691049175470808931841056146322338217465637321248226383092103297701648054726243842374862411453093812206564914032751086643394517512161526545361333111314042436854805106765843493523836959653428071768775328348234345557366719731392746273629108210679280784718035329131176778924659089938635459327894523777674406192240337638674004021330343297496902028328145933418826817683893072003634795623117103101291953169794607632737589253530772552375943788434504067715555779056450443016640119462580972216729758615026968443146952034614932291105970676243268515992834709891284706740862008587135016260312071903172086094081298321581077282076353186624611278245537208532365305775956430072517744315051539600905168603220349163222640885248852433158051534849622434848299380905070483482449327453732624567755879089187190803662058009594743150052402532709746995318770724376825907419939632265984147498193609285223945039707165443156421328157688908058783183404917434556270520223564846495196112460268313970975069382648706613264507665074611512677522748621598642530711298441182622661057163515069260029861704945425047491378115154139941550671256271197133252763631939606902895650288268608362241082050562430701794976171121233066073310059947366875
です。
おわりに
平方根型を実装して、フィボナッチ数列の一般項を使用して 10000 番目の値を求めました。
ポイントは、
・平方根部分は無理に計算しないので誤差が出ない。
・先頭から順に求めないので計算がとても速い。
ところです。
参考文献
・計算機プログラムの構造と解釈