Scheme - データによる抽象の構築(区間得算術演算子(中央値とパーセント相対許容誤差、積の相対許容誤差の近似式))

計算機プログラムの構造と解釈
Gerald Jay Sussman, Julie Sussman
Harold Abelson (原著) 和田 英一 (翻訳)

開発環境

• OS X Lion - Apple(OS)
• Emacs、BBEdit - Bare Bones Software, Inc. (Text Editor)
• プログラミング言語: MIT/GNU Scheme

計算機プログラムの構造と解釈(Gerald Jay Sussman(原著)、Julie Sussman(原著)、Harold Abelson(原著)、和田 英一(翻訳)、ピアソンエデュケーション)の2(データによる抽象の構築)、2.1(データ抽象入門)、2.1.4(拡張問題: 区間算術演算)の問題 2.12、問題 2.13を解いてみる。

その他参考書籍

• Instructor's Manual to Accompany Structure & Interpretation of Computer Programs

問題 2.12

コード

sample.scm

(define (make-center-percent c p)
  (let ((w (* c (/ p 100))))
    (make-interval (- c w) (+ c w))))

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

(define (percent i)
  (let ((c (center i))
        (w (width i)))
    (* (/ w c) 100)))

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

(define (make-interval a b) (cons a b))

(define (upper-bound x) (cdr x))

(define (lower-bound x) (car x))

(define (print-interval i)
  (newline)
  (display "(")
  (display (lower-bound i))
  (display ", ")
  (display (upper-bound i))
  (display ")"))

(define a (make-center-percent 100 5))

入出力結果(Terminal, REPL(Read, Eval, Print, Loop))

1 ]=> (print-interval a)

(95, 105)
;Unspecified return value

1 ]=> (center a)

;Value: 100

1 ]=> (percent a)

;Value: 5

問題 2.13

区間i1、i2のそれぞれの中央値をc1、c2、パーセント相対許容誤差をp1、p2とする。(問題を単純化するために、区間i1、i2の下限はともに正とする。)

e1 = p1/100
e2 = p2/100
とおく。p1とp2を十分小さいと仮定し、e1 * e2、e1 * e1、e2 * e2を0として近似する。
区間i1とi2の積(i3)の
下限
(c1 - c1 * e1) * (c2 - c2 * e2)
= c1 * c2 - c1 * c2 * e2 - c1 * c2 * e1
= c1*c2*(1 - e1 - e2)
上限
(c1 + c1 * e1) * (c2 + c2 * e2)
= c1 * c2 + c1 * c2 * e2 + c1 * c2 * e1
= c1*c2*(1 + e1 + e2)
区間i3の中央値
(c1*c2*(1 + e1 + e2) + c1*c2*(1 - e1 - e2)) / 2
= c1*c2
区間i3のパーセント相対許容誤差をp3とし、e3 = p3/100とおくとe3は、
(c1*c2*(1 + e1 + e2) - c1*c2) / (c1*c2)
= e1 + e2
よってパーセント相対許容誤差が小さいと仮定できるときは、二つの区間の積の相対許容誤差は、因子の許容誤差を使って
p1 + p2
という簡単な式で求めることができる。

確認。

コード

sample.scm

(define (mul-interval x y)
  (let ((x1 (lower-bound x))
        (y1 (upper-bound x))
        (x2 (lower-bound y))
        (y2 (upper-bound y))
        (negative? (lambda (n) (< n 0))))
    (let ((x1-sign (negative? x1))
          (y1-sign (negative? y1))
          (x2-sign (negative? x2))
          (y2-sign (negative? y2)))
      (cond ((and x1-sign (not y1-sign) x2-sign (not y2-sign))
             (let ((p1 (* x1 y2))
                   (p2 (* x2 y1))
                   (p3 (* x1 y1))
                   (p4 (* y2 y2)))
                   (make-interval (min p1 p2) (max p3 p4))))
            ((and (not x1-sign) (not x2-sign))
             (make-interval (* x1 x2) (* y1 y2)))
            ((and (not x1-sign) (not y2-sign))
             (make-interval (* y1 x2) (* y1 y2)))
            ((not x1-sign) (make-interval (* y1 x2) (* x1 y2)))
            ((and x1-sign (not y1-sign) (not x2-sign))
             (make-interval (* x1 y2) (* y1 y2)))
            ((and x1-sign (not y1-sign))
             (make-interval (* y1 x2) (* x1 x2)))
            (y2-sign (make-interval (* y1 y2) (* x1 x2)))
            ((not x2-sign) (make-interval (* x1 y2) (* y1 x2)))
            (else (make-interval (* x1 y2) (* x1 x2)))))))

(define (make-center-percent c p)
  (let ((w (* c (/ p 100))))
    (make-interval (- c w) (+ c w))))

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

(define (percent i)
  (let ((c (center i))
        (w (width i)))
    (* (/ w c) 100)))

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

(define (make-interval a b) (cons a b))

(define (upper-bound x) (cdr x))

(define (lower-bound x) (car x))

(define p1 0.01)
(define p2 0.1)

(define i1 (make-center-percent 10 p1))
(define i2 (make-center-percent 100 p2))

(define i (mul-interval i1 i2))
(define p (percent i))
(define p-test (+ p1 p2))

入出力結果(Terminal, REPL(Read, Eval, Print, Loop))

1 ]=> p

;Value: .109999988999992

1 ]=> p-test

;Value: .11