Interval Arithmetic Programming Reference

Chapter 1

Using Interval Arithmetic With f95

1.1 f95 INTERVAL Type and Interval Arithmetic Support

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 narrow-width 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 pre-processors, C++ classes, or Fortran 90 modules. The goals of intrinsic compiler support for interval data types in f95 are:

• Reliability
• Speed
• Ease-of-use

Interval support in the Sun WorkShop 6 release of f95 is a significant extension to Fortran.

1.2 f95 Interval Support Goal: Implementation Quality

The goal of intrinsic INTERVAL support in f95 is to stimulate development of commercial interval solver libraries and applications by providing program developers with:

• Quality interval code
• Narrow-width interval results
• Rapidly executing interval code
• An easy to use interval software development environment that includes interval-specific language support and compiler features

Support and features are components of implementation quality. Throughout this book, various quality of implementation opportunities are described. Additional suggestions from users are welcome.

1.2.1 Quality Interval Code

As a consequence of evaluating any interval expression, a valid interval-supporting compiler must produce an interval that contains the set of all possible results. The requirement to contain the set of all possible results is called the containment constraint of interval arithmetic. The failure to satisfy the containment constraint is a containment failure. A silent failure (with no warning or documentation) to satisfy the interval containment constraint is a fatal error in any interval computing system. By satisfying this single constraint, intervals provide unprecedented computing quality.

Given the containment constraint is satisfied, implementation quality is determined by the location of a point in the two-dimensional 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 interval-system 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 WorkShop 6 tools for performance profiling can be used to tune interval programs. However, in f95, no interval-specific 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 interval-specific code development and debugging tools are quality of implementation opportunities.

1.2.2 Narrow-Width Interval Results

All the normal language and compiler quality of implementation opportunities exist for intervals, including rapid execution and ease-of-use.

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 floating-point 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:

• Individual intervals are sharp approximations of constants.
• Individual interval arithmetic operators produce sharp results.
• Intrinsic mathematical functions usually produce sharp results.

1.2.3 Rapidly Executing Interval Code

By providing compiler optimization and hardware instruction support, INTERVAL operations are not necessarily slower than their REAL floating-point counterparts. In f95, the following can be said about the speed of intrinsic interval operators and mathematical functions:

• Arithmetic operations are reasonably fast.
• The speed of default INTERVAL mathematical functions is generally less than 2 times that of their DOUBLE PRECISION counterparts. KIND = 4 intrinsic interval math functions are provided, but are not tuned for speed (unlike their KIND = 8 counterparts). KIND = 16 mathematical functions are not provided in this release. However, other INTERVAL KIND = 16 functions are supported.
• The following intrinsic INTERVAL array functions are optimized for performance:
  
  
  • SUM
    
  • PRODUCT
    
  • DOT_PRODUCT
    
  • MATMUL

1.2.4 Easy to Use Development Environment

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 Fortran become immediately available. Sun WorkShop 6 f95 includes the following interval Fortran extensions:

• INTERVAL data types
• INTERVAL arithmetic operations and intrinsic mathematical functions form a closed mathematical system. (This means that valid results are produced for any possible operator-operand combination, including division by zero and other indeterminate forms involving zero and infinities.)
• Three classes of interval relational operators:
  
  
  • Certainly
    
  • Possibly
    
  • Set
• Intrinsic INTERVAL-specific operators, such as .IX. (intersection) and .IH. (interval hull)
• INTERVAL-specific functions, such as INF, SUP, and WID
• INTERVAL input/output, including single-number input/output
• Expression-context-dependent INTERVAL constants
• Mixed-mode interval expression processing

For examples and more information on these and other intrinsic interval functions, see CODE EXAMPLE 1-11 through CODE EXAMPLE 1-14 and Section 2.9.4.4 Intrinsic Functions.

Chapter 2 describes these and other interval features.

1.3 Writing Interval Code for f95

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 cen-m.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

1.3.1 Command-Line Options

Including the following command-line macro in the f95 command line invokes recognition of INTERVAL data types as intrinsic and controls INTERVAL expression processing:

• Compiler support for widest-need interval expression processing is invoked by including:
  -xia or -xia=widestneed
• Compiler support for strict interval expression processing is invoked by including:
  -xia=strict

For intrinsic INTERVAL data types to be recognized by the compiler, either -xia or -xinterval must be entered in the f95 command line.

All command-line options that interact with intervals are described in Section 2.3.3 Interval Command-Line Options. Widest-need and strict expression processing are described in Section 2.3 INTERVAL Arithmetic Expressions.

The simplest command-line invocation of f95 with interval support is shown in CODE EXAMPLE 1-1.

1.3.2 Hello Interval World

Unless explicitly stated otherwise, all code examples are compiled using the -xia command-line option. The -xia command-line option is required to use the interval extensions to f95.

CODE EXAMPLE 1-1 is the interval equivalent of "hello world."

CODE EXAMPLE 1-1   Hello Interval World

math% cat ce1-1.f95 PRINT *, "[2, 3] + [4, 5] = ", [2, 3] + [4, 5] ! line 1 END math% f95 -xia ce1-1.f95 math% a.out [2, 3] + [4, 5] = [6.0,8.0]

CODE EXAMPLE 1-1 uses list-directed output to print the labeled sum of the intervals [2, 3] and [4, 5].

1.3.3 Interval Declaration and Initialization

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 1-2 uses INTERVAL variables and initialization to perform the same operation as CODE EXAMPLE 1-1.

CODE EXAMPLE 1-2   Hello Interval World with INTERVAL Variables

math% cat ce1-2.f95 INTERVAL :: X = [2, 3], Y = [4, 5] ! Line 1 PRINT *, "[2, 3] + [4, 5] = ", X+Y ! Line 2 END math% f95 -xia ce1-2.f95 math% a.out [2, 3] + [4, 5] = [6.0,8.0]

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 list-directed output to print the labeled interval sum of X and Y.

1.3.4 INTERVAL Input/Output

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. Because reading and interactively entering interval data can be tedious, a single-number interval format is introduced. The single-number 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

2.345 = [2.344, 2.346],

2.34500 = [2.34499, 2.34501],

and

23 = [22, 24].

Symbolically,

[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."

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 input and representing degenerate literal INTERVAL constants in Fortran code. Thus, type [0.1] to indicate the input value is an exact decimal number, even though 0.1 is not machine representable.

For example, during input to a program, [0.1,0.1] = [0.1] represents the point, 0.1, while using single-number input/output, 0.1 represents the interval

0.1 + [-1, +1]_uld = [0, 0.2].

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 1-ulp (unit-in-the-last-place of the mantissa) is constructed.

---

Note – A uld and an ulp are different. A uld refers to implicitly constructing an interval using the single number input/output format to add and subtract one unit to and from the last displayed digit. An ulp is the smallest possible increment or decrement that can be made to an internal machine number.

---

The simplest way to read and print INTERVAL data items is with list-directed input and output.

CODE EXAMPLE 1-3 is a simple tool to help users become familiar with interval arithmetic and single-number INTERVAL input/output using list-directed READ and PRINT statements. Complete support for formatted INTERVAL input/output is provided, as described in Section 2.9.2 Input and Output.

---

Caution – The interval containment constraint requires that directed rounding be used both during input and output. With single-number input followed immediately by single-number output, a decimal digit of accuracy can appear to be lost. In fact, the width of the input interval is increased by at most 1-ulp, when the input value is not machine representable. See Section 1.3.5 Single-Number Input/Output and CODE EXAMPLE 1-6.

CODE EXAMPLE 1-3   Interval Input/Output  

math% cat ce1-3.f95 INTERVAL :: X, Y INTEGER :: IOS = 0 PRINT *, "Press Control/D to terminate!" WRITE(*, 1, ADVANCE = 'NO') READ(*, *, IOSTAT = IOS) X, Y DO WHILE (IOS >= 0) PRINT *, " For X =", X, ", and Y =", Y PRINT *, "X+Y =", X+Y PRINT *, "X-Y =", X-Y PRINT *, "X*Y =", X*Y PRINT *, "X/Y =", X/Y PRINT *, "X**Y =", X**Y WRITE(*, 1, ADVANCE = 'NO') READ(*, *, IOSTAT=IOS) X, Y END DO 1 FORMAT(" X, Y = ? ") END math% f95 -xia ce1-3.f95 math% a.out Press Control/D to terminate! X, Y = ? [1,2] [3,4] For X = [1.0,2.0] , and Y = [3.0,4.0] X+Y = [4.0,6.0] X-Y = [-3.0,-1.0] X*Y = [3.0,8.0] X/Y = [0.25,0.66666666666666675] X**Y = [1.0,16.0] X, Y = ? [1,2] -inf For X = [1.0,2.0] , and Y = [-Inf,-1.7976931348623157E+308] X+Y = [-Inf,-1.7976931348623155E+308] X-Y = [1.7976931348623157E+308,Inf] X*Y = [-Inf,-1.7976931348623157E+308] X/Y = [-1.1125369292536012E-308,0.0E+0] X**Y = [0.0E+0,Inf] X, Y = ? <Control-D>

---

1.3.5 Single-Number Input/Output

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 1-4 and CODE EXAMPLE 1-5 shows the interval output of a program that generates different random width INTERVAL data.

---

Note – Only program output is shown in CODE EXAMPLE 1-4 and CODE EXAMPLE 1-5. The code that generates the output is included with the examples located in the /opt/SUNWspro/examples/intervalmath/docExamples directory.

CODE EXAMPLE 1-4   [inf, sup] Interval Output 

math% f95 -xia ce1-4.f95 math% a.out Press Control/D to terminate! Enter number of intervals, KTPV (4,8,16) and 1 for single-number output: 5,4,0 [ 0.2017321E-029, 0.2017343E-029] [ 0.2176913E-022, 0.2179092E-022] [-0.3602303E-006,-0.3602302E-006] [-0.3816341E+038,-0.3816302E+038] [-0.1011276E-039,-0.1011261E-039] Enter number of intervals, KTPV (4,8,16) and 1 for single-number output: 5,8,0 [ -0.3945547546440221E+035, -0.3945543600894656E+035] [ 0.5054960140922359E-270, 0.5054960140927415E-270] [ -0.2461623589326215E-043, -0.2461623343163864E-043] [ -0.2128913523672577E+204, -0.2128913523672576E+204] [ -0.3765492464030608E-072, -0.3765492464030606E-072] Enter number of intervals, KTPV (4,8,16) and 1 for single-number output: 5,16,0 [ 0.199050353252318620256245071374058E+055, 0.199050353252320610759742664557447E+055] [ -0.277386431989417915223682516437493E+203, -0.277386431989417915195943874118822E+203] [ 0.132585288598265472316856821380503E+410, 0.132585288598265472316856822706356E+410] [ 0.955714436647437881071727891682804E+351, 0.955714436647437881071727891683760E+351] [ -0.224211897768824210398306994401732E+196, -0.224211897768824210398306994177519E+196] Enter number of intervals, KTPV (4,8,16) and 1 for single-number output: <Control-D>

---

Compare the output readability in CODE EXAMPLE 1-4 with CODE EXAMPLE 1-5.

CODE EXAMPLE 1-5   Single-Number Output. 

math% a.out Press Control/D to terminate! Enter number of intervals, KTPV (4,8,16) and 1 for single-number output: 5,4,1 0.20173 E-029 0.218 E-022 -0.3602303E-006 -0.38163 E+038 -0.10112 E-039 Enter number of intervals, KTPV (4,8,16) and 1 for single-number output: 5,8,1 -0.394554 E+035 0.505496014092 E-270 -0.2461623 E-043 -0.2128913523672577E+204 -0.3765492464030607E-072 Enter number of intervals, KTPV (4,8,16) and 1 for single-number output: 5,16,1 0.19905035325232 E+055 -0.2773864319894179152 E+203 0.132585288598265472316856822 E+410 0.955714436647437881071727891683 E+351 -0.224211897768824210398306994 E+196 Enter number of intervals, KTPV (4,8,16) and 1 for single-number output: <Control-D>

In the single-number display format, trailing zeros are significant. See Section 2.9.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 1-6   Character Input with Internal Data Conversion 

math% cat ce1-6.f95 INTERVAL :: X INTEGER :: IOS = 0 CHARACTER*30 BUFFER PRINT *, "Press Control/D to terminate!" WRITE(*, 1, ADVANCE='NO') READ(*, '(A12)', IOSTAT=IOS) BUFFER DO WHILE (IOS >= 0) PRINT *, ' Your input was: ', BUFFER READ(BUFFER, '(Y12.16)') X PRINT *, "Resulting stored interval is:", X PRINT '(A, Y12.2)', ' Single number interval output is:', X WRITE(*, 1, ADVANCE='NO') READ(*, '(A12)', IOSTAT=IOS) BUFFER END DO 1 FORMAT(" X = ? ") END math% f95 -xia ce1-6.f95 math% a.out Press Control/D to terminate! X = ? 1.37 Your input was: 1.37 Resulting stored interval is: [1.3599999999999998,1.3800000000000002] Single number interval output is: 1.3 X = ? 1.444 Your input was: 1.444 Resulting stored interval is: [1.4429999999999998,1.4450000000000001] Single number interval output is: 1.44 X = ? <Control-D>

CODE EXAMPLE 1-6 notes:

• Single numbers in square brackets represent degenerate intervals.
• When a non-machine representable number is read using single-number input, conversion from decimal to binary (radix conversion) and the containment constraint force the number's interval width to be increased by 1-ulp (unit in the last place of the mantissa). When this result is displayed using single-number output, it can appear that a decimal digit of accuracy has been lost. This is not so. To echo single-number interval inputs, use character input together with internal READ statement data conversion, as shown in CODE EXAMPLE 1-6.

  ---

  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.9.2.1 External Representations and illustrated in CODE EXAMPLE 2-37.

  ---

1.3.6 Interval Statements and Expressions

The f95 compiler contains the following INTERVAL-specific statements, expressions, and extensions:

• The INTERVAL data type, related instructions, and statements described in TABLE 1-1 are supported.
• All intrinsic functions that accept real arguments have corresponding interval versions. (It is a known error that the same is not true for integers. See Section 1.5.1.1 Integer Overflow.)
• A number of intrinsic INTERVAL-specific functions and operators have been added, including INTERVAL-specific relational operators and set-theoretic functions. For a complete list of intrinsic INTERVAL functions and INTERVAL operators, see Section 2.9.3 Intrinsic INTERVAL Functions and Section 2.9.4 Mathematical Functions.

1.3.7

TABLE 1-1   INTERVAL Specific Statements and Expressions

Statement/expression	Description
INTERVAL INTERVAL(4) INTERVAL(8) INTERVAL(16)	Default INTERVAL type declaration KIND=4 INTERVAL KIND=8 INTERVAL KIND=16 INTERVAL
[a,b] See Note 1	Literal INTERVAL constant: [a,b]
[a] See Note 2	[a,a]
INTERVAL A PARAMETER A=[c,d]	Named constant: A
V = expr See Note 3	Value assignment
FORMAT(E, EN, ES, F, G, VE, VF, VG, VEN, VES, Y) See Note 4	E, EN, ES, F, G, VE, VF, VG, VEN, VES, Y edit descriptors
(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. Mixed-mode INTERVAL expressions are not permitted under the -xia=strict command line option. Under the -xia or -xia=widestneed option, mixed-mode expressions are correctly evaluated using widest-need expression processing. Before expression evaluation under widest-need, all integer and floating-point 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 single-number interval format. For a complete description of all edit descriptors that can process intervals, see Section 2.9.2 Input and Output.

Default Kind Type Parameter Value (KTPV)

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:

• Intervals are often used to perform numerically intense computations, as have been performed on CDC and Cray machines.
• When evaluating a single arithmetic expression, the width of intervals necessarily grows because of accumulated rounding errors, dependence, and cancellation. Extra precision can help to reduce the effect of accumulated rounding errors.

1.3.8 Value Assignment V = expr

The INTERVAL assignment statement assigns the value of an interval expression, denoted by the placeholder expr, to an INTERVAL variable, array element, or array, V. The syntax is:

V = expr

where V must have an INTERVAL type, and expr denotes any non-COMPLEX numeric expression. Under widest-need 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.

1.3.9 Mixed-Type Expression Evaluation

Gracefully handling mixed-type INTERVAL expressions is an important ease-of-use feature, because it facilitates writing transparent (easy to understand) mathematical expressions.

Mixed-type 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 mixed-mode expressions using either widest-need or strict expression processing.

1.3.9.1 Widest-Need and Strict Expression Processing

Computing narrow-width 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 mixed-KTPV expressions, shown in CODE EXAMPLE 1-7, dynamically increasing the KTPV of INTERVAL variables can also decrease the width of INTERVAL expression results.

CODE EXAMPLE 1-7   Mixed Precision with Widest-Need

math% cat ce1-7.f95 INTERVAL(4) :: X = [1, 2], Y = [3, 4] INTERVAL :: Z1, Z2 ! Widest-need Code Z1 = X*Y ! Line 3 ! Equivalent Strict Code Z2 = INTERVAL(X, KIND=8)*INTERVAL(Y, KIND=8) ! Line 4 IF (Z1 .SEQ. Z2) PRINT *, 'Check.' END math% f95 -xia ce1-7.f95 math% a.out Check.

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 widest-need expression processing, see Section 2.3 INTERVAL Arithmetic Expressions.

For syntax and semantics of the intrinsic INTERVAL constructor functions, see Section 2.8 Extending Intrinsic INTERVAL Operators.

1.3.9.2 Mixed-Mode (Type and KTPV) Expressions

If the widest-need principle is used with both KTPVs and data types, mixed-mode (type and KTPV) INTERVAL expressions can be safely and predictably evaluated. For example, in CODE EXAMPLE 1-8, the expression for Y1 in line 3 is an interval expression, because X and Y1 are INTERVAL variables.

CODE EXAMPLE 1-8   Mixed Types

math% cat ce1-8.f95 INTERVAL(16) :: X = [0.1, 0.3] INTERVAL(4) :: Y1, Y2 ! Widest-need code Y1 = X + 0.1 ! Line 3 ! Equivalent strict code Y2 = INTERVAL(X + [0.1_16], KIND=4) ! Line 4 IF (Y1 == Y2) PRINT *, "Check" END math% f95 -xia ce1-8.f95 math% a.out Check

with Widest-Need

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 mixed-type nor mixed-KTPV expressions are permitted.

The logical steps in widest-need expression processing are:

1. Scan the entire statement, including the left-hand 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.

4. Evaluate the expression.

5. Convert the result to a lower KTPV if needed to match the left-hand side's KTPV.

6. Assign the resulting value to the left-hand side.

These steps guarantee that mixed-mode INTERVAL expression processing satisfies the containment constraint and efficiently produces reasonably narrow interval results.

Mixed-mode INTERVAL expression evaluation using widest-need expression processing is supported by default with the -xia command-line 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.

1.3.10 Arithmetic Expressions

Writing arithmetic expressions that contain INTERVAL data items is simple and straightforward. Except for INTERVAL literal constants and intrinsic INTERVAL-specific functions, INTERVAL expressions look like REAL arithmetic expressions. In particular, with widest-need expression processing, REAL and INTEGER variables and literal constants can be freely used anywhere in an INTERVAL expression, such as in CODE EXAMPLE 1-9.

CODE EXAMPLE 1-9   Simple INTERVAL Expression Example

math% cat ce1-9.f95 INTEGER :: N = 3 REAL :: A = 5.0 INTERVAL :: X X = 0.1*A/N ! Line 5 PRINT *, "0.1*A/N = ", X END math% f95 -xia ce1-9.f95 math% a.out 0.1*A/N = [0.16666666666666662,0.16666666666666672]

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 WorkShop 6 f95 performs sharp data conversions. For example, the internal approximation of [0.1] is 1-ulp 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 widest-need expression processing, which is required to satisfy the containment constraint when evaluating mixed-mode INTERVAL expressions. See Section 1.3.9 Mixed-Type 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, or array. For more information on the widest-need processing mode, see Section 2.3.1 Mixed-Mode INTERVAL Expressions.

Because the interval system implemented in Sun WorkShop 6 f95 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.

---

Note – Not all mathematically equivalent INTERVAL expressions produce intervals having the same width. Additionally, it is often not possible to compute a sharp result by simply evaluating a single INTERVAL expression. In general, interval result width depends on the value of INTERVAL arguments and the form of the expression.

---

1.3.11 Interval Order Relations

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:

• certainly
• possibly
• set

For a certainly-relation to be true, every element of the operand intervals must satisfy the relation. A possibly-relation is true if it is satisfied by any elements of the operand intervals. The set-relations 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 two-letter 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:

• .CLT. certainly less than
• .PLT. possibly less than
• .SLT. set less than

See Section 2.4 Intrinsic Operators for the syntax and semantics of all INTERVAL operators.

The following program demonstrates the use of a set-equality test.

CODE EXAMPLE 1-10   Set-Equality Test

math% cat ce1-10.f95 INTERVAL :: X = [2, 3], Y = [4, 5] ! Line 1 IF(X+Y .SEQ. [6, 8]) PRINT *, "Check." ! Line 2 END math% f95 -xia ce1-10.f95 math% a.out Check.

Line 2 uses the set-equality test to verify that X+Y is equal to the interval [6, 8].

An equivalent line 2 is:

  IF(X+Y == [6, 8]) PRINT *, "Check." ! line 2 

UseCODE EXAMPLE 1-11 and CODE EXAMPLE 1-12 to explore the result of INTERVAL-specific relational operators.

CODE EXAMPLE 1-11   Interval Relational Operators  

math% cat ce1-11.f95 INTERVAL :: X, Y INTEGER :: IOS = 0 PRINT *, "Press Control/D to terminate!" WRITE(*, 1, ADVANCE='NO') READ(*, *, IOSTAT=IOS) X, Y DO WHILE (IOS >= 0) PRINT *, " For X =", X, ", and Y =", Y PRINT *, 'X .CEQ. Y, X .PEQ. Y, X .SEQ. Y =', & X .CEQ. Y, X .PEQ. Y, X .SEQ. Y PRINT *, 'X .CNE. Y, X .PNE. Y, X .SNE. Y =', & X .CNE. Y, X .PNE. Y, X .SNE. Y PRINT *, 'X .CLE. Y, X .PLE. Y, X .SLE. Y =', & X .CLE. Y, X .PLE. Y, X .SLE. Y PRINT *, 'X .CLT. Y, X .PLT. Y, X .SLT. Y =', & X .CLT. Y, X .PLT. Y, X .SLT. Y PRINT *, 'X .CGE. Y, X .PGE. Y, X .SGE. Y =', & X .CGE. Y, X .PGE. Y, X .SGE. Y PRINT *, 'X .CGT. Y, X .PGT. Y, X .SGT. Y =', & X .CGT. Y, X .PGT. Y, X .SGT. Y WRITE(*, 1, ADVANCE='NO') READ(*, *, IOSTAT=IOS) X, Y END DO 1 FORMAT( " X, Y = ") END math% f95 -xia ce1-11.f95 math% a.out Press Control/D to terminate! X, Y = [2] [3] For X = [2.0,2.0] , and Y = [3.0,3.0] X .CEQ. Y, X .PEQ. Y, X .SEQ. Y = F F F X .CNE. Y, X .PNE. Y, X .SNE. Y = T T T X .CLE. Y, X .PLE. Y, X .SLE. Y = T T T X .CLT. Y, X .PLT. Y, X .SLT. Y = T T T X .CGE. Y, X .PGE. Y, X .SGE. Y = F F F X .CGT. Y, X .PGT. Y, X .SGT. Y = F F F X, Y = 2 3 For X = [1.0,3.0] , and Y = [2.0,4.0] X .CEQ. Y, X .PEQ. Y, X .SEQ. Y = F T F X .CNE. Y, X .PNE. Y, X .SNE. Y = F T T X .CLE. Y, X .PLE. Y, X .SLE. Y = F T T X .CLT. Y, X .PLT. Y, X .SLT. Y = F T T X .CGE. Y, X .PGE. Y, X .SGE. Y = F T F X .CGT. Y, X .PGT. Y, X .SGT. Y = F T F X, Y = <Control-D>

CODE EXAMPLE 1-12 demonstrates the use of the INTERVAL-specific operators listed in TABLE 1-2.

TABLE 1-2   Interval-Specific Operators

Operator	Name	Mathematical Symbol
.IH.	Interval Hull
.IX.	Intersection
.DJ.	Disjoint
.IN.	Element
.INT.	Interior	See Section 2.7.3 Interior: (X .INT. Y).
.PSB.	Proper Subset
.PSP.	Proper Superset
.SB.	Subset
.SP.	Superset

CODE EXAMPLE 1-12   Set Operators  

math% cat ce1-12.f95 INTERVAL :: X, Y INTEGER :: IOS = 0 REAL(8) :: R = 1.5 PRINT *, "Press Control/D to terminate!" WRITE(*, 1, ADVANCE='NO') READ(*, *, IOSTAT=IOS) X, Y DO WHILE (IOS >= 0) PRINT *, " For X =", X, ", and Y =", Y PRINT *, 'X .IH. Y =', X .IH. Y PRINT *, 'X .IX. Y =', X .IX. Y PRINT *, 'X .DJ. Y =', X .DJ. Y PRINT *, 'R .IN. Y =', R .IN. Y PRINT *, 'X .INT. Y =', X .INT. Y PRINT *, 'X .PSB. Y =', X .PSB. Y PRINT *, 'X .PSP. Y =', X .PSP. Y PRINT *, 'X .SP. Y =', X .SP. Y PRINT *, 'X .SB. Y =', X .SB. Y WRITE(*, 1, ADVANCE='NO') READ(*, *, IOSTAT=IOS) X, Y END DO 1 FORMAT(" X, Y = ? ") END math% f95 -xia ce1-12.f95 math% a.out Press Control/D to terminate! X, Y = ? [1] [2] For X = [1.0,1.0] , and Y = [2.0,2.0] X .IH. Y = [1.0,2.0] X .IX. Y = [EMPTY] X .DJ. Y = T R .IN. Y = F X .INT. Y = F X .PSB. Y = F X .PSP. Y = F X .SP. Y = F X .SB. Y = F X, Y = ? [1,2] [1,3] For X = [1.0,2.0] , and Y = [1.0,3.0] X .IH. Y = [1.0,3.0] X .IX. Y = [1.0,2.0] X .DJ. Y = F R .IN. Y = T X .INT. Y = F X .PSB. Y = T X .PSP. Y = F X .SP. Y = F X .SB. Y = T X, Y = ? <Control-D>

1.3.12 Intrinsic INTERVAL-Specific Functions

A variety of intrinsic INTERVAL-specific functions are provided. See Section 2.9.4.4 Intrinsic Functions. Use CODE EXAMPLE 1-13 to explore how intrinsic INTERVAL functions behave.

CODE EXAMPLE 1-13   Intrinsic INTERVAL-Specific Functions  

math% cat ce1-13.f95 INTERVAL :: X, Y PRINT *, "Press Control/D to terminate!" WRITE(*, 1, ADVANCE='NO') READ(*, *, IOSTAT=IOS) X DO WHILE (IOS >= 0) PRINT *, " For X =", X PRINT *, 'MID(X)= ', MID(X) PRINT *, 'MIG(X)= ', MIG(X) PRINT *, 'MAG(X)= ', MAG(X) PRINT *, 'WID(X)= ', WID(X) PRINT *, 'NDIGITS(X)= ', NDIGITS(X) WRITE(*, 1, ADVANCE='NO') READ(*, *, IOSTAT=IOS) X END DO 1 FORMAT(" X = ?") END math% f95 -xia ce1-13.f95 math% a.out Press Control/D to terminate! X = ?[1.23456,1.234567890] For X = [1.2345599999999998,1.2345678900000002] MID(X)= 1.234563945 MIG(X)= 1.2345599999999998 MAG(X)= 1.2345678900000001 WID(X)= 7.890000000232433E-6 NDIGITS(X)= 6 X = ?[1,10] For X = [1.0,10.0] MID(X)= 5.5 MIG(X)= 1.0 MAG(X)= 10.0 WID(X)= 9.0 NDIGITS(X)= 1 X = ? <Control-D>

1.3.13 Interval Versions of Standard Intrinsic Functions

Every Fortran intrinsic function that accepts REAL arguments has an interval version. See Section 2.9.4.4 Intrinsic Functions. Use CODE EXAMPLE 1-14 to explore how some intrinsic functions behave.

CODE EXAMPLE 1-14   Interval Versions of Standard Intrinsic Functions  

math% cat ce1-14.f95 INTERVAL :: X, Y INTEGER :: IOS = 0 PRINT *, "Press Control/D to terminate!" WRITE(*, 1, ADVANCE='NO') READ(*, *, IOSTAT=IOS) X DO WHILE (ios >= 0) PRINT *, "For X =", X PRINT *, 'ABS(X) = ', ABS(X) PRINT *, 'LOG(X) = ', LOG(X) PRINT *, 'SQRT(X)= ', SQRT(X) PRINT *, 'SIN(X) = ', SIN(X) PRINT *, 'ACOS(X)= ', ACOS(X) WRITE(*, 1, ADVANCE='NO') READ(*, *, IOSTAT=IOS) X END DO 1 FORMAT(" X = ?") END math% f95 -xia ce1-14.f95 math% a.out Press Control/D to terminate! X = ?[1.1,1.2] For X = [1.0999999999999998,1.2000000000000002] ABS(X) = [1.0999999999999998,1.2000000000000002] LOG(X) = [0.095310179804324726,0.18232155679395479] SQRT(X)= [1.0488088481701514,1.0954451150103324] SIN(X) = [0.89120736006143519,0.93203908596722652] ACOS(X)= [EMPTY] X = ?[-0.5,0.5] For X = [-0.5,0.5] ABS(X) = [0.0E+0,0.5] LOG(X) = [-Inf,-0.69314718055994528] SQRT(X)= [0.0E+0,0.70710678118654758] SIN(X) = [-0.47942553860420307,0.47942553860420307] ACOS(X)= [1.0471975511965976,2.0943951023931958] X = ? <Control-D>

1.4 Code Development Tools

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

sun-dp-comments@Sun.COM

Include the following text in the Subject line of the email message:

WORKSHOP "6.0 mm/dd/yy" Interval

where mm/dd/yy is the month, day, and year.

1.4.1 Debugging Support

In Sun WorkShop 6, interval data types are supported by dbx to the following extent:

• The values of individual INTERVAL variables can be printed using the print command.
• The value of all INTERVAL variables can be printed using the dump command.
• New values can be assigned to INTERVAL variables using the assign command.
• Printing the value of INTERVAL expressions is not supported.
• There is no provision to visualize INTERVAL data arrays.
• All generic functionality that is not data type specific should work.

For additional details on dbx functionality, see Debugging a Program With dbx.

1.4.2 Global Program Checking

Global program checking (GPC) in Sun WorkShop 6 Fortran 95 detects one interval-specific error: INTERVAL type mismatches in user-supplied routine calls.

CODE EXAMPLE 1-15   INTERVAL Type Mismatch

math% cat ce1-15.f95 INTERVAL X X = [-1.0,+2.9] PRINT *,X CALL SUB(X) END SUBROUTINE SUB(Y) INTEGER Y(2) PRINT *,Y END math% f95 -xia ce1-15.f95 -Xlist --- See ce1-15.lst --- Global Call-Chain Considerata =============================== 1) <503> At line 4, MainPgm() calls SUB(fileline 6): MainPgm() sends argument 1 as type "Interval(16)," but SUB() expects type "Integer(4)" 2) <507> At line 4, MainPgm() calls SUB(fileline 6): MainPgm() sends argument 1 as a "Scalar," but SUB() expects a "1-D Array"

1.4.3 Interval Functionality Provided in Sun Fortran Libraries

The following libraries contain intrinsic INTERVAL routines.

TABLE 1-3   Interval Libraries

Library	Name	Needed Options
intrinsic INTERVAL array functions	libifai	None
intrinsic INTERVAL library	libsunimath	None

1.4.4 Porting Code and Binary Files

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 listed in the Interval Arithmetic README.

The representation of intervals in binary files will change as compilers supporting narrower interval systems are made available.

1.4.5 Parallelization

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.

1.5 Error Detection

The following code samples list interval-specific error messages. Each code sample includes the error message and the sample code that produced the error.

CODE EXAMPLE 1-16   Invalid Endpoints

math% cat ce1-16.f95 INTERVAL :: I = [2., 1.] END math% f95 -xia ce1-16.f95 INTERVAL :: I = [2., 1.] ^ "ce1-14.f95", Line = 1, Column = 24: ERROR: The left endpoint of the interval constant must be less than or equal to the right endpoint. f90: COMPILE TIME 0.150000 SECONDS f90: MAXIMUM FIELD LENGTH 4117346 DECIMAL WORDS f90: 2 SOURCE LINES f90: 1 ERRORS, 0 WARNINGS, 0 OTHER MESSAGES, 0 ANSI

CODE EXAMPLE 1-17   Equivalence of Intervals and Non-Intervals

math% cat ce1-17.f95 INTERVAL :: I REAL :: R EQUIVALENCE (I, R) END math% f95 -xia ce1-17.f95 EQUIVALENCE (I, R) ^ "ce1-15.f95", Line = 3, Column = 14: ERROR: Equivalence of INTERVAL object "I" and REAL object "R" is not allowed. f90: COMPILE TIME 0.160000 SECONDS f90: MAXIMUM FIELD LENGTH 4117346 DECIMAL WORDS f90: 4 SOURCE LINES f90: 1 ERRORS, 0 WARNINGS, 0 OTHER MESSAGES, 0 ANSI

CODE EXAMPLE 1-18   Equivalence of INTERVAL Objects with Different KTPVs

math% cat ce1-18.f95 INTERVAL(4) :: I1 INTERVAL(8) :: I2 EQUIVALENCE (I1, I2) END math% f95 -xia ce1-18.f95 EQUIVALENCE (I1, I2) ^ "ce1-16.f95", Line = 3, Column = 14: ERROR: Equivalence of the interval objects "I1" and "I2" with the different kind type parameters is not allowed. f90: COMPILE TIME 0.190000 SECONDS f90: MAXIMUM FIELD LENGTH 4117346 DECIMAL WORDS f90: 4 SOURCE LINES f90: 1 ERRORS, 0 WARNINGS, 0 OTHER MESSAGES, 0 ANSI

CODE EXAMPLE 1-19   Assigning a REAL Expression to an INTERVAL Variable in Strict Mode

math% cat ce1-19.f95 INTERVAL :: X REAL :: R X = R END math% f95 -xia=strict ce1-19.f95 X = R ^ "ce1-17.f95", Line = 3, Column = 3: ERROR: Assignment of a REAL expression to a INTERVAL variable is not allowed. f90: COMPILE TIME 0.350000 SECONDS f90: MAXIMUM FIELD LENGTH 4117346 DECIMAL WORDS f90: 4 SOURCE LINES f90: 1 ERRORS, 0 WARNINGS, 0 OTHER MESSAGES, 0 ANSI

CODE EXAMPLE 1-20   Assigning an INTERVAL Expression to INTERVAL Variable in Strict Mode

math% cat ce1-20.f95 INTERVAL :: X INTERVAL(16) :: y X = Y END math% f95 -xia=strict ce1-20.f95 X = Y ^ "ce1-18.f95", Line = 3, Column = 3: ERROR: Assignment of an interval expression to an interval variable is not allowed when they have different kind type parameter values. f90: COMPILE TIME 0.170000 SECONDS f90: MAXIMUM FIELD LENGTH 4117346 DECIMAL WORDS f90: 4 SOURCE LINES f90: 1 ERRORS, 0 WARNINGS, 0 OTHER MESSAGES, 0 ANSI

1.5.1 Known Containment Failures

Whenever an interval containment failure can occur, a compile-time warning should be issued. An integer expression outside the scope of widest-need expression processing is the only known situation in which such a warning is necessary.

1.5.1.1 Integer Overflow

Numerical inaccuracies are normally associated with REAL rather than INTEGER expressions. In one respect, INTEGER expressions are more dangerous 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 floating-point 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 possibly-opposite-signed 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.

The following example shows what can happen if the scope of widest-need expression processing is not extended to all intrinsic INTEGER operations and functions, including the ** operation with an INTEGER exponent.

CODE EXAMPLE 1-21   INTEGER Overflow Containment Failure

math% cat ce1-21.f95 INTERVAL :: X = [2], Y = [2] INTEGER :: I = HUGE(0) X = X**(I+1) Y = Y*(Y**I) IF(X .DJ. Y) PRINT *, "X and Y are disjoint." END math% f95 -xia ce1-21.f95 math% a.out X and Y are disjoint.

This code demonstrates a silent containment failure. It is a known error because the scope of widest need expression processing does not presently extend to the integer exponent of the ** operation. For information on the power operator, see Section 2.5 Power Operators X**N and X**Y.

---

Caution – This error has not been fixed in the Sun WorkShop 6 Fortran 95 release, and no warning messages are issued.

---