J3/97-164
Date: April 24, 1997
To: J3
From: Dick Hendrickson
Subject: Preliminary Module Proposal for Interval Arithmetic
Reference: x3j3/97-105, x3j3/97-141, x3j3/97-136, x3j3/97-126, x3j3/97-113
This is a proposal to provide Interval Arithmetic as an optional
intrinsic module in Fortran 2000. I believe it includes all of the
functionality in the reference papers, although not all of the details are
completely spelled out. I did not intend to leave out anything
significant.
The major addition is a new function "accuracy" which has 2 or 3
arguments and returns the accuracy when the first argument is applied
to the second and third. For example accuracy("+",2.0,3.0) gives the
machine accuracy when 2 and 3 are added together. It is obviously
hardware specific but I believe we can
provide a reasonable overapproximation that is portable. On specific
machines vendors can provide better approximations.
J3 provides both the definitions of the necessary module interfaces and
functions as well as a defined source code for the module. The intent
is that a vendor can "support" the module merely by making the user
find the defined source code somewhere on the net and compiling it.
Vendors have the option to provide direct optimized support in the
compiler, effectively ignoring all or some of the defined module code
(but not the interface definitions).
Since most new machines are becoming IEEE like, at least in the
sense of having directed roundings, many of the proposed functions
have an exact IEEE analog. We probably want to design them to fit
very well with an IEEE machine and not too badly on a nonIEEE
machine.
Rationale
Some people regard Interval Arithmetic is an important feature, but
some do not believe there is enough demand to add the complexity to
the language. It's my belief that an INTRINSIC module lets us have
the advantages, for those who believe there is an advantage, without
the high costs, for those who believe there is a high cost.
The mental model I have is based on the way other things are
currently implemented in existing practice.
Consider the SQRT function. On some machines it executes quickly
and accurately as a single inlined instruction, on some it's a few
inlined instructions, on some it's implemented as a call to a hand
crafted assembly language routine using a special calling sequence
with special conventions about register usage, on some it's a call to a
hand crafted assembly language routine using normal calling
conventions, on some it's a call to an optimized C routine, on some it's
a call to an unoptimized routine. Some processors "execute"
SQRT(9.0) in exactly 0 microseconds at "run-time", some don't. All
of these are standard conforming. The difference is the importance the
hardware/software designers attach to speed and accuracy in the SQRT
routine.
Consider the BLAS, especially the SAXPY routine. It's not standard
Fortran, but it's a de facto intrinsic routine for vendors who compete in
Dongarra's LINPACK derby. It's available in standard Fortran.
Vendors do the same thing with SAXPY that they do with SQRT (at
least as far as optimization strategy goes, most of them actually
compute different answers). Some have highly optimized assembly
language versions available, others generate perfect inline code, others
make the user find the code on the net and compile it along with the
rest of his code. In the official LINPACK test some vendors inline the
vanilla Fortran, propagate the arguments (noting that almost all of the
interesting ones are 1), eliminate the resulting dead code, recognize
the remaining loop as a SAXPY and either generate perfect inline
code or call an optimized external assemble language routine. All
standard conforming and all choices that vendors make in response to
market place pressure.
Consider the varying string module. As far as I know it's only been
"implemented" as Fortran source code available over the web. No
vendors have felt enough pressure to do an optimized version. Not
surprising since character operations tend to be inefficient on Fortran
oriented machines and character operations are usually not the time-
dominant part of important Fortran programs. But, let a real character
customer with $100,000,000 in his pockets walk in and watch what
happens.
I believe Interval Arithmetic as an intrinsic module can fall into this
pattern.
Module approach
Interval becomes a defined type with lots of overloaded operators and
functions that operate on it. The type has private components.
J3 standardizes the definition and interfaces for interval arithmetic. It
also provides a module which implements this definition BUT
carefully words the standard to say that the module code need not be
used in the actual implementation. It must be clear that a standard
conforming implementation can be less sharp or more sharp in the
intervals it calculates than it would be if it merely compiled the
module. Vendors must be free to treat the INTERFACES as intrinsic
and get different computational results. Sort of like the definition of
ANINT (or whatever it was) that someone from SUN asked about a
year or two ago. This is different from how we normally standardize
things, but we are going to have to do something like this even if we
make intervals a "first class" member of the standard. There just isn't
enough common language practice to say "go and do intervals",
whereas everybody understand what "go and multiply" means.
To be successful as a module we must do this without adding any new
syntax to the language. In effect J3 must agree at the outset to add
whatever (with luck, small) syntaxs are needed to make intervals
work. The additions will have to be reasonably general and likely to
be useful for other things or the module approach will not work.
The following areas need to be considered.
I/O
Nothing much needs to be done, the work the I/O group is doing is
sufficient. The drawback is that the interval format conversions will
have names like "DT?", rather than "?". I, personally, don't think that
is a major problem, all of the good single letters have already been
used for other things anyhow. This can be integrated with the IEEE
I/O inquiry functions for the module implementation. If the processor
rn time library does not do accurate I/O conversions (perhaps in the
IEEE sense?) then the module routines will have to generate or read
character strings and do the conversions by hand.
Overloaded operators and rank
The module will define overload functions for all of the operators and
functions that make sense. For the binary operators there are 1 + 7 + 7
+ 7 = 22 function definitions (to cover all combinations of scalar and
array rank). The functions will all be identical except for rank. It
would be nice to have some syntactic help if specifying them only
once.
New (comparison) operators and their precedence
There are so many new comparison operators needed (to cover partial
overlap, complete containment, complete disjointedness...) that I don't
think it is reasonable to add a dozen or so new dot operators.
Expressions using more than one operator are likely to be so complex
that the user will use parenthesis to group things. In current practice
most expressions that have both an .AND. and .OR. operator have
parens to separate the terms because humans can't remember which
comes first. No need for new operators and new precedence. Use
functions (ala LLT) for the new comparisons.
Interval function
I believe the INTERVAL function as described in Baker's paper 105 is
fine for converting non-interval variables into interval variables. It
also works for constants, but the notation isn't perfect. A function is a
bulky way to express a constant.
In my view there are 2 kinds of numbers. Those like pi or the speed of
light that are know to more than machine precision and those like the
probability of rain which are not. The interval function as currently
proposed covers the second class. The user can specify an uncertainty
when he converts from non-interval to interval. The problem with
writing something like INTERVAL(3.141592653) is that the rules of
Fortran define the argument as a single precision number and then the
precision gets lost. I propose also allowing a character string as the
first argument. INTERVAL("3.141592653") passes the string intact
to the module routine which can decode it to whatever precision is
needed. Naturally, a native mode compiler could do all of the
conversion at compile-time. If the INTERVAL function is invoked
with literal real constants accuracy will probably be lost unless the
constant is exactly representable.
As an alternative we could consider treating "(<" and ">)" as
overloadable intrinsic operators. The syntax would be to pass their
argument(s) AS CHARACTER STRINGS to the resolved function.
This would allow (<3.141592653>) or (<1.0, 3.0>) as ways to write an
interval constant. The major drawback seems to be that usually only
one overload could be visible in any scoping unit, since all of the
useful arguments types are likely to be real.
Other intrinsic functions
All of the other interval functions from paper 105 can be implemented
as module functions.
whats_my_interval
To let the user be able to find out what's going on I propose an inquiry
function that defines the level of processor support, somewhat like the
IEEE inquiry functions. It ranges from "no support", the user must
arrange to have the defined Fortran source code compiled and linked
to "native support" where the compiler is "truely aware" of just exactly
what is going on in interval arithmetic. My aim would be to return a
number that defines the level, not a whole series of logical values that
define particular types of support. Safe portable code would test the
return value and abort if it wasn't good enough. As a first try:
-1 no direct processor support, everything comes from the defined
module. The module is compiled "as is" with no tuning for the
particular processor. Normal optimizations are performed just as for
any other compiled code.
0 no direct processor support, everything comes from the defined
module but with user supplied ACCURACY and FUZZ functions
tuned to the particular machine, otherwise the module is compiled as
above.
1 minimal processor support, everything comes from the defined
module but with processor supplied ACCURACY and FUZZ functions
tuned to the hardware. The module is compiled as above.
For the remaining values the processor has some level of awareness.
2 accurate I/O conversions
3 most math functions are accurate
4 all math functions are accurate
5 native vendor implementation, not using the defined module
functions, for almost everything
6 native vendor implementation.
The intent is that higher numbers require the lower ones to be true.
As a refinement we'd return named_constants of some module specific
type rather than numbers. I used numbers now to avoid spelling
arguments (we all agree how to spell "4" :-) )
special cases, x**2, etc.
There seems to be a consensus within the interval group that X**2 is
not the same operation as X*X. We will need to make sure that the
chapter 7 rules don't allow the general algebraic operations on interval
arguments. I think this happens by magic because of the parsing
overloading rules, but we need to make sure. Note that the module
interface can define allowed transformations and a native mode
compiler can perform them, (for example X**3 is the same as
X*X*X) it's just that some of them are different from what's allowed
with reals.
use, intrinsic vs use
I'd like to find a way to allow a knowledgeable user or any system
administrator be able to augment the defined module with their own
versions of some or all of the routines. For example on a machine that
is known to have an IEEE SQRT instruction the defined SQRT routine
could probably be both sped up and sharpened up. As currently
defined I think the USE,.INTRINSIC statement forces use of the
intrinsic module and the USE statement first looks for a user supplied
version of the module. To be usefully portable we want the USE to
select the INTRINSIC version if it's compiled in native mode, then to
look for a user supplied version, then to look for the defined module
version.
Round up/down functions
The module will have to define a whole set of rounding functions.
Things like
round_up_add(x1,x2,x3,...)
which forms the sum x1 + x2 + x3... using directed rounding away
from zero. On an IEEE machine this would translate into a simple
series of adds.
Accuracy
We need a function that gives an estimate of the accuracy of an
operation. It is obviously hardware specific. The module will have to
contain a version that does a portable estimate. The user could
compute (probably in a separate run) the accuracy of various
operations for various operands (that's what programs like
PARANOIA do) and feed that information into the module.
Alternatively, a vendor who provides minimal support (level 1) could
provide an accuracy function and nothing else. The function would
have a character argument and one or 2 real arguments.
ACCURACY ("+", 3.0, 4.0) gives the accuracy of the + when applied
to 3 and 4. A reasonable version might be max (nearest(3.0,1.0)-3.0,
nearest(4.0,1.0), nearest(3.0+4.0)-(3.0+4.0), epsilon(1.0)). On an
IEEE machine this function is easier to compute.
ACCURACY ("SQRT", 9.0) gives the accuracy of the SQRT routine
for arguments near 9.0. Vendors probably have some sort of mini-max
estimate of this number.
ACCURACY ("F20.10, X) gives the accuracy of an I/O conversion for
the value of X.
This might also be a class of numeric inquiry functions we should add
to Fortran 2000 as a general aid.
Fuzz
A safety factor number used in defining results. Processor dependent
and user definable for the pessimist.
Combining all of this the interval addition routine could look like:
function add_two_intervals (x,y) result (z)
type (interval), intent(in) :: x,y
type interval :: z
z%high = round_up_add(x%high, y%high) + accuracy("+", &
x%high, y%high)*(1.0+fuzz())
z%low = round_down_add(x%low, y%low) - accuracy("+", x%low, &
y%low)*(1.0+fuzz())
! also need code to detect overflow and do halting if needed.
! I don't know of a portable way to detect overflows, this is overkill
If (x%high > huge(1.0)/2.1 .or. y%high > huge(1.0)/2.1) call
interval_overflow_routine()
end function add_two_intervals