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

[オンラインSICP読書女子会 #12 (1.3.4)](http://ladiespp.connpass.com/event/32991/)

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
```

### 実装:

```scm
; [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.  

実装:

```scm
; [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)))
```

実行結果:


```scm
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
```