Implementation of the log1p (x) function (continued)

At that time, a GSL log1p (x) code implementation was provided. However, I didn't want to understand why it was like that. I suddenly understood it when I read a post on shuimu recently.

The Code on GSL is as follows:

[Cpp] view plaincopyprint? Double gsl_log1p (const double x)
{
Volatile double y, z;
Y = 1 + x;
Z = y-1;
Return log (y)-(z-x)/y;/* cancels errors with IEEE arithmetic */
}

Double gsl_log1p (const double x)
{
Volatile double y, z;
Y = 1 + x;
Z = y-1;
Return log (y)-(z-x)/y;/* cancels errors with IEEE arithmetic */
}

We know that floating point numbers on computers only have limited precision.

Therefore, the value assignment statement "y: = 1 + x" is not exactly equal to 1 + x when x is calculated in hours.

Or we can think that what we actually calculate is

Y: 1 + x'

X' is the result of the loss of the Number of valid digits of x in the decimal place when the value of 1 + x is calculated.

Therefore, directly calculating log (1 + x) actually calculates log (1 + x ').

Since x is very similar to x', the following formula is similar:

 

 
 

The above formula is changed to the following:

 

This is the code used on GSL (x' corresponds to z in the Code ).

The following statement cannot be written in a C program:

Z = 1 + X-1;

Otherwise, the compiler will intelligently optimize it to z = x, and it does not seem to be able to turn off this optimization through the compiler's compilation options.

In addition, there is a very famous article, What Every Computer Scientist shocould Know About Floating-Point Arithmetic, which also discusses this computing. Another calculation method is provided.

 

This algorithm is based on the fact that when x is very close to x,

 

The value of x (x'-x) is very small, so the error of this algorithm is very small.