Interval Arithmetic

December 21, 2010

We represent an interval as a pair and state the operations as given:

(define (plus x y)
  (cons (+ (car x) (car y)) (+ (cdr x) (cdr y))))

(define (minus x y)
  (cons (- (car x) (cdr y)) (- (cdr x) (car y))))

(define (times x y)
  (cons (min (* (car x) (car y)) (* (car x) (cdr y))
             (* (cdr x) (car y)) (* (cdr x) (cdr y)))
        (max (* (car x) (car y)) (* (car x) (cdr y))
             (* (cdr x) (car y)) (* (cdr x) (cdr y)))))

(define (divide x y)
  (if (< (car y) 0 (cdr y))
      (error 'divide "divide by zero")
      (cons (min (/ (car x) (car y)) (/ (car x) (cdr y))
                 (/ (cdr x) (car y)) (/ (cdr x) (cdr y)))
            (max (/ (car x) (car y)) (/ (car x) (cdr y))
                 (/ (cdr x) (car y)) (/ (cdr x) (cdr y))))))

To represent the interval as a pair containing its central point and a tolerance (width), we provide two conversions back-and-forth between the two representations:

(define (ends->center x)
  (cons (/ (+ (car x) (cdr x)) 2)
        (/ (- (cdr x) (car x)) 2)))

(define (center->ends x)
  (cons (- (car x) (cdr x)) (+ (car x) (cdr x))))

Here are some examples, which are reproduced at http://programmingpraxis.codepad.org/JZQwpzdb.

> (define x (cons 1 2))
> (define y (cons 3 4))
> (plus x y)
(4 . 6)
> (minus x y)
(-3 . -1)
> (times x y)
(3 . 8)
> (divide x y)
(1/4 . 2/3)
> (ends->center x)
(3/2 1/2)
> (center->ends (cons 3/2 1/2))
(1 . 2)
> (divide x x)
(1/2 . 2)

For the other three exercises, note that the problem is that normal algebra doesn’t necessarily work on intervals. You can see the basic problem by dividing an interval by itself: if x is the interval [1, 2], x/x is not 1, as we are accustomed to in algebra, but the interval [1/2, 2] that includes 1 as it’s geometric mean. Eva is correct in pointing out that par2 is better than par1 because it calculates each interval once, instead of applying one interval to another. The last exercise is the subject of several books and academic papers, which we won’t even try to summarize here; if you’re interested, a good place to start is Vladik Kreinovich’s web site.

About these ads

Pages: 1 2

10 Responses to “Interval Arithmetic”

  1. [...] Praxis – Interval Arithmetic By Remco Niemeijer In today’s Programming Praxis exercise, our goal is to implement some functions to do interval arithmetic. [...]

  2. My Haskell solution (see http://bonsaicode.wordpress.com/2010/12/21/programming-praxis-interval-arithmetic/ for a version with comments):

    plus :: (Num a, Num b) => (a, b) -> (a, b) -> (a, b)
    plus (a,b) (c,d) = (a+c, b+d)
    
    minus :: (Num a, Num b) => (a, b) -> (b, a) -> (a, b)
    minus (a,b) (c,d) = (a-d, b-c)
    
    times :: (Num a, Ord a) => (a, a) -> (a, a) -> (a, a)
    times (a,b) (c,d) = (min a b * min c d, max a b * max c d)
    
    divide :: (Fractional a, Ord a) => (a, a) -> (a, a) -> (a, a)
    divide (a,b) (c,d) = if c < 0 && d > 0 then error "divide by 0"
        else (min a b / max c d, max a b / min c d)
    
    toCenter :: Fractional a => (a, a) -> (a, a)
    toCenter (a,b) = ((a+b) / 2, (b-a) / 2)
    
    fromCenter :: Num a => (a, a) -> (a, a)
    fromCenter (a,b) = (a-b, a+b)
    
  3. Lautaro Pecile said
    from collection import namedtuples
    
    class Interval (namedtuple('Interval', 'a b')):
        __slots__ = ()
        
        def __add__ (self, other):
            return Interval(self.a + other.a, self.b + other.b)
    
        def __sub__ (self, other):
            return Interval(self.a - other.a, self.b - other.b)
        
        def __mul__ (self, other):
            l = [x*y for x in self for y in other]
            return Interval(min(l), max(l))
        
        def __div__ (self, other):
            l = [x/float(y) for x in self for y in other]
            return Interval(min(l), max(l))
        
        @property
        def toCenter(self):
            a = (self.a + self.b) / 2.0
            b = (self.b - self.a) / 2.0
            return Interval(a, b)
        
        @property
        def fromCenter(self):
            a = (self.a - self.b)
            b = (self.a + self.b)
            return Interval(a, b)
    
  4. Graham said

    My Python submission is available on codepad.org.
    It requires the fractions module (new in Python 2.6), which adds support for rational number arithmetic.
    My work is slightly hackish and could use better variable names, but it’ll do for now.

  5. Interval =
      Struct.new(:a, :b) do
        def +(other)
          self.class.new(a + other.a, b + other.b)
        end
    
        def -(other)
          self.class.new(a - other.a, b - other.b)
        end
    
        def *(other)
          products =
            [
              [ a, other.a ],
              [ a, other.b ],
              [ b, other.a ],
              [ b, other.b ],
            ].map { |(c, d)| c * d }
          self.class.new(products.min, products.max)
        end
    
        def /(other)
          quotients =
            [
              [ a, other.a ],
              [ a, other.b ],
              [ b, other.a ],
              [ b, other.b ],
            ].map { |(c, d)| Rational(c, d) }
          self.class.new(quotients.min, quotients.max)
        rescue ZeroDivisionError
          raise ZeroDivisionError, "divided by #{other}", caller
        end
      end
    
    # Examples:
    a = Interval.new(5, 7)
    b = Interval.new(0, 3)
    
    a + b # => #<struct Interval a=5, b=10>
    a - b # => #<struct Interval a=5, b=4>
    a * b # => #<struct Interval a=0, b=21>
    b / a # => #<struct Interval a=(0/1), b=(3/5)>
    
    ## ZeroDivisionError
    a / b rescue $! # => #<ZeroDivisionError: divided by #<struct Interval a=0, b=3>>
    
  6. Axio said

    I guess something like this is good enough:

    (defclass interval-class ()
      ((min :initarg :min)
       (max :initarg :max)))

    (defmethod print-object ((iv interval-class) stream)
      (format stream “[~a,~a]” (slot-value iv ‘min)
              (slot-value iv ‘max)))

    (defvar it1 (make-instance ‘interval-class :min 20 :max 30))
    (defvar it2 (make-instance ‘interval-class :min 10 :max 20))
    (defvar it3 (make-instance ‘interval-class :min -10 :max 20))

    (defmacro it-op (name min max &optional guard)
      `(defmethod ,name ((i1 interval-class)
                         (i2 interval-class))
         (macrolet ((a () `(slot-value i1 ‘min))
                    (b () `(slot-value i1 ‘max))
                    (c () `(slot-value i2 ‘min))
                    (d () `(slot-value i2 ‘max)))
           (unless ,guard
             (make-instance ‘interval-class :min ,min :max ,max)))))

    (it-op it-add (+ (a) (c)) (+ (b) (d)))
    (it-op it-sub (- (a) (c)) (- (b) (d)))
    (it-op it-mul
           (min (* (a) (c)) (* (a) (d))
                (* (b) (c)) (* (b) (d)))
           (max (* (a) (c)) (* (a) (d))
                (* (b) (c)) (* (b) (d))))
    (it-op it-div
           (min (/ (a) (c)) (/ (a) (d))
                (/ (b) (c)) (/ (b) (d)))
           (max (/ (a) (c)) (/ (a) (d))
                (/ (b) (c)) (* (b) (d)))
           (or (and (<= (a) 0) (<= 0 (b)))
               (and (<= (c) 0) (<= 0 (d)))))

    Next time, I’ll implement a reader to add intervals to the syntax :)

  7. Axio said

    With less bugs :)
    Please disregard the above untested code :)

    (defclass interval-class ()
      ((min :reader i-min :initarg :min)
       (max :reader i-max :initarg :max)))

    (set-macro-character #\[
                         (lambda (stream char)
                           (let* ((lst (read-delimited-list #\] stream t)))
                             (if (= 2 (length lst))
                                 (make-instance ‘interval-class :min (apply #’min lst) :max (apply #’max lst))
                                 (error “Invalid interval”)))))
    (set-macro-character #\] (get-macro-character #\)))

    (defmethod print-object ((iv interval-class) stream)
      (format stream “[~a,~a]” (i-min iv) (i-max iv)))

    (defmethod plus ((i1 interval-class) (i2 interval-class))
      (let ((a (+ (i-min i1) (i-min i2)))
            (b (+ (i-max i1) (i-max i2))))
        (make-instance ‘interval-class :min (min a b) :max (max a b))))

    (defmethod minus ((i1 interval-class) (i2 interval-class))
      (let ((a (- (i-min i1) (i-max i2)))
            (b (- (i-max i1) (i-min i2))))
        (make-instance ‘interval-class :min (min a b) :max (max a b))))

    (defmacro m-d (opname op &optional guard)
      `(defmethod ,opname ((i1 interval-class) (i2 interval-class))
         (unless ,guard
           (let* ((b0 (loop for i in (list #’i-min #’i-max)
                         append (loop for j in (list #’i-min #’i-max)
                                   collect (,op (funcall i i1) (funcall j i2)))))
                  (b1 (apply #’min b0))
                  (b2 (apply #’max b0)))
             (make-instance ‘interval-class :min b1 :max b2)))))

    (m-d mul *)
    (m-d div / (and (<= (i-min i2) 0) (<= 0 (i-max i2))))

    (defvar x [1 2])
    (defvar y [4 3]) ;; reverse order \o/
    (list (plus x y) (minus x y) (mul x y) (div x y))

  8. slabounty said

    A ruby version using the Rational class. The to_s is has a bunch of extra code to make it look “pretty” if the denominator is 1. It might be better to monkey patch this in the Rational class itself.

    class Interval
    
        attr_reader :a, :b
        def initialize(a, b)
            @a = Rational(a/1)
            @b = Rational(b/1)
        end
    
        def add(i)
            return Interval.new(@a + i.a, @b + i.b)
        end
    
        def subtract(i)
            return Interval.new(@a - i.b, @b - i.a)
        end
    
        def multiply(i)
            return Interval.new(
                [@a*i.a, @a*i.b, @b*i.a, @b*i.b].min,
                [@a*i.a, @a*i.b, @b*i.a, @b*i.b].max
            )
        end
    
        def divide(i)
            begin
                return Interval.new(
                    [@a/i.a, @a/i.b, @b/i.a, @b/i.b].min,
                    [@a/i.a, @a/i.b, @b/i.a, @b/i.b].max
                )
            rescue ZeroDivisionError
                return "Divide by zero"
            end
        end
    
        def ends_center
            return Interval.new(((@a + @b) / 2.0), ((@b - @a) / 2.0))
        end
    
        def center_ends
            return Interval.new(@a - @b, @a + @b)
        end
    
        def to_s
            "#{a.numerator}" + ((a.denominator !=  1) ? "/#{a.denominator}" : "" ) + ", " +
            "#{b.numerator}" + ((b.denominator !=  1) ? "/#{b.denominator}" : "" )
        end
    end
    
    x = Interval.new(1, 2)
    y = Interval.new(3, 4)
    puts "x = #{x}"
    puts "y = #{y}"
    puts "interval x+y = #{x.add(y)}"
    puts "interval x+y = #{x.subtract(y)}"
    puts "interval x*y = #{x.multiply(y)}"
    puts "interval x/y = #{x.divide(y)}"
    puts "ends_center x = #{x.ends_center}"
    puts "center_ends x = #{Interval.new(1.5, 0.5).center_ends}"
    puts "interval x/x = #{x.divide(x)}"
    
  9. Khanh Nguyen said
    open System
    
    let plus (x:float*float) (y:float*float) =
        (fst x + fst y, snd x + snd y)
    
    let minus (x:float*float) (y:float*float) =
        (fst x - fst y, snd x - snd y)
    
    let times (x:float*float) (y:float*float) = 
        match x, y with
        | (a,b), (c,d) -> (List.min [a*c; a*d; b*c; b*d], List.max [a*c; a*d; b*c; b*d])
    
    let div (x:float*float) (y:float*float) = 
        match x, y with
        | (a,b), (c,d) -> (List.min [a/c; a/d; b/c; b/d], List.max [a/c; a/d; b/c; b/d])
    
  10. As Chun Kin Lee pointed out on my blog, my original attempt at removing duplication from times and divide doesn’t work on ranges with negative numbers, so you actually have to get the minimum/maximum of the four possible combinations.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 609 other followers

%d bloggers like this: