/internet

/internet

[tprince]

The question of possible over/underflow of intermediate results is more relevant for floating point than
for integer multiplication.  However, I believe most people consider that the programmer is obligated
to use parentheses to enforce safe grouping, if there is a known such grouping. If the program was
originally tested on Intel, this is not likely to have been considered!  Reassociation is commonly
required for adequate performance of a series of additions, and there of course accuracy is more
likely to be an issue than over/underflow.

I try to make a practice of using appropriate parentheses in an expression such as

(a-b) + (c-d) + (e-f)

where that may improve accuracy (given that the differences are known to be relatively small) as well
as suggesting an effective pipelined grouping.  Possibly this would deserve an FAQ treatment if
egcs/gnu begins to do reassociations.

Treatment of reassociation is by no means uniform.  I just accepted Arnaud Desitter's suggestion of a
usage of Kahan summation to improve the portability of the check-summing in Livermore Fortran
Kernel.  It turns out that the Irix (MipsPro7.2) compiler requires the option -OPT:fold_reassociate=OFF
in conjunction with optimization, in order for this to work.   They treat neglect of parentheses and order
of assignments, as well as changing the sense of comparison, as a normal consequence of full
optimization.  Of course, in that case, the intended data dependencies require 4 times as many cycles
as the "optimized" code, but that's still faster than REAL*16 on some machines.

This optimization doesn't matter on Intel in extended precision mode, but one of the compilers I tested
apparently sets double precision mode.  I tried declaring "real(selected_real_kind(18)) sum" to make
a strong suggestion that the compiler should use extended precision or REAL*16, but certain
compilers rejected this entirely, and g77 has not adopted this syntax.
____________________________________________

>>For FP, we would like the ability to reassociate some expressions.  Take
(a * b * c * d) * e

>>Right now we'll genrate

t1 = a * b;
t2 = t1 * c;
t3 = t2 * d;
t4 = t3 * e;

>>Note the dependency of each insn on the previous insn.  This can be a major
performance penalty -- especially on targets which have dual FP units or where
a fpmul isn't incredibly fast (data dependency stalls at each step).

t1 = a * b;
t2 = c * d;
t3 = t1 * t2;
t4 = t3 * e;


>>Is a much better (and safe as far as I know) sequence.  The first two insns
are totally independent, which at the minimum reduces one of the 3 stall
conditions due to data dependency.  For a target with a pipelined FPU or
dual FPUs the second sequence sequence will be significantly faster.


>>For integer, we need to know where the parens are to preserve integer overflow
semantics in languages like Ada for similar transformations.


jeff<<

-
Dr. Timothy C. Prince
Consulting Engineer
Solar Turbines, a Caterpillar Company
alternate e-mail: tprince@computer.org

           To:                                              INTERNET - IBMMAIL
                                                            N0520484 - IBMMAIL

Re: /internet

[Stephen L Moshier]

There has been general agreement that gcc is not supposed to make
optimizations that could change the value of an expression.
Associative law optimizations of floating-point expressions certainly
can change the value and they should continue to be disallowed.  Users
should not have to fear that gcc is going to take arbitrary liberties
and rewrite their programs for them without very explicit permission.

Re: /internet

[Joe Buck]

> There has been general agreement that gcc is not supposed to make
> optimizations that could change the value of an expression.
> Associative law optimizations of floating-point expressions certainly
> can change the value and they should continue to be disallowed.  Users
> should not have to fear that gcc is going to take arbitrary liberties
> and rewrite their programs for them without very explicit permission.

If I understand correctly, Jeff's proposal was to respect ordering when
the user provides parentheses, but to leave the compiler free to do
reordering when the user does not provide parentheses.

That is, for

	y = a * b * c * d;

the compiler could do the multiplications in any order, but for

	y = a * (b * (c * d));

the compiler is forced to follow the specified order.

I believe that this is consistent with what the relevant language standards
specify.

Re: /internet

[Jeffrey A Law]

  In message < 199812160320.TAA11586@atrus.synopsys.com >you write:
  > If I understand correctly, Jeff's proposal was to respect ordering when
  > the user provides parentheses, but to leave the compiler free to do
  > reordering when the user does not provide parentheses.
  > 
  > That is, for
  > 
  > 	y = a * b * c * d;
  > 
  > the compiler could do the multiplications in any order, but for
  > 
  > 	y = a * (b * (c * d));
  > 
  > the compiler is forced to follow the specified order.
  > 
  > I believe that this is consistent with what the relevant language standards
  > specify.
You understand my suggestions correctly.

But I'm neither a language lawyer, nor an ieee/fp expert, so you have to take
my suggestions with some care.  I certainly wouldn't actually make such changes
without being 100% they're safe according to the language specs and don't
violate the principle of least suprise too often.  For that I would have to
lean on folks in on this list to guide me through the pitfalls of each
language and ieee/fp.

jeff

Re: /internet

[Matthias Urlichs]

Hi,

Joe Buck:
> That is, for
> 
> 	y = a * b * c * d;
> 
> the compiler could do the multiplications in any order, but for

No it cannot. C specifies whether operations are left- or
right-associative, meaning it's well-defined which is executed first.
Think "-" vs. "=". ;-)

-- 
Matthias Urlichs  |  noris network GmbH   |   smurf@noris.de  |  ICQ: 20193661

Re: /internet

[Stephen L Moshier]

> If I understand correctly, Jeff's proposal was to respect ordering when
> the user provides parentheses, but to leave the compiler free to do
> reordering when the user does not provide parentheses.

That would lead immediately to getting different answers on different
machines, or even different answers after slight code changes on the
same machine.  The order of evaluation would be changed to suit each
particular machine's scheduling algorithm.  It should not be permitted
as the normal behavior.

Since reordering can produce infinitely large numerical changes,
I have argued against allowing associative law optimization even
in fast-math mode.  It ought to be well isolated with its own
command switch, -fvectorize or some such.

Re: /internet

[Tim Hollebeek]

Jeffrey A Law writes ...
> 
> You understand my suggestions correctly.
> 
> But I'm neither a language lawyer, nor an ieee/fp expert, so you have to take
> my suggestions with some care.  I certainly wouldn't actually make such changes
> without being 100% they're safe according to the language specs and don't
> violate the principle of least suprise too often.  For that I would have to
> lean on folks in on this list to guide me through the pitfalls of each
> language and ieee/fp.

I believe the example you were talking about yesterday violates IEEE,
since the reordering you suggested can change the value of the result
(by causing an intermediate (c*d) to overflow, when left to right
evaluation would not cause an overflow):

#include <math.h>

int main() {
    int z = 400;

    double a = scalb(1, -z);
    double b = scalb(1, -z);
    double c = scalb(1, z);
    double d = scalb(1, 2*z);

    double t1, t2, t3;

    t1 = a*b;
    t2 = t1*c;
    t3 = t2*d;

    printf("--> %g\n", t3);

    t1 = a * b;
    t2 = c * d;
    t3 = t1 * t2;

    printf("--> %g\n", t3);
}

---output---
--> 2.58225e+120
--> inf

---------------------------------------------------------------------------
Tim Hollebeek                           | "Everything above is a true
email: tim@wfn-shop.princeton.edu       |  statement, for sufficiently
URL: http://wfn-shop.princeton.edu/~tim |  false values of true."

Re: /internet

[Joe Buck]

> I believe the example you were talking about yesterday violates IEEE,
> since the reordering you suggested can change the value of the result
> (by causing an intermediate (c*d) to overflow, when left to right
> evaluation would not cause an overflow):

Does IEEE require left-to-right evaluation of y = a * b * c * d;
when there are no parentheses present?

Re: /internet

[Craig Burley]

>> 	y = a * b * c * d;
>> 
>> the compiler could do the multiplications in any order, but for
>
>No it cannot. C specifies whether operations are left- or
>right-associative, meaning it's well-defined which is executed first.
>Think "-" vs. "=". ;-)

Read the pertinent standards, then re-think.  You're wrong.

        tq vm, (burley)

Re: /internet

[Craig Burley]

>There has been general agreement that gcc is not supposed to make
>optimizations that could change the value of an expression.
>Associative law optimizations of floating-point expressions certainly
>can change the value and they should continue to be disallowed.  Users
>should not have to fear that gcc is going to take arbitrary liberties
>and rewrite their programs for them without very explicit permission.

Even if this was true -- and IMO it isn't -- it doesn't apply to well-
designed numerical-processing languages like Fortran, in which it
is quite clear that A*B*C*D can be evaluated as (D*A)*(C*B), for
example.

Therefore, we should make the *back end* work properly in a language-
independent sense.  If people want the gcc to produce slower code
than g77, that's okay -- they can have gcc explicitly group the
operations for consumption by the back end, something Fortran doesn't
generally need to do.

        tq vm, (burley)

Re: /internet

[Stephen L Moshier]

> Does IEEE require left-to-right evaluation of y = a * b * c * d;
> when there are no parentheses present?

I think that IEEE does not actually say anything about it.  C9X, however,
does contain some remarks such as the following, which might have come
from NCEG.


       [#14] EXAMPLE 5 Rearrangement for floating-point expressions
       is often restricted because of limitations in  precision  as
       well  as  range.   The implementation cannot generally apply
       the  mathematical  associative   rules   for   addition   or
       multiplication,   nor  the  distributive  rule,  because  of
       roundoff  error,  even  in  the  absence  of  overflow   and
       underflow.    Likewise,   implementations  cannot  generally
       replace decimal constants in order to rearrange expressions.
       In  the  following  fragment,  rearrangements  suggested  by
       mathematical rules for real numbers are often not valid (see
       F.8).

               double x, y, z;
               /* ... */
               x = (x * y) * z;  // not equivalent to x *= y * z;
               z = (x - y) + y ; // not equivalent to z = x;
               z = x + x * y;    // not equivalent to z = x * (1.0 + y);
               y = x / 5.0;      // not equivalent to y = x * 0.2;





       5.1.2.3                  Environment                 5.1.2.3
\f



       16           Committee Draft  --  August 3, 1998   WG14/N843


       [#15]  EXAMPLE 6 To  illustrate  the  grouping  behavior  of
       expressions, in the following fragment

               int a, b;
               /* ... */
               a = a + 32760 + b + 5;

       the expression statement behaves exactly the same as

               a = (((a + 32760) + b) + 5);

       due to the associativity and precedence of these  operators.
       Thus,  the result of the sum (a + 32760) is next added to b,
       and that result is then added to  5  which  results  in  the
       value  assigned  to  a.   On  a  machine  in which overflows
       produce an explicit trap and in which the  range  of  values
       representable   by   an   int   is   [-32768,  +32767],  the
       implementation cannot rewrite this expression as

               a = ((a + b) + 32765);

       since if the values for a and b were,  respectively,  -32754
       and  -15,  the  sum  a + b  would  produce  a trap while the
       original expression would not; nor  can  the  expression  be
       rewritten either as

               a = ((a + 32765) + b);
       or
               a = (a + (b + 32765));

       since  the values for a and b might have been, respectively,
       4 and -8 or -17 and 12.  However,  on  a  machine  in  which
       overflow  silently  generates  some value and where positive
       and  negative  overflows  cancel,   the   above   expression
       statement  can  be rewritten by the implementation in any of
       the above ways because the same result will occur.

Re: /internet

[Joe Buck]

> > Does IEEE require left-to-right evaluation of y = a * b * c * d;
> > when there are no parentheses present?
> 
> I think that IEEE does not actually say anything about it.  C9X, however,
> does contain some remarks such as the following, which might have come
> from NCEG.

>                int a, b;
>                /* ... */
>                a = a + 32760 + b + 5;
> 
>        the expression statement behaves exactly the same as
> 
>                a = (((a + 32760) + b) + 5);
>
>        due to the associativity and precedence of these  operators.
>        Thus,  the result of the sum (a + 32760) is next added to b,
>        and that result is then added to  5  which  results  in  the
>        value  assigned  to  a.   On  a  machine  in which overflows
>        produce an explicit trap and in which the  range  of  values
>        representable   by   an   int   is   [-32768,  +32767],  the
>        implementation cannot rewrite this expression as
> 
>                a = ((a + b) + 32765);

Amazing.  These guys are trying to turn C into Ada.

If this is the rule, then a*b*c*d can't be rearranged in C9X.  However,
it appears that it can be rearranged in Fortran (where the user must
use parentheses to force the order of evaluation) as well as in the
current ANSI/ISO C.

Re: /internet

[Chip Salzenberg]

According to Joe Buck:
> NCEG:
> >        Thus,  the result of the sum (a + 32760) is next added to b,
> >        and that result is then added to  5  which  results  in  the
> >        value  assigned  to  a.   On  a  machine  in which overflows
> >        produce an explicit trap and in which the  range  of  values
> >        representable   by   an   int   is   [-32768,  +32767],  the
> >        implementation cannot rewrite this expression as
> > 
> >                a = ((a + b) + 32765);
> 
> Amazing.  These guys are trying to turn C into Ada.

Hey, it only applies to machines where integer overflows trap.  Would
it be a problem to say that EGCS targets should never include such
environments?  Integer overflow traps are almost unheard of in C.
-- 
Chip Salzenberg      - a.k.a. -      <chip@perlsupport.com>
      "When do you work?"   "Whenever I'm not busy."

Re: /internet

[Jeffrey A Law]

  In message < 19981216162248.Z22090@perlsupport.com >you write:
  > According to Joe Buck:
  > > NCEG:
  > > >        Thus,  the result of the sum (a + 32760) is next added to b,
  > > >        and that result is then added to  5  which  results  in  the
  > > >        value  assigned  to  a.   On  a  machine  in which overflows
  > > >        produce an explicit trap and in which the  range  of  values
  > > >        representable   by   an   int   is   [-32768,  +32767],  the
  > > >        implementation cannot rewrite this expression as
  > > > 
  > > >                a = ((a + b) + 32765);
  > > 
  > > Amazing.  These guys are trying to turn C into Ada.
  > 
  > Hey, it only applies to machines where integer overflows trap.  Would
  > it be a problem to say that EGCS targets should never include such
  > environments?  Integer overflow traps are almost unheard of in C.
We'd have to introduce some new code to handle such an architecture if
one was to ever come along.

Right now GCC assumes integer instructions do not trap on overflow and wrap
in the expected manner.  As long as gcc is presented with that kind of
architecture it can reassociate integer operations in useful ways.

Some of gcc's existing targets have instructions which trap on overflow, but
they also include corresponding instructions which do not trap.  GCC uses
the non trapping variants.



jeff

Re: /internet

[Jeffrey A Law]

  In message < 199812161958.LAA29146@atrus.synopsys.com >you write:
  > > > Does IEEE require left-to-right evaluation of y = a * b * c * d;
  > > > when there are no parentheses present?
  > > 
  > > I think that IEEE does not actually say anything about it.  C9X, however,
  > > does contain some remarks such as the following, which might have come
  > > from NCEG.
  > 
  > >                int a, b;
  > >                /* ... */
  > >                a = a + 32760 + b + 5;
  > > 
  > >        the expression statement behaves exactly the same as
  > > 
  > >                a = (((a + 32760) + b) + 5);
  > >
  > >        due to the associativity and precedence of these  operators.
  > >        Thus,  the result of the sum (a + 32760) is next added to b,
  > >        and that result is then added to  5  which  results  in  the
  > >        value  assigned  to  a.   On  a  machine  in which overflows
  > >        produce an explicit trap and in which the  range  of  values
  > >        representable   by   an   int   is   [-32768,  +32767],  the
  > >        implementation cannot rewrite this expression as
  > > 
  > >                a = ((a + b) + 32765);
  > 
  > Amazing.  These guys are trying to turn C into Ada.
  > 
  > If this is the rule, then a*b*c*d can't be rearranged in C9X.  However,
  > it appears that it can be rearranged in Fortran (where the user must
  > use parentheses to force the order of evaluation) as well as in the
  > current ANSI/ISO C.
Someone is starting to mix integer vs FP.

For integer the reassociation rules are much simpler.  Our current target
set does not trap on overflow and overflows wrap in the expected manner.

For such targets aggressive integer reassociation is a valid transformation
(as long as the language does not mandate overflow checking at this level).

The rules for FP are different becuase it's not 100% certain that the
results after reassociation will be the same as before reassociation.  Though
I believe in the case reassociating a series of multiplies we are safe.

Again, I have no intention of reassociating a + 5 + b + c for FP because of
overflow concerns.  The same restrictions are not necessary for multiplies
though as far as I can tell.

I challenge anyone to come up with a case where a reassociation of
a * b * c * d  produces different results than ((a * b) * c) * d.


jeff

Re: /internet

[Per Bothner]

> The rules for FP are different becuase it's not 100% certain that the
> results after reassociation will be the same as before reassociation.  Though
> I believe in the case reassociating a series of multiplies we are safe.

How about:
	x = "semi-Infinity"; /* something that when squared will overflow */
	(x * x) * (1/x) /* +Infinity */
	x * (x * (1/x)) == x * 1.0 (approximately) == x;

	--Per Bothner
Cygnus Solutions     bothner@cygnus.com     http://www.cygnus.com/~bothner

Re: /internet

[Joern Rennecke]

> We'd have to introduce some new code to handle such an architecture if
> one was to ever come along.

Well, if the machine can support C at all, it can support unsigned integers.
So the worst case would be to implement all signed arithmetic with
unsigned arithmetic.

Re: /internet

[Joern Rennecke]

> The rules for FP are different becuase it's not 100% certain that the
> results after reassociation will be the same as before reassociation.  Though
> I believe in the case reassociating a series of multiplies we are safe.

Only if the exponent freely overflows and underflows.  However, often you'll
get Inf for an overflow and 0 for an underflow.
And if you come close to an underflow, you might get de-normalized numbers.

> Again, I have no intention of reassociating a + 5 + b + c for FP because of
> overflow concerns.  The same restrictions are not necessary for multiplies
> though as far as I can tell.
> 
> I challenge anyone to come up with a case where a reassociation of
> a * b * c * d  produces different results than ((a * b) * c) * d.

(DBL_MAX * (DBL_MIN * DBL_MIN)) * 4 is different from
((DBL_MAX * DBL_MIN) * DBL_MIN) * 4

Re: /internet

[Sylvain Pion]

On Wed, Dec 16, 1998 at 05:34:47PM -0700, Jeffrey A Law wrote:
> [...] Though
> I believe in the case reassociating a series of multiplies we are safe.
>
> I challenge anyone to come up with a case where a reassociation of
> a * b * c * d  produces different results than ((a * b) * c) * d.

If you want that, you'd better first _prove_ that it works in every case.
But you won't be able to do so because it's false.

It might be true when the following both conditions are met:
- no intermediate overflow/underflow appears (nor denormalized numbers).
- the rounding mode is set to nearest.

I'd be interested in reading the proof in this case (if it's true).

Case if you have an intermediate overflow is easy:
DBL_MAX * DLB_MAX * DBL_MIN  gives you either Inf or DBL_MAX.
Case when rounding is set to +Infinity:
0.1 * 0.1 * -1.0  doesn't give you the same result, try it.

-- 
Sylvain

Re: /internet

[Horst von Brand]

Joern Rennecke <amylaar@cygnus.co.uk> said:

[...]

> Well, if the machine can support C at all, it can support unsigned integers.
> So the worst case would be to implement all signed arithmetic with
> unsigned arithmetic.

Who says unsigned arithmetic isn't implemented by splitting each word in
half and computing the result from there to avoid overflows?

(Yes, I love to spoilsport :)
-- 
Dr. Horst H. von Brand                       mailto:vonbrand@inf.utfsm.cl
Departamento de Informatica                     Fono: +56 32 654431
Universidad Tecnica Federico Santa Maria              +56 32 654239
Casilla 110-V, Valparaiso, Chile                Fax:  +56 32 797513

Re: /internet

[Craig Burley]

>If this is the rule, then a*b*c*d can't be rearranged in C9X.  However,
>it appears that it can be rearranged in Fortran (where the user must
>use parentheses to force the order of evaluation) as well as in the
>current ANSI/ISO C.

Thus preserving Fortran as the language of choice for scientific/
numeric programming, at least over anything stemming from C, since
the Fortran code can be implemented faster than the equivalent C code.

So it strikes me as strange that they would impose these rules, given
how they (same people, I *think*) have been trying so hard to wedge
in new gimmicks to make C *potentially* as fast as Fortran, like
`restrict'.  (Which is much harder to use than writing ((a*b)*c)*d,
AFAICT.)

I wouldn't be surprised to see these rules dropped, if there's still
time, so we probably shouldn't take them as the gospel according to C9X
until it's reached the stage where substantial changes are no longer made
to the draft standard (though I have no idea where the process is at
this point, myself).

        tq vm, (burley)

Re: /internet

[Jeffrey A Law]

  In message < 199812170055.QAA14233@cygnus.com >you write:
  > > The rules for FP are different becuase it's not 100% certain that the
  > > results after reassociation will be the same as before reassociation.  Th
  > ough
  > > I believe in the case reassociating a series of multiplies we are safe.
  > 
  > How about:
  > 	x = "semi-Infinity"; /* something that when squared will overflow */
  > 	(x * x) * (1/x) /* +Infinity */
  > 	x * (x * (1/x)) == x * 1.0 (approximately) == x;
Ah.  Shit.  That's what I was missing.  

Bummer.  This may severely restrict what we can do in this space for floating
point since we can't use the as-if rule.  Oh well...  

jeff

Re: /internet

[Jeffrey A Law]

  In message < 9812162034.AA17456@nile.gnat.com >you write:
  > Another useful difference is that in Ada evaluations of
  > functions in "pure" packages are guaranteed to be free of 
  > side-effects and may reuse results of earlier invocations 
  > with the same parameters. GCC takes full advantage of this 
  > by moving them out of loops etc. Of course all elementary 
  > fpt functions (including complex ones) are defined in pure 
  > packages and this leads to nice optimizations without sacrifying 
  > one bit of accuracy.
Right.  You just set CONST_CALL_P and all the right things will happen;  most
of the optimizers know about pure functions.   It was really lame of C not to
have such a concept from day 1.

So to make use of that feature in C we're stuck with attributes.  Though maybe
c9x has a reasonable syntax for declaring functions as pure functions.


jeff

Re: /internet

[Harvey J. Stein]

"Geert Bosch" <bosch@gnat.com> writes:

 > Those guys should copy these features instead of ones they make up!

I especially like the part about how it can rearrange & allow overflow
& underflow to cancel, but if it traps it can't rearrange at all.
Makes for nice portable behavior.

-- 
Harvey J. Stein
BFM Financial Research
hjstein@bfr.co.il

Re: /internet

[Geert Bosch]

On Wed, 16 Dec 98 11:58:08 PST, Joe Buck wrote:

  Amazing.  These guys are trying to turn C into Ada.

If that's true, they are doing it all wrong! :-)
In Ada95 rules have be precisely chosen such to allow 
useful optimizations like these. See my previous message 
where I explain how checks (and as a result operations) 
may be removed. 

Another useful difference is that in Ada evaluations of
functions in "pure" packages are guaranteed to be free of 
side-effects and may reuse results of earlier invocations 
with the same parameters. GCC takes full advantage of this 
by moving them out of loops etc. Of course all elementary 
fpt functions (including complex ones) are defined in pure 
packages and this leads to nice optimizations without sacrifying 
one bit of accuracy.

Those guys should copy these features instead of ones they make up!

Regards,
   Geert

/internet

[tprince]

>T: I think what you are getting at is that it's usually acceptable for the
>results to be calculated in the declared precision; extra precision is usually
>desirable, but unpredictable combinations of extra precision and no extra
>precision may be disastrous.  See Kahan's writings about the quadratic
>formula.  Your proposal would make an improvement here.

>C:That feedback is helpful, and does seem to reflect what I was trying to
>say originally.  (I haven't seen Kahan's writings, or at least very little
>of them, at this point.)

T:  Look at http://http.cs.berkeley.edu/~wkahan/ieee754status/ieee754.ps
I figured out how to do his quadratic algorithm in C with volatiles on the
R8K and PowerPC (neither of which I use any more) but it needed function
calls to do it in g77.  If someone gets the fused MACs going for hppa2.0,
those issues will come up there.  It really is possible to find these fused
MAC's introducing NaN's in a program which works correctly otherwise,
and I've not been able to track them down in a program bigger than that
quadratic formula thing.  The R10K has a gentler form of fused MAC
where the behavior has been made to be generally the same as with
individual IEEE-754 compliant operations.  Kahan didn't analyze what
might happen with unfavorable combinations of two rounding modes
on an Intel-like processor; I'm suspecting it would not be good news.
Has anyone looked into this?

>>C: I think a substantial portion of the audience asking for REAL*16 is
>*non-Intel*.  SPARC and Alpha people come to mind.  I agree that those
>who want enough extra precision to more reliably compute 64-bit results
>from 64-bit inputs would likely prefer the faster, native support
>provided by REAL*10 on Intel, and ideally "we" (g77/egcs/whatever) would
>be able to provide REAL*10 somewhat faster than REAL*16 on other machines
>as well, even though, unlike on Intels, the REAL*10 would be emulated.

T:  There are 2 major varieties of REAL*16.  The one which HP (and, I believe,
Sun Lahey and DEC) use is the more accurate and slower one, which conforms
nominally to IEEE P854 and has roughly the same exponent range as the
Intel REAL*10.  SGI and IBM use a faster version, which is facilitated by the
fused Multiply-accumulate instructions, which has roughly 6 fewer bits of
precision, a range less than that of double precision, and doesn't conform to
IEEE P854.

T; In both the HP and SGI libraries, the math functions give up
accuracy so as not to lose as much speed, so it is possible in either case to
wind up with little more accuracy than you would get with a carefully
implemented REAL*10.  I don't know about the other vendors' libraries.  On
the pentiums, some of the math functions inherently take advantage of the
full precision (log() but not log10(), sqrt(), sin()/cos(), tan(), atan()), while a
few require more of the style of programming found in non-Intel math
libraries, but with asm() mixed in, putting the proper usage of clobbers to
the test.

>>C:I don't think aligned spills happen reliably at all on any *released*
>version of egcs or gcc yet (well, except maybe for old versions of
>gcc patched with those big g77 patches that *seemed* to do most of the
>aligned-double thing).  But it looks like egcs 1.2 or 1.3 will align
>doubles on the stack, covering spills, at or near a rock-solid level
>of reliability.

T: Treatment of spills in general seems to be one area where gnu has some
room for improvement, in comparison to commercial compilers, particularly
for Intel.  I'm sure amazed that Lahey lost track of their alignments for lf95,
but they seem otherwise to be able to avoid spill performance problems.

>>C: if the
>compiler decides to call library (or inline) functions for constructs
>not explicitly, in the code, involving such calls, and those functions
>are not 80-bit, the result might indeed be similar to spilling to 64-bit
>values in that the programmer doesn't expect a sudden loss of precision
>there.

>>C:I'm thinking, for example, of complex divides, which g77 implements
>by avoiding the back end's version and going straight for c_div (or
>whatever) in libF77, to support a larger domain of inputs with
>greater accuracy.

T:  There, of course, the straightforward use of extended precision takes
care of the situation more effectively, where special-case coding is needed
otherwise.  But that can be done by using conditional compilation inside
c_div, according to whether the target architecture has long double of
greater precision and range than double.

>>C:Though, in this example, the loss of precision is a bit easier to
>predict: it currently happens for complex divides.  Someday, though,
>we might decide to have it apply to complex multiplies, and/or it
>might be desirable to have the compiler choose, based on less visible
>data (than the source code) to do a call rather than in-line the code.
>It's important to preserve the precision in such cases.

T:  It's more a question of avoiding unexpected exceptions.  The overhead
of the function call is not a serious matter for c_div, but it could be for
multiplication.  I looked up some of the implementations when you
brought this up over a year ago, and the only one I found which takes
special precautions on complex multiplication was VAX/VMS.  It's
needed more on VAX floating point, as even with the precautions, the
range of working operands is less than with IEEE floating point and no
special precautions.

Dr. Timothy C. Prince
Consulting Engineer
Solar Turbines, a Caterpillar Company
alternate e-mail: tprince@computer.org

           To:                                              INTERNET - IBMMAIL
                                                            N3356140 - IBMMAIL