Interval Arithmetic in Mathematica

Jerry B. Keiper
                                                   Wolfram Research, Inc.
                                                        keiper@wri.com

Introduction

Mathematica, like all computer-algebra systems, has several types of arithmetic.  In addition to exact arithmetic (i.e., integers and rational numbers) and machine-precision floating-point arithmetic, it has high-precision floating-point arithmetic and interval
arithmetic.

Range Arithmetic

The high-precision floating-point arithmetic has two parameters that can be set to control its behavior, but the default behavior is essentially a form of interval arithmetic in which the intervals are assumed to be "short" relative to the magnitude of the numbers represented.  This form of arithmetic is sometimes referred to as significance
arithmetic or range arithmetic.  The length of the interval that is implicitly represented by a high-precision floating-point number is
                                                                 10-a
where a is the "accuracy" of the number, i.e., the number of digits to the right of the decimal point.  Internally each high-precision floating-point number has associated with it a machine float that represents the logarithm of the length of the interval; the length is not explicitly calculated nor are the digits actually counted.

Interval Arithmetic

There are problems, however, where range arithmetic is not sufficient and you really do want interval arithmetic.  The interval arithmetic in Mathematica version 2.2 is genuine interval arithmetic, complete with outward rounding and multi-intervals.  Its deficiencies are that complex intervals have not yet been implemented and the outward rounding is done a posteriori rather than as directed rounding in hardware, which many machines do not support.

Range Arithmetic

The "accuracy" of a number is given by the function Accuracy[ ], but by default its result is rounded to the nearest integer.  This rounding can be disabled with the command

SetPrecision[Round, False];

The length of the implicit interval associated with a real number can now be found with the function

len[x_Real] := 10^(-Accuracy[x])

Let's find the lengths of the intervals associated with 30-digit representations of the numbers 3 and 4.

a = N[3, 30]; len[a]

b = N[4, 30]; len[b]

Now let's find the lengths of the intervals resulting from various arithmetic operations on these two numbers.

{len[a+b], len[a-b], len[a b], len[a/b]}

Note that from the point of view of rigorous interval arithmetic, the lengths of the intervals implicit in range arithmetic are a bit sloppy.  This is because the purpose of range arithmetic is merely to give a good estimate of the number of correct digits in the result.  To reduce memory requirements, the logarithm of the length of the interval is stored only as a float.  To increase the speed  of the length calculations, rational minimax approximations are used.  Because interval arithmetic (and hence also range arithmetic) tends to be rather pessimistic and because range arithmetic is not intended to be a substitute for interval arithmetic, sloppy length calculations are not problematic.  (Note that a much more serious problem with attempting to view range arithmetic as interval arithmetic is that the length calculations assume that the interval lengths are "short" relative to the magnitude of the numbers they represent.  For even moderate precision, say 10 digits or so, this works quite nicely, but for 

For many calculations significance arithmetic is all that is needed.  If we subtract two numbers  we can lose many digits.  The number of digits displayed in a number corresponds to the logarithm of the relative length of the interval.

N[Pi, 30] - 31415926535897932384626/10000000000000000000000

For very well-conditioned functions we can gain many digits.

Erf[N[20, 20]]

The following example is from S. M. Rump, "Algorithms for verified inclusions---Theory and practice" in Reliability  in Computing, R. E. Moore, ed., Academic Press, San Diego,
CA, 1988, pp 109-126.  The problem is to evaluate

                      c = 333.75 b6 + a2 (11 a2 b2 - b6 - 121 b4 - 2) + 5.5 b8 + a / (2 b)
                    
where a = 77617 and b = 33096.  With fixed-precision floating-point arithmetic one does not know how many digits of the result are correct.  In fact, unless rather high-precision arithmetic is used the answer will be completely wrong.  In fact in Rump's calculations, single precision, double precision, and even extended precision gave the same wrong result, even the sign was wrong.  We start by defining a function that evaluates c starting with n digits.

c[n_] := Block[{a = N[77617, n], b = N[33096, n]},
             N[333+3/4, n] b^6 + a^2 (11 a^2 b^2 - b^6
             - 121 b^4 - 2) + N[5 + 1/2, n] b^8 + a/(2 b)
             ]

{c[10], c[20], c[30], c[40], c[50]}

10 digits is less than $MachinePrecision so the first calculation was done using machine numbers and is completely wrong.

SetPrecision[%[[3]], 10]

The first three digits of the result using 30 digits were in fact correct, but none of the digits was known to be correct so none was displayed.  We can let Mathematica worry about  how much precision is required to get a specified precision in the result.

n = 20; While[Precision[nc = c[n]] < 10, n += 5]; nc

Interval Arithmetic

There are problems where range arithmetic is not sufficient and you really do want interval arithmetic.  The interval arithmetic in Mathematica Version 2.2 is genuine interval arithmetic, complete with outward rounding and multi-intervals.  Its only deficiency is that the outward rounding is done a posteriori rather than as directed rounding in hardware, which many machines do not support. Thus the intervals tend to grow slightly more rapidly than is absolutely necessary.  Also, complex intervals have not yet been implemented.

a = Interval[{2, 3}]

b = Interval[{-2, 1}]

c = a + b

a/b

a > b

d = IntervalUnion[a, b]

IntervalIntersection[c, d]

Sin[a]

Sin[b]

Interval[Pi]

d = Interval[2.3]

d[[1,2]] - d[[1,1]]

Interval Rootfinding

The following examples are not intended to be complete, robust algorithms.  They are merely intended to illustrate some of the possibilities for interval algorithms in Mathematica.

Bisection

In a bisection search a given interval is recursively bisected and each half is tested for roots. In interval arithmetic, testing for roots can be done by testing whether 0 is in the image of the interval.  IntervalMemberQ[ ] works nicely for this.

intervalbisection::rec = "MaxRecursion exceeded.";

split[f_, x_, int_Interval, eps_, n_] :=
    Block[{a = int[[1,1]], b = int[[1,2]], c},
        If[!IntervalMemberQ[f /. x -> int, 0], Return[{}]];
        If[b-a < eps, Return[int]];
        If[n == 0, Message[intervalbisection::rec];
            Return[int]];
        c = (a+b)/2;
        split[f, x, #, eps, n-1]& /@
            {Interval[{a,c}], Interval[{c,b}]}
    ];

Options[intervalbisection] = {MaxRecursion -> 7};

intervalbisection[f_, x_, intab_, eps_, opts___] :=
    Block[{int, n},
        n = MaxRecursion /. {opts} /.
            Options[intervalbisection];
        int = Interval /@ (List @@ intab);
        IntervalUnion @@
            Flatten[split[f, x, #, eps, n]& /@ int]
        ];

intervalbisection[Sin[x], x, Interval[{2., 20.}], .1]

intervalbisection[Sin[x], x, %, .1]

intervalbisection[Sin[x], x, Interval[{2., 20.}], .1,
     MaxRecursion -> 10]

intervalbisection[Sin[1/x], x, Interval[{-1., 1.}], .01]

Newton's Method

With Newton's method the idea is to pick a point in the given interval and evaluate the function at that point. Evaluating the derivative of the function on the interval allows us to eliminate parts of the original interval (assuming that the derivative is bounded).

intervalnewton::rec = "MaxRecursion exceeded.";

intnewt[f_, df_, x_, {a_, b_}, eps_, n_] :=
    Block[{xmid, int = Interval[{a, b}]},
        If[b-a < eps, Return[int]];
        If[n == 0, Message[intervalnewton::rec];
            Return[int]];
        xmid = Interval[(a+b)/2];
        int = IntervalIntersection[int,
        xmid - N[f /. x -> xmid]/N[df /. x -> int]];
        (intnewt[f, df, x, #, eps, n-1])& /@
            (List @@ int)
        ];

Options[intervalnewton] = {MaxRecursion -> 7};

intervalnewton[f_, x_, int_Interval, eps_, opts___] :=
    Block[{df, n},
        n = MaxRecursion/.{opts}/.Options[intervalnewton];
        df = D[f, x];
        IntervalUnion @@ Select[
            Flatten[(intnewt[f, df, x, #, eps, n])& /@
               (List @@ int)],
            IntervalMemberQ[N[f /. x -> #], 0]&]
        ];

intervalnewton[Sin[x], x, Interval[{2., 20.}], .1]

intervalnewton[Sin[x], x, %, 0.00001]

intervalnewton[Sin[x]^2, x, Interval[{2., 20.}], .1]

intervalnewton[Sin[1/x], x, Interval[{-1., 1.}], .01]

Multidimensional Bisection

In principle, multidimensional bisection is no more difficult than one-dimensional bisection.  In practice, however, combinatorial explosion becomes a problem.  Here is a simple implementation of the multidimensional bisection method.

bis[int_Interval, eps_] :=
    Block[{a = int[[1,1]], b = int[[1,2]], c},
        If[b-a < eps, Return[{int}]];
        c = (a+b)/2;
        {Interval[{a, c}], Interval[{c, b}]}
        ];

split[f_, x_, int_, eps_] :=
    Block[{intlist},
        If[Or @@ (!IntervalMemberQ[#, 0]& /@
            (f /. Apply[Rule, Transpose[{x, int}], 1])),
            Return[{}]];
        intlist = Flatten[Outer[List, Sequence @@
            (bis[#, eps]& /@ int)], Length[int]-1];
        If[Length[intlist] == 1, Return[int]];
        Sequence @@ (split[f,x,#,eps]& /@ intlist)
        ];

recombine[
   {d1___List, {a___Interval, b1_Interval, c___Interval},
    d2___List, {a___Interval, b2_Interval, c___Interval},
    d3___List}] :=
  recombine[{d1,{a,IntervalUnion[b1,b2],c},d2,d3}]

intervalbisection[f_, x_, intab_, eps_] :=
    Block[{int},
        int = (Interval /@ (List @@ #))& /@ intab;
        int = Flatten[Outer[List, Sequence @@ int],
              Length[int]-1];
        recombine[Select[split[f,x,#,eps]& /@ int,
              (Length[#]>0)&]][[1]]
        ];

intervalbisection[{x^2 + y^2 - 2, x^2 - y^2 - 1}, {x, y},
        {Interval[{-3., 3.}], Interval[{-3., 3.}]}, .001]