計算機プログラムの構造と解釈 第二版 P53 問題2.14

Lemの正しさを示すには実際やってみることです。

実装

#!/usr/local/bin/gosh
;; -*- coding: utf-8 -*-

(use ggc.debug.trace)
(use math.mt-random)

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

(define upper-bound cdr)
(define lower-bound car)

;;和
(define (add-interval x y)
  (make-interval (+ (lower-bound x) (lower-bound y))
                 (+ (upper-bound x) (upper-bound y))))
;;差
(define (sub-interval x y)
  (make-interval (- (lower-bound x) (upper-bound y))
                 (- (upper-bound x) (lower-bound y))))

;;積
(define (mul-interval x y)
  (let ((xl (lower-bound x))
        (xu (upper-bound x))
        (yl (lower-bound y))
        (yu (upper-bound y)))
    (cond ((>= xl 0)
           (cond ((>= yl 0) (make-interval (* xl yl) (* xu yu)))
                 ((<= yu 0) (make-interval (* xu yl) (* xl yu)))
                 (else
                   (make-interval (* xu yl) (* xu yu)))))
          ((<= xu 0)
           (cond ((>= yl 0) (make-interval (* xl yu) (* xu yl)))
                 ((<= yu 0) (make-interval (* xu yu) (* xl yl)))
                 (else
                   (make-interval (* xl yu) (* xl yl)))))
          (else
            (cond ((>= yl 0) (make-interval (* xl yu) (* xu yu)))
                  ((<= yu 0) (make-interval (* xu yl) (* xl yl)))
                  (else
                    (make-interval (min (* xl yu) (* xu yl)) (max (* xl yl) (* xu yu)))))))))

;;商
(define (straddle? x)
  (if (>= 0 (* (upper-bound x) (lower-bound x)))
    #t
  #f))

(define (div-interval x y)
 (if (straddle? y)
    (error "straddle 0" y))
 (mul-interval x
               (make-interval (/ 1.0 (upper-bound y))
                              (/ 1.0 (lower-bound y)))))

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


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

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


;;2-12Answer
(define (make-center-percent c p)
  (let ((w (/ (* c (/ p 100)) 2)))
    (make-center-width c w)))

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



;;2-13
(define (part1 r1 r2)
  (div-interval (mul-interval r1 r2)
                (add-interval r1 r2)))

(define (part2 r1 r2)
  (let ((one (make-interval 1 1)))
    (div-interval one
                  (add-interval (div-interval one r1)
                                (div-interval one r2)))))
                                                                            
(define interval1 (make-interval 2 6))
(define interval2 (make-interval 5 7))

(define A (make-interval 99.5 100.5))
(define B (make-center-percent 5 10))




;; main
(define (main args)

  (newline)
  (display "interval1: ")
  (display interval1)(newline)
  (display "interval2: ")
  (display interval2)(newline)
  (display "A: ")
  (display A)(newline)
  (display "B: ")
  (display B)(newline)
  (newline)


  (display "Lem曰く、二つの計算法で異なる結果になる。")(newline)
  (display "(part1 interval1 interval2): ")
  (display (part1 interval1 interval2))
  (newline)
  (display "(part2 interval1 interval2): ")
  (display (part2 interval1 interval2))
  (newline)
  (newline)


  (display "教科書に書いてあるとおりにやってみる")(newline)
  (display "(div-interval A A): ")
  (display (div-interval A A))
  (newline)

  (display "(div-interval A B): ")
  (display (div-interval A B))
  (newline)

  (display "(percent (div-interval A A)): ")
  (display (percent (div-interval A A)))
  (newline)

  (display "(percent (div-interval A B)): ")
  (display (percent (div-interval A B)))

  (newline)
  (display "許容誤差が足された値になっております。")(newline)
  (newline)

  (display "和とか差とか、積をやってみるとどうなる?")(newline)
  (display "(percent (mul-interval A A)): ")
  (display (percent (mul-interval A A)))
  (newline)

  (display "(percent (mul-interval A B)): ")
  (display (percent (mul-interval A B)))
  (newline)

  (display "(percent (add-interval A A)): ")
  (display (percent (add-interval A A)))
  (newline)

  (display "(percent (add-interval A B)): ")
  (display (percent (add-interval A B)))
  (newline)

  (display "(percent (sub-interval A A)): ")
  (display (percent (sub-interval A A)))
  (newline)

  (display "(percent (sub-interval A B)): ")
  (display (percent (sub-interval A B)))
  (newline)

  (display "積算も許容誤差が足された値になっております。")(newline)

  (newline)
0)

結果

interval1: (2 . 6)
interval2: (5 . 7)
A: (99.5 . 100.5)
B: (19/4 . 21/4)

Lem曰く、二つの計算法で異なる結果になる。
(part1 interval1 interval2): (0.7692307692307693 . 6.0)
(part2 interval1 interval2): (1.4285714285714286 . 3.230769230769231)

教科書に書いてあるとおりにやってみる
(div-interval A A): (0.990049751243781 . 1.0100502512562815)
(div-interval A B): (18.952380952380953 . 21.157894736842103)
(percent (div-interval A A)): 1.9999500012499796
(percent (div-interval A B)): 10.997250687328155
許容誤差が足された値になっております。

和とか差とか、積をやってみるとどうなる?
(percent (mul-interval A A)): 1.999950001249969
(percent (mul-interval A B)): 10.997250687328169
(percent (add-interval A A)): 1.0
(percent (add-interval A B)): 1.4285714285714286
(percent (sub-interval A A)): +inf.0
(percent (sub-interval A B)): 1.5789473684210527
積算も許容誤差が足された値になっております。


足し算だと、許容誤差は計算しにくくなる。
普通に考えて、大きい方の値に引きずられる。