SICP 読書ノート#6 - 2.1.4 区間算術演算 (pp.52-54)

「§2.1.4 拡張問題：区間算術演算」から。

なんというか、この章はひたすら数学。。。

問題 2.7

(define (make-interval a b) (cons a b))
(define (upper-bound i) (max (car i) (cdr i)))
(define (lower-bound i) (min (car i) (cdr i)))

問題 2.8

減算手続きの実装。

すでに定義している add-interval、mul-interval から導き出せばかんたん。

(define (sub-interval x y)
  (add-interval x (mul-interval (make-interval -1 -1) y)))

テスト。

racket@> (sub-interval (make-interval 4 8) (make-interval 1 7))
'(-3 . 7)

問題 2.9

任意の区間iの中央値を $c(i)$ 、幅を $w(i)$ と定義すると、区間xと区間yの加算は、

${ \displaystyle x + y = (c(x) + w(x)) + (c(y) + w(y)) = c(x) + c(y) + (w(x) + w(y)) }$

となるため、結果の幅は $w(x) + w(y)$ となり、幅だけの関数となる。減算も同様。

一方、乗算では、

${ \displaystyle x \times y = (c(x) + w(x)) \times (c(y) + w(y)) \\ = c(x) \times c(y) + (c(x) \times w(y) + c(y) \times w(x) + w(x) \times w(y)) }$

となるため、結果の幅は $c(x) \times w(y) + c(y) \times w(x) + w(x) \times w(y)$ となり、幅だけの関数とはならない。除算も同様。

問題 2.10

正の区間か、負の区間か、0をまたぐ区間かをチェックする手続き。

(define (positive-interval? i)
  (and (> (upper-bound i) 0) (> (lower-bound i) 0)))

(define (negative-interval? i)
  (and (< (upper-bound i) 0) (< (lower-bound i) 0)))

(define (zero-interval? i)
  (and (not (positive-interval? i)) (not (negative-interval? i))))

除算手続き。除数が0をまたぐ区間の場合は例外を発生させる。

(define (div-interval x y)
  (if (zero-interval? y)
      (error "divied by zero interval")
      (mul-interval x 
                    (make-interval (/ 1.0 (upper-bound y))
                                   (/ 1.0 (lower-bound y))))))

問題 2.11

2つの引数がそれぞれ正の区間、負の区間、0をまたぐ区間で9パターンに分けられる。

そのパターンのうち、双方とも0をまたぐ区間はlower,upperの乗算の組み合わせをチェックする必要がある。

(define (mul-interval x y)
  (cond ((positive-interval? x)
         (cond ((positive-interval? y)
                (make-interval (* (lower-bound x) (lower-bound y))
                               (* (upper-bound x) (upper-bound y))))
               ((negative-interval? y)
                (make-interval (* (upper-bound x) (lower-bound y))
                               (* (lower-bound x) (upper-bound y))))
               ((zero-interval? y)
                (make-interval (* (upper-bound x) (lower-bound y))
                               (* (upper-bound x) (upper-bound y))))))
        ((negative-interval? x)
         (cond ((positive-interval? y)
                (make-interval (* (lower-bound x) (upper-bound y))
                               (* (upper-bound x) (lower-bound y))))
               ((negative-interval? y)
                (make-interval (* (upper-bound x) (upper-bound y))
                               (* (lower-bound x) (lower-bound y))))
               ((zero-interval? y)
                (make-interval (* (lower-bound x) (upper-bound y))
                               (* (lower-bound x) (lower-bound y))))))
        ((zero-interval? x)
         (cond ((positive-interval? y)
                (make-interval (* (lower-bound x) (upper-bound y))
                               (* (upper-bound x) (upper-bound y))))
               ((negative-interval? y)
                (make-interval (* (upper-bound x) (lower-bound y))
                               (* (lower-bound x) (lower-bound y))))
               ((zero-interval? y)
                (make-interval (min (* (lower-bound x) (upper-bound y))
                                    (* (upper-bound x) (lower-bound y)))
                               (max (* (lower-bound x) (lower-bound y))
                                    (* (upper-bound x) (upper-bound y)))))))
        ))

元々の mul-interval を mul-interval-orig に名前を変えてテスト。

(define (equal-interval x y)
  (and (= (lower-bound x) (lower-bound y))
       (= (upper-bound x) (upper-bound y))))

(let ((z (make-interval -3 3))  ;; 0をまたぐ区間
      (p (make-interval 2 4))   ;; 正の区間
      (n (make-interval -1 6))) ;; 負の区間
  (display "test01 => ")
  (display (equal-interval (mul-interval-orig p p) (mul-interval p p)))
  (newline)
  (display "test02 => ")
  (display (equal-interval (mul-interval-orig p n) (mul-interval p n)))
  (newline)
  (display "test03 => ")
  (display (equal-interval (mul-interval-orig p z) (mul-interval p z)))
  (newline)
  (display "test04 => ")
  (display (equal-interval (mul-interval-orig n p) (mul-interval n p)))
  (newline)
  (display "test05 => ")
  (display (equal-interval (mul-interval-orig n n) (mul-interval n n)))
  (newline)
  (display "test06 => ")
  (display (equal-interval (mul-interval-orig n z) (mul-interval n z)))
  (newline)
  (display "test07 => ")
  (display (equal-interval (mul-interval-orig z p) (mul-interval z p)))
  (newline)
  (display "test08 => ")
  (display (equal-interval (mul-interval-orig z n) (mul-interval z n)))
  (newline)
  (display "test09 => ")
  (display (equal-interval (mul-interval-orig z z) (mul-interval z z)))
  (newline))

問題 2.12

最小・最大ではなく中央値と幅から区間を生成する手続きがあったとして、

(define (center i)
  (/ (+ (lower-bound i) (upper-bound i)) 2))

(define (width i)
  (/ (- (upper-bound i) (lower-bound i)) 2))

(define (make-center-width c w)
  (make-interval (- c w) (+ c w)))

中央値とパーセント誤差から区間を生成する手続きと、区間からパーセント誤差を取り出す手続きを生成する。

(define (make-center-percent c p)
  (let ((w (* c (/ p 100.0))))
    (make-center-width c w)))

(define (percent i)
  (/ (* (width i) 100.0) (center i)))

問題 2.13

問題 2.12で解いた任意の区間iのパーセント誤差を取り出す手続きを数式で表現すると

${ \displaystyle p(i) = \frac{w(i)}{c(i)} \times 100 }$

となる。

よって、区間xと区間yの積のパーセント誤差は

${ \displaystyle p(x \times y) = \frac{w(x \times y)}{c(x \times y)} \times 100 }$

となる。

一方、区間xと区間yの積は、問題2.9より、

${ \displaystyle x \times y = c(x) \times c(y) + c(y) \times w(x) + c(x) \times w(y) + w(x) \times w(y) \\ }$

$w(x) \times w(y)$ が中央値に対して無視できるほど小さいとすると、

${ \displaystyle x \times y \fallingdotseq c(x) \times c(y) + c(y) \times w(x) + c(x) \times w(y) }$

となるから、中央値からの誤差をTとすると、

${ \displaystyle T = \frac{c(y) \times w(x) + c(x) \times w(y)}{c(x) \times c(y)} = \frac{w(x)}{c(x)} + \frac{w(y)}{c(y)} }$

となるため、区間xと区間yの積の近似パーセント誤差の手続き(※)は

${ \displaystyle p'(x \times y) = \left( \frac{w(x)}{c(x)} + \frac{w(y)}{c(y)} \right) \times 100 }$

と表現できる。

(※)の手続きを実装すると、

(define (percent-mul-approx x y)
  (* (+ (/ (width x) (center x)) (/ (width y) (center y))) 100.0))

実行すると、

racket@> (define i1 (make-center-width 1000 1))
racket@> (define i2 (make-center-width 2000 1))

racket@> (percent (mul-interval i1 i2))
0.1499999250000375
racket@> (percent-mul-approx i1 i2)
0.15

細かい誤差については捨てられていることがわかる。

問題 2.14、2.15

par1, par2 について、

(define (par1 r1 r2)
  (div-interval (mul-interval r1 r2)
                (add-interval r1 r2)))

(define (par2 r1 r2)
  (let ((one (make-interval 1 1))) 
    (div-interval one
                  (add-interval (div-interval one r1)
                                (div-interval one r2)))))

適当な区間を代入してみる。par2は正しい模様。

racket@> (define r1 (make-center-width 3 1))
racket@> (define r2 (make-center-width 5 1))
racket@> r1
'(2 . 4)
racket@> r2
'(4 . 6)

racket@> (par1 r1 r2)
'(0.8 . 4.0)

racket@> (par2 r1 r2)
'(1.3333333333333333 . 2.4000000000000004)

これは、加算、減算および区間 $[1, 1$ ] と任意の区間の乗算、除算では誤差が蓄積されないが、不確かな数同士の乗算、除算は計算する度に誤差が蓄積されていくからだと思う。

問題 2.16

全くわからなかったので、電気系の仕事をしている友人に聞いてみたら「最大・最小の組み合わせをチェックしてみたら」といったヒントをもらったので、 div-intervalをそのように書き直してみた。

(define (div-interval2 x y)
  (let ((l1 (/ (lower-bound x) (lower-bound y)))
        (l2 (/ (lower-bound x) (upper-bound y)))
        (u1 (/ (upper-bound x) (lower-bound y)))
        (u2 (/ (upper-bound x) (upper-bound y))))
    (make-interval (max l1 l2) (min u1 u2))))

par1、par2をdiv-interval2に置き換えた手続きをpar3、par4とする

(define (par3 r1 r2)
  (div-interval2 (mul-interval r1 r2)
                 (add-interval r1 r2)))

(define (par4 r1 r2)
  (let ((one (make-interval 1 1))) 
    (div-interval2 one
                   (add-interval (div-interval2 one r1)
                                 (div-interval2 one r2)))))

評価してみる。

racket@> (par3 r1 r2)
'(4/3 . 12/5)
racket@> (par4 r1 r2)
'(12/5 . 4/3)

惜しい。div-interval2を少し手直ししてみる。

(define (div-interval2 x y)
  (let ((l1 (/ (lower-bound x) (lower-bound y)))
        (l2 (/ (lower-bound x) (upper-bound y)))
        (u1 (/ (upper-bound x) (lower-bound y)))
        (u2 (/ (upper-bound x) (upper-bound y))))
    (make-interval (min (max l1 l2) (min u1 u2))
                   (max (max l1 l2) (min u1 u2)))))

評価すると、

racket@> (par3 r1 r2)
'(4/3 . 12/5)
racket@> (par4 r1 r2)
'(4/3 . 12/5)

できたっぽいが、どうなんだろう？

正直よくわかっていないけど、いったんこれでパスで。

次は「§2.2 階層データ構造と閉包性」から。

@uents blog

Code wins arguments.