(j3.2006) Interval arithmetic

[Van Snyder]

On Tue, 2009-03-17 at 20:32 -0700, Malcolm Cohen wrote:
> These hammers are *WAY* too big for this alleged "problem", which does
> not (IMO) even arise given the current wording of F2003 in particular 
> with the IEEE modules.  Those modules *already* give the precise
> answer, viz the one required by IEEE 754, for (qua the right
> declarations)
>    CALL ieee_set_rounding(ieee_up)
>     X%HI = A%HI + B%HI
>    CALL ieee_set_rounding(ieee_down)
>    X%LO = A%LO + B%LO
> which is the kind of code we would be seeing for the user-defined 
> implementation of interval-arithmetic "+".

True but not relevant.

> In that

excessively trivial contrived

> example there is no intermediate result for interp F95/whatever to be
> applied to,...

but in a real computation, such as interval divide followed by stacking
of exterior intervals resulting from denominator intervals containing
zero, real problems do arise.

> All us implementors know we're not supposed to move computations
> across a statement that changes the rounding mode.

If you actually look into an interval arithmetic package, you'll see
that they store their results into volatile variables.  This is
correspondence from Christian Keil c.keil at tu-harburg.de on the P1788
mailing list:

! That's exactly what interval packages (PROFIL for example) do right now
! to enforce (among others) rounding of intermediate results... Hopefully
! the standard will allow us sometime to get rid of these.

So, no, although NAG probably does know not to move intermediate results
across statements that change the rounding mode, there are evidently
plenty who don't, and correcting this "quality of implementation" defect
is usually quite low on their priority lists.  Apparently, unless
NASTRAN or SPEC say it's a problem, the standard(s) needs to say it's
prohibited.

If a variable carried an attribute specifying the direction of rounding
to apply when a higher-precision intermediate result is rounded to its
canonical representation, as would be the case with the components of a
type the compiler understood, or an explicitly specified attribute, it
wouldn't be necessary to prohibit moving high-precision intermediate
results across rounding-mode changes.

I opposed intrinsic interval arithmetic for f2003 partly because I
thought it could be done by user procedures with IEEE support.  The
people who have actually tried to implement it tell me it can be done
only by taking extreme defensive measures that cripple performance,
which wouldn't be the case if a compiler understood, or could be told,
what to do.

I got reinvolved with interval arithmetic when at an IFIP WG 2.5
(numerical software) meeting there was discussion of a proposal to send
a letter to IEEE P754 to beg them to add interval arithmetic.  I said
that was naive, that somebody actually had to write the words to put
into the standard, and that sending a one-page letter to beg P754 to do
it would almost certainly not be fruitful, especially if there was
nobody on P754 with both expertise and interest.  I got appointed, along
with Bo Einarsson, to cast a paper by Ulrich Kulisch into IEEE
standardese (I'll send this to anybody who wants it).  That draft, which
P754 decided not to incorporate into their standard at the last minute,
resulted in the formation of P1788.  Dan Zuras (a standards maven with
experience in SCSI and 754), our old friend Baker Kearfott, John Pryce,
and Nathalie Revol are leaders of P1788, with active participation by
Ulrich Kulisch, Siegfried Rump, Christian Keil, Arnold Neumaier, Nate
Hayes, Gilles Chabert, Michael Hack, Wolff von Gudenberg, Michael
Schulte, Chenyi Hu, David Lester, Gabriel Dos Rios, Guillaume Melquiond,
Jim Demmel, Vladik Kreinovich, ....

My concerns about lack of control of rounding of intermediate results
arises from correspondence of that committee of experienced experts, who
have had to work around this problem.  It's not idle speculation.