(j3.2006) Interval arithmetic

[Malcolm Cohen]

Van Snyder wrote:
> 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.
>   
If true, then it is relevant and incidentally realistic as in real.
>   
>> In that
>>     
>
> excessively trivial contrived
>   
Not trivial, not contrived.  Realistic.  *You* are the one making claims 
of impossibility, *YOU* are the one who needs to prove your case!  You 
have yet to come up with a single concrete example of code.
>   
>> 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.
>   
No they don't.  Each interval operation is a defined operation i.e. a 
function reference.  Each function reference works *EXACTLY LIKE MY 
"CONTRIVED" EXAMPLE*.  Sure, I chose the simplest example.  So what?  
It's a realistic concrete example nonetheless.

Seeing as how the rounding mode *ONLY CHANGES DUE TO EXECUTION OF CALL 
STATEMENTS* it is obviously IMPOSSIBLE for a single expression to do 
rounding mode stuff.  Especially since on return from any function 
referenced by the expression the rounding mode is reset.
>   
>> 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.
Some packages existed and worked before VOLATILE existed, so no, when I 
actually looked they didn't.

There certainly is a place for volatile when intermediate calculations 
*inside the interval package itself* have to be correctly rounded, at 
least on x87.  I agree with Bill that it doesn't make much sense on 
anything else that is current.  I also think that volatile is much less 
necessary, even on x87, if the compiler understands what 
IEEE_SET_ROUNDING actually does; not knowing this falls into the realms 
of compiler bugginess, probably to be expected in early 
implementations.  And if the compiler doesn't even offer the IEEE 
modules (and many probably didn't when some of these codes were written) 
and the rounding mode is being twiddled by some unknown procedure, it's 
hardly surprising that they didn't work reliably without volatile.
>   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.
>   
The standard already requires the correct result.  SPEC and its ilk 
trump the standard though.  Maybe we need a standard-mandated option:
   
-do_what_the_standard_says_I_want_the_right_result_I_am_not_a_SPEC_benchmark
> If a variable carried an attribute
Why do you hate programming languages (or maybe just programmers)?  
That's the most complicated and difficult-to-understand way of doing 
these things.
>  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.
>   
This is just so wrong.  For a start, as Bill said the big performance 
hit on most current hardware is doing the rounding mode changes.  Doing 
a round-trip to L1 cache is not going to be totally insignificant, but 
it's not the super-killer you might imagine.  (A local variable with 
volatile usually doesn't need any special treatment that a dummy 
argument or module variable or common variable might require - just the 
ordinary store/load does the job and even that could in theory be 
omitted when there is no change in format.)  Certainly not even close to 
being as significant as doing at least twice as many arithmetic ops, 
interspersed with expensive rounding mode changes and sign tests (yes, 
more pipeline bubbles).  And unless the processor does sufficiently 
agressive function inlining, procedure call overhead to boot.

Secondly, if the hardware has static rounding modes built-in to the 
instructions, instead of a single dynamic rounding mode, there's no 
problem moving the computations around.  The DEC Alpha did this... 
unfortunately, most of the world decided it preferred x86 ISA 
compatibility to high performance.

Certainly the hopes and desires of the 754r committee were that static 
rounding modes were to be encouraged for performance reasons, even 
though they couldn't get rid of the dynamic ones for compatibility 
reasons.  Whether such designs will escape purpose-built one-off 
machines and become commodity hardware remains to be seen, but I don't 
expect it for another decade at least.

> 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.
>   
Well, if that's what they say they're just wrong.  In a Fortran 
calculation with rounding mode changes, there are NO intermediate 
results that can be propagated across any rounding mode change.  That's 
a straightforward consequence of the Fortran+754 standards, and I find 
it is written up quite plainly in c14.  That some grad student couldn't 
work things out or was working prior to F2003 (not exactly a lot of 
F2003 compilers around even today, and we still have F95 compilers that 
lack the IEEE modules) or ran into a compiler bug doesn't change that.  
(Hey, you don't think all those high-flying academics do their own grunt 
work, do you?)

And since when does being an interval expert make someone a hardware 
expert or a compiler expert?  You are certainly arguing from 
inappropriate authority here.  I don't care whose "idle speculation" it 
is, but calling the idea that some wording change to the standard will 
remove compiler bugs, badly written user packages, and inefficient 
hardware "idle speculation" is to praise it to the skies!

Finally, the C standard with the right #pragma definitely allows none of 
the intermediate result hanky-panky you are claiming is the problem, and 
we can call C functions from Fortran quite easily.  So if these guys 
didn't already do their implementation in C, or they found that it ran 
just as slow (or even slower) in C, that shows that the wild goose they 
are chasing is the wrong red herring.

Cheers,
-- 
...........................Malcolm Cohen, Nihon NAG, Tokyo.