[SICP] 問題 1.07 : good-enough?の改良

平方根の計算で使ったgood-enough?テストは, 非常に小さい数の平方根をとる時には効果的ではない. また, 実際の計算機では, 算術演算は殆んどの場合, 限られた精度で実行される. それでわれわれのテストは非常に大きい数にも不適切である. 小さい数, 大きい数の場合, どのようにテストが失敗するかの例を使ってこのことを説明せよ. good-enough?を実装するもう一つの戦略は, ある繰返しから次へのguessの変化に注目し, 変化が予測値に比べ非常に小さくなった時に止めるのである. こういう終了テストを使う平方根手続きを設計せよ. これは小さい数, 大きい数に対してうまく働くか.

問題文にあるように, 小さい数と大きい数の平方根を求め, 組み込みの平方根を求める手続きの結果と比較します。

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (my-sqrt 10000.01)
100.00005025489767
> (sqrt 10000.01)
100.0000499999875
> (my-sqrt 0.0001)
0.03230844833048122
> (sqrt 0.0001)
0.01
> 

大きな数の場合は、組み込みの手続きと同等の結果が得られていますが, 小さな数の場合は正しい結果といえません.

ループの回数を調べるためにtraceを使います. 手続きsqrt-iterの呼び出され方を調べてみます. プログラムは次のようになります. １行目でtraceのためのライブラリを呼び出し, 最後の行でトレースする手続きを指定しています。

(require (lib "racket/trace.ss"))

(define (square x)
  (* x x))

(define (average x y)
  (/ (+ x y) 2))

(define (good-enough? guess x)
  (< (abs (- (square guess) x)) 0.001))

(define (improve guess x)
  (average guess (/ x guess)))

(define (sqrt-iter guess x)
  (if (good-enough? guess x)
      guess
      (sqrt-iter (improve guess x)
                 x)))

(define (my-sqrt x)
  (sqrt-iter 1.0 x))

(trace sqrt-iter)

平方根を求めると次のような結果が得られました.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (my-sqrt 10000.01)
>(sqrt-iter 1.0 10000.01)
>(sqrt-iter 5000.505 10000.01)
>(sqrt-iter 2501.252400010099 10000.01)
>(sqrt-iter 1252.6252005857104 10000.01)
>(sqrt-iter 630.3042212483191 10000.01)
>(sqrt-iter 323.08479587589 10000.01)
>(sqrt-iter 177.0182267724306 10000.01)
>(sqrt-iter 116.75482057222084 10000.01)
>(sqrt-iter 101.20223735102081 10000.01)
>(sqrt-iter 100.00719042723898 10000.01)
>(sqrt-iter 100.00005025489767 10000.01)
<100.00005025489767
100.00005025489767
> (my-sqrt 0.0001)
>(sqrt-iter 1.0 0.0001)
>(sqrt-iter 0.50005 0.0001)
>(sqrt-iter 0.2501249900009999 0.0001)
>(sqrt-iter 0.12526239505846617 0.0001)
>(sqrt-iter 0.06303035962394366 0.0001)
>(sqrt-iter 0.03230844833048122 0.0001)
<0.03230844833048122
0.03230844833048122
> (my-sqrt 0.000001)
>(sqrt-iter 1.0 1e-006)
>(sqrt-iter 0.5000005 1e-006)
>(sqrt-iter 0.250001249999 1e-006)
>(sqrt-iter 0.12500262498950004 1e-006)
>(sqrt-iter 0.06250531241075212 1e-006)
>(sqrt-iter 0.031260655525445276 1e-006)
<0.031260655525445276
0.031260655525445276
> 

実行結果から, 小さい数の場合は手続きの呼び出し回数に差がないことが分かります. 手続きgood-enough?は, 予測値guessの二乗と被開平数xの差が0.001より小さい場合に「十分よい」と判定しています. 被開平数が基準となる0.001より小さい場合, good-enough?の判定に与える影響が少なくなります. 例えば, 被開平数を0.000001, guessを0.03としgood-enough?を計算してみます.

  (< (abs (- (square guess) x)) 0.001)
= (< (abs (- (square 0.03) 0.000001)) 0.001)
= (< (abs (- (square 0.03) 0.000001)) 0.001)
= (< (abs (- 0.0009 0.000001)) 0.001)
= (< (abs 0.000899) 0.001)
= (< 0.000899 0.001)
= #t

この結果から, 被開平数が0.0001より小さいと予測値を2乗した結果だけで「十分よい」と判断し, 約0.03を求める平方根として返してしまいます.

問題文にあるように予測値の変化に着目し, 変化が少なくなった時に「十分に良い」と判定するように手続きを修正します. i番目の予測値の値をg_iとすると, 予測値の変化が次の式を満たせば予測値の変化が小さくなったと言えます.

| (gi+1 - gi) / gi | < 0.001

この式を使って手続きgood-enough?を改良します。 上の式でg_iはguessに相当し, g_i+1は手続きimproveを用いて計算出来ます. プログラムの全体は次のようになります.

(require (lib "racket/trace.ss"))

(define (square x)
  (* x x))

(define (average x y)
  (/ (+ x y) 2))

(define (improve guess x)
  (average guess (/ x guess)))

(define (good-enough? guess x)
  (< (abs (/ (- (improve guess x) guess)
             guess))
     0.001))

(define (sqrt-iter guess x)
  (if (good-enough? guess x)
      guess
      (sqrt-iter (improve guess x)
                 x)))

(define (my-sqrt x)
  (sqrt-iter 1.0 x))

(trace sqrt-iter)

実行結果は次のようになりました.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (my-sqrt 0.0001)
>(sqrt-iter 1.0 0.0001)
>(sqrt-iter 0.50005 0.0001)
>(sqrt-iter 0.2501249900009999 0.0001)
>(sqrt-iter 0.12526239505846617 0.0001)
>(sqrt-iter 0.06303035962394366 0.0001)
>(sqrt-iter 0.03230844833048122 0.0001)
>(sqrt-iter 0.017701806998329746 0.0001)
>(sqrt-iter 0.011675473894984788 0.0001)
>(sqrt-iter 0.010120218365353947 0.0001)
>(sqrt-iter 0.010000714038711746 0.0001)
<0.010000714038711746
0.010000714038711746
> (my-sqrt 0.000001)
>(sqrt-iter 1.0 1e-006)
>(sqrt-iter 0.5000005 1e-006)
>(sqrt-iter 0.250001249999 1e-006)
>(sqrt-iter 0.12500262498950004 1e-006)
>(sqrt-iter 0.06250531241075212 1e-006)
>(sqrt-iter 0.031260655525445276 1e-006)
>(sqrt-iter 0.015646322308953218 1e-006)
>(sqrt-iter 0.007855117545897352 1e-006)
>(sqrt-iter 0.003991211544161647 1e-006)
>(sqrt-iter 0.0021208810160681787 1e-006)
>(sqrt-iter 0.0012961915927068783 1e-006)
>(sqrt-iter 0.0010338412392442034 1e-006)
>(sqrt-iter 0.0010005538710539446 1e-006)
<0.0010005538710539446
0.0010005538710539446
> (my-sqrt 0.00000001)
>(sqrt-iter 1.0 1e-008)
>(sqrt-iter 0.500000005 1e-008)
>(sqrt-iter 0.25000001249999987 1e-008)
>(sqrt-iter 0.12500002624999892 1e-008)
>(sqrt-iter 0.06250005312499106 1e-008)
>(sqrt-iter 0.03125010656242753 1e-008)
>(sqrt-iter 0.015625213280668165 1e-008)
>(sqrt-iter 0.007812926635966154 1e-008)
>(sqrt-iter 0.003907103283034967 1e-008)
>(sqrt-iter 0.001954831361974762 1e-008)
>(sqrt-iter 0.0009799734463768973 1e-008)
>(sqrt-iter 0.0004950889022510404 1e-008)
>(sqrt-iter 0.0002576436474057565 1e-008)
>(sqrt-iter 0.00014822847335384218 1e-008)
>(sqrt-iter 0.00010784594750729777 1e-008)
>(sqrt-iter 0.00010028540197249 1e-008)
>(sqrt-iter 0.00010000040611237676 1e-008)
<0.00010000040611237676
0.00010000040611237676
> 

期待通りの結果が得られました.