LoginSignup
0
0

More than 5 years have passed since last update.

@hio-harp

オンラインSICP読書女子会 #12 (1.3.4)

Last updated at Posted at 2016-06-08

オンラインSICP読書女子会 #12 (1.3.4)

1.3.4 (練習問題1.45~1.46)

ex-1.45.

問題文からまとめ:

  • $y = \sqrt{x}$ → $y^2 = x$ → $y = \frac{x}{y}$
    • 不動点 → 収束しない
    • 平均緩和 → 収束する
  • $y = \sqrt[3]{x}$ → $y^3 = x$ → $y = \frac{x}{y^2}$
    • 不動点 → ?
    • 平均緩和 → 収束する
  • $y = \sqrt[4]{x}$ → $y^4 = x$ → $y = \frac{x}{y^3}$
    • 不動点 → ?
    • 平均緩和 → 収束せず?
    • 2重平均緩和 → 収束する

確認結果:

  • $\sqrt[2]{5.0}$ → 1重平均緩和 → $2.23606797749997 ^ 2 = 5.000000000000843$
  • $\sqrt[3]{5.0}$ → 1重平均緩和 → $1.70997404839059 ^ 3 = 4.999983348151809$
  • $\sqrt[4]{5.0}$ → 2重平均緩和 → $1.50948899881336 ^ 4 = 5.191822173316605$
  • $\sqrt[5]{5.0}$ → 2重平均緩和 → $1.98711426317909 ^ 5 = 30.98233912369552$
  • $\sqrt[6]{5.0}$ → 7重平均緩和 → $1.03477395094990 ^ 6 = 1.227645353466118$
  • $\sqrt[7]{5.0}$ → 6重平均緩和 → $3.08298443886412 ^ 7 = 2647.276647682332$
  • $\sqrt[8]{5.0}$ → 2重平均緩和 → $1.06820209772746 ^ 8 = 1.695225248112375$
  • $\sqrt[9]{5.0}$ → 6重平均緩和 → $0.85594475635858 ^ 9 = 0.246610532823243$
  • $\sqrt[10]{5.0}$ → NG

2乗根から5乗根の４つはそれっぽさあるけれど, 5乗根は検算合わないからだめっぽい…
6乗根から先はどうにもだめっぽ…＞＜。
close-enough? の判定と検算時のn乗での誤差拡大具合か時々すごい値も….

最初, n重平均緩和を算出するコードで 0重平均緩和 (つまり平均緩和かけない) の時に $f$ の適用自体が抜けてて, それが巡って $\sqrt[3]{x}$ の1重平均緩和使用時も発散しててなやんだり. (修正済み)

sqrt の n重平均緩和の確認:

; (f (lambda (y) (/ x y)))
; x = 5.0, initial = 1.0.

; n重平均緩和 (筆算)
; 0重平均緩和 = 5.0 (guessed まま, 処理なし)
; 1重平均緩和 = 3.0 ; next = x / guessed = 5.0 / 5.0 = 1.0, avg(5.0, 1.0) = 3.0
; 2重平均緩和 = 2.3 ; next = 5.0 / 3.0 = 1.666, avg(3.0, 1.666) = 2.333
; 3重平均緩和 = 2.238 ; next = (5.0 / 2.333) = 2.143, avg(2.333, 2.148) = 2.238
; 4重平均緩和 = 2.236 ; next = (5.0 / 2.238) = 2.234, avg(2.238, 2.234) = 2.236
; 5重平均緩和 = 2.236 ; next = (5.0 / 2.236) = 2.236, avg(2.236, 2.236) = 2.236 (高い精度で変化はしている)

; n重平均緩和 (ループ)
n-0 5.0
n-1 3.0
n-2 2.3333333333333335
n-3 2.238095238095238
n-4 2.2360688956433634
n-5 2.236067977499978

; n重平均緩和 (手動展開)
c-0 5.0
c-1 3.0
c-2 2.3333333333333335
c-3 2.238095238095238
c-4 2.2360688956433634
c-5 2.236067977499978

$\sqrt{5}$ の試行結果:

; nth-root:2 x:5.0 n-avg:1
OK #5 = 2.236067977499978
#==> 2.236067977499978
#==> [power:2] 5.000000000000843

; nth-root:3 x:5.0 n-avg:1
OK #17 = 1.709974048390598
#==> 1.709974048390598
#==> [power:3] 4.999983348151809

; nth-root:4 x:5.0 n-avg:2
OK #700 = 1.5094889988133677
#==> 1.5094889988133677
#==> [power:4] 5.191822173316605

; nth-root:5 x:5.0 n-avg:2
OK #41 = 1.9871142631790968
#==> 1.9871142631790968
#==> [power:5] 30.982339123695528

; nth-root:6 x:5.0 n-avg:7
OK #1117 = 1.034773950949901
#==> 1.034773950949901
#==> [power:6] 1.2276453534661187

; nth-root:7 x:5.0 n-avg:6
OK #5498 = 3.0829844388641208
#==> 3.0829844388641208
#==> [power:7] 2647.2766476823326

; nth-root:8 x:5.0 n-avg:2
OK #6746 = 1.0682020977274655
#==> 1.0682020977274655
#==> [power:8] 1.6952252481123755

; nth-root:9 x:5.0 n-avg:6
OK #1010 = 0.85594475635858
#==> 0.85594475635858
#==> [power:9] 0.24661053282324305

; nth-root:10 x:5.0 n-avg:12
NG: NEED MORE #10000 = 16.022456665477986
#==> 16.022456665477986
#==> [power:10] 1.115041564677067e12

実装:

; [ex-1.45.scm]

; {{{ configurations.

; 不動点探索の完了判定差分.
(define tolerance 0.00001)
;(define tolerance 0.0001)


; 不動点探索の繰り返し上限.
;(define MAX-RECURSION 20)
(define MAX-RECURSION 10000)


; n乗根.
;(define MAX-NTH-ROOT 4)
(define MAX-NTH-ROOT 10)


; n重平均緩和.
(define MAX-AVG-DAMP 12)


;(define (print-guess i guess) (print "guess #" i " = " guess))
(define (print-guess i guess) #t)
; }}} configurations.


(define (ex-1.45)
    (define x 5.0)
    (let
        (
            (f (lambda (y) (/ x y)))
        )
        (print "c-0 " (calc-next-simple f 1.0))
        (print "n-0 " (n-fold-average-damp 0 f 1.0))
        (print "c-1 " (calc-next-avg-damping f 1.0))
        (print "n-1 " (n-fold-average-damp 1 f 1.0))
        (print "c-2 " (calc-next-2nd-avg-damping f 1.0))
        (print "n-2 " (n-fold-average-damp 2 f 1.0))
        (print "c-3 " (calc-next-3rd-avg-damping f 1.0))
        (print "n-3 " (n-fold-average-damp 3 f 1.0))
        (print "c-4 " (calc-next-4th-avg-damping f 1.0))
        (print "n-4 " (n-fold-average-damp 4 f 1.0))
        (print "c-5 " (calc-next-5th-avg-damping f 1.0))
        (print "n-5 " (n-fold-average-damp 5 f 1.0))
        #t)
    (newline)
    (print "(sqrt-1.45 5.0 0)")
    (print "#==> " (sqrt-1.45 5.0 0))
    (print "#==> " (power (sqrt-1.45 5.0 0) 2))
    (print "(sqrt-1.45 5.0 1)")
    (print "#==> " (sqrt-1.45 5.0 1))
    (print "#==> [power:" 2 "] " (power (sqrt-1.45 5.0 1) 2))

    (for-each
        2
        MAX-NTH-ROOT
        (lambda (i) (test-nth-root-1.45 x i)))
    ;(test-nth-root-1.45 x 2)
    ;(test-nth-root-1.45 x 3)
    ;(test-nth-root-1.45 x 4)
    ;(test-nth-root-1.45 x 5)
    ;(test-nth-root-1.45 x 6)
    #t)


; [first..last]
(define (for-each first last f)
    (define (loop i acc-fake)
        (if
            (<= i last)
            (loop (+ i 1) (f i))
            #t))
    (loop first #t))


(define (test-nth-root-1.45 x root-nth)
    (define (f y) (/ x (power y (- root-nth 1))))
    (define first-guess 1.0)
    (define (nth-root-n-avg-damp n-avg-damp)
        (print "; nth-root:" root-nth " x:" x " n-avg:" n-avg-damp)
        (let
            (
                (result
                    (trace-fixed-point-1.45
                        f
                        first-guess
                        (lambda (f a) (n-fold-average-damp n-avg-damp f a))
                    )
                )
            )
            (print "#==> " result)
            (print "#==> [power:" root-nth "] " (power result root-nth))
        )
    )
    (newline)
    (print "[nth-root:" root-nth " x:" x "]")
    (for-each
        0
        MAX-AVG-DAMP
        (lambda (i) (nth-root-n-avg-damp i)))
)


; n重平均緩和.
(define (n-fold-average-damp n f guess)
    (define (loop i x)
        (if
            (= i 0)
            x
            (loop (- i 1) (average x (f x)))))
    (if
        (= n 0)
        (f guess)
        (loop n guess)))


(define (calc-next-simple f guess)
    (f guess)) ; 0th avg-damping.

(define (calc-next-avg-damping f guess)
    (average (f guess) guess))

(define (calc-next-2nd-avg-damping f guess)
    (calc-next-avg-damping f (calc-next-avg-damping f guess)))

(define (calc-next-3rd-avg-damping f guess)
    (calc-next-avg-damping f (calc-next-2nd-avg-damping f guess)))

(define (calc-next-4th-avg-damping f guess)
    (calc-next-avg-damping f (calc-next-3rd-avg-damping f guess)))

(define (calc-next-5th-avg-damping f guess)
    (calc-next-avg-damping f (calc-next-4th-avg-damping f guess)))


; トレース出力付き不動点探索.
(define (trace-fixed-point-1.45 f first-guess calc-next)
    (define (close-enough? v1 v2)
        (<
            (abs (- v1 v2))
            tolerance))
    (define (try guess i prev-guess-1 prev-guess-2 prev-guess-3 prev-guess-4)
        (let ((next (calc-next f guess)))
            ;(print "guess #" i " = " next)
            (print-guess i next)
            (cond
                ((close-enough? guess next)
                    (print "OK #" i " = " next)
                    next)
                ((= i MAX-RECURSION)
                    (print "NG: NEED MORE #" i " = " next)
                    next)
                ((infinite? next)
                    (print "NG: INF #" i)
                    next)
                (
                    (or
                        (= next prev-guess-1)
                        (= next prev-guess-2)
                        ;(close-enough? next prev-guess-1)
                        ;(close-enough? next prev-guess-2)
                        ;(close-enough? next prev-guess-3)
                        ;(close-enough? next prev-guess-4)
                    )
                    (print "NG: LOOPED #" i " = " next)
                    next)
                (else (try next (+ i 1) guess prev-guess-1 prev-guess-2 prev-guess-3)))))
    ;(print "guess #" 0 " = " first-guess)
    (print-guess 0 first-guess)
    (try first-guess 1 first-guess first-guess first-guess first-guess))


(define (average a b)
    (/ (+ a b) 2))


(define (power x n)
    (define (loop i y)
        (if
            (= i 0)
            y
            (loop (- i 1) (* y x))))
    (loop n 1))


(define (sqrt-1.45 x n-avg-damp)
    (trace-fixed-point-1.45
        (lambda (y) (/ x y)) ; f.
        1.0                  ; first-guess.
        ;calc-next-simple     ; calc-next.
        ;calc-next-avg-damping ; calc-next.
        (lambda (f a) (n-fold-average-damp n-avg-damp f a))
    )
)


(define (cbrt-1.45-1avg x)
    (trace-fixed-point-1.45
        (lambda (y) (/ x (square y))) ; f.
        1.0                    ; first-guess.
        ;calc-next-avg-damping ; calc-next.
        (lambda (f a) (n-fold-average-damp 1 f a))
    )
)

ex-1.46.

実装:

; [ex-1.46.scm]


(define (ex-1.46)
    (print "(sqrt-1.46 2.0)")
    (print ";==> " (sqrt-1.46 2.0))
    (print "(sqrt-1.46 5.0)")
    (print ";==> " (sqrt-1.46 5.0))
    (newline)
    (print "(fixed-point-1.46 cos 1.0) ; 0.7390822985224023")
    (print ";==> " (fixed-point-1.46 cos 1.0))
    #t)


(define (iterative-improve good-enough? improve)
    (define (loop guess)
        (if
            (good-enough? guess)
            guess
            (loop (improve guess))))
    loop)

; sqrt (sec-1.1.7)
(define (sqrt-1.46 x)
    (define (average a b) (/ (+ a b) 2))
    (define (solver)
        (iterative-improve
            ; (good-enough? guess)
            (lambda (guess)
                (<
                    (abs (- (square guess) x))
                    0.001))
            ; (improve guess)
            (lambda (guess)
                (average guess (/ x guess)))))
    ((solver) 1.0))


; fixed-point (sec-1.3.3)
(define (fixed-point-1.46 f first-guess)
    (define tolerance 0.00001)
    ; 1つ前の guess がわからないので,
    ; 1つ次の値を生成してそれに対して評価するように.
    ; 1つ次の値の生成 (improve 相当) が反復のたびに
    ; 2回よばれるのがすこし残念・ω・`
    ; 得られた解から最終的な解にするのにさらに追加でもう一回とか・ω・`
    (define (solver)
        (iterative-improve
            ; (good-enough? guess)
            (lambda (guess)
                (<
                    (abs (- (f guess) guess))
                    tolerance))
            ; (improve guess)
            (lambda (guess)
                (f guess))))
    (f ((solver) first-guess)))

実行結果:

gosh> (ex-1.46)
(sqrt-1.46 2.0)
;==> 1.4142156862745097
(sqrt-1.46 5.0)
;==> 2.2360688956433634

(fixed-point-1.46 cos 1.0) ; 0.7390822985224023
;==> 0.7390822985224024
#t
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0