X3J3/96156
Date: October 12, 1996
To: X3J3
From: R. Baker Kearfott
Subject: Interval Arithmetic  The Data Type and LowLevel Operations
INTERVAL ARITHMETIC
THE DATA TYPE AND LOWLEVEL OPERATIONS
The following is a synthesis of extensive email discussions
between Keith Bierman, George Corliss, David Hough, Andrew
Pitonyak, Michael Schulte, William Walster, Wolfgang Walter, and
me. This represents substantial refinement over corresponding
portions of paper 1996066. I am now presenting it for X3J3's
consideration.
R. Baker Kearfott
================================================================
Name and Structure
____ ___ _________
The INTERVAL type is a numeric type. Its values are closed and
bounded real intervals which are defined by an ordered pair of
real values. The first real value is the lower bound (or
infimum), the second real value is the upper bound (or supremum).
The lower bound shall be less than or equal to the upper bound.
The continuum of real numbers between and including these two
bounds (or endpoints) is said to be contained in the interval.
Note: Thus intervals contain real numbers that are not otherwise
machine representable.
The actual internal structure of an interval data type is
implementationdependent, and the bit structure is not directly
accessible to the user. That is, the user may change the values
of the endpoints of an interval only through setting
an entire interval, through the IVAL constructor described
below, or through an implicit conversion to interval.
The values of the lower bound and upper bound of an interval can
be obtained with the functions INF and SUP defined below.
The precision of each interval data type shall correspond to the
precision of a real data type. The kind type parameter of each
supported interval type shall be the same as the kind type
parameter of the corresponding real type.
Note: Although the interval data type is opaque, a common model
of intervals is that of two reals. For example, for the interval
data type corresponding to the DOUBLE PRECISION IEEE binary type,
the lower bound can be viewed as an IEEE DOUBLE PRECISION type
and the upper bound can be viewed as an IEEE DOUBLE PRECISION
type.
Note: It is recommended that the default INTERVAL kind correspond
to a REAL kind with at least 64 bits. ("DOUBLE PRECISION" on
many machines). It is the consensus of experts that 32bit
interval arithmetic is of limited use.
Interval Constants
________ _________
Both where literal constants are admitted in expressions and in
input or output data, INTERVAL's shall be represented by a single
REAL or INTEGER or a pair of REAL's, INTEGER's, or combinations
thereof, beginning with "(<", separated by "," if there are two
numbers, and ending in ">)". For example
(<1, 2>), (<1E0, 3>), (<1>), and (<.1234D5>)
are all valid INTERVAL constants. An INTERVAL constant specified
by a single number is the same as an INTERVAL constant specified
by two numbers, both of whose endpoints are equal to the single
number. When such a decimal constant is converted to its internal
representation, the internal representation shall contain the
decimal constant, regardless of how many digits are specified by
the decimal constant. For example, upon execution of the
statement:
X = (<0.31415926535897932384626433832795028D+01>)
the interval X shall contain the smallestwidth machine interval
that contains the number 3.1415926535897932384626433832795028.
Interval constants shall admit kind type parameters, such as in
the construction (<1,2>)_2. An interval constant with no kind
type parameter shall be of default type.
Note: Thus, on machines in which the interval data type X
appearing in the example above contained endpoints with accuracy
that corresponded to less than 35 decimal digits, the interval X
would contain the mathematical number PI.
Note: The statement
X = 0.31415926535897932384626433832795028D+01,
is allowed, with mixedmode operations as defined below. The
result of this statement is an interval in X that contains
3.1415926535897932384626433832795028, just as as
with the statement
X = (<0.31415926535897932384626433832795028D+01>)
Arithmetic Operations
__________ __________
The four basic operations +, , *, and / are defined to contain
the ranges of the corresponding operations on real numbers.
Specifically, let X = [xl,xu] and Y = [yl,yu] be intervals, where
xl, xu, yl, and yu represent the lower and upper bounds of X and Y,
respectively. Then:
X + Y shall contain the exact value [xl + yl, xu + yu],
X  Y shall contain the exact value [xl  yu, xu  yl],
X * Y shall contain the exact value
[min{xl*yl, xl*yu, xu*yl, xu*yu}, max{xl*yl, xl*yu, xu*yl, xu*yu}],
1 / X shall contain the exact value [1/xu, 1/xl] if xu < 0 or xl > 0,
and shall signal 'denominator contains zero' otherwise.
X / Y shall contain the exact value
[min{xl/yu, xl/yl, xu/yu, xu/yl}, max{xl/yu, xl/yl, xu/yu, xu/yl}]
if yu < 0 or yl > 0, and shall signal 'denominator contains zero'
otherwise.
Note: Particular processors may support an extended interval
data type, in which division by intervals that contain zero
yields a meaningful result. For example, a meaningful
interpretation of [1,2]/[1,2] is the settheoretic union of the
two intervals (infinity,1] and [1/2,infinity).
Note: Using floating point arithmetic, the operations on the
righthand sides may first be computed, then the lower bound may
be rounded down to a number known to be less than or equal to the
exact mathematical result, and the upper bound may be rounded up
to a number known to be greater than or equal to the exact
mathematical result. If a processor provides directed roundings
upwards (towards plus infinity) and downwards (towards minus
infinity), then the operation and the rounding can be performed
in one step, e.g. if the processor conforms to the IEEE 754
standard. The excess interval width caused by this outward
rounding is called ROUNDOUT ERROR.
Note: There is an alternate implementation of interval
multiplication that also gives the range of the real operator "*"
over the intervals X and Y. This alternative involves nine cases
determined by the algebraic signs of the endpoints of X and Y;
see page 12 of R. E. Moore, "Methods and Applications of Interval
Computations," SIAM, Philadelphia, 1979. The average number of
multiplications required for this alternative is less than above,
but one or more comparisons are required.
Implemented in software, the relative efficiencies of the
alternative above and the ninecase alternative are
architecturedependent, although the ninecase alternative is
often preferred in lowlevel implementations designed for
efficiency.
Note: The only processor requirement is that the computed
intervals contain the exact mathematical range of the
corresponding point operations. In an ideal implementation (not
required), the result of the operations is the smallestwidth
machine interval that contains the exact mathematical range.
Note: IEEE arithmetic can be used to perform ideal (minimum width)
interval operations. For example, take
[xl,xu] + [yl,yu] = [xl+yl,xu+yu]
in exact interval arithmetic. The IEEE 754 standard defines a
downwardly (upwardly) rounded operation as producing the same
result as would be obtained by computing the exact result, then
rounding it to the nearest floating point number less (greater)
than or equal to the exact result. Thus, if the result xl+yl is
rounded down and xu+yu is rounded up according to the IEEE
specifications, the result is an ideal interval addition.
Mixed Mode Operations and Conversions
_____ ____ __________ ___ ___________
Explicit conversion to interval is performed with the
constructor IVAL as follows:
Z = IVAL(R,S[,kind]) Z < [R,S] or Z = IVAL(R[,kind]) Z < [R,R],
where Z is an interval data type, and R and S are real or
integer data types. The optional argument "kind" is the kind type
parameter of the target interval. If "kind" is omitted, the
result of IVAL shall be of default interval kind. If both
R and S are present, then R shall be less than or equal to S.
The interval value stored in Z shall be an enclosure for the
specified interval, with an ideal enclosure equal to a machine
interval of minimum width that contains the exact mathematical
interval in the specification. When the argument(s) of IVAL
are literal constants, the specification of IVAL shall be the same
as if "(<" and ")>" were used.
Example: IVAL(1.1,2.2) and (<1.1,2.2>) yield the same INTERVAL's.
Note: The specification of IVAL for literal constants may be best
done with special handling by the compiler. This is because,
if IVAL were an ordinary Fortran function, the literal
constants may first be converted to their internal
representations. For example, for IVAL(1.1,2.2),
1.1 and 2.2 could first be converted to internal
representations, and IVAL(1.1,2.2) would then not
necessarily contain the actual decimal interval
[1.1,2.2].
Note: An IVAL with correct specifications can be designed
without special compiler handling. In particular,
it can round EACH lower bound down and EACH upper bound
up, if conversion accuracies are known.
However, such an IVAL would then not be optimally
accurate, since IVAL(1,1) or IVAL(R), R a floating point
variable, would be intervals of nonzero width. The width
of the returned interval would necessarily correspond to
the minimum precision of any real data type supported by the
processor. Note that, although the specification of IVAL
is the same for literal constant arguments, the actual
value returned may be different from the internal
representation generated by a literal constant.
Mixedmode INTERVAL/REAL and mixed mode INTERVAL/INTEGER
operations shall be defined. The implicit conversions shall
be such that, at each conversion to interval, the resulting
interval shall contain the exact mathematical value specified
by the floating point or integer variable. Conversion of literal
constants to interval shall be as if the literal constant were
converted directly to interval, with no intervening conversion
to an internal floating point representation.
Examples:
1. If X is an interval variable, then the expression
2*X**2 + 3.1D0*X + 0.1
shall be evaluated in such a way that the interval result is the
same as if the literal constants 2, 3.1D0 and 0.1 were first
replaced by interval variables of the same type as X, such that
each interval variable contains the exact mathematical
values represented by these constants. This does not preclude
2, 3.1D0, and 0.1 from first being converted to internal
REAL representations; however, the ultimate
conversion to interval would then need to encompass the
conversion errors of the original conversion to internal
REAL representation.
2. If X is an interval variable and A, B, and C are real or integer
variables, then the expression
A*X**2 + B*X + C
shall be evaluated in such a way that the interval result is the
same as if the real variables A, B, and C were first converted
to interval variables of the same type as X, such that each
interval contains the exact mathematical result represented by
these real or integer variables.
3. If X is an interval and A is either a real or integer variable,
then execution of
X = A
shall cause an interval to be stored in X that contains the exact
value represented in A. (Note: Most implementations will give
xl = xu = A.)
4. If X is an interval, then execution of
X = 2.7
shall cause an interval to be stored in X that contains the exact
value 2.7.
5. If X is an interval, then execution of
X = 2.7 + 1.1 + 3.3 + 4.9 + R
will be an interval that may not necessarily contain the
exact value 12+R. This is because the 2.7, 1.1, 3.3,4.9, and R
may first be added as single precision floating point numbers,
then converted. This is analogous to 0 being stored in the
real variable R, with the statement.
R = 1/7
Users should be made aware of this possibility, and should
be encouraged to avoid it.
6. Similarly, if R is a real and X is an interval, then the
statements
R = 0.1
X = R
may result in the degenerate interval X < (), and
X will not contain the decimal constant 0.1. Good programming
practice would be to avoid such constructs.
7. In general, if X is an interval, then the assignment
X = <>
does not require that an interval containing the exact mathematical
value of <> be stored in X. Ideally,
whenever the compiler performs this assignment without
containment, the compiler should issue a warning.
Conversion between intervals of different types is also
possible, either implicitly or with the function IVAL. Conversion
with IVAL is with the expression IVAL(X,KIND), where X is the
interval to be converted and KIND is the kind type parameter of the
target interval. The result IVAL(X,TYPE) shall be an interval of
type KIND that contains or is equal to the exact interval
represented by X.
Implicit conversion between intervals of one kind and another
within expressions obeys the same rules as conversion between
reals and integers and intervals. In expressions that contain
more than one kind of interval, the conversion shall be to the
kind corresponding to highest precision.
Example: Suppose the kind type parameter 1 represents "single
precision" intervals, while the kind type parameter 2 represents
"double precision" intervals, and consider the following
program.
PROGRAM CONVERT_INTERVALS
INTERVAL(KIND=2) X=(<0.1_2>)
INTERVAL(KIND=1) Y
Y = X
END PROGRAM CONVERT_INTERVALS
Execution of this program will cause a machine interval that
contains the mathematical number 0.1 to be stored into the
single precision interval variable Y.
Note: Mixedmode INTERVAL/COMPLEX is not defined.
Implicit conversion from interval to other data types, except to
other kinds of interval, shall not be allowed.
Note: The functions INF and SUP defined below may be used to
convert an INTERVAL to another data type. Similarly, MID may
also be viewed as producing a real approximation to an
interval, while IVAL converts from real to interval.
In some computations, the implicit or default conversions will
not be adequate to ensure enclosure of the intended quantity.
Example:
3.14159 is a sixdigit representation of the
mathematical constant pi. However, suppose X is an interval
variable of a type with precision equal to approximately 15
decimal digits. Then execution of the statement
X = 3.14159
may cause an interval to be stored into X that is contained in
the mathematical interval [3.141589999999999,3.141590000000001].
Yet this interval does not contain the actual value of pi.
This is the reason for the following special conversion function
CONVERT_DECIMAL_DIGITS.
The conversion function CONVERT_DECIMAL_DIGITS shall be of the form
CONVERT_DECIMAL_DIGITS(, [, ])
where is a character string representing a valid
integer or real literal constant, where is a
nonnegative integer constant or variable, and where is a
valid interval kind type parameter. The return value of
CONVERT_DECIMAL_DIGITS shall be an interval of kind . This
interval shall contain the mathematical interval [l,u], where
l and u are constructed from by subtracting and
adding 5 to the +1 digit of , respectively.
Thus, after being rounded to digits, l and u will be equal,
but after being rounded to +1 digits, they will not.
Note: Thus, the function CONVERT_DECIMAL_DIGITS provides an
internal mechanism for constructing intervals from a single
number in precisely the same way as the SE and SF Edit
Descriptors (to be described later in I/O).
Example: Suppose a kind parameter of 2 corresponds to real and
interval data types with approximately 15 decimal digits of
precision. Then CONVERT('3.14159',6,2) will be an interval that
contains the mathematical interval [3.141585, 3.14195].
CONVERT('3.14159',20,2) will be a machine interval (kind
parameter 2) of nonzero width that contains the decimal
constant 3.14159. This is precisely the same result that
will be obtained using single number I/O.
The function CONVERT_WITHIN_BOUNDS shall be available to convert
REAL values with known relative error bounds, as follows.
CONVERT_WITHIN_BOUNDS(R,EPS[,KIND])
where R is a REAL value to be converted, EPS is a real value
that represents relative error, and the optional argument KIND
represents the target kind type parameter, shall return an
interval that contains the mathematical interval
[R*(1EPS),R*(1+EPS)]. This holds true even if R or EPS or both
are literal constants.
Note: Design of CONVERT_WITHIN_BOUNDS has similar considerations
as the design of IVAL; see the notes below IVAL.
Example: Suppose it is known that the real value R has been
computed to within an accuracy of one part in a
thousand, and R happens to be equal to 2. Then the
call
CONVERT_WITHIN_BOUNDS(R,0.001)
will return an INTERVAL of default kind that contains
the mathematical interval [0.999,1.001].
New Infix Operators
___________________
The following infix operators shall be a part of standard interval
support.
Syntax function
______ ________
Z = X.IS.Y Z < intersection of X and Y, that is,
[max{xl,yl},min{xu,yu}] if
max{xl,yl} < = min{xu,yu} and
signals an "intersection of disjoint
intervals" otherwise.
Z = X.CH.Y Z < [min{xl,yl}, max{xu,yu}]
("interval hull" of X and Y. The mnemonic is
"convex hull")
X.SB.Y .TRUE. if X is a subset of Y
( i.e. if xl >= yl .AND. xu <= yu )
X.PRSB.y .TRUE. if X is a proper subset of Y
( i.e. if X.SB.Y .AND. (xl > yl .OR. xu < yu )
X.SP.Y .TRUE. if and only if Y.SB.X is true
( i.e. if xl <= yl .AND. xu >= yu )
X.PRSP.Y .TRUE. if and only if Y.PSB.X is true
( i.e. if Y.SB.X .AND. (yl > xl .OR. yu < xu )
X.DJ.Y .TRUE. if X and Y are disjoint sets
( i.e. if xl > yu or xu < yl )
R.IN.X .TRUE. if the REAL value R is contained in the
interval X (i.e. if xl <= R <= xu)
Note: Intervals are closed. So, if R.IN.X, then R may be equal
to one of the endpoints of X.
Interval Versions of Relational Operators
_________________________________________
The following relational operators shall be extended to interval
operations, in the "certainly true" sense. That is, the result
is .TRUE. if and only if it is true for each pair of real values
taken from the corresponding interval values.
Syntax function
______ ________
X.LT.Y .TRUE. if xu < yl
X.GT.Y .TRUE. if xl > yu
X.LE.Y .TRUE. if xu <= yl
X.GE.Y .TRUE. if xl >= yu
Another set of relational operators, the POSSIBLY TRUE relationals,
shall be defined as follows.
Syntax function
______ ________
X.PLT.Y .TRUE. if xl < yu (i.e. if .NOT.(X.GE.Y) )
X.PGT.Y .TRUE. if xu > yl (i.e. if .NOT.(X.LE.Y) )
X.PLE.Y .TRUE. if xl <= yu (i.e. if .NOT.(X.GT.Y) )
X.PGE.Y .TRUE. if xu >= yl (i.e. if .NOT.(X.LT.Y) )
Finally, equality and inequality of intervals are defined by viewing
the intervals as sets.
Syntax function
______ ________
X.EQ.Y .TRUE. if xl=yl and xu=yu
X.NE.Y .TRUE. if .NOT. (X.EQ.Y)
Operator Precedence
________ __________
The precedence of these operators is as follows:
Category of Operators Precedence
operation
_____________________________________________________________________
_____________________________________________________________________
Extension Highest
Numeric ** .
Numeric *, /, .IS. .
Numeric unary + or  .
Numeric binary + or , .CH. .
Character // .
Relational .EQ., .NE., .LT., .PLT., .LE., .PLE., .GT., .PGT., .
.GE., .PGE., .SB., .PRSB., .SP., .PRSP., .DJ.,.IN. .
Logical .NOT. .
Logical .AND. .
Logical .OR. .
Logical .EQV. or .NEQV. .
Extension Lowest
_____________________________________________________________________
Special Interval Functions
_______ ________ _________
The following utility functions shall be provided for conversion from
INTERVAL to REAL, etc.
Syntax function attainable accuracy
______ ________ ___________________
R = INF(X) Lower bound of X (the value in the lower
storage unit of the interval
datum X)
R = SUP(X) Upper bound of X (the value in the upper
storage unit of the interval
datum X)
R = MID(X) Midpoint of X (a floating point approximation,
always greater than or equal to
the value returned by INF and
less than or equal to the value
returned by SUP)
R = WID(X) R < xu  xl (the value shall be rounded up,
"Width" to be greater than or equal to
the actual value)
R = MAG(X) R < max { xl, xu }
"Magnitude"
 min { xl, xu } if .NOT.(0.IN.X),
R = MIG(X) R <
 0 otherwise.
"Mignitude"

Z = ABS(X) Z <  [min{x}, max{x}]
 x.IN.X x.IN.X
Range of absolute value
Z = MAX(X,Y) Z < [max {xl,yl}, max {xu,yu}]
Range of maximum
MAX shall be extended analogously
for more than two
arguments.
Z = MIN(X,Y) Z < [min {xl,yl}, min {xu,yu}]
Range of minimum
MIN shall be extended analogously
for more than two
arguments.
N = NDIGITS(X) Number of leading decimal digits that are the same in
xl and xu. n digits shall be counted as the same if
rounding xl to the nearest decimal number with n
significant digits gives the same result as rounding
xu to the nearest decimal number with n significant
digits. If xl=xu, , then NDIGITS(X) shall return
PRECISION(R)+1, where R is any real of kind the same as
the kind of the INTERVAL X (that is, whose precision
is the same as the precision of X).
Note: On many machines, INF, SUP, MAG, MIG, ABS, MAX, and MIN
can be exact, if the target is of a type that corresponds to the
input. This is because these functions merely involve storing
one of the endpoints of the interval into the target variable.
Similarly, the conversion function IVAL can be exact on such
machines if it specifies conversion from REAL data of
corresponding type.
Except for NDIGITS, all of these special interval functions
shall be elemental.
Note regarding NDIGITS: For example, if X = [0.1996,0.2004],
then three leading decimal digits of this function are the same,
and NDIGITS(X) is equal to 3. This is because, if .1996 and .2004
are each rounded to the nearest decimal number with three
significant digits, they both round to .200, yet they round to
different fourdigit decimal numbers. This value can
be computed as INT(LOG10(2.0041.996). In many cases,
the value can be computed as INT(LOG10(xuxl) + e LOG10(b)),
where e = MAX(EXPONENT(xl),EXPONENT(xu)) and b=RADIX(X).
Note: three interval functions, MAG, MIG, and ABS, correspond to
the point intrinsic ABS. The specification of ABS is as the
range of the absolute value function, consistent with the general
principle that the results of interval functions shall contain
the ranges of corresponding point intrinsics. Although "MAG(X)"
is written X in much of the interval literature, it is more
natural in various applications to have ABS(X) denote the range
of the absolute value function.
If WID(X) is not exact, then its value shall be upwardly rounded.
Note: WID(X) often appears in convergence criteria of the form
WID(X) < EPS. The criterion is certain to be satisfied if the
computed value WID(X), used in the comparison, is greater than
or equal to the exact value.
Optimization of Interval Expressions
____________ __ ________ ___________
In an interval expression, any transformations by the optimizer
that are permissible for REAL's shall also be permissible for
INTERVAL's. For example 1/x/y may be replaced by 1/(x*y).
Note: Different code optimizations may yield significantly
different interval values, but each will be an enclosure of the
range of the expression over the input intervals.
Note: If the program requires the precise form of an interval
expression, then either parentheses may be used or
optimization may be turned off through compiler options.