[Bug-gsl] enforcing the order of numerical operations

[Eugene Loh]

I'm playing with GSL and trying to get it to "behave correctly" with
aggressive optimizations turned on in Sun Studio compilers.

Generally, reorganizing arithmetic operations according to traditional 
mathemetical semantics is fine.  In a handful of spots, however, careful 
treatment of arithmetic ordering is important for handling limitations 
of machine arithmetic.  What I would like is to express those careful 
orderings specially so that the compiler will have liberal license to 
optimize the bulk of the code.

One practice I've seen in GSL is to use the volatile keyword to enforce 
certain orderings.  For example, in sys/log1p.c, I find:

     double gsl_log1p (const double x)
     {
       volatile double y;
       y = 1 + x;
       return log(y) - ((y-1)-x)/y ;
     }

While conventional mathematical semantics would allow us to return

    log(y)-((y-1)-x)/y = log(y)-(x-x)/y = log(y)

the code stores 1+x into a volatile variable to force the intermediate 
computation.

This trick seems to work, but I assume there must also be an expensive 
memory flush in there.

There are other places in GSL where careful orderings must be observed. 
 I am somewhat motivated to look for them, but I wanted to get an idea 
from you what sort of fixes you'd be open to putting back into the code. 
 Is this "volatile" trick acceptable?

Another possibility is using a dummy routine.  E.g., something like

    double d_dummy(double x) { return x; }

Then, you could write stuff like

     double gsl_log1p (const double x)
     {
       double y;
       y = 1 + x;
       return log(y) - ((d_dummy(y)-1)-x)/y ;
     }

This would save a memory flush and in some cases improve readability of 
the code.

Anyhow, I'm looking for guidance as to what you'd support in putbacks. 
Do you favor the "volatile" approach?  Are you open to the dummy() function?