Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!

implementing LOG1P

zp3
Beginner
216 Views

Hi,

As everyone knows, log(1+dx) suffers heavily on cancellation errors when dx gets sufficient small. The solution is, as in most other programming languages already implemented, the use of the log1p() function. One common way to implement a workaround DIY is the rewriting of log(1+dx) as:

log1p=log(1+dx)*dx/((1-dx)-1)

The problem with this now is, that most compilers probably optimize it simply back to log(1+dx). Therefore I'd like to ask if there is a safe way to tell the compiler to just not touch this line of code?

Thanks for any advice!

0 Kudos
2 Replies
mecej4
Black Belt
216 Views

The formula that you listed has a sign error. Rather than use a trick, why not use the C-library function log1p(), or write code such as

if(abs(dx) < 1d-8)then
   log1p = dx*(1d0 - 0.5d0*dx)
else
   log1p = log(1d0 + dx)
endif

See also https://www.johndcook.com/blog/2010/06/07/math-library-functions-that-seem-unnecessary/ .

TimP
Black Belt
216 Views
Most compilers have options to respect parentheses, such as ifort -standard-semantics .
Reply