C H A P T E R 1 
Using Interval Arithmetic With f95 
Interval arithmetic is a system for computing with intervals of numbers. Because interval arithmetic always produces intervals that contain the set of all possible result values, interval algorithms have been developed to perform surprisingly difficult computations. For more information on interval applications, see the Interval Arithmetic Readme.
Since the inception of interval arithmetic, interval algorithms that produce narrowwidth results have been developed, and the syntax and semantics for interval language support have been designed. However, relatively little progress has been made in providing commercially available and supported interval compilers. With one exception (M77 Minnesota FORTRAN 1977 Standards Version Edition 1), interval systems have been based on preprocessors, C++ classes, or Fortran 90 modules. The goals of intrinsic compiler support for interval data types in f95 are:
Sun Studio Fortran 95 interval support is a significant extension to Fortran.
The goal of intrinsic INTERVAL support in f95 is to stimulate development of commercial interval solver libraries and applications by providing program developers with:
Support and features are components of implementation quality. Not all possible quality of implementation features have been implemented. Throughout this book, various unimplemented quality of implementation opportunities are described. Additional suggestions from users are welcome.
As a consequence of evaluating any interval expression, a valid intervalsupporting compiler must produce an interval that contains the set of all possible results. The set of all possible results is called the containment set (cset) of the given expression. The requirement to enclose an expression's cset is the containment constraint of interval arithmetic. The failure to satisfy the containment constraint is a containment failure. A silent containment failure (with no warning or documentation) is a fatal error in any interval computing system. By satisfying this single constraint, intervals provide otherwise unprecedented computing quality.
Given the containment constraint is satisfied, implementation quality is determined by the location of a point in the twodimensional plane whose axes are runtime and interval width. On both axes, small is better. How to trade runtime for interval width depends on the application. Both runtime and interval width are obvious measures of intervalsystem quality. Because interval width and runtime are always available, measuring the accuracy of both interval algorithms and implementation systems is no more difficult than measuring their speed.
The Sun Studio tools for performance profiling can be used to tune interval programs. However, in f95, no intervalspecific tools exist to help isolate where an algorithm may gain unnecessary interval width. As described in Section 1.4, Code Development Tools, some interval dbx and global program checking (GPC) support are provided. Adding additional intervalspecific code development and debugging tools are quality of implementation opportunities.
All the normal language and compiler quality of implementation opportunities exist for intervals, including rapid execution and easeofuse.
Valid interval implementation systems include a new additional quality of implementation opportunity: Minimize the width of computed intervals while always satisfying the containment constraint.
If an interval's width is as narrow as possible, it is said to be sharp. For a given floatingpoint precision, an interval result is sharp if its width is as narrow as possible.
The following can be said about the width of intervals produced by the f95 compiler:
By providing compiler optimization and hardware instruction support, INTERVAL operations are not necessarily slower than their REAL floatingpoint counterparts. In f95, the following can be said about the speed of intrinsic interval operators and mathematical functions:
The intrinsic INTERVAL data type in Fortran facilitates interval code development, testing, and execution. To make interval code transparent (easy to write and read), interval syntax and semantics have been added to Fortran. User acceptance will ultimately determine which interval features are added to standard Fortran.
By introducing intervals as an intrinsic data type to Fortran, all of the applicable syntax and semantics of the Fortran language become immediately available. Sun Studio Fortran 95 includes the following intervalspecific Fortran extensions:
For examples and more information on these and other intrinsic interval functions, see CODE EXAMPLE 111 through CODE EXAMPLE 114 and Section 2.10.4.5, Intrinsic Functions.
Chapter 2 contains detailed descriptions of these and other interval features.
The examples in this section are designed to help new interval programmers to understand the basics and to quickly begin writing useful interval code. Modifying and experimenting with the examples is strongly recommended.
All code examples in this book are contained in the directory:
/opt/SUNWspro/examples/intervalmath/docExamples
The name of each file is cenm.f95, where n is the chapter in which the example occurs, and m is the number of the example. Additional interval examples are contained in the directory:
/opt/SUNWspro/examples/intervalmath/general
The following f95 commandline macro is the simplest way to invoke recognition of INTERVAL data types as intrinsic and to control INTERVAL expression processing:
For intrinsic INTERVAL data types to be recognized by the compiler, either xia or xinterval must be entered in the f95 command line.
All commandline options that interact with intervals are described in Section 2.3.3, Interval CommandLine Options. Widestneed and strict expression processing are described in Section 2.3, INTERVAL Arithmetic Expressions.
The simplest commandline invocation of f95 with interval support is shown in CODE EXAMPLE 11.
Unless explicitly stated otherwise, all code examples are compiled using the xia commandline option. The xia or xinterval commandline option is required to use the interval extensions to f95.
CODE EXAMPLE 11 is the interval equivalent of "hello world."
math% cat ce11.f95 PRINT *, "[2, 3] + [4, 5] = ", [2, 3] + [4, 5] ! line 1 END math% f95 xia ce11.f95 math% a.out [2, 3] + [4, 5] = [6.0,8.0] 
CODE EXAMPLE 11 uses listdirected output to print the labeled sum of the intervals [2, 3] and [4, 5].
The INTERVAL declaration statement performs the same functions for INTERVAL data items as the REAL, INTEGER, and COMPLEX declarations do for their respective data items. The default INTERVAL kind type parameter value (KTPV) is twice the default INTEGER KTPV. This permits any default INTEGER to be exactly represented using a degenerate default INTERVAL. See Section 1.3.7, Default Kind Type Parameter Value (KTPV) for more information.
CODE EXAMPLE 12 uses INTERVAL variables and initialization to perform the same operation as CODE EXAMPLE 11.
In line 1, the variables, X and Y are declared to be default type INTERVAL variables and are initialized to [2, 3] and [4, 5], respectively. Line 2 uses listdirected output to print the labeled interval sum of X and Y.
Full support for reading and writing intervals is provided. Reading and writing INTERVAL and COMPLEX data items are similar. Intervals use square brackets, instead of parentheses as delimiters.
In f95 the input conversion process constructs a sharp interval that contains the input decimal value. If the value is machine representable, the internal machine approximation is degenerate. If the value is not machine representable, an interval having width of 1ulp (unitinthelastplace of the mantissa, pronounced "ulp") is constructed.
The simplest way to read and print INTERVAL data items is with listdirected input and output.
CODE EXAMPLE 13 is a simple tool to help users become familiar with interval arithmetic and singlenumber INTERVAL input/output using listdirected READ and PRINT statements. Complete support for formatted INTERVAL input/output is provided, as described in Section 2.10.2, Input and Output.
Note  The interval containment constraint requires that directed rounding be used during input and output. With singlenumber input followed immediately by singlenumber output, a decimal digit of accuracy can appear to be lost. In fact, the width of the input interval is increased by at most 1ulp, when the input value is not machine representable. See Section 1.3.5, SingleNumber Input/Output and CODE EXAMPLE 16 
Note  The empty interval is supported in f95. The empty interval can be entered as "[empty]". Infinite interval endpoints are also supported, as described in Section 2.10.2.1, External Representations and illustrated in CODE EXAMPLE 237. 
One of the most frustrating aspects of reading interval output is comparing interval infima and suprema to count the number of digits that agree. For example, CODE EXAMPLE 14 and CODE EXAMPLE 15 shows the interval output of a program that generates different random width INTERVAL data.
Note  Only program output is shown in CODE EXAMPLE 14 and CODE EXAMPLE 15. The code that generates the output is included with the examples locatedat http://developer.sun.com/prodtech/cc/reference/codesamples/ 
Compare the output readability in CODE EXAMPLE 14 with CODE EXAMPLE 15.
Because reading and interactively entering interval data can be tedious, a singlenumber interval format is introduced. The singlenumber convention is that any number not contained in brackets is interpreted as an interval whose lower and upper bounds are constructed by subtracting and adding 1 unit to the last displayed digit.
Thus during interval input and output,
[2.34499, 2.34501] = 2.34500 + [1, +1]_{uld }
where [1, +1]_{uld} means that the interval [1, +1] is added to the last digit of the preceding number. The subscript, uld, is a mnemonic for "unit in the last digit."
Note  The single number input/output representation is not used to represent INTERVAL literal constants in f95 code. 
To represent a degenerate interval, a single number can be enclosed in square brackets. For example,
[2.345] = [2.345, 2.345] = 2.345000000000.....
This convention is used both for singlenumber input/output and to represent degenerate literal INTERVAL constants in Fortran code. Thus, type [0.1] to enter an exact decimal number, even though 0.1 is not machine representable.
During input to a program, both [0.1,0.1] and [0.1] represents the point, 0.1. However, the singlenumber input/output value 0.1 represents the interval
0.1 + [1, +1]_{uld} = [0, 0.2].
In the singlenumber display format, trailing zeros are significant. See Section 2.10.2, Input and Output for more information.
Intervals can always be entered and displayed using the traditional [inf, sup] display format. In addition, a single number in square brackets denotes a point. For example, on input, [0.1] is interpreted as the number 1/10. To guarantee containment, directed rounding is used to construct an internal approximation that is known to contain the number 1/10.
CODE EXAMPLE 16 notes:
The f95 compiler contains the following INTERVALspecific statements, expressions, and extensions:
In f95 the default INTEGER KTPV is KIND(0) = 4. To represent any default INTEGER with a degenerate default INTERVAL requires the default INTERVAL KTPV, KIND([0]), to be 2*KIND(0) = 8. Choosing 8 for the default INTERVAL KTPV is also done because:
TABLE 11 notes:
1. The letters a and b are placeholders for literal decimal constants, such as 0.1 and 0.2.
2. A single decimal constant contained in square brackets denotes a degenerate INTERVAL constant. The same convention is used in input/output.
3. Let expr stand for any Fortran arithmetic expression, whether or not it contains items of type INTERVAL. An assignment statement, V = expr, evaluates the expression, expr, and assigns the resulting value to V. Mixedmode INTERVAL expressions are not permitted under the xia=strict command line option. Under the xia or xia=widestneed option, mixedmode expressions are correctly evaluated using widestneed expression processing. Before expression evaluation under widestneed, all integer and floatingpoint data items are promoted to containing intervals with the largest KTPV found anywhere in the expression, including, V. For details, see Section 2.3.2, Value Assignment.
4. Interval input/output support is designed to provide flexibility, readability, and ease of code development. The most important new edit descriptor is Y, which is used to read and display intervals using the singlenumber interval format. For a complete description of all edit descriptors that can process intervals, see Section 2.10.2, Input and Output.
The INTERVAL assignment statement assigns the value of an interval expression, denoted by the placeholder expr, to an INTERVAL variable, array element, array, array section, or structure component V. The syntax is:
where V must have an INTERVAL type, and expr denotes any nonCOMPLEX numeric expression. Under widestneed expression processing, the expression expr need not be an INTERVAL expression. Under strict expression processing, expr must be an INTERVAL expression with the same KTPV as V.
Gracefully handling mixedtype INTERVAL expressions is an important easeofuse feature, because it facilitates writing transparent (easy to understand) mathematical expressions.
Mixedtype INTERVAL expressions are supported to make writing and reading interval code no more difficult than it is for REAL code. The interval containment constraint is satisfied in mixedmode expressions using either widestneed or strict expression processing.
Computing narrowwidth interval results is facilitated if the width of INTERVAL constants is dynamically defined by expression context, as described in Section 2.3, INTERVAL Arithmetic Expressions. In mixedKTPV expressions, shown in CODE EXAMPLE 17, dynamically increasing the KTPV of INTERVAL variables can also decrease the width of INTERVAL expression results.
In line 3, KTPV_{max} = KIND(Z) = 8. This value is used to promote the KTPV of X and Y to 8 before computing their product and storing the result in Z1.
These steps are shown explicitly in the equivalent strict code in line 4.
The process of scanning a statement to determine the maximum KTPV and performing the necessary promotions, is called widestneed expression processing, see Section 2.3, INTERVAL Arithmetic Expressions.
For syntax and semantics of the intrinsic INTERVAL constructor functions, see Section 2.9, Extending Intrinsic INTERVAL Operators.
If the widestneed principle is used with both KTPVs and data types, mixedmode (type and KTPV) INTERVAL expressions can be safely and predictably evaluated. For example, in CODE EXAMPLE 18, the expression for Y1 in line 3 is an interval expression, because X and Y1 are INTERVAL variables.
To guarantee containment, a containing interval must be used in place of a real approximation to the constant 0.1. However, KTPV_{max} = 16, because KIND(X) = 16. Therefore, the INTERVAL constant [0.1_16], a sharp KTPV = 16 interval containing the exact value, 1/10, is used to update X. Finally, the result is converted to a KTPV = 4 containing interval and assigned to Y1. Line 4 contains the equivalent strict code. Under strict expression processing, neither mixedtype nor mixedKTPV expressions are permitted.
The logical steps in widestneed expression processing are:
1. Scan the entire statement, including the lefthand side, for any INTERVAL data items.
The presence of any INTERVAL constants, variables, or intrinsic functions, makes the expression's type INTERVAL.
2. Scan the INTERVAL expressions for KTPV_{max}, based on the KTPV of each INTERVAL, REAL, INTEGER, constant, or variable.
Note  Integers are converted to intervals with twice their KTPV so all integer values can be exactly represented. 
3. Promote all variables and constants to intervals with KTPV_{max}.
5. Convert the result to a lower KTPV if needed to match the lefthand side's KTPV.
6. Assign the resulting value to the lefthand side.
These steps guarantee that mixedmode INTERVAL expression processing satisfies the containment constraint and efficiently produces reasonably narrow interval results.
Mixedmode INTERVAL expression evaluation using widestneed expression processing is supported by default with the xia commandline flag. Using xia=strict eliminates any automatic type conversions to intervals and any automatic KTPV increases of INTERVAL variables. In strict mode, all interval type and precision conversions must be explicitly coded.
Writing arithmetic expressions that contain INTERVAL data items is simple and straightforward. Except for INTERVAL literal constants and intrinsic INTERVALspecific functions, INTERVAL expressions look like REAL arithmetic expressions. In particular, with widestneed expression processing, REAL and INTEGER variables and literal constants can be freely used anywhere in an INTERVAL expression, such as in CODE EXAMPLE 19.
Because X, the variable to which the assignment is made in line 5, is an INTERVAL, the following steps are taken before evaluating the expression 0.1*A/N:
1. The literal constant 0.1 is converted to the default INTERVAL variable containing the degenerate interval [0.1].
While not required in a valid interval system implementation, Sun Studio Fortran 95 performs sharp data conversions. For example, the internal approximation of [0.1] is 1ulp wide.
2. The REAL variable A is converted to the degenerate interval [5].
3. The INTEGER variable N is converted to the degenerate interval [3].
The expression [0.1] × [5]/[3] is evaluated using interval arithmetic. The above steps are part of widestneed expression processing, which is required to satisfy the containment constraint when evaluating mixedmode INTERVAL expressions. See Section 1.3.9, MixedType Expression Evaluation.
An INTERVAL assignment statement must satisfy one requirement: the variable to which the assignment is made must be an INTERVAL variable, array element, array, array section, or structure component. For more information on the widestneed processing mode, see Section 2.3.1, MixedMode INTERVAL Expressions.
Because the interval system implemented in Sun Studio Fortran 95 is closed, if any INTERVAL expression fails to produce a valid interval result, it is a compiler error that should be reported. See Section 1.4, Code Development Tools for information on how to report a suspected error and Section 1.5.1, Known Containment Failures for a list of known errors.
Ordering intervals is more complicated than ordering points. Testing whether 2 is less than 3 is unambiguous. With intervals, while the interval [2,3] is certainly less than the interval [4,5], what should be said about [2,3] and [3,4]?
Three different classes of INTERVAL relational operators are implemented:
For a certainlyrelation to be true, every element of the operand intervals must satisfy the relation. A possiblyrelation is true if it is satisfied by any elements of the operand intervals. The setrelations treat intervals as sets. The three classes of INTERVAL relational operators converge to the normal relational operators on points if both operand intervals are degenerate.
To distinguish the three operator classes, the normal twoletter Fortran relation mnemonics are prefixed with the letters C, P, or S. In f95 the set operators .SEQ. and .SNE. are the only operators for which the point defaults (.EQ. or == and .NE. or /=) are supported. In all other cases, the relational operator class must be explicitly identified, as for example in:
See Section 2.4, Intrinsic Operators for the syntax and semantics of all INTERVAL operators.
The following program demonstrates the use of a setequality test.
math% cat ce110.f95 INTERVAL :: X = [2, 3], Y = [4, 5] ! Line 1 IF(X+Y .SEQ. [6, 8]) PRINT *, "Check." ! Line 2 END math% f95 xia ce110.f95 math% a.out Check. 
Line 2 uses the setequality test to verify that X+Y is equal to the interval [6, 8].
IF(X+Y == [6, 8]) PRINT *, "Check." ! line 2
Use CODE EXAMPLE 111 and CODE EXAMPLE 112 to explore the result of INTERVALspecific relational operators.
CODE EXAMPLE 112 demonstrates the use of the INTERVALspecific operators ed in TABLE 12.
.IH. 

.IX. 

.DJ. 


.IN. 

.INT. 

.PSB. 

.PSP. 

.SB. 

.SP. 
A variety of intrinsic INTERVALspecific functions are provided. See Section 2.10.4.5, Intrinsic Functions. Use CODE EXAMPLE 113 to explore how intrinsic INTERVAL functions behave.
Every Fortran intrinsic function that accepts REAL arguments has an interval version. See Section 2.10.4.5, Intrinsic Functions. Use CODE EXAMPLE 114 to explore how some intrinsic functions behave.
Information on interval code development tools is available online. See the Interval Arithmetic Readme for a list of interval web sites and other online resources.
To report a suspected interval error, send email to
Include the following text in the Subject line of the email message:
FORTEDEV "7.0 mm/dd/yy" Interval
where mm/dd/yy is the month, day, and year.
In Sun Studio, interval data types are supported by dbx to the following extent:
For additional details on dbx functionality, see Debugging a Program With dbx.
Global program checking (GPC) in Sun Studio Fortran 95 detects one intervalspecific error: INTERVAL type mismatches in usersupplied routine calls. CODE EXAMPLE 115 shows an example of GPC detecting an INTERVAL type mismatch.
The following libraries contain intrinsic INTERVAL routines.
There is limited legacy interval Fortran code with which to contend. Until language syntax and semantics are standardized, different providers of interval compiler support will inevitably diverge. The standardization process will be facilitated if users provide feedback regarding the most favored INTERVAL syntax and semantics. Comments can be sent to the email alias ed in the Interval Arithmetic Readme.
The representation of intervals in binary files will change as compilers supporting narrower interval systems are made available.
In this release, the autopar compiler option has no effect on loops containing interval arithmetic operations. These loops are not automatically parallelized. The explicitpar compiler option must be used to parallelize loops marked with explicit parallelization directives.
The following code samples list intervalspecific error messages. Each code sample includes the error message and the sample code that produced the error.
Whenever an interval containment failure can occur, a compiletime warning should be issued. There are no know containment failures under widestneed expression processing. In xia=strict mode, it is possible to violate the containment constraint with an interval ** (integer expression) operation if the integer expression overflows.
Numerical inaccuracies are normally associated with REAL rather than INTEGER expressions. In one respect, INTEGER expressions are harder to detect than REAL expressions. When REAL expressions overflow, an exception is raised and an IEEE infinity is generated. The exception is a warning that overflow has occurred. Infinities tend to propagate in floatingpoint computations, thereby alerting users of a potential problem. It is also possible to trap on overflow.
When INTEGER expressions overflow, they silently wrap around to some possiblyoppositesigned value. Moreover, the only practical way to detect integer overflow is to perform the inverse operation and test for equality on every integer operation. Integer constant expressions are safe because they are evaluated during compilation where overflow is detected and signalled with a warning message.
Under xia=widestneed expression processing when the second operand of the ** operator is an integer expression that overflows, the returned interval is guaranteed to contain the correct result. However, the same is not true under xia=strict processing, because it is not possible to promote integers to intervals prior to evaluating the given expression without widestneed expression processing, The same is true if the second operand of the ** operator is the INTERVAL type conversion routine.
CODE EXAMPLE 121 shows that widestneed expression processing is extended to all intrinsic INTEGER operations and functions inside integer expressions in the second operand of the ** operator. This is not true under xia=strict mode.
This code demonstrates a silent containment failure in xia=strict mode and the correct interval results in xia=widestneed mode. For information on the power operator, see Section 2.5, Power Operators X**N and X**Y.
Copyright © 2004, Sun Microsystems, Inc. All Rights Reserved.