区間演算パッケージ

計算機プログラムの構造と解釈

「計算機プログラムの構造と解釈」（第2版）問題2.16（p54）

一般に等価な式が異なる答えになる理由を説明せよ。この欠点の無い区間算術演算パッケージを工夫することはできるか。あるいはこれは不可能な仕事か。（注意：この問題は非常に難しい。）

区間算術演算パッケージというのは[3.2, 5.6]みたいな幅（誤差）を持つ数どうしで計算をするためのプログラムで、教科書のこの節で例題として説明されています。

教科書の実装は

1. そのような数を上端と下端を持つ区間とし、上端と下端をまとめてデータとして持っておく。
2. 上端と下端にアクセスする方法を定義する。
3. 区間どうしの基本的な演算、すなわち足し算、引き算、掛け算、割り算を定義する。
4. 基本的な演算の組合せによってより複雑な演算も実現する

となっています。ここで基本的な演算の実装は「区間の端点どうしから計算（四則演算）して、その最大値および最小値から新しい区間を作る」となっています。しかしこのような実装だと、例えばXY/(X+Y)と1/(1/X + 1/Y)という等価な式なのに、計算結果が異なるいう問題が教科書中では指摘されます。なぜそのようなことが起きるのか？それを解決する手段はあるか？という問題です。

問題にわざわざ難しいと書いてある通り自分で考えてもよくわからなかったので検索してみたところ、http://oss.timedia.co.jp/show/SICP/ex-2.16というページを見つけました。

なるほど、理由としては式の上で同じ変数は同じ値を指しているはずなのに、今の実装では端点どうしの計算の最大値・最小値だから、そうした前提が崩れ問題になるということですね。基本的な演算を直接定義してしまい、それの組合せを使うという方法自体に問題がありそう。だいたい、その方針ではsinやcosとかも使えないですね。

この欠点の無い区間算術演算パッケージを工夫することはできるか？という問いに対しては、演算は高階関数で扱いなさい、という解。見て、あ、これは1章でやったことだぞ、と思い出しました。

高階関数で定義された式に区間の端点の値を可能な組合せ全部代入して得られた結果の最大・最小を演算結果の区間とするので、注意するべきは

1. （特別な目的が無ければ）結果をさらに別の演算に使わない
2. 演算する式を表現する関数において、与えられた区間において不連続な点がある場合は使えない

というところでしょうか。2番目の問題点は問題2.10

経験あるシステムプログラマ、Ben BitdiddleはAlyssaの肩越しに眺め、零を跨る区間で割った時、どうなるかよくわからないと評した。この状態が生じたことを調べ、起きたらエラーとするようにAlyssaのプログラムを修正せよ。

で指摘されています。0に十分近い所で割ると値が大きくなり、区間の端点どうしの計算が意味を持たなくなるということです。同じことが例えばタンジェントのπ/2などについても言える。どこに不連続点があるかはすぐに求められるわけではないので、これをどうにかしようと思ったら工夫が必要ということですね。

リンク貼ったページでは実装も示されていますが、私も実装してみました。

(define (make-interval x y)
  (cons (min x y) (max x y)))
(define (lower-bound i) (car i))
(define (upper-bound i) (cdr i))
(define (make-center-width c w)
  (make-interval (- c w) (+ c w)))
(define (center i)
  (/ (+ (lower-bound i) (upper-bound i)) 2))
(define (width i)
  (/ (- (upper-bound i) (lower-bound i)) 2))
(define (make-center-percent c p)
  (make-center-width c (* c (/ p 100))))
(define (percent i)
  (/ (width i) (center i)))

(define (test-values intervals)
  (if (null? intervals)
      (list '())
      (let ((rest (test-values (cdr intervals))))
	(append (map (lambda (x) (cons (lower-bound (car intervals)) x))
		     rest)
		(map (lambda (x) (cons (upper-bound (car intervals)) x))
		     rest)))))

(define (calc-interval formula . intervals)
  (let* ((test (test-values intervals))
	 (gain (map (lambda (v) (apply formula v)) test)))
      (make-interval (apply max gain) (apply min gain))))

(define (par1 i1 i2)
  (define (formula r1 r2)
    (/ (* r1 r2) (+ r1 r2)))
  (calc-interval formula i1 i2))

(define (par2 i1 i2)
  (define (formula r1 r2)
    (/ (+ (/ r1) (/ r2))))
  (calc-interval formula i1 i2))

(define (pr-interval i)
  (print (center i) " +/- " (percent i)))

(define i1 (make-center-percent 100 0.8))
(define i2 (make-center-percent 50 1.12))

(pr-interval (par1 i1 i2)) ;=> 33.33325747524561 +/- 1.0133377027691373
(pr-interval (par2 i1 i2)) ;=> 33.33325747524561 +/- 1.0133377027691373

多変数関数はlambdaをlambdaでくるむやり方よりも直裁的に表現した方がすぐれたインターフェイスだろう、という発想と、演算に使う区間の端点の組合せは冪集合を求める方法と同じ、という点を利用しています。それでapplyもうまく使えてるかなという感じ。applyの便利さが書いててわかります。