SICP 2.1.4 Extended Exercise: Interval Arithmetic - プログラミング再入門

ノート

追加演習？幅を持った値の算術。

Alyssa P. Hackerちゃんが誤差範囲（？）を持った数値の演算を支援するシステムを作成していると言う。
例えば並列に繋いだ二つの抵抗器の合成の抵抗値の計算。抵抗器は抵抗値に対して一定の誤差範囲が定義されている。

Alyssaのアイデアは最低値と最高値のペアに対して算術演算を定義するもの。
掛け算の説明の「reciprocal」は逆数か。0除算のチェックが無いのが気になるが。

Exercise 2.7

selectorを定義せよ。杓子定規に書くと：

(define (lower-bound i)
  (let ((a (car i))
        (b (cdr i)))
    (if (< a b) a b)))
(define (upper-bound i)
  (let ((a (car i))
        (b (cdr i)))
    (if (> a b) a b)))

minとmaxが使えるのでもう少しスマートに：

(define (lower-bound i)
  (min (car i) (cdr i)))
(define (upper-bound i)
  (max (car i) (cdr i)))

動作確認

> (lower-bound (make-interval 1.1 0.9))
0.9
> (upper-bound (make-interval 1.1 0.9))
1.1
>

Exercise 2.8

引かれる数のlower boundから引く数のupper boundが最小
引かれる数のupper boundから引く数のlower boundが最大

二つの範囲が重なっていない条件では考えやすいが、重なっていても小さい方から大きい方を引いてもこの関係は変わらない。

(define (sub-interval i j)
  (make-interval (- (lower-bound i) (upper-bound j))
                 (- (upper-bound i) (lower-bound j))))

動作確認

> (sub-interval (make-interval 1.2 1.4) (make-interval 1.5 1.7))
(-0.5 . -0.10000000000000009)
> (sub-interval (make-interval 1.2 1.4) (make-interval 1.3 1.5))
(-0.30000000000000004 . 0.09999999999999987)
> (sub-interval (make-interval 1.2 1.4) (make-interval 1.1 1.3))
(-0.10000000000000009 . 0.2999999999999998)
> (sub-interval (make-interval 1.2 1.4) (make-interval 0.9 1.1))
(0.09999999999999987 . 0.4999999999999999)
> (- 1.4 0.9)
0.4999999999999999
>

半端な小数になるのは仕方ない模様。

Exercise 2.9

widthをupper-boundとlower-boundの差の半分と定義する。

二つの値AとBに対して、Aの中央値をa、widthをwa、Bの中央値をb、widthをwbとした時、すなわち
A=(a-wa, a+wa)
B=(b-wb, b+wb)
とすると、
A+B={a+b+(wa+wb), a+b-(wa+wb)}
A-B={a-b-(wa+wb), a-b+(wa+wb)}
なのでA+Bのwidthはwa+wbである。

掛け算の場合、例えば
A=(0.9, 1.1), wa=0.1
B=(1.9 2.1), wb=0.1
でA×Bを考えると、A×B=(1.71, 2.31), width=0.3で上記の条件は成り立たない。

Exercise 2.10

気になっていた0除算のチェック。

lower	upper	判定
負	負	0を含まないのでOK
負	正	跨いでいるのでNG
正	負	存在しない
正	正	0を含まないのでOK

lower-boundが正、あるいはupper-boundが負なら0を跨いでいない。
エラーを報告してプログラムを停止する関数errorはR5RSには無いので、ここではR6RSで動作させる事にする。

#!r6rs
(import (rnrs))
; 中略
(define (div-interval x y)
  (define (spans-zero? x)
    (not (or (positive? (lower-bound x)) 
             (negative? (upper-bound x)))))
  (if (spans-zero? y)
      (error "Divisor spans zero" y)
      (mul-interval x 
                    (make-interval (/ 1.0 (upper-bound y))
                                   (/ 1.0 (lower-bound y))))))

動作確認。

> (div-interval (make-interval 0.9 1.1) (make-interval 0.4 0.6))
(1.5 . 2.75)
> (div-interval (make-interval 0.9 1.1) (make-interval -0.1 0.1))
. . Divisor spans zero: (-0.1 . 0.1)
> (div-interval (make-interval 0.9 1.1) (make-interval -0.6 -0.4))
(-2.75 . -1.5)
>

ちゃんとエラーを報告して止まった。

Exercise 2.11

intervalの値には、上限と下限が

両方とも正
正負に跨がる
両方とも負

の３通りがあり、二つのintervalでは9通りの組み合わせとなる。

正正×正正＝小×小〜大×大
正正×負正＝大×小〜大×大
正正×負負＝大×小〜小×大
負正×正正＝小×大〜大×大
負正×負正＝min(小×大、大×小）〜max(小×小、大×大)
負正×負負＝大×小〜小×小
負負×正正＝小×大〜大×小
負負×負正＝小×大〜小×小
負負×負負＝大×大〜小×小

この場合分けをすると、一つの場合を除いて掛け算が２回で済む。従来では掛け算を必ず４回行っていた。
プログラムが長くなるので、lower-bound、upper-boundをそれぞれlow、upに改名して、前問のspans-zero?を独立の関数として取り出して利用する。

(define (low i)
  (let ((a (car i))
        (b (cdr i)))
    (if (< a b) a b)))
(define (up i)
  (let ((a (car i))
        (b (cdr i)))
    (if (> a b) a b)))
(define (positive-interval? x)
  (and (positive? (low x)) (positive? (low x))))
(define (spans-zero? x)
  (not (or (positive? (low x)) 
           (negative? (up x)))))
(define (mul-interval a b)
  (cond ((positive-interval? a)
         (cond ((positive-interval? b)
                (make-interval (* (low a) (low b)) (* (up a) (up b))))
               ((spans-zero? b)
                (make-interval (* (up a) (low b)) (* (up a) (up b))))
               (else
                (make-interval (* (up a) (low b)) (* (low a) (up b))))))
        ((spans-zero? a)
         (cond ((positive-interval? b)
                (make-interval (* (low a) (up b)) (* (up a) (up b))))
               ((spans-zero? b)
                (make-interval (min (* (low a) (up b)) (* (up a) (low b))) (max (* (low a) (low b)) (* (up a) (up b)))))
               (else
                (make-interval (* (up a) (low b)) (* (low a) (low b))))))
        (else
         (cond ((positive-interval? b)
                (make-interval (* (low a) (up b)) (* (up a) (low b))))
               ((spans-zero? b)
                (make-interval (* (low a) (up b)) (* (low a) (low b))))
               (else
                (make-interval (* (up a) (up b)) (* (low a) (low b))))))))

動作確認

> (mul-interval (make-interval 1 3) (make-interval 2 4))
(2 . 12)
> (mul-interval (make-interval 1 3) (make-interval -1 1))
(-3 . 3)
> (mul-interval (make-interval 1 3) (make-interval -4 -2))
(-12 . -2)
> (mul-interval (make-interval -1 1) (make-interval 1 3))
(-3 . 3)
> (mul-interval (make-interval -1 1) (make-interval -1 1))
(-1 . 1)
> (mul-interval (make-interval -1 1) (make-interval -4 -2))
(-4 . 4)
> (mul-interval (make-interval -3 -1) (make-interval 2 4))
(-12 . -2)
> (mul-interval (make-interval -3 -1) (make-interval -1 1))
(-3 . 3)
> (mul-interval (make-interval -3 -1) (make-interval -4 -2))
(2 . 12)
>

元のmul-intervalをmul-interval-orgとして定義し直して検算。

> (mul-interval-org (make-interval 1 3) (make-interval 2 4))
(2 . 12)
> (mul-interval-org (make-interval 1 3) (make-interval -1 1))
(-3 . 3)
> (mul-interval-org (make-interval 1 3) (make-interval -4 -2))
(-12 . -2)
> (mul-interval-org (make-interval -1 1) (make-interval 1 3))
(-3 . 3)
> (mul-interval-org (make-interval -1 1) (make-interval -1 1))
(-1 . 1)
> (mul-interval-org (make-interval -1 1) (make-interval -4 -2))
(-4 . 4)
> (mul-interval-org (make-interval -3 -1) (make-interval 2 4))
(-12 . -2)
> (mul-interval-org (make-interval -3 -1) (make-interval -1 1))
(-3 . 3)
> (mul-interval-org (make-interval -3 -1) (make-interval -4 -2))
(2 . 12)
>

Exercise 2.12

最低値、最高値の組み合わせではなく中央値と誤差率（％）での表現に変更。
データ構造そのものは変更せず追加のconstructorとselectorで対応する。

(define (make-center-percent center tolerance)
  (let ((width (* center tolerance (/ 1 100))))
    (make-interval (- center width) (+ center width))))

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

(define (percent i)
  (let ((w (- (upper-bound i) (lower-bound i))))
    (* (/ (/ w 2) (center i)) 100)))

検算。前述の掛け算をやってみる。

> (make-center-percent 2 50)
(1 . 3)
> (center (make-center-percent 2 50))
2
> (percent (make-center-percent 2 50))
50
> (make-center-percent 3 (/ 100 3))
(2 . 4)
> (center (make-center-percent 3 (/ 100 3)))
3
> (percent (make-center-percent 3 (/ 100 3)))
33 1/3
> (mul-interval (make-center-percent 2 50) (make-center-percent 3 (/ 100 3)))
(2 . 12)
> (mul-interval (make-interval 1 3) (make-interval 2 4))
(2 . 12)
>

Exercise 2.13

中央値に対して幅が小さい正の値の場合、掛け算の式を簡略化する事が出来る。a、bのトレランスをt_{a、t_b賭したとき、

を展開すると、+、-によって４つのパターンに分けられる。}

a、b、に対してt_a、t_bが十分に小さいとき、abt_at_bは無視出来る程小さくなる。
また、t_a、t_bは正数を前提と出来るので、上記４つの式の第２項は±ab(t_a+t_b)の範囲に入ってしまう。

そうすると、
$=ab\pm ab(t_a + t_b)$
と簡略化出来、centerはab、widthはt_a+t_bとみなせる。

Exercise 2.14

計算してみる。

> (define a (make-center-percent 1000 1))
> (define b (make-center-percent 1500 2))
> (center (par1 a b))
600.5617438064144
> (center (par2 a b))
599.9855963126561
> (percent (par1 a b))
4.597193908142459
> (percent (par2 a b))
1.400072020165655
>

Lemの言う通りcenterとtolerance共に結構な差が出る。
intervalの計算では同じあたいで割っても正確に１とはならない。

> (center (div-interval a a))
1.000200020002
> (percent (div-interval a a))
1.9998000199979968
>

なので、数式上a/a=1と置き換えて計算するのと、a/aを実際に計算させるのでは当然結果が異なる。色々やってみろと書いてはあるが、例えばa+a-aでもpercentは変わってしまうし、a*a/aだとcenterもpercentも随分変わってしまう。どちらも代数的には元のaと同じになる筈。

> (define a (make-center-percent 1000 1))
> (define aa (sub-interval (add-interval a a) a))
> (center aa)
1000
> (percent aa)
3
> (define aaa (div-interval (mul-interval a a) a))
> (center aaa)
1000.4000400040004
> (percent aaa)
2.9992002399280224
>

Exercise 2.15

Exercise 2.14で見た通り、intervalの計算を行う毎に誤差は大きくなって行く。
一方oneの要に誤差を含んでいないintervalの場合計算で誤差は大きくならない。

> (define one (make-interval 1 1))
> (center (add-interval a one))
1001
> (percent (add-interval a one))
1000/1001
> (center (mul-interval a one))
1000
> (percent (mul-interval a one))
1
> (center (sub-interval a one))
999
> (percent (sub-interval a one))
1 1/999
> (center (div-interval a one))
1000.0
> (percent (div-interval a one))
1.0
>

per1の計算ではdiv-interval、mul-interval、add-intervalをそれぞれ１回ずつ、合計３回行っている。
per2の計算は誤差を含んだ演算はadd-intervalの１回だけ。後はoneとの演算なので誤差は大きくならない。
Evaの言っている事は正しい。

Exercise 2.16

例えば二つの値aとbの中央値と誤差がたまたま同じだったとして、演算の段階でa/aなのかa/bなのかが区別出来ない限り、代数的に等価な二つの計算式の結果を同じにする事は無理。その計算がa/aだったらきっかり１を返して、そうでない場合は１に誤差をつけて返す必要がある。