(j3.2006) Interval arithmetic

Van Snyder Van.Snyder
Wed Mar 18 23:42:09 EDT 2009


On Wed, 2009-03-18 at 18:47 -0700, Malcolm Cohen wrote:
> You have yet to come up with a single concrete example of code.

  pure real function add_up ( x, y )
    use, intrinsic :: ieee_arithmetic
    real, intent(in) :: x, y
    call ieee_set_rounding_mode ( ieee_up )
    add_up = x + y                             ! Instant one
  end function add_up

The answer to question two in interp f95/000104 (which hasn't been
refuted in f08) allows the processor to compute x+y in double precision
at instant one, notwithstanding f08:7.1.9, so rounding it to double
precision does not violate Clause 14 or 754.  If you want to quibble "no
processor would do that with today's Intel chips" just change the kinds
of the dummy arguments to double precision.  The standard doesn't
prohibit inter-statement optimization, and users want it, so processors
do it.  The standard doesn't prohibit returning function results in
registers, so many processors do it.  Therefore, the processor is not
compelled by the standard to store x+y into a default real variable.
The answer to question four allows the processor to return the function
result as double precision even though the text says default real, so
there's no requirement to round it to default real, and no violation of
Clause 14 or P754.  There's nothing in C14 that says "respect 7.1.9;
ignore f95/000104."

etc. similar with defined operators such as .addup. .subdown. ...

  real :: x

  ...

  x = (y .addup. z) .mulup. w                  ! Instant two

  ! The standard doesn't prohibit passing arguments to functions
  ! in registers, and many do it, so a processor might find it
  ! convenient to copy the result of (y.addup.z) to the register
  ! it will use for the first argument of .mulup.  After all, that
  ! will get rid of a store from double precision to single in
  ! memory (or maybe cache) and a load from single to double.

  ...

  ! Code wanders in the wilderness for a while.

  ! More computations involving (y.addup.z), which the processor
  ! discovers is conveniently still in a register as double precision.
  ! So far, no violations of f95, f2003, 15580 or 754.

  ...

  call ieee_set_rounding_mode ( ieee_down )
                                              ! Instant three
  ! Optimizer runs out of registers, so it invisibly does

  temp = register containing y.addup.z ! rounds from double to single

  ! There's nothing here that mentions temp or the register
  ! containing y.addup.z, so the program can't demand

  call ieee_get_rounding_mode ( foo )
  call ieee_set_rounding_mode ( ieee_up )
  temp = the anonymous register containing y.addup.z
  call ieee_set_rounding_mode ( foo )

  ! and the processor doesn't know enough any more to do it on its own

  a = b + c ! with downward rounding, which needed the register
            ! containing y.addup.z

  ...

  q = (y .addup. z) .divup. denom           ! y.addup.z is in temp here
                                            ! but because of downward
                                            ! rounding at instant three,
                                            ! the "correct" result is
                                            ! greater than the value in
                                            ! temp.

  ...

  print *, "The bridge won't fall down; fortunately Bill wasn't on it."

There's nothing in 14.4, 14.11.21, or IEEE P754 that addresses the
separation of instants one and three. The standard doesn't prohibit
inter-statement optimization, so most processors assume it's allowed,
and most customers want it.  At instant one, X+Y is computed with upward
rounding, giving a double-precision result.  The answer to question two
in interp f95/000104 allows this.  The rounding done at this instant
conforms to 14 and p754.

At instant two the value of (y.addup.z) is in a register represented as
a double-precision value, even though the result of add_up is default
real.  f95/000104 allows this.

At instant three, the processor rounds the value of (y.addup.z) from
double precision to default real, conforming to p754.  Unfortunately,
there's no program text here, so clause 14 has nothing to say.

Sure, the processor shouldn't have carried the double precision result
across the change of rounding mode at "end function add_up" (and the NAG
compiler almost surely doesn't) but there's nothing in the standard that
prohibits it, and interp f95/000104 explicitly allows a double-precision
result, with no caveats about not doing it if the rounding mode is to
change.

Vendors are much less responsive to quality of implementation complaints
(except from NASTRAN and SPEC) than to violations of the standard, so
correcting the "bug" (which isn't a violation of the standard) of
carrying a double-precision result that will ultimate be a default real
result across a change in rounding mode has very low priority.

> Some packages existed and worked before VOLATILE existed, so no, when
> I actually looked they didn't.

Well look again for yourself if you don't believe Christian Keil or any
of the other members of P1788.  The packages existed, but were found
_not_ to work.  The cause of the problem was found to be exactly as I
explained.  That's why implementors started shoving results through
volatile variables.

> There certainly is a place for volatile when intermediate calculations
> *inside the interval package itself* have to be correctly rounded, at
> least on x87.

This shouldn't be necessary, and it took a long time to discover that
omitting it was a source of error.  We really shouldn't expect everybody
to understand the subtle interaction between f95/000104, rounding mode,
and volatile.  The standard should require it to work correctly as the
text + 7.1.9 appear to require.  I think interp f95/000104 means it
doesn't.  If it does, please tell me where.  Just shouting "P754 +
Clause 14" doesn't convince me.

Is it too much to ask for at least a note in 14.4 that advises a
processor not to carry a result of higher precision than implied by the
kind demanded by 7.1.9 across a change of rounding mode?  If we really
want to enshrine f95/000104 in the standard, it should be a normative
requirement.  It would be better not to be so glib about ignoring the
requirements of 7.1.9.  Then the problem would go away.





More information about the J3 mailing list