[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?
`

**[Bug-gsl] enforcing the order of numerical operations**,
*Eugene Loh* **<=**