ウォンツテック

そでやまのーと

SICPの問題を解こう 22

SICP

問題 1.45

n乗根を求めるのに y -> x/y^(n-1)のfixed-pointを計算していき差が少なくなった値を近似値とする方法において、y -> x/y^(n-1)にaverage dampingを何回行えばそれが収束するかという問題。
解析的に気になったのでSICP問題を解いているブログ界隈で探したけど記述があまりなかったので考察してみました。

解析

k乗根のfixed-pointを求める変換式 y -> x / y^(k-1) においてaverage damping
をm重回行った変換を考え、その変換式をn回行った後の値をy_nと置く。
この時y_1は最初に与える値first-guessである。
y_nが収束するとして収束値をyと置き、y_nにおけるyとの誤差を\delta_nと置く。

y_nはaverage dampingをm回行った変換式であるので以下の関係が成り立つ
\begin{eqnarray} y_{n+1}&=&\frac{y_n + \frac{y_n + ... \frac{y_n + \frac{x}{y_n^{k-1}}}{2}}{2}}{2}\\&=&\frac{(2^{m-1} + 2^{m-2} + ... + 1) y_n^k + x}{2^m y_n^{k-1}}\\&=&(1-\frac{1}{2^m})y_n + \frac{x}{2^m y_n^{k-1}}\end{eqnarray}

ここでy_nが収束するとするならば、nが十分大きい時\delta_nは収束値yに比べて
十分小さくなるので\delta_nの2次以上の項は無視をした以下の近似式が成り立つ。
\begin{eqnarray}\frac{x}{2^m y_n^{k-1}} &=& \frac{x}{2^m}(y + \delta_n)^{1-k}\\&\simeq& \frac{x}{2^m}(y^{1-k} + (1-k)y^{-k}\delta_n) \\ &=& \frac{y + (1-k)\delta_n}{2^m} \end{eqnarray}
※最後の等式でy^k = xを使用した。

これよりy_{n+1}は以下のように近似出来る
\begin{eqnarray} y_{n+1} &\simeq& (1-\frac1{2^m})y_n + \frac{y + (1-k)\delta_n}{2^m} \\&=& y_n - \frac{y_n}{2^m} + \frac{y + \delta_n}{2^m} - \frac{k}{2^m}\delta_n \\&=& y_n - \frac{k}{2^m}\delta_n \end{eqnarray}
※最後の等式でy_n = y + \delta_nを使用。

また、y_{n+1} = y + \delta_{n+1}より
\begin{eqnarray} y_{n+1} - y_n &=& y + \delta_{n+1} - (y + \delta_n) \\&=& \delta_{n+1} - \delta_n \\&\simeq& - \frac{k}{2^m}\delta_n \end{eqnarray}
※先ほどの近似式を利用

となるので\delta_nについての漸化式
 \delta_{n+1} = (1-\frac{k}{2^m})\delta_n
が得られる。y_nが収束するためには\delta_nが0に収束する必要があり、その条件は上記漸化式より
 -1 < 1-\frac{k}{2^m} < 1
となる。この右側の不等号は\frac{k}{2^m} > 0より常に成り立つので左側の不等号を考えると

\frac{k}{2^m} < 2
すなわち
2^{m+1} > k
m > log_2{k} - 1
という条件が得られる。

考察終わり。(mは整数なので m = \lfloor log_2{k} \rfloorが求める整数値ですね）

コード

(define (repeated f count)                                                      
  (if (= count 1)                                                               
      (lambda (x) (f x))                                                        
      (lambda (x) (f ((repeated f (- count 1)) x)))))                           
                                                                                
(define (fixed-point f guess)                                                   
  (let ((next (f guess)))                                                       
    (if (close-enough? guess next)                                              
        next                                                                    
        (fixed-point f next))))                                                 
(define (close-enough? a b)                                                     
  (define tolerance 0.0001)                                                     
  (print a " " b)                                                               
  (< (abs (- a b)) tolerance))                                                  
                                                                                
(define (average-damp f)                                                        
  (lambda (x) (/ (+ x (f x)) 2)))                                               
                                                                                
(define (expr x n)                                                              
  (if (= n 1)                                                                   
      x                                                                         
      (* x (expr x (- n 1)))))                                                  
                                                                                
(define (nth-root x n k)                                                        
  (fixed-point ((repeated average-damp k)                                       
                (lambda (y) (/ x (expr y (- n 1))))) 2.0))                      
                                                                                
(print (nth-root 10 2 1))                                                       
(print (nth-root 10 3 1))                                                       
;(print (nth-root 10 4 1)) ;don't converge                                      
(print (nth-root 10 4 2))                                                       
;(print (nth-root 10 8 2)) ;don't converge                                      
(print (nth-root 10 8 3))

問題 1.46

終了条件のチェックと次の改良値を算出する手続きの抽象化

(define (interative-improve f-enough f-improve)                                 
  (lambda (guess x) (let ((next (f-improve guess x)))                           
                        (if (f-enough guess next)                               
                             next                                               
                             ((interative-improve f-enough f-improve)           
                              next x)))))                                       
(define tolerance 0.0001)                                                       
(define (enough a b)                                                            
  (print a " " b)                                                               
  (< (abs (- a b)) tolerance))                                                  
(define (improve-sqrt guess x)                                                  
  (average guess (/ x guess)))                                                  
(define (average a b)                                                           
  (/ (+ a b) 2))                                                                
                                                                                
(define (fixed-point f first-guess x)                                           
  ((interative-improve enough f) first-guess x))                                
                                                                                
(define (sqrt first-guess x)                                                    
  ((interative-improve enough improve-sqrt) first-guess x))                     
                                                                                
(print (sqrt 1.0 2))                                                            
(print (fixed-point (lambda (y x) (/ (+ y (/ x y)) 2)) 1.0 2))

defineで引数を全て記述する場合はなんの迷いもなく書けるけど、lambdaをこんな感じで使おうとすると必ず失敗してちょっとづつ修正する事になる。。慣れるしか無さそう。
やっと1章終了。1ヶ月に数日思い出したようにやってただけだし、数学的に問題を解析してたから長いことかかった。。後4章終わらせるのに2年くらいかかるかも。