#lang racket

(define (add-interval x y)
  (make-interval (+ (lower-bound x) (lower-bound y))
                 (+ (upper-bound x) (upper-bound y))))

(define (mul-interval x y)
  (define p1 (* (lower-bound x) (lower-bound y)))
  (define p2 (* (lower-bound x) (upper-bound y)))
  (define p3 (* (upper-bound x) (lower-bound y)))
  (define p4 (* (upper-bound x) (upper-bound y)))
  (make-interval (min p1 p2 p3 p4)
                 (max p1 p2 p3 p4)))

(define (div-interval x y)
  (mul-interval x
                (make-interval (/ 1.0 (upper-bound y))
                               (/ 1.0 (lower-bound y)))))

;;;;;;;;;
;; 2.7 ;;
;;;;;;;;;

(define (make-interval a b)
  (if (<= a b)
    (cons a b)
    (error "lower bound must not be greater than upper bound" a b)))

(define (lower-bound i)
  (car i))

(define (upper-bound i)
  (cdr i))

;;;;;;;;;
;; 2.8 ;;
;;;;;;;;;

(define (sub-interval x y)
  (make-interval (- (lower-bound x) (upper-bound y))
                 (- (upper-bound x) (lower-bound y))))

;;;;;;;;;
;; 2.9 ;;
;;;;;;;;;

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

(define a (make-interval 0.2 0.4))
(define b (make-interval 1.1 1.5))

(width a)
;; 0.1
(width b)
;; 0.19999999999999996
(width (add-interval a b))
;; 0.29999999999999993
(width (sub-interval a b))
;; 0.3
(width (mul-interval a b))
;; 0.19000000000000003
(width (div-interval a b))
;; 0.11515151515151516

;; In general, the width of the sum or difference of intervals is the sum of the
;; widths.  However the widths of the product or quotient of intervals depends on
;; the upper and lower bounds of the intervals and not just on their widths.

(define b2 (make-interval 0.1 0.5)) ; same width as b
(width (mul-interval a b2))
;; 0.09
(width (div-interval a b2))
;; 1.8

;;;;;;;;;;
;; 2.10 ;;
;;;;;;;;;;

(define (contains-zero? i)
  (<= (* (upper-bound i)
         (lower-bound i))
      0))

(define (new-div-interval x y)
  (if (contains-zero? y)
    (error "divisor interval cannot contain zero"
           (lower-bound y)
           (upper-bound y))
    (mul-interval x
                  (make-interval (/ 1.0 (upper-bound y))
                                 (/ 1.0 (lower-bound y))))))

;;;;;;;;;;
;; 2.11 ;;
;;;;;;;;;;

;; The nine cases are:
;; x positive and y positive,
;; x positive and y negative,
;; x negative and y positive,
;; x negative and y negative,
;; x contains zero and y positive,
;; x contains zero and y negative,
;; y contains zero and x positive,
;; y contains zero and x negative,
;; x and y both contain zero; only this last case needs four multiplications.

(define (new-mul-interval x y)
  (define xl (lower-bound x))
  (define xu (upper-bound x))
  (define yl (lower-bound y))
  (define yu (upper-bound y))
  (cond [(and (>= xl 0) (>= yl 0)) (make-interval (* xl yl) (* xu yu))]
        [(and (>= xl 0) (<= yu 0)) (make-interval (* xu yl) (* xl yu))]
        [(and (<= xu 0) (>= yl 0)) (make-interval (* xl yu) (* xu yl))]
        [(and (<= xu 0) (<= yu 0)) (make-interval (* xu yu) (* xl yl))]
        [(>= yl 0) (make-interval (* xl yu) (* xu yu))]
        [(<= yu 0) (make-interval (* xu yl) (* xl yl))]
        [(>= xl 0) (make-interval (* xu yl) (* xu yu))]
        [(<= xu 0) (make-interval (* xl yu) (* xl yl))]
        [else
          (define p1 (* xl yl))
          (define p2 (* xl yu))
          (define p3 (* xu yl))
          (define p4 (* xu yu))
          (make-interval (min p1 p2 p3 p4)
                         (max p1 p2 p3 p4))]))

(define x1 (make-interval 1 2))
(define x2 (make-interval -2 -1))
(define x3 (make-interval -1 2))

(define y1 (make-interval 3 5))
(define y2 (make-interval -5 -3))
(define y3 (make-interval -3 5))

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

(define xs (list x1 x1 x1 x2 x2 x2 x3 x3 x3))
(define ys (list y1 y2 y3 y1 y2 y3 y1 y2 y3))

(for/and ([x xs]
          [y ys])
         (equal-interval? (mul-interval x y)
                          (new-mul-interval x y)))
;; #t

;;;;;;;;;;
;; 2.12 ;;
;;;;;;;;;;

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

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

(define (make-center-percent c p)
  (make-center-width c (* c p)))

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

;;;;;;;;;;
;; 2.13 ;;
;;;;;;;;;;

;; x = interval with center c1 and percent p1
;; y = interval with center c2 and percent p2

;; x has lower-bound c1 - c1 * p1 and upper-bound c1 + c1 * p1 
;; y has lower-bound c2 - c2 * p2 and upper-bound c2 + c2 * p2 

;; assuming all numbers are positive, the product
;; has lower-bound c1 * (1 - p1) * c2 * (1 - p2) ~= c1 * c2 * (1 - p1 - p2)
;; and upper-bound c1 * (1 + p1) * c2 * (1 + p2) ~= c1 * c2 * (1 + p1 + p2)
;; where we drop the insignificant second-order terms involving p1 * p2

;; so the product has center c1 * c2 and width c1 * c2 * (p1 + p2)

;; so the tolerance of the product is p1 + p2

;;;;;;;;;;
;; 2.14 ;;
;;;;;;;;;;

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

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

(define r1 (make-center-percent 4 0.1))
(define r2 (make-center-percent 10 0.2))

(define circ1 (par1 r1 r2))
(center circ1)
;; 3.1539108494533226
(percent circ1)
;; 0.4432000000000001

(define circ2 (par2 r1 r2))
(center circ2)
;; 2.8511354079058036
(percent circ2)
;; 0.12920353982300886

;;;;;;;;;;
;; 2.15 ;;
;;;;;;;;;;

;; I think so.  It stands to reason that more uncertain inputs should increase the
;; uncertainty of the output.

;;;;;;;;;;
;; 2.16 ;;
;;;;;;;;;;

;; Simplifying an algebraic expression before doing the interval arithmetic will
;; usually decrease the uncertainty in the final answer, like in 2.15.

;; One way around this issue would be to choose a normal form for each equivalence
;; class of algebraic expressions and always use that form to compute the answer.
;; I'm not sure if that's possible though.