[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Bug-gsl] enforcing the order of numerical operations

From: Eugene Loh
Subject: [Bug-gsl] enforcing the order of numerical operations
Date: Wed, 25 Jul 2007 22:59:14 -0700
User-agent: Mozilla/5.0 (X11; U; SunOS sun4u; en-US; rv:1.7) Gecko/20041221

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?

reply via email to

[Prev in Thread] Current Thread [Next in Thread]