Goldberg log1p vs gsl_log1p

Question

I am looking for a simple portable implementation of log1p. I came across two implementations.

The first appears here as Theorem 4 http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html ,

Implementation of the above

double log1p(double p) { volatile double y = p; return ( (1 + y) == 1 ) ? y : y * ( log( 1 + y) / ( ( 1 + y) - 1 ) ); } 

Second in GSL http://fossies.org/dox/gsl-1.16/log1p_8c_source.html

 double gsl_log1p (const double x) { volatile double y, z; y = 1 + x; z = y - 1; return log(y) - (zx)/y ; /* cancels errors with IEEE arithmetic */ } 

Is there a reason to prefer one over the other?

Answer 1

I tested these two approaches using a log() implementation with a maximum error of <0.53 ulps, compared to a multi-point arithmetic library. Using this implementation of log() as the building block for the two variants of log1p() , I found that the maximum error in the Goldberg version is <1.51 ulps, while the maximum error in the GSL variant was <0.67 ulps. This indicates that the latter will be significantly more accurate.

In terms of handling special cases, the Goldberg variant showed one inconsistency in which it returns NaN for input + infinity, while the correct result is + infinity. For special cases with the GSL implementation, there were three inconsistencies: inputs -1 and + of infinity delivered NaN, and the correct results should be -infection and + infinity, respectively. In addition, for entering -0, this code returns +0, while the correct result is -0.

It is difficult to evaluate performance without knowledge of resource allocation. As others noted in the comments, the Goldberg version is potentially faster when many arguments are close to zero, as it misses the expensive log() call for such arguments.

Answer 2

There is no right answer between them. The GNU Science Library is fairly reliable, well used, and actively supported in all later versions of gcc. You are unlikely to encounter many surprises with its use. As for any other code that you type, there is absolutely no reason not to use it after you have checked its logic and are satisfied with its level / method of error checking. GSL's slight drawback is that it is another library that you should carry with you, and depending on how widespread use of your code can become more difficult for other users on other platforms. It concerns its size.

The best piece of code is the one that most closely matches the requirements of your project.